基于多特征融合和水平集的碧根果品质检测

刘 哲 邹小波 宋余庆 王 明 苏 骏

(1.江苏大学计算机科学与通信工程学院， 镇江 212013； 2.江苏大学食品与生物工程学院， 镇江 212013)

0 引言

碧根果也叫薄壳山核桃，具有良好的营养价值与经济价值[1]。但碧根果在生产加工过程中易酸败，并产生刺鼻的气味, 降低了其商品价值和保健功效。一旦误食酸败的碧根果将会对人体健康产生多方面的危害，甚至致癌[2]。

目前，国内外检测坚果的技术主要集中在机器视觉、X射线、光谱分析和其他检测方法。研究的对象大多为整幅图像，并未针对目标有效区域。近年来，将水平集方法运用到几何活动轮廓模型中，可分割具有拓扑结构变化的目标，其在农业检测领域中占据着越来越重要的地位[3]。文献[4]以杂草﹑小麦和苹果为研究对象，提出了基于水平集和先验信息的计算机视觉(Computer vision, CV)模型，可分割复杂环境下的农业图像，为图像后续处理提供了新思路。文献[5]提出了一种无需重新初始化水平集函数、基于目标区域停止准则的水平集分割算法，由此实现完整和破损巴西坚果二分类，准确率高达99.55%。考虑到单一特征识别能力的局限性，研究者开始尝试多特征融合，实现更高精度的检测。文献[6]提取行人图像的色调饱和度值(Hue saturation value, HSV)和LAB颜色特征以及SILTP纹理特征和方向梯度直方图(Histogram of oriented gradient, HOG)特征,融合多种特征得到行人图像对的整体与局部相似度度量函数，并结合产生相似度函数,最后使用交替方向乘子优化算法更新出最优的测度矩阵，实现行人再识别,效果较好。文献[7]提出一种将颜色、形状、纹理特征与稀疏表示相融合的害虫识别方法。该方法利用已标注的训练样本构造不同特征下的训练样本矩阵，通过求解样本的最优稀疏系数实现害虫图像识别，该方法的害虫识别率获得较大的提高。文献[8]利用颜色特征与纹理特征融合的方式表示害虫图像样本，并通过有监督的字典训练方式构建害虫图像特征的完备字典，利用改进的SRC分类模型对害虫图像进行稀疏表示，并实现自动分类，分类准确率达到94.49%。文献[9]在不去除背景条件下，利用颜色和HOG特征融合训练获得的稀疏表示害虫识别模型，有效地克服了环境、光照和害虫姿态变化等引起的问题，获得了较高的农业测报害虫识别率和较低的误检率。

本文将提取的碧根果灰度直方图、灰度共生矩阵等多特征进行融合分析，结合水平集算法和多特征融合的优势，对碧根果酸败与否及酸败情况进行研究，并将其分为正常、轻度酸败、中度酸败和重度酸败4个等级，以期为坚果品质无损检测提供新的方案。

1 材料和方法

1.1 试验样本

本试验由江苏沃顿食品有限公司提供正常、酸败碧根果样本。样品横向直径为(2.5±0.3)cm，纵向直径为(4.6±0.2)cm，质量为(9.4±4)g。精挑选正常果、酸败碧根果各20个，并对样本进行编号，如表1所示：A1～A20(酸败的碧根果标签序号)，B1～B20(正常的碧根果标签序号)。采集碧根果横断面的T2(自旋-自旋或横向弛豫时间)加权核磁共振图像。图像分辨率为512像素×512像素，每个样本采集5幅图像，总计200幅。

1.2 核磁共振采集设备

MRI(磁共振成像)设备采用上海纽迈公司的中型尺寸核磁共振成像分析仪，其磁体类型为永磁体，磁场强度为(0.5±0.08)T，仪器主频率为21.3 MHz；探头线圈直径为60 mm，磁体温度为32℃。采集碧根果横断面的T2加权MRI图像。试验参数如下：切片厚度为2 mm，切片间距为0.3 mm，重复时间(TR)为1 600 ms，回波时间(TE)为46 ms。

由于采集管空间有限，为保证采集图像的质量，每组采集4个碧根果，每个样本果采集5幅横切面图像，每组采集耗时5～8 min，图像分辨率为512像素×512像素，总计200幅。

1.3 图像预处理

碧根果核磁共振图像分辨率虽然很高，但是碧根果完整轮廓区域面积占比非常小。因此本文对碧根果轮廓作最小外接矩形处理，以突出目标域，如图1所示。由于每个坚果的切片大小不同，得到的最小外接矩形存在差异，导致处理过的图像分辨率各不相同。因此本文采用平滑效果好且运算时间快的双线性插值实现统一的分辨率(190像素×180像素)。

1.4 碧根果图像酸败分析

通常情况下，带正电的原子核不会有自旋现象，但有些原子核(比如1H、19F和31P等)却可以自发地产生自旋运动。这类原子核自旋轴无规律，在一定的外来磁场等因素作用下将产生核磁共振现象。原子核从激化的状态恢复到平衡排列状态的过程叫弛豫过程。弛豫时间有两种：T1(自旋-点阵或纵向驰豫时间)和T2。由于T2加权成像表现在组织的T2越长，恢复越慢，信号就越强，成像部分就越亮；组织的T2越短，恢复越快，信号就越弱，成像部分就越暗。因此本文采用信号较强的氢原子核质子1H，在T2加权下采集碧根果核磁共振图像[10]。

酸败后的果仁油脂含量将减少，氢氧基酸含量也将随之减少，所以酸败后的果仁将呈现低灰度级别。为了直观地分析由于酸败引起的灰度变化情况，本文对酸败果和正常果各100幅图像统计了每幅切片图像的灰度折线图，考虑到背景为黑色(即灰度为0)的点较多且无意义，所以选择统计灰度根据图2在灰度区间1～100的曲线密集度可知，酸败果明显多于正常果，并且酸败果在灰度100的波峰处有多条曲线重叠即多幅切片图像相似。总体而言，正常果灰度主要集中在灰度为175的波峰附近，酸败碧根果灰度主要集中在灰度150～175的波峰附近，即正常果整体灰度较酸败果稍微偏大。

图1 MRI图像预处理过程Fig.1 MRI image preprocessing

区间为[1，255]的像素点，得到的试验结果如图2所示。

图3 部分样品果的疑似酸败域和正常域的灰度直方图Fig.3 Gray histogram of suspected rancid domains and normal domains in some samples

图2 部分样品果的疑似酸败域和正常域的灰度折线图Fig.2 Gray line chart of suspected rancid domains and normal domains in some samples

为了进一步分析图像中酸败域与正常域的区别，手动勾勒区块，并统计灰度直方图。红色框为手动勾选的酸败区块，蓝色框为手动勾选一个正常区块，若无酸败域，则红色框选中的仅作为对比区块不代表酸败区域。由此测定红蓝区块的灰度直方图，可得酸败果对比图如图3所示。

由图3可知，酸败域的灰度集中在50和100(波峰)处，而正常域灰度集中在150和175(波峰)处，碧根果果仁中酸败域灰度明显低于正常域，若能够自动分割出最大疑似酸败区块，将有利于最终的酸败分类识别。

1.5 水平集分割算法

文献[11-13]主要依赖外部能量项不断推进曲线向目标边缘附近靠近。其中基于图像梯度信息的边界指示函数对噪声敏感且难以分割边界模糊或灰度不均匀的图像。本文在DRLSE模型基础上，针对噪声敏感、易陷入局部极小值以及初始轮廓敏感的问题，提出一种改进边缘指示函数的自适应DRLSE方法，主要是改进面积项系数，使其可以自适应地向曲线内或曲线外演变，避免初始轮廓位置的影响；边缘指示函数中引入图像局部区域像素拟合值避免了噪声导致曲线演化到局部极值处的问题。

能量泛函定义为

ε(φ)=μRp(φ)+εext(φ)

(1)

式中Rp(φ)——内部能量项，也称距离规则项

μ——大于0的常数

εext(φ)——外部能量项，由长度项和面积项组成

φ——水平集函数

在曲线演化过程中，外部能量项εext(φ)不断推进曲线向目标边缘附近靠近，而距离正则项Fp(φ)则用于消除曲线与符号距离函数的偏差。距离正则项作为一个惩罚项，定义为

(2)

其中

P1——势阱函数

Ω——定义域空间

外部能量函数为

εext(φ)=λLg(φ)+aAg(φ)

(3)

式中Lg(φ)——曲线长度项

Ag(φ)——曲线内部区域面积，用于加快曲线的演化速度

λ、a——定值常数

当初始化轮廓曲线完全包围目标域，此时a必须为正数，以保证曲线向内收缩至目标边缘；反之，初始化轮廓曲线在目标域内，此时a必须为负数，以保证曲线向外扩张至目标边缘；但是当初始化轮廓曲线与目标域相交则无法正确分割，并且长度项和面积项中都引入了边缘指示函数g，导致曲线在演化过程中易陷入局部极值且对噪声敏感，边缘指示函数为

(4)

式中Gσ——标准差为σ的高斯核函数

I——待分割图像

本文改进了面积项系数，使曲线能够自适应地向内或向外演化，公式为

(5)

式中S(I)——曲线演化速率

A、b——调整S(I)的常数

TOTSU——图像卷积后的最大类间方差法所求分割阈值

本文改进的边缘指示函数为

(6)

其中

M=A[Kσ(1-f1(x))2Hε+Kσ(1-f2(x))2(1-Hε)]+b

(7)

(8)

(9)

式中Kσ——标准差为σ的高斯核函数

Hε——正则化的Heaviside函数

f1(x)和f2(x)用来拟合零水平集轮廓曲线内、外的图像局部区域灰度。

根据式(1)可得

(10)

δ(φ)——一阶变分

由此可得水平集函数的迭代公式为

(11)

式中 Δt——时间步长

D(·)——式(1)的离散化求解形式

i、j——像素值横纵坐标

k——当前时刻

改进边缘指示函数自适应的DRLSE方法步骤如下：①读入碧根果MRI图像。②初始化参数μ、λ、a、b、c、ε、σ与水平集函数Φ。③根据式(5)计算S(I)。④根据式(7)～(9)计算f1(x)、f2(x)、M。⑤根据式(11)更新水平集函数φ。⑥检测更新后的φ是否收敛，若收敛，则φ停止演化，所在位置记为目标边缘完成分割，否则转步骤③继续演化。⑦输出最终分割结果。

1.6 图像特征提取

特征一般分为颜色特征、形状特征和纹理特征等，由于本文研究对象是碧根果的核磁共振图像，因此本文主要关注灰度直方图统计特征与纹理特征[14-15]。在纹理特征中比较常用的特征量为：灰度共生矩阵、Tamura纹理特征、局部二值模式特征。本文级联各特征向量以获取最终的特征描述子。

1.6.1灰度直方图统计特征

灰度直方图描述图像中灰度出现的频率。如果一幅图像有L个灰度级(对于8位灰度图像，L=256)，则灰度直方图的数学表达式为

(12)

式中NG——灰度为G的像素数

N——图像总像素数

灰度直方图提供了图像外观最简单可视的全局描述，具有旋转、比例及位移不变性。经常使用的7个灰度直方图统计特征[16-17]为均值、均方差、偏斜度、信息熵、能量、峭度、平滑度。

1.6.2灰度共生矩阵

灰度共生矩阵就是从图像坐标(x,y)处的灰度为i的像素出发，统计距离为d、角度为θ、灰度为j的像素同时出现的概率。一般距离d为1、2、3，角度为0°、45°、90°、135°中任意一个。图像的纹理特征不同，其共生矩阵也会明显不同。一般采用4个最常用的特征来提取影像的纹理特征:角二阶矩、对比度、相关性、熵。

1.6.3Tamura特征

文献[18]从心理学的角度提出了一种纹理特征的描述方法(Tamura特征)，其中定义了6个用于描述纹理的特征量，分别为粗糙度(Coarseness)、对比度(Contrast)、方向度(Directionality)、线像度(Linelikeness)、规则度(Regularity)和粗略度(Roughness)。其中对比度、方向度和粗糙度3个定量分析指标对图像纹理的分析比较重要。

(1)粗糙度

反映图像灰度变化的剧烈程度，纹理基元尺寸越大则感觉越粗糙，其定义如下：计算图像中以坐标(x,y)为中心、大小为2l×2l的活动窗口,窗口中每个像素元的灰度均值Ak(x,y)为

(13)

式中l——参与均值计算的像素值

对于每个像元，分别计算其在水平及垂直方向上互不重叠的活动窗口之间的平均灰度差Ek,v、Ek,h,公式为

Ek,v=|Ak(x,y+2l-1)-Ak(x,y-2l-1)|

(14)

Ek,h=|Ak(x+2l-1,y)-Ak(x-2l-1,y)|

(15)

对于每一个像元，寻找最优的尺寸Skbest=2l，使其平均灰度达到最大Ekbest，即

Ekbest=argmax(Ek,v,Ek,h)

(16)

粗糙度Fcrs定义为M′×N′图像中各像元最优尺寸Skbest的平均值，即

(17)

(2)对比度

反映灰度图像中最暗和最亮灰度的层级，其差异范围决定对比度，此外灰度级的动态范围、边缘的尖锐程度、模式的重复周期等也会影响对比度。定义对比度Fcon为

(18)

其中

a4=μ4/σ4

(19)

式中a4——图像灰度统计量的峰态

μ4——四阶矩均值

(3)方向度

描述纹理沿某方向集中的强度。方向度的计算基于每个像素点的梯度向量，该向量的模和方向分别定义为

(20)

(21)

式中Δh、Δv——(x,y)处水平与垂直方向灰度差分

量化各像素点的方向角θ，统计其相应方向角上边缘像素数量Nθ(i)，挑选出边缘强度(模)大于预定阈值的边缘像素点，设其数量为Nθ(φ)，则相应方向上的直方图HD(φ)为

(φ=0,1,…,n-1)

(22)

则方向度定义为

(23)

式中np——峰的数量

wp——围绕该角度峰值与谷值之间的范围

r——与角度θ量化水平相关的归一化因子

φp——第p个峰值的位置

1.6.4局部二值模式

局部二值模式(Local binary patterns, LBP)最早是由芬兰Oulu大学机器视觉小组的研究者OJALA等[19-20]提出,是一种典型的局部纹理描述子，该方法在计算机视觉和模式识别领域引起了广泛的关注，目前已经成为纹理分类和人脸识别领域主要的特征提取方法之一。

LBP算子最初是定义在某中心像素及其周围大小为3×3的矩形邻域系统上，如图4所示：将中心像素的每个邻域像素值以该中心像素的灰度为阈值进行二值量化，大于或等于中心像素的像素值则编码为1，小于则编码为0，按照给定的编码方向依次对不同的邻域点赋予不同的权重，并将量化后得到的二值编码序列转换成一个无符号十进制数值，然后将该值作为对应像素点的LBP特征模式值。计算公式为

(24)

其中

式中gi——邻域像素点的像素值

gc——中心像素点的像素值

PS——采样点个数

图4 局部二值模式方法示意图Fig.4 Schematic of local binary mode method

2 试验结果与分析

本文试验运行环境配置：Core(TM)i3-3320 CPU @3.30 GHz,4 GB RAM，64位Windows 7操作系统，Matlab R2016a。

2.1 改进的自适应DRLSE方法分割结果

2.1.1碧根果完整轮廓分割

为验证本文改进算法相比于原始DRLSE算法具有初始轮廓不敏感且演化速度更快的优点，将本文算法与原始DRLSE算法进行试验对比。试验样本选择酸败的碧根果核磁共振图像A3-1，试验效果如图5所示，其中蓝色线为初始轮廓曲线，红线为最终分割结果。

图5 改进水平集分割图Fig.5 Improved level set segmentation graphs

由图5可知，原始DRLSE受面积项系数a的影响只能单向向内收缩或向外膨胀，当初始曲线与碧根果轮廓边缘相交时便无法完成正常分割，只有初始轮廓曲线包围目标域才勉强完成分割，效果相比于本文模型方法较差。本文所提模型无论初始曲线在目标内部或外部甚至相交都能很好地完成分割，如图5d所示。

2.1.2疑似酸败域分割

在1.4节分析了碧根果图像酸败特性后，本文在分割出碧根果完整轮廓基础上结合OTSU算法分割出疑似酸败域，如图6所示。

图6 疑似酸败域分割图Fig.6 Suspected rancid domain segmentation graphs

由图6可知，碧根果存在边缘模糊的情况，OTSU算法无法正常分割，导致异或处理后出现了很多边缘散点，如图6c所示，因此本文结合形态学开操作，去除散点如图6d所示。

为了证明提取“疑似酸败域”的方法可行性，试验样本选择之前的酸败果A2-3、A3-3、A8-3和正常果B6-1、B16-3、B20-5，对200幅MRI图像分析可知，提取其中面积占比较大的前4幅即可。倘若面积过小(像素点之和少于16)则舍去，试验结果分别如图7、8所示。

图7 酸败果处理示意图Fig.7 Diagram of processing rancid pecan

图8 正常果处理示意图Fig.8 Diagram of processing normal pecan

2.2 碧根果酸败识别

采用支持向量机[21]完成碧根果酸败分类识别。分别提取出疑似酸败域的灰度直方图统计特征、灰度共生矩阵特征、Tamura特征和局部二值模式特征。将提取的特征数据作标记0或1(0表示正常果，1表示酸败果)，随机选择其中160个数据作为训练集，余下40个数据作为测试集，一共3组，每组计算20次独立试验下分类准确率的平均值如表2所示。

由表2可知，Tamura特征分类准确率最低，灰度直方图统计特征较好，可达到88.53%，当融合其他所有特征时，酸败特征描述能力得到增强，高达96.15%。这也说明了本文提出的灰度直方图统计特征在融合其他纹理特征后，能够较好地满足碧根果酸败与否的检测。

表2 不同特征下碧根果酸败分类准确率Tab.2 Accuracy of pecan rancidity classification under different characteristics %

2.3 碧根果酸败程度判别

将碧根果酸败程度划分为4个等级：正常、轻度、中度、重度。试验中碧根果酸败的程度取决于酸败特征面积占整个碧根果有效域面积的比值，面积比越大说明酸败的程度越严重。在实际测量酸败程度上，将碧根果的果仁取出，将具有酸败的部分果仁(品尝具有哈喇味，外表发黑、泛油，可闻到异臭味均视为酸败)取出称量，测出其所占整个碧根果的质量比值并作为试验参考值。在已判别是否酸败的基础上，根据标签找出该碧根果的5幅切片图像，测出其酸败的面积，整个碧根果酸败程度的计算公式为

(25)

式中Pi——单个碧根果中第i幅切片图像被判别为酸败果的酸败特征面积占碧根果完整轮廓面积百分比；若其被判别为正常果，则Pi=0

根据式(25)，将酸败面积比大于0且小于等于30%视为轻度酸败；大于30%且小于等于50%视为中度酸败，大于50%视为重度。经过3次分类并测量酸败面积比所得结果的识别率如表3所示。

表3 碧根果酸败程度识别率Tab.3 Recognition rate of pecan rancidity degree %

由表3可知，正常和重度酸败果的识别率很高，这是因为“是否酸败”检测准确率高达96.15%，所以这个前提下再判别酸败程度会更准确。而轻度酸败的识别率相对低一些，主要由于疑似酸败域分割产生的偏差以及正常果被判别为酸败果时，酸败面积比偏高，识别率稍微降低。中度酸败识别率为89.47%，也存在此情况但稍高于轻度酸败果，主要因为存在酸败果的切片较正常果多。酸败程度的总体识别率平均值为90.81%。

3 结论

(1)基于核磁共振成像技术采集碧根果T2加权的图像，通过对比酸败、正常区块比灰度直方图，可直观地发现，碧根果酸败域和正常域的灰度存在明显差异。

(2)采用改进边缘指示函数的自适应DRLSE算法，结合OTSU算法和形态学方法，提取碧根果疑似酸败域。融合灰度直方图统计特征、GLCM、Tamura和LBP特征判别碧根果是否酸败。试验结果表明，融合后的多特征在酸败特性上描述能力强于单一特征，准确率高达96.15%。在此基础上，通过被判别为酸败果的疑似酸败域面积占比粗略计算碧根果酸败程度，碧根果酸败程度的总体识别率平均值为90.81%。

(3)本文方法能够自适应地分割出碧根果疑似酸败域，并建立了多特征融合描述子，提高了分类识别率，试验效果较好，为坚果类的内部品质无损检测提供了新思路。