CN105701805A - 一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法 - Google Patents

Abstract

本发明涉及一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法。采用摄像机标定法对CCD数码摄像机进行标定，选取刚从屠宰场购买的猪肉眼肌作为实验样本，对其横截面进行拍照采样，用改进的样本块修复方法对得到的样本图像进行预处理，恢复图像反光区域原有信息。结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割，提取出猪肉大理石花纹。从得到的大理石花纹图像中提取出脂肪数量指标、脂肪分布指标和脂肪纹理指标等291个特征值，根据特征值和化学方法检测结果建立猪肉肌内脂肪含量预测模型，最后通过该模型对猪肉肌内脂肪含量进行预测。利用本发明能很好地预测出猪肉肌内脂肪含量，使得猪肉的营养检测具备客观性、准确性和高效性。

Description

一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法

技术领域

本发明涉及营养检测、模式识别和计算机视觉技术领域，尤其涉及一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法。

背景技术

肌内脂肪是影响猪肉品质的重要因素，肉的风味和多汁性随肌内脂肪含量的增加而持续改善，对于猪肉嫩度主要通过切断肌纤维束间的交联结构，有利于咀嚼过程中肌纤维的断裂，这样就改变了肉质的感官品质。所以在保证较高生长速度和瘦肉率的前提下也要维持适当的体脂和肌内脂肪水平。2％-3％的肌内脂肪含量对猪肉的食用特性比较理想。

传统的肌内脂肪含量检测一般使用索氏抽提法，可参考国标GB9695.1—1988《肉与肉制品—游离脂肪含量的测定》，使用乙醚、石油醚或正己烷抽提单块肌肉待测样中脂肪。该方法步骤繁琐，属破坏性检测，在抽提时使用的萃取剂如乙醚或石油醚等对人体也有危害。

大理石花纹是肌内脂肪的一种外在表象，刘强等人在文章《大理石花纹评分与淮南猪背最长肌感官特性的关系研究》中提到，大理石花纹评分越高，则肌内脂肪含量越高，因此通过测定猪肉背最长肌大理石花纹的丰富程度，可以间接得出肌内脂肪的含量。

目前，我国肉品大理石花纹分级标准主要采用人工测量和感官评定相结合的方法，这种方法是主观的、定性的，在实际生产中效率低、误差大、公平和公正性难以保证。

参考国内外猪肉分级技术研究的最新进展，釆用计算机视觉技术对猪肉进行分级被认为是最理想的解决方案。2005年彭增起等人在发明《一种猪肉颜色质量等级的划分方法》中根据色度学上推荐的颜色空间和色度计测定的肉色参数L、a、b值进行颜色分级，并做成猪肉颜色等级图谱；2009年成芳等人在发明《基于机器视觉的猪肉外观品质实时检测分级方法与装置》中对眼肌肉完成颜色、大理石花纹和嫩度以及综合品质的等级评定；2013年陈坤杰等人在发明《基于多光谱图像纹理特征的检测牛肉嫩度的装置及方法中》通过提取猪肉大理石花纹的纹理特征对猪肉的嫩度进行预测。上述发明中介绍的图像处理算法和检测装置没有考虑到样本表面反光对后续工作的影响，只是应用中值滤波等方法进行简单的预处理，同时选取的特征值有的是关于脂肪数量的，有的是关于脂肪分布情况的，有的是关于纹理情况的，但都没有对这些指标进行综合考虑。

发明内容

本发明的目的是提供一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法，为解决上述技术问题，本发明的主要技术内容如下：

一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法，包括以下步骤：

(1)标定摄像机，得到图像中单位像素对应于真实场景中的实际长度；

(2)选取猪肉第5-6根肋骨处背最长肌(俗称为眼肌)作为样本，将样本放置于黑色背景布上，并用白色LED灯对样本进行补光，利用数码照相机对样本横截面进行俯拍，获取猪肉眼肌样本图像，图像大小为1000*1000像素；

(3)利用改进的样本块修复方法对得到的样本图像进行预处理，恢复图像反光区域原有信息；

(4)结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割，提取出猪肉大理石花纹；

(5)从得到的大理石花纹图像中提取出脂肪数量指标、脂肪分布指标和脂肪纹理指标等特征值；

(6)根据特征值和化学方法检测结果建立逐步回归预测模型和Fisher判别模型，而后对模型进行检验和结果比较，最后通过该模型对猪肉肌内脂肪含量进行预测。

上述步骤(3)中利用样本块修复方法对得到的样本图像进行预处理包括以下步骤：

(a)针对图像反光区域特征，即反光区域中心像素点有较高的亮度(Value，V)值和较高的饱和度(Saturation，S)值，对图像逐像素点进行阈值分割以得到反光中心区域，分割公式如下：

$F(i,j)=\left\{\begin{array}{cc}255& V>thresh1,S>thresh2\\ 0& else\end{array}\right.$

其中的F(i，j)代表坐标为(i，j)的像素点所对应的颜色值，255对应的是白色，0对应的是黑色，thresh1和thresh2分别代表了某一固定的阈值，阈值大小视实际拍摄条件而定；

(b)对图像进行2次全方向膨胀，使得图像中相互靠近的反光中心区域相互合并，然后对膨胀后的图像进行2次全方向腐蚀，使得面积变大的反光中心区域还原为原始大小；

(c)对于面积较大(大于9个像素点)的反光中心区域，把这些区域的外接矩形作为图像的反光区域，而对于面积很小(小于9个像素点)的区域则进行简单的平滑处理；

(d)利用改进的样本块修复方法对反光区域进行图像复原，改进方法如下：

首先修补过程中每次只针对一个待修复区域进行修补，然后再复原其它修复块，直到全部修补完成为止；其次在搜索最佳匹配块的过程中采用二次匹配法进行搜索，第一次匹配是粗略匹配，通过隔行隔列粗略搜索最佳匹配块的中心点位置，每次扫描的数据量降为原图四分之一；第二次匹配是精确匹配，在第一次匹配搜索到的中心点位置的八邻域点内进行遍历搜索，以得到最后的最佳匹配块进行修复。

上述步骤(4)中结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割包括以下步骤：

(a)去除图像背景：将去除反光后的图像转换为灰度图像，根据眼肌图像的特征，转换的灰度值公式为：Gray＝G*0.6+B*0.4，其中Gray为图像灰度值，G为原始图像绿色通道值，B为原始图像蓝色通道值，运用大津法对得到的灰度图像进行整体阈值分割去除图像背景，得到图像P1；

(b)提提取眼肌中的脂肪部分：把P1均分为400(20*20)份并分别进行阈值分割，每一份小图像根据其图像特征为其选取不同的自适应阈值分割算法，选取的流程如下：(迭代法分割出的大理石花纹面积为S1，最大熵法分割出的大理石花纹面积为S2)

当图像块中包含背景区域时，选用迭代法对图像块进行自适应阈值分割；

当S2≤S1＜1.3*S2时，选用迭代法对图像块进行自适应阈值分割；

当S1≤S2＜1.3*S1时，选用最大熵法对图像块进行自适应阈值分割；

当S2≥1.3*S1时，选用迭代法对图像块进行自适应阈值分割；

当S1≥1.3*S2时，选用最大熵法对图像块进行自适应阈值分割；

每个小块分割完成后得到图像P2；

(c)去除肌间脂肪：对P2进行轮廓提取，选出周长最大的一个轮廓，对该轮廓进行填充，对填充后的区域进行形态学膨胀处理得到图像P3，对图像P2和P3进行差值运算，得到图像P4；

(d)对P4进行形态学膨胀处理得到图像P5，对P5进行轮廓提取，找出轮廓像素点个数小于10的轮廓区域，这些区域为噪声点，去除这些区域得到图像P6，对P5和P6进行与计算得到图像P7，此时的图像即为提取出的眼肌大理石花纹。

上述步骤(5)中从得到的大理石花纹图像中提取出脂肪数量指标特征值共8个，所述特征的详细表述如下：

F1有效眼肌面积：即为图像P7中的大理石花纹面积；

F2脂肪颗粒个数：夹杂于眼肌中的脂肪斑点，每一个联通部分为一个脂肪颗粒；

F3大颗粒脂肪密度：单位面积上拥有的大颗粒脂肪数量(个/cm²)；

F4中颗粒脂肪密度：单位面积上拥有的中颗粒脂肪数量(个/cm²)；

F5小颗粒脂肪密度：单位面积上拥有的小颗粒脂肪数量(个/cm²)；

F6大颗粒脂肪面积比：大颗粒脂肪总面积与有效眼肌面积的比例(％)；

F7中颗粒脂肪面积比：中颗粒脂肪总面积与有效眼肌面积的比例(％)；

F8小颗粒脂肪面积比：小颗粒脂肪总面积与有效眼肌面积的比例(％)；

其中脂肪颗粒按面积(A)划分成三类：小颗粒(L1)，A1≤3.72mm²(30个像素)；中颗粒(L2)，3.72＜A2≤14.88mm²(120个像素)；大颗粒(L3)，A3＞14.88mm²。

上述步骤(5)中从得到的大理石花纹图像中提取出脂肪分布指标特征值共3个，所述特征的详细表述如下：

F1脂肪颗粒均匀度：指的是脂肪颗粒分布变异系数，设有效眼肌图像有n行，每行的脂肪像素占该行总像素比例为w_i(i＝1,2,…,n)，平均值为，则脂肪颗粒均匀度C的计算公式如下：

$C=\frac{\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(w_i-\stackrel{\‾}{w})}^2}{n-1}}}{\stackrel{\‾}{w}}$

F2计盒维数；

F3信息维数；

上述步骤(5)中从得到的大理石花纹图像中提取出脂肪纹理指标特征值，本发明用的是灰度共生矩阵特征值，分别选取像素间距d＝1,2,3,4,5共五组距离、方向角度θ＝0°、45°、90°、135°共四个方向的14个特征参数，则特征值共有280个(5*4*14)，所述特征参数的详细表述如下：

设f(x，y)为一幅二维数字图像，灰度级别为n＝256，其灰度共生矩阵为p(i，j，d，θ)，其中i、j分别表示矩阵中的行与列，d表示像素间距，θ表示方向角度；

F1角二阶矩：

$W_1=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{p}^2(i,j,d,\mathrm{\θ})$

F2对比度：

$W_2=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\[{(i-j)}^2\×p^2(i,j,d,\mathrm{\θ})\]$

F3相关：

$W_3=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\[i\×j\×p^2(i,j,d,\mathrm{\θ})-u_1\×u_2\]/(d_1\×d_2)$

$\begin{array}{cc}{u}_1=\underset{i=1}{\overset{n}{\Σ}}i\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& u_2=\underset{j=1}{\overset{n}{\Σ}}j\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

$\begin{array}{cc}{d}_1^2=\underset{i=1}{\overset{n}{\Σ}}{(i-u_1)}^2\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& d_2^2=\underset{j=1}{\overset{n}{\Σ}}{(j-u_1)}^2\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

F4熵：

$W_4=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\×{\log }_{10}p(i,j,d,\mathrm{\θ})$

F5方差：

$W_5=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{(i-m)}^{2}p(i,j,d,\mathrm{\θ})$

其中，m为p(i，j，d，θ)的均值；

F6均值和：

$W_6=\underset{k=0}{\overset{2n}{\Σ}}k\×\underset{i=0}{\overset{n}{\Σ}}\underset{j=0}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})$

F7方差和：

$W_7=\underset{k=2}{\overset{2n}{\Σ}}{(k-W_6)}^2\×\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})$

F8逆差矩：

$W_8=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})/\[1+{(i-j)}^2\]$

F9差的方差：

$W_9=\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}{\[k-\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}k\×\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]}^2\×\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})$

F10和熵：

$W_{10}=-\underset{k=2}{\overset{2n}{\mathrm{\Σ}}}\[\underset{i=0}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]\×\log \[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]$

F11差熵：

$W_{11}=-\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}\[\underset{i=0}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]\×\log \[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]$

F12聚类阴影：

$W_{12}=-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\[(i-u_1)+(j-u_2)\]}^3\×p(i,j,d,\mathrm{\θ})$

$\begin{array}{cc}{u}_1=\underset{i=1}{\overset{n}{\Σ}}i\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& u_2=\underset{j=1}{\overset{n}{\Σ}}j\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

F13显著聚类：

$W_{13}=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[(i-u_1)+(j-u_2)\]}^4\×p(i,j,d,\mathrm{\θ})$

$\begin{array}{cc}{u}_1=\underset{i=1}{\overset{n}{\Σ}}i\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& u_2=\underset{j=1}{\overset{n}{\Σ}}j\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

F14最大概率：

$W_{14}=\underset{i,j}{MAX}\[p(i,j,d,\mathrm{\θ})\]$

上述步骤(6)中的逐步回归预测模型的具体实现步骤如下：

(a)根据猪肉眼肌实际肌内脂肪含量对猪肉眼肌样本进行打分，分为预测集和测试集；

(b)根据预测集样本的291个特征值建立逐步回归模型；

(c)将测试集样本代入建好的逐步回归预测模型；

(d)将逐步回归预测结果与实际肌内脂肪含量进行比较，得出预测准确率；

上述步骤(6)中的Fisher判别预测模型的具体实现步骤如下：

(a)根据猪肉眼肌实际肌内脂肪含量对猪肉眼肌样本进行打分，分为预测集和测试集；

(b)根据预测集样本的291个特征值建立Fisher判别预测模型；

(c)将测试集样本代入建好的Fisher判别预测模型；

(d)将Fisher判别预测结果与实际肌内脂肪含量进行比较，得出预测准确率。

本发明的优点

1、本发明采用改进的样本块图像修复方法，可以很好的复原图像反光区域原有信息，并在修复效率上优于原算法。

2、本发明结合最大熵法与迭代法对眼肌图像进行图像分割，可以精确提取出猪肉眼肌大理石花纹。

3、本发明提取了脂肪数量指标、脂肪分布指标和脂肪纹理指标等291个特征值，用逐步回归与Fisher判别对猪肉肌内脂肪含量进行预测，其中逐步回归模型对猪肉肌内脂肪含量的预测相对误差在5％以内，Fisher判别模型对猪肉肌内脂肪含量的预测相对误差在4％以内，都达到了很高的准确率。

附图说明

图1为本发明的方法流程图；

图2为反光区域及其邻域示意图；

图3为存在反光区域的猪肉眼肌图像；

图4为改进样本块算法对反光区域的复原效果；

图5为拍摄所得猪肉眼肌图像；

图6为提取出的猪肉眼肌大理石花纹。

具体实施方式

为更进一步阐述本发明为达成预定发明目的所采取的技术手段及功效，以下结合附图及较佳实施例，对依据本发明提出的其具体实施方式、结构、特征及其功效，详细说明如后。

如图1所示，为一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法，包括以下步骤：

(1)标定摄像机，得到图像中单位像素对应于真实场景中的实际长度；

(2)选取猪肉第5-6根肋骨处背最长肌(俗称为眼肌)作为样本，将样本放置于黑色背景布上，并用白色LED灯对样本进行补光，利用数码照相机对样本横截面进行俯拍，获取猪肉眼肌样本图像，图像大小为1000*1000像素；

(3)利用改进的样本块修复方法对得到的样本图像进行预处理，恢复图像反光区域原有信息；

(4)结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割，提取出猪肉大理石花纹；

(5)从得到的大理石花纹图像中提取出脂肪数量指标、脂肪分布指标和脂肪纹理指标等特征值；

(6)根据特征值和化学方法检测结果建立逐步回归预测模型和Fisher判别模型，而后对模型进行检验和结果比较，最后通过该模型对猪肉肌内脂肪含量进行预测。

上述步骤(1)中的摄像机标定包括以下步骤：

(a)选取一个标定用的棋盘，该棋盘为6*8黑白正方格棋盘，每一正方格的边长为3cm；

(b)通过摄像头对放置在固定位置的棋盘进行10个不同角度的拍摄，共得到10幅棋盘图像；

(c)以OpenCV计算机视觉库为基础，调用cvFindChessboardCorners()函数找到其中的角点，并标记出来，得到10幅图像的角点，然后调用cvCalibrateCamera2()函数即可得到摄像头的内外参数；

(d)图像中物体的实际尺寸即可通过以下公式计算得出：

$L=\frac{l*d}{f}$

其中l为物体在图像中所占的像素个数，d为摄像头到物体的距离，f即为内参数矩阵中的第一个参数值，物体的实际尺寸L与距离d单位相同。

上述步骤(3)中利用样本块修复方法对得到的样本图像进行预处理。基于样本块的图像修复是从缺损区域的边缘开始逐步向内进行修复，所以修复顺序的确定将和最终得到的修复效果、质量紧密相连。由传统的纹理合成图像修复技术知，待修复的样本区域内包含有效信息越多，则具有更高的修复优先级。在借鉴传统纹理合成技术的同时，引入基于偏微分方程思想的等照度线强度约束修复顺序，能够在保持纹理信息的同时更好地保持图像的线性结构。算法实现过程如下：

(1)区域建模

对反光区域的建模如图2所示，图中Ω表示待修复的反光区域，表示反光区域的边界线，ψ_p是以位于反光区域边界线上的p点为中心的待修复样本块，n_p作为p点的法线矢量，表示p点的等照度线方向，φ表示信息完好区域。

(2)计算待修补块的优先权

对于以轮廓上的点p为中心的目标块，定义其优先权为：

P(p)＝C(p)D(p)

其中：

$C\left(p\right)=\frac{\underset{q\∈{\mathrm{\ψ}}_p\∩\mathrm{\Φ}}{\Σ}M\left(q\right)}{\left|{\mathrm{\ψ}}_p\right|}$

$D\left(p\right)=\frac{|\▿I_p^{\⊥}\·n_p|}{\mathrm{\α}}$

C(p)称为置信度项，用于衡量目标块ψ_p内可靠信息的数量。初始时，M(q)＝0，否则M(q)＝1。D(p)称为数据项，确保了邻域内线性结构强度越大的p点优先得到修复，表示p点的等照度线的方向和强度，n_p表示待修补区域轮廓的单位外法线方向，α为归一化因子。

(3)搜索最佳匹配块

确定具有最佳优先权的目标块ψ_P后，就要在已知区域内搜索此目标块的最佳匹配块ψ_Q，目标块与最佳匹配块有如下匹配准则：

${\mathrm{\ψ}}_Q=\arg \underset{{\mathrm{\ψ}}_q\∈\mathrm{\φ}}{min}d({\mathrm{\ψ}}_P,{\mathrm{\ψ}}_q)$

其中，d(ψ_P,ψ_q)表示目标块ψ_P和样本块ψ_q中对应已知像素的颜色差的平方和，可由下式表示：

d(ψ_P,ψ_q)＝Σ[(P_r-q_r)²+(P_g-q_g)²+(P_b-q_b)²]

其中P、q分别对应块ψ_P和ψ_q中的已知像素点。

(4)模块拷贝修复及优先权信息更新

在找到最佳匹配块ψ_Q后，将块ψ_Q中的像素对应复制到目标块ψ_P中的未知像素点，该目标块内未知像素点转变为已知像素点，因此这些点的置信度需要重新更新。

(5)重复以上过程，直至目标区域修复完毕。

具体包括以下步骤：

(a)针对图像反光区域特征，即反光区域中心像素点有较高的亮度(Value，V)值和较高的饱和度(Saturation，S)值，对图像逐像素点进行阈值分割以得到反光中心区域，分割公式如下：

$F(i,j)=\left\{\begin{array}{cc}255& V>thresh1,S>thresh2\\ 0& else\end{array}\right.$

其中的F(i，j)代表坐标为(i，j)的像素点所对应的颜色值，255对应的是白色，0对应的是黑色，thresh1和thresh2分别代表了某一固定的阈值，阈值大小视实际拍摄条件而定；

(b)对图像进行2次全方向膨胀，使得图像中相互靠近的反光中心区域相互合并，然后对膨胀后的图像进行2次全方向腐蚀，使得面积变大的反光中心区域还原为原始大小；

(c)对于面积较大(大于9个像素点)的反光中心区域，把这些区域的外接矩形作为图像的反光区域，而对于面积很小(小于9个像素点)的区域则进行简单的平滑处理；

(d)利用改进的样本块修复方法对反光区域进行图像复原，改进方法如下：

首先修补过程中每次只针对一个待修复区域进行修补，然后再复原其它修复块，直到全部修补完成为止；其次在搜索最佳匹配块的过程中采用二次匹配法进行搜索，第一次匹配是粗略匹配，通过隔行隔列粗略搜索最佳匹配块的中心点位置，每次扫描的数据量降为原图四分之一；第二次匹配是精确匹配，在第一次匹配搜索到的中心点位置的八邻域点内进行遍历搜索，以得到最后的最佳匹配块进行修复。

上述步骤(4)中结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割；其中迭代法是基于逼近的思想，其算法如下：

(1)求出图象的最大灰度值和最小灰度值，分别记为Z_max和Z_min，令初始阈值T₀＝(Z_max+Z_min)/2；

(2)根据阈值T_K将图象分割为前景和背景，分别求出两者的平均灰度值Z_O和Z_B；

(3)求出新阈值T_K+1＝(Z_O+Z_B)/2；

(4)若T_K＝T_K+1，则所得即为阈值；否则转(2)，迭代计算。

迭代所得的阈值分割的图像效果良好。基于迭代的阈值能区分出图像的前景和背景的主要区域所在；

其中最大熵阀值法的原理如下：

熵是图像统计特性的一种表现形式，反映了图像包含信息量的大小，分割后图像熵值越大，说明分割后从原图得到信息量越大，分割图像细节越丰富，因而总体分割效果也越好。

香农熵公式为：H(S)＝-P₁lnP₁-P₀lnP₀，其中P₁和P₀分别表示分割图像二值输出为1和0的概率，选取图像直方图的谷值中熵最大的值作为分割值。

具体包括以下步骤：

(a)去除图像背景：将去除反光后的图像转换为灰度图像，根据眼肌图像的特征，转换的灰度值公式为：Gray＝G*0.6+B*0.4，其中Gray为图像灰度值，G为原始图像绿色通道值，B为原始图像蓝色通道值，运用大津法对得到的灰度图像进行整体阈值分割去除图像背景，得到图像P1；

(b)提取眼肌中的脂肪部分：把P1均分为400(20*20)份并分别进行阈值分割，每一份小图像根据其图像特征为其选取不同的自适应阈值分割算法，选取的流程如下：(迭代法分割出的大理石花纹面积为S1，最大熵法分割出的大理石花纹面积为S2)

当图像块中包含背景区域时，选用迭代法对图像块进行自适应阈值分割；

当S2≤S1＜1.3*S2时，选用迭代法对图像块进行自适应阈值分割；

当S1≤S2＜1.3*S1时，选用最大熵法对图像块进行自适应阈值分割；

当S2≥1.3*S1时，选用迭代法对图像块进行自适应阈值分割；

当S1≥1.3*S2时，选用最大熵法对图像块进行自适应阈值分割；

每个小块分割完成后得到图像P2；

(c)去除肌间脂肪：对P2进行轮廓提取，选出周长最大的一个轮廓，对该轮廓进行填充，对填充后的区域进行形态学膨胀处理得到图像P3，对图像P2和P3进行差值运算，得到图像P4；

(d)对P4进行形态学膨胀处理得到图像P5，对P5进行轮廓提取，找出轮廓像素点个数小于10的轮廓区域，这些区域为噪声点，去除这些区域得到图像P6，对P5和P6进行与计算得到图像P7，此时的图像即为提取出的眼肌大理石花纹。

上述步骤(5)中从得到的大理石花纹图像中提取出脂肪数量指标特征值共8个，所述特征的详细表述如下：

F1有效眼肌面积：即为图像P7中的大理石花纹面积；

F2脂肪颗粒个数：夹杂于眼肌中的脂肪斑点，每一个联通部分为一个脂肪颗粒；

F3大颗粒脂肪密度：单位面积上拥有的大颗粒脂肪数量(个/cm²)；

F4中颗粒脂肪密度：单位面积上拥有的中颗粒脂肪数量(个/cm²)；

F5小颗粒脂肪密度：单位面积上拥有的小颗粒脂肪数量(个/cm²)；

F6大颗粒脂肪面积比：大颗粒脂肪总面积与有效眼肌面积的比例(％)；

F7中颗粒脂肪面积比：中颗粒脂肪总面积与有效眼肌面积的比例(％)；

F8小颗粒脂肪面积比：小颗粒脂肪总面积与有效眼肌面积的比例(％)；

其中脂肪颗粒按面积(A)划分成三类：小颗粒(L1)，A1≤3.72mm²(30个像素)；中颗粒(L2)，3.72＜A2≤14.88mm²(120个像素)；大颗粒(L3)，A3＞14.88mm²。

上述步骤(5)中从得到的大理石花纹图像中提取出脂肪分布指标特征值共3个，所述特征的详细表述如下：

F1脂肪颗粒均匀度：指的是脂肪颗粒分布变异系数，设有效眼肌图像有n行，每行的脂肪像素占该行总像素比例为w_i(i＝1,2,…,n)，平均值为，则脂肪颗粒均匀度C的计算公式如下：

$C=\frac{\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(w_i-\stackrel{\‾}{w})}^2}{n-1}}}{\stackrel{\‾}{w}}$

F2计盒维数：计盒维数的计算，直观的理解，就是数格子数。设有一个平面集合F，首先用一个矩形将其覆盖住，再以r为边长，将此矩形均匀地划分为若干小矩形(称为盒子)，然后开始计数，只要网格中含有F的任意一部分，均可计数，逐个地数出这样地网格数目N(r)。改变标度r，N(r)与1/r将满足下列关系：

$N\left(r\right)\∝{(1/r)}^{D_B}$

则F为分形集，D_B为其计盒维数。

F3信息维数：对于一幅经过二值化处理的牛肉大理石花纹图像，可用1或0来表示一个脂肪信息元素，一个覆盖中包含的脂肪信息量，可以通过计算一个覆盖中包含了多少个脂肪像素1或0来确定。根据信息维的定义，设N为脂肪信息元素总量，N_i为每个覆盖中所含脂肪信息元素的数量，则每个覆盖中脂肪信息分布概率为P_i＝N_i/N，脂肪信息I_i＝-P_ilnP_i，总脂肪信息量为：

$I\left(r\right)=\underset{i=1}{\overset{N\left(r\right)}{\Σ}}{P}_i{\mathrm{lnP}}_i$

其中N(r)为非空盒个数。

上述步骤(5)中从得到的大理石花纹图像中提取出脂肪纹理指标特征值，本发明用的是灰度共生矩阵特征值，灰度共生矩阵是对图像上保持某距离的两象素分别具有某灰度的状况进行统计得到的。

取图像(N×N)中任意一点(x，y)及偏离它的另一点(x+a，y+b)，设该点对的灰度值为(g1，g2)。令点(x，y)在整个画面上移动，则会得到各种(g1，g2)值，设灰度值的级数为k，则(g1，g2)的组合共有k的平方种。对于整个画面，统计出每一种(g1，g2)值出现的次数，然后排列成一个方阵，再用(g1，g2)出现的总次数将它们归一化为出现的概率P(g1，g2)，这样的方阵称为灰度共生矩阵。距离差分值(a，b)取不同的数值组合，可以得到不同情况下的联合概率矩阵。(a，b)取值要根据纹理周期分布的特性来选择，对于较细的纹理，选取(1，0)、(1，1)、(2，0)等小的差分值。

当a＝1，b＝0时，像素对是水平的，即0度扫描；当a＝0，b＝1时，像素对是垂直的，即90度扫描；当a＝1，b＝1时，像素对是右对角线的，即45度扫描；当a＝-1，b＝1时，像素对是左对角线，即135度扫描。

这样，两个象素灰度级同时发生的概率，就将(x，y)的空间坐标转化为“灰度对”(g1，g2)的描述，形成了灰度共生矩阵。

分别选取像素间距d＝1,2,3,4,5共五组距离、方向角度θ＝0°、45°、90°、135°共四个方向的14个特征参数，则特征值共有280个(5*4*14)，所述特征参数的详细表述如下：

设f(x，y)为一幅二维数字图像，灰度级别为n＝256，其灰度共生矩阵为p(i，j，d，θ)，其中i、j分别表示矩阵中的行与列，d表示像素间距，θ表示方向角度；

F1角二阶矩：

$W_1=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{p}^2(i,j,d,\mathrm{\θ})$

F2对比度：

$W_2=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\[{(i-j)}^2\×p^2(i,j,d,\mathrm{\θ})\]$

F3相关：

$W_3=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\[i\×j\×p^2(i,j,d,\mathrm{\θ})-u_1\×u_2\]/(d_1\×d_2)$

$\begin{array}{cc}{u}_1=\underset{i=1}{\overset{n}{\Σ}}i\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& u_2=\underset{j=1}{\overset{n}{\Σ}}j\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

$\begin{array}{cc}{d}_1^2=\underset{i=1}{\overset{n}{\Σ}}{(i-u_1)}^2\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& d_2^2=\underset{j=1}{\overset{n}{\Σ}}{(j-u_1)}^2\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

F4熵：

$W_4=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\×{\log }_{10}p(i,j,d,\mathrm{\θ})$

F5方差：

$W_5=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{(i-m)}^{2}p(i,j,d,\mathrm{\θ})$

其中，m为p(i，j，d，θ)的均值；

F6均值和：

$W_6=\underset{k=0}{\overset{2n}{\Σ}}k\×\underset{i=0}{\overset{n}{\Σ}}\underset{j=0}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})$

F7方差和：

$W_7=\underset{k=2}{\overset{2n}{\Σ}}{(k-W_6)}^2\×\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})$

F8逆差矩：

$W_8=\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})/\[1+{(i-j)}^2\]$

F9差的方差：

$W_9=\underset{k=0}{\overset{n-1}{\Σ}}\[k-{\underset{k=0}{\overset{n-1}{\Σ}}k\×\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\]}^2\×\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})$

F10和熵：

$W_{10}=-\underset{k=2}{\overset{2n}{\mathrm{\Σ}}}\[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]\×\log \[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]$

F11差熵：

$W_{11}=-\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}\[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]\×\log \[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}p(i,j,d,\mathrm{\θ})\]$

F12聚类阴影：

$W_{12}=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[(i-u_1)+(j-u_2)\]}^3\×p(i,j,d,\mathrm{\θ})$

$\begin{array}{cc}{u}_1=\underset{i=1}{\overset{n}{\Σ}}i\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& u_2=\underset{j=1}{\overset{n}{\Σ}}j\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

F13显著聚类：

$W_{13}=-\underset{i=1}{\overset{n}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[(i-u_1)+(j-u_2)\]}^4\×p(i,j,d,\mathrm{\θ})$

$\begin{array}{cc}{u}_1=\underset{i=1}{\overset{n}{\Σ}}i\underset{j=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})& u_2=\underset{j=1}{\overset{n}{\Σ}}j\underset{i=1}{\overset{n}{\Σ}}p(i,j,d,\mathrm{\θ})\end{array}$

F14最大概率：

$W_{14}=\underset{i,j}{MAX}\[p(i,j,d,\mathrm{\θ})\]$

上述步骤(6)中的逐步回归预测模型的具体实现步骤如下：

(a)根据猪肉眼肌实际肌内脂肪含量对猪肉眼肌样本进行打分，分为预测集和测试集；

(b)根据预测集样本的291个特征值建立逐步回归模型；

(c)将测试集样本代入建好的逐步回归预测模型；

(d)将逐步回归预测结果与实际肌内脂肪含量进行比较，得出预测准确率；

上述步骤(6)中的Fisher判别预测模型的具体实现步骤如下：

(a)根据猪肉眼肌实际肌内脂肪含量对猪肉眼肌样本进行打分，分为预测集和测试集；

(b)根据预测集样本的291个特征值建立Fisher判别预测模型；

(c)将测试集样本代入建好的Fisher判别预测模型；

(d)将Fisher判别预测结果与实际肌内脂肪含量进行比较，得出预测准确率。

Claims (8)

1.一种基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于包括以下步骤：

(1)标定摄像机，得到图像中单位像素对应于真实场景中的实际长度；

(2)选取猪肉第5-6根肋骨处背最长肌(俗称为眼肌)作为样本，将样本放置于黑色背景布上，并用白色LED灯对样本进行补光，利用数码照相机对样本横截面进行俯拍，获取猪肉眼肌样本图像，图像大小为1000*1000像素；

(3)利用改进的样本块修复方法对得到的样本图像进行预处理，恢复图像反光区域原有信息；

(4)结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割，提取出猪肉大理石花纹；

(5)从得到的大理石花纹图像中提取出脂肪数量指标、脂肪分布指标和脂肪纹理指标等特征值；

(6)根据特征值和化学方法检测结果建立逐步回归预测模型和Fisher判别模型，而后对模型进行检验和结果比较，最后通过该模型对猪肉肌内脂肪含量进行预测。

2.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述步骤(3)中利用样本块修复方法对得到的样本图像进行预处理包括以下步骤：

(2a)针对图像反光区域特征，即反光区域中心像素点有较高的亮度(Value，V)值和较高的饱和度(Saturation，S)值，对图像逐像素点进行阈值分割以得到反光中心区域，分割公式如下：

其中的F(i，j)代表坐标为(i，j)的像素点所对应的颜色值，255对应的是白色，0对应的是黑色，thresh1和thresh2分别代表了某一固定的阈值，阈值大小视实际拍摄条件而定；

(2b)对图像进行2次全方向膨胀，使得图像中相互靠近的反光中心区域相互合并，然后对膨胀后的图像进行2次全方向腐蚀，使得面积变大的反光中心区域还原为原始大小；

(2c)对于面积较大(大于9个像素点)的反光中心区域，把这些区域的外接矩形作为图像的反光区域，而对于面积很小(小于9个像素点)的区域则进行简单的平滑处理；

(2d)利用改进的样本块修复方法对反光区域进行图像复原，改进方法如下：

首先修补过程中每次只针对一个待修复区域进行修补，然后再复原其它修复块，直到全部修补完成为止；其次在搜索最佳匹配块的过程中采用二次匹配法进行搜索，第一次匹配是粗略匹配，通过隔行隔列粗略搜索最佳匹配块的中心点位置，每次扫描的数据量降为原图四分之一；第二次匹配是精确匹配，在第一次匹配搜索到的中心点位置的八邻域点内进行遍历搜索，以得到最后的最佳匹配块进行修复。

3.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述步骤(4)中结合最大熵法与迭代法对预处理后的眼肌图像进行图像分割包括以下步骤：

(3a)去除图像背景：将去除反光后的图像转换为灰度图像，根据眼肌图像的特征，转换的灰度值公式为：Gray＝G*0.6+B*0.4，其中Gray为图像灰度值，G为原始图像绿色通道值，B为原始图像蓝色通道值，运用大津法对得到的灰度图像进行整体阈值分割去除图像背景，得到图像P1；

(3b)提取眼肌中的脂肪部分：把P1均分为400(20*20)份并分别进行阈值分割，每一份小图像根据其图像特征为其选取不同的自适应阈值分割算法，选取的流程如下：(迭代法分割出的大理石花纹面积为S1，最大熵法分割出的大理石花纹面积为S2)

当图像块中包含背景区域时，选用迭代法对图像块进行自适应阈值分割；

当S2≤S1＜1.3*S2时，选用迭代法对图像块进行自适应阈值分割；

当S1≤S2＜1.3*S1时，选用最大熵法对图像块进行自适应阈值分割；

当S2≥1.3*S1时，选用迭代法对图像块进行自适应阈值分割；

当S1≥1.3*S2时，选用最大熵法对图像块进行自适应阈值分割；

每个小块分割完成后得到图像P2；

(3c)去除肌间脂肪：对P2进行轮廓提取，选出周长最大的一个轮廓，对该轮廓进行填充，对填充后的区域进行形态学膨胀处理得到图像P3，对图像P2和P3进行差值运算，得到图像P4；

(3d)对P4进行形态学膨胀处理得到图像P5，对P5进行轮廓提取，找出轮廓像素点个数小于10的轮廓区域，这些区域为噪声点，去除这些区域得到图像P6，对P5和P6进行与计算得到图像P7，此时的图像即为提取出的眼肌大理石花纹。

4.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述上述步骤(5)中从得到的大理石花纹图像中提取出脂肪数量指标特征值共8个，所述特征的详细表述如下：

F1有效眼肌面积：即为图像P7中的大理石花纹面积；

F2脂肪颗粒个数：夹杂于眼肌中的脂肪斑点，每一个联通部分为一个脂肪颗粒；

F3大颗粒脂肪密度：单位面积上拥有的大颗粒脂肪数量(个/cm²)；

F4中颗粒脂肪密度：单位面积上拥有的中颗粒脂肪数量(个/cm²)；

F5小颗粒脂肪密度：单位面积上拥有的小颗粒脂肪数量(个/cm²)；

F6大颗粒脂肪面积比：大颗粒脂肪总面积与有效眼肌面积的比例(％)；

F7中颗粒脂肪面积比：中颗粒脂肪总面积与有效眼肌面积的比例(％)；

F8小颗粒脂肪面积比：小颗粒脂肪总面积与有效眼肌面积的比例(％)；

其中脂肪颗粒按面积(A)划分成三类：小颗粒(L1)，A1≤3.72mm²(30个像素)；中颗粒(L2)，3.72＜A2≤14.88mm²(120个像素)；大颗粒(L3)，A3＞14.88mm²。

5.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述步骤(5)中从得到的大理石花纹图像中提取出脂肪分布指标特征值共3个，所述特征的详细表述如下：

F1脂肪颗粒均匀度：指的是脂肪颗粒分布变异系数，设有效眼肌图像有n行，每行的脂肪像素占该行总像素比例为w_i(i＝1,2,…,n)，平均值为w，则脂肪颗粒均匀度C的计算公式如下：

F2计盒维数；

F3信息维数。

6.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述步骤(5)中从得到的大理石花纹图像中提取出脂肪纹理指标特征值，本发明用的是灰度共生矩阵特征值，分别选取像素间距d＝1,2,3,4,5共五组距离、方向角度θ＝0°、45°、90°、135°共四个方向的14个特征参数，则特征值共有280个(5*4*14)，所述特征参数的详细表述如下：

设f(x，y)为一幅二维数字图像，灰度级别为n＝256，其灰度共生矩阵为p(i，j，d，θ)，其中i、j分别表示矩阵中的行与列，d表示像素间距，θ表示方向角度；

F1角二阶矩：

F2对比度：

F3相关：

F4熵：

F5方差：

其中，m为p(i，j，d，θ)的均值；

F6均值和：

F7方差和：

F8逆差矩：

F9差的方差：

F10和熵：

F11差熵：

F12聚类阴影：

F13显著聚类：

F14最大概率：

。

7.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述步骤(6)中的逐步回归预测模型的具体实现步骤如下：

(7a)根据猪肉眼肌实际肌内脂肪含量对猪肉眼肌样本进行打分，分为预测集和测试集；

(7b)根据预测集样本的291个特征值建立逐步回归模型；

(7c)将测试集样本代入建好的逐步回归预测模型；

(7d)将逐步回归预测结果与实际肌内脂肪含量进行比较，得出预测准确率。

8.根据权利要求1所述的基于计算机视觉的猪肉肌内脂肪含量无损检测方法，其特征在于所述步骤(6)中的Fisher判别预测模型的具体实现步骤如下：

(8a)根据猪肉眼肌实际肌内脂肪含量对猪肉眼肌样本进行打分，分为预测集和测试集；

(8b)根据预测集样本的291个特征值建立Fisher判别预测模型；

(8c)将测试集样本代入建好的Fisher判别预测模型；

(8d)将Fisher判别预测结果与实际肌内脂肪含量进行比较，得出预测准确率。