基于水色图像的水质评价

背景
有经验的渔业生产的从业者可以通过观察水质变化调控水质，来维持养殖水体生态系统中的浮游植物、微生物、浮游动物等的动态平衡，然而这些判断是通过经验和肉眼观察得出的，存在主观性引起的观察性偏差，使观察结果的可比性、可重复性降低，不易推广使用。

数字图像处理技术为计算机监控技术在水产养殖业的应用提供了更大的空间。

在水质在线监测方面，数字图像处理技术是基于计算机视觉的，以专家经验为基础，对水色进行优劣分级，实现对水色的准确快速判别。

目标
利用给出的已经分类的数据，利用图像处理技术实现水质自动评价。

分析
通过拍摄得到的图像数据维度太大，不易分析，需要提取图像的特征，提取反映图像本质的关键指标，以达到自动进行图像识别或者分类的目的。

显然，图像特征提取是图像识别和分类的关键步骤，提取的结果直接关系到图像识别和分类的好坏。

图像特征
主要包括颜色特征、纹理特征、形状特征和空间关系特征等。

与几何特征相比，颜色特征更为稳健，对物体的大小和方向不敏感，表现出较强的鲁棒性。

由于本案例中的水色图像是均匀的，所以主要关注颜色特征。

颜色特征是一种全局特征，描述图像或者图像对应景物的表明性质。

一般颜色特征是基于像素点的特征，所以图像区域的像素点都有自己的贡献。

在利用图像颜色特征进行分析的研究中，实现方法上已经有啦很多研究成果，主要采用直方图法和颜色矩方法。

其中直方图法反映颜色分布，即哪些颜色及其出现概率，优点是适合描述难以自动分割的图像和不需要考虑物体空间位置的图像；缺点是无法描述颜色的局部分布及每种色彩所处的空间位置，也就是无法描述图像中的某一具体的对象和物体。

其中颜色矩法的数字基础为图像中的任何颜色分布均可以用它道德矩来表示。

根据概率论，随机变量的概率分布可以由其各阶矩唯一描述和表示。

一个图像的颜色分布完全可以看做是一种概率分布，图像就可以由各阶矩描述。

颜色矩包含各个颜色的一阶矩、二阶矩和三阶矩，一个RGB图像有R、G、B三个通道，共9个分量。

直方图法产生的特征维数一般会大于颜色矩产生的，为了避免过多变量干扰分类结果，使用颜色矩法。

数据预处理
图像切割
由于图像较大，我们提取每张图片最中间的100x100像素构成子图象集。

特征提取
我们采用颜色矩来提取水样特征。

这种方法的数学基础在于图像中任何的颜色分布均可以用它的矩来表示。

此外，由于颜色分布信息主要集中在低阶矩中，因此仅采用颜色的一阶矩（mean）、二阶矩（variance）和三阶矩（skewness）就足以表达图像的颜色分布。

与颜色直方图相比，该方法的另一个好处在于无需对特征进行向量化。

因此，图像的颜色矩一共只需要9个分量（3个颜色分量，每个分量上3个低阶矩），与其他的颜色特征相比是非常简洁的。

在实际应用中为避免低次矩较弱的分辨能力，颜色矩常和其它特征结合使用，而且一般在使用其它特征前起到过滤缩小范围（narrow down）的作用。

其中Ei是第i个颜色通道的一阶颜色矩，对于RGB颜色空间的图像，i=1,2,3，pij是第i个像素点的第j个颜色通道的颜色值。

在颜色矩的提取中，提取每个文件名中的类别和序号，同时针对所有图片进行同样操作。

#导入库文件
#import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
import os
import pandas as pd
#查看
#img = mpimg.imread('eeeee/chapter9/chapter9/test/data/images/1_1.jpg')#读取一张图片#l,w,t=img.shape #提取长，宽，通道数
#print(l, w, t)
#plt.imshow(img)
#寻找中心100x100像素点
#l_max =img.shape[0]//2+50
#l_min=img.shape[0]//2-50
#w_max =img.shape[1]//2+50
#w_min =img.shape[1]//2-50
#计算颜色矩特征
def img2vector(filename):
returnvect=np.zeros((1,9))
fr=mpimg.imread(filename)
l_max =fr.shape[0]//2+50
l_min=fr.shape[0]//2-50
w_max =fr.shape[1]//2+50
w_min =fr.shape[1]//2-50
water=fr[l_min:l_max, w_min:w_max, :].reshape(1, 10000, 3)
for i in range(3):
this = water[:,:,i]/255.
print this
returnvect[0,i]=np.mean(this) #0,1,2存储一阶颜色矩
returnvect[0,3+i]=np.sqrt(np.mean(np.square(this-returnvect[0,i])))#3,4,5存储二阶颜色矩
returnvect[0,6+i]=np.cbrt(np.mean(np.power(this-returnvect[0,i], 3)))#6,7,8存储三阶颜色矩
print returnvect
return returnvect
#计算每个图片的特征
trainfilelist=os.listdir('eeeee/chapter9/chapter9/test/data/images')#读取目录下文件列表
m=len(trainfilelist) #计算文件数目
labels=np.zeros((1,m))
train=np.zeros((1,m))
#trainingMat=[]
#print(trainfilelist)
trainingMat=np.zeros((m,9))
for i in range(m):
filenamestr=trainfilelist[i] #获取当前文件名，例1_1.jpg
filestr=filenamestr.split('.')[0] #按照.划分，取前一部分
classnumstr=int(filestr.split('_')[0])#按照_划分，后一部分为该类图片中的序列
picture_num = int(filestr.split('_')[1])
labels[0,i]=classnumstr #前一部分为该图片的标签
train[0,i]=picture_num
trainingMat[i,:]=img2vector('eeeee/chapter9/chapter9/test/data/images/%s' % filenamestr) #构成数组
#保存
d=np.concatenate((labels.T,train.T,trainingMat),axis=1)#连接数组
dataframe=pd.DataFrame(d,columns=['label','num','R_1','G_1','B_2','R_2','G_2','B_2','R_3','G_ 3','B_3'])
dataframe.to_csv('eeeee/chapter9/chapter9/test/moment.csv',header=None,index=None)#保存文件
# 把数据分为两部分：训练数据、测试数据
# 读入数据
Data <- read.csv("./data/moment.csv")
# 数据命名
colnames(Data) <- c("class", "id", "R1", "G1", "B1", "R2",
"G2", "B2", "R3", "G3", "B3")
# 数据分割
set.seed(1234) # 设置随机种子
# 定义序列ind，随机抽取1和2,1的个数占80%，2的个数占20%
ind <- sample(2, nrow(Data), replace = TRUE, prob = c(0.8, 0.2))
trainData <- Data[ind == 1, ] # 训练数据
testData <- Data[ind == 2, ] # 测试数据
# 数据存储
write.csv(trainData, "./tmp/trainData.csv", s = FALSE)
write.csv(testData, "./tmp/testData.csv", s = FALSE)
# 读取数据
trainData <- read.csv("./data/trainData.csv")
testData <- read.csv("./data/testData.csv")
# 将class列转换为factor类型
trainData <- transform(trainData, class = as.factor(class))
testData <- transform(testData, class = as.factor(class))
# 支持向量机分类模型构建
library(e1071) # 加载e1071包
# 利用svm建立支持向量机分类模型
svm.model <- svm(class ~ ., trainData[, -2])
summary(svm.model)
# 建立混淆矩阵
confusion <- table(trainData$class, predict(svm.model, trainData, type = "class")) accuracy <- sum(diag(confusion)) * 100 / sum(confusion)
# 保存输出结果
output_trainData <- cbind(trainData, predict(svm.model, trainData, type = "class"))
colnames(output_trainData) <- c("class", "id", "R1", "G1", "B1", "R2", "G2",
"B2", "R3", "G3", "B3", "OUTPUT")
write.csv(output_trainData, "./tmp/output_trainData.csv", s = FALSE)
# 保存支持向量机模型
save(svm.model, file = "./tmp/svm.model.RData")
# 读取数据
testData <- read.csv("./data/testData.csv")
# 读取模型
load("./tmp/svm.model.RData")
# 建立混淆矩阵
confusion <- table(testData$class, predict(svm.model, testData, type = "class")) accuracy <- sum(diag(confusion)) * 100 / sum(confusion)
# 保存输出结果
output_testData <- cbind(testData, predict(svm.model, testData, type = "class")) colnames(output_testData) <- c("class", "id", "R1", "G1", "B1", "R2",
"G2", "B2", "R3", "G3", "B3", "OUTPUT") write.csv(output_testData, "./tmp/output_testData.csv", s = FALSE)。