CN102799860A - 一种用于显微图像全息识别的方法 - Google Patents

Abstract

一种用于显微图像全息识别的方法，其特征是，包括以下过程：用RGB颜色模型表示彩色显微图像；建立IUV颜色模型过程；检测图像边缘过程；提取显微图像几何特征的过程；颜色识别过程；提取显微图像颜色特征的过程。本发明是以颜色模型理论基础，以识别微粒曲面和微粒全息特性为目标，研究微粒图像全息识别的方法，该全息识别方法是在三维流形中进行的，识别过程中各个颜色通道保持独立，没有产生降维，其识别结果是微粒的曲面片，包含微粒的颜色特征信息和几何特征信息，因此，本发明不仅解决了彩色图像的识别问题，提高了微粒识别精度和速度，为微粒图像识别提供了坚实的理论基础。

Description

一种用于显微图像全息识别的方法

技术领域

本发明涉及一种图像处理方法，具体地说是一种用于显微图像全息识别的方法，属于图像处理技术领域。

背景技术

目前，微粒显微图像识别技术存在的主要问题就是不能真正识别彩色微粒、识别精度低，因此，如何能突破现有识别方法，真正实现微粒显微图像的彩色识别、提高识别精度和速度，让微粒图像识别系统真正发挥其应有的作用，是图像识别研究者必须担负的责任和紧迫的任务。

彩色图像的边缘检测是显微图像几何特征的识别的第一步。

活动轮廓模型(Active contour models)是20世纪80年代后期发展起来的一种图像分割方法，早期的活动轮廓模型又称为Snake模型。在活动轮廓模型或蛇模型中最基本的理念是为了检测图像中的目标，曲线演化的约束条件受制于一幅给定的图像。例如，一个被曲线包围检测目标，曲线朝其内部移动，并最终停止在目标的边界^[27]。活动轮廓模型在一个给定的图像中检测对象，即基于曲线演化功能和水平集方法。我们的模型检测到的对象不一定是由梯度定义的边界。我们使能看作为最小的区分问题的特例的能量减到最小。在水平集公式中,这个问题就会成为一个“平均曲率流”-主动轮廓演变的，它将停在预期的边界。虽然经典的活动轮廓模型速度的停止不依赖于图像的梯度，但实际上与特定的图像分割有关，近几年来活动轮廓法及水平集方法由于本身的优异性得到了快速的发展。

设闭合曲线C将定义域为Ω的图像u₀(x,y)划分为两部分区域：曲线内部区域Ω₁和曲线外部区域Ω₂，活动轮廓模型的能量泛函定义如下式(1.1)所示：

E^CV(c₁,c₂,C)＝λ∫_Ω1|u₀(x,y)-c₁|²dxdy+λ₂∫_Ω2|u₀(x,y)-c₂|²dxdy+μL_Length(C)(1.1)

其中，λ₁＞0,λ₂＞0,μ≥0为各项的权重系数，c₁代表闭合曲线C内部区域的灰度均值，c₂分别是闭合曲线C外部区域的灰度均值，L_Length(C)表示闭合曲线的长度。

在水平集方法中可以用C＝{(x,y)∈Ω:φ(x,y)＝0}表示闭合曲线C，其中φ为轮廓内为正轮廓外为负的符号距离函数。Chan和Vese根据欧拉—拉格朗日方程和梯度下降流法推导出水平集函数的演化方程如式(1.2)所示：

$\left\{\begin{array}{c}{c}_1\left(\mathrm{\φ}\right)=\frac{{\∫}_{\mathrm{\Ω}}{u}_0(x,y)H_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x.y)\right)\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}{H}_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}},c_2\left(\mathrm{\φ}\right)=\frac{{\∫}_{\mathrm{\Ω}}{u}_0(x,y)(1-H_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}(1-H_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}\\ \frac{\∂\mathrm{\φ}}{\∂t}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)[\mathrm{\μ}\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)-{\mathrm{\λ}}_1{(u_0-c_1)}^2+{\mathrm{\λ}}_2{(u_0-c_2)}^2]\end{array}\right.---\left(1.2\right)$

其中，H(φ)和δ(φ)是Heaviside函数和Dirac函数的规则，表示为式(1.3)所示：

$H\left(\mathrm{\φ}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\π}}\arctan \left(\frac{\mathrm{\φ}}{\mathrm{\ϵ}}\right))$ $\mathrm{\δ}\left(\mathrm{\φ}\right)=\frac{1}{\mathrm{\π}}\frac{\mathrm{\ϵ}}{{\mathrm{\ϵ}}^2+{\mathrm{\φ}}^2}---\left(1.3\right)$

另外，水平集的演化需要一定的稳定性，所以每次更新后需要进行重新初始化水平集函数为符号距离函数。重新初始化通常采用如式(1.4)迭代方程：

$\left\{\begin{array}{c}{\mathrm{\ψ}}_t=\mathrm{sign}\left(\mathrm{\φ}\left(t\right)\right)(1-|\▿\mathrm{\ψ}\left|\right)\\ \mathrm{\ψ}(0,\·)=\mathrm{\φ}(t,\·)\end{array}\right.---\left(1.4\right)$

但是，由于每次演化都要重新初始化过程，所以大大增加了该模型的计算量，因此具有成本高、效率底和耗时间的特点，从而限制了模型的实际应用。

发明内容

为克服上述现有技术存在的不足，本发明所要解决的技术问题在于提供了一种用于显微图像全息识别的方法，其不仅解决了彩色图像的识别问题，而且提高了微粒识别精度和速度，为微粒图像识别提供了坚实的理论基础。

本发明解决其技术问题所采取的技术方案是：一种用于显微图像全息识别的方法，其特征是，包括以下过程：

用RGB颜色模型表示彩色显微图像；

建立IUV颜色模型过程；

检测图像边缘过程；

提取显微图像几何特征的过程；

颜色识别过程；

提取显微图像颜色特征的过程。

进一步地，所述建立IUV颜色模型过程为将彩色显微图像从RGB颜色模型转换为IUV颜色模型的过程。

进一步地，所述IUV颜色模型与RGB颜色模型的转换转换关系为：

[I U V 1]＝[R G B 1]T

其中， $T=\left[\begin{array}{cccc}\sqrt{3}/3& -\sqrt{2}/2& -\sqrt{6}/6& 0\\ \sqrt{3}/3& \sqrt{2}/2& -\sqrt{6}/6& 0\\ \sqrt{3}/3& 0& \sqrt{6}/3& 0\\ 0& 0& 0& 1\end{array}\right].$

进一步地，所述IUV颜色模型的IUV分量的分布为：

$\left\{\begin{array}{c}I=\sqrt{3}/3R+\sqrt{3}/3G+\sqrt{3}/3B\\ U=-\sqrt{2}/2R+\sqrt{2}/2G\\ V=-\sqrt{6}/6R-\sqrt{6}/6G+\sqrt{6}/3\end{array}\right.$

其中，IUV分别为显微图像的使用色调、饱和度和透明度。

进一步地，所述检测图像边缘过程为提取IUV颜色模型的I分量并使用改进活动轮廓法完成微粒边缘检测的过程。

进一步地，所述使用改进活动轮廓法完成微粒边缘检测的过程包括以下步骤：

第1步用水平集函数初始化，所述水平集函数为

$\mathrm{\φ}(x,t=0)=\left\{\begin{array}{c}-\mathrm{\ρ}(x\∈{\mathrm{\Ω}}_1)\\ 0(x\∈\∂{\mathrm{\Ω}}_1)\\ \mathrm{\ρ}(x\∈{\mathrm{\Ω}}_2)\end{array}\right.$

其中，ρ＞0是常数，轮廓线的内部区域为Ω₁，轮廓线的外部区域为Ω₂，

为区域的边界；

第2步根据迭代方程式更新C₁和C₂，所述迭代方程式为

$\left\{\begin{array}{c}{\mathrm{\ψ}}_t=\mathrm{sign}\left(\mathrm{\φ}\left(t\right)\right)(1-|\▿\mathrm{\ψ}\left|\right)\\ \mathrm{\ψ}(0,\·)=\mathrm{\φ}(t,\·)\end{array}\right.;$

第3步根据

通过离散化式来演化水平集函数，所述离散化式为

$\frac{{\mathrm{\φ}}_{i,j}^{n+1}-{\mathrm{\φ}}_{i,j}^n}{\mathrm{\Δt}}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left({\mathrm{\φ}}_{i,j}^n\right)[\mathrm{\μ}{K}_{i,j}^n-{\mathrm{\λ}}_1{({\mathrm{\μ}}_{0,i,j}-c_1\left({\mathrm{\φ}}^n\right))}^2+{\mathrm{\λ}}_2{({\mathrm{\μ}}_{0,i,j}-c_2\left({\mathrm{\φ}}^n\right))}^2]+$

$\mathrm{\μ}({\mathrm{\φ}}_{i+1,j}^n+{\mathrm{\φ}}_{i-1,j}^n+{\mathrm{\φ}}_{i,j+1}^n+{\mathrm{\φ}}_{i,j-1}^n-4{\mathrm{\φ}}_{i,j}^n-K_{i,j}^n);$

第4步判断水平集函数的演化是否达到稳定状态，如果没有，则转向第2步，进行下一次迭代，否则算法结束。

进一步地，所述几何特征包括：表征颗粒大小的参数包括颗粒周长、颗粒面积、等效面积直径、周长等效直径和等效粒径，以及表征颗粒的形状的参数有球型度、长宽比、庞大率和表面率。

进一步地，所述根颜色识别过程为根据由IUV颜色模型I分量创建的曲面以及UV分量的分布，使用最接近点法完成颜色特征识别的过程。

进一步地，所述最接近点法为曲面上的很多微分算子可以在高维的空间中由相应的笛卡尔微分算子来替代，用最接近点cp(x)来计算数量。

进一步地，所述颜色特征包括：全局颜色特征和空间颜色特征，所述全局颜色特征包括颜色直方图、累计直方图、颜色矩和颜色熵；所述颜色空间特征包括颜色聚合向量和颜色相关图。

本发明的有益效果如下：本发明是以颜色模型理论基础，以识别微粒曲面和微粒全息特性为目标，研究微粒图像全息识别的方法，该全息识别方法是在三维流形中进行的，识别过程中各个颜色通道保持独立，没有产生降维，其识别结果是微粒的曲面片，包含微粒的颜色特征信息和几何特征信息，因此，本发明不仅解决了彩色图像的识别问题，提高了微粒识别精度和速度，为微粒图像识别提供了坚实的理论基础。另外本发明也为图像识别技术在颗粒检测和细胞检测中的应用提供了坚实的基础。

附图说明

下面结合附图对本发明进一步说明：

图1是本发明的原理框架图；

图2（a）是本发明所述RGB颜色空间模型的示意图；

图2（b）是本发明所述IUV颜色空间模型的示意图。

具体实施方式

如图1所示，本发明的一种用于显微图像全息识别的方法，其过程为：彩色显微图像根据IUV颜色模型，完成RGB色空间到IUV色彩空间的转换，提取I分量使用改进活动轮廓法完成微粒的边缘检测，实现微粒的几何特征提取；创建I曲面，使用最接近点法识别颜色特征，实现颜色特征的提取，几何特征提取和颜色特征的提取实现显微图像的全息识别。

其中，构建颜色模型是关键的第一步，目前，显微图像无论采取那种分割方法都与颜色的坐标选择和颜色的特征有关，但是所使用某一种分割方法不一定对所有的颜色特征都有效，同样也不可能有单一的颜色坐标适用于所有的分割方法。我们若想获得最佳的图像边缘就需要采用适用这种分割方法的颜色空间。为了更好的研究显微图像颜色特征的提取，构造IUV颜色模型，表示的使用色调、饱和度和透明度的一种颜色空间，称为IUV(intensity，saturation,verification)颜色空间为例来说明。

下面对本发明的以下几个主要部分进行详细说明。

一、建立IUV颜色模型

图2（a）是本发明所述RGB颜色空间模型的示意图，图2（b）是本发明所述IUV颜色空间模型的示意图。如图2（a）所示，取立方体点O（0,0,0）和点B（1,1,1）的对角线矢量坐标（1,0,1）设为点A，

与XOZ平面的夹角为α，OA与X轴夹角为β，β＝45°。

$\mathrm{tg\α}=\frac{\mathrm{AB}}{\mathrm{OA}}=\frac{1}{\sqrt{2}}---\left(1.5\right)$

由式(1.2)解得α＝35.2644°。

在RGB色彩空间里进行两步复合变换，第一步向量绕Z轴旋转α角，使向量

与向量重合，第二步向量

绕Y轴旋转β角，使向量

与X轴重合。设向量

绕Z轴旋转角为θ，则变换矩阵为T_Z，如式(1.6)所示：

$T_Z=\left[\begin{array}{cccc}\cos \mathrm{\θ}& \sin \mathrm{\θ}& 0& 0\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]---\left(1.6\right)$

其中θ＝-α时则有，如式(1.7)所示：

$T_Z=\left[\begin{array}{cccc}\cos \mathrm{\α}& -\sin \mathrm{\α}& 0& 0\\ \sin \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}\sqrt{2}/2& -\sqrt{2}/2& 0& 0\\ \sqrt{2}/2& \sqrt{2}/2& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]---\left(1.7\right)$

设向量

绕Y轴旋转角为β，则变换矩阵为Ty，如式(1.8)所示：

$\mathrm{Ty}=\left[\begin{array}{cccc}\cos \mathrm{\θ}& 0& -\sin \mathrm{\θ}& 0\\ 0& 1& 0& 0\\ \sin \mathrm{\θ}& 0& \cos \mathrm{\θ}& 0\\ 0& 0& 0& 1\end{array}\right]---\left(1.8\right)$

其中θ＝β则有，如式(1.9)所示：

$\mathrm{Ty}=\left[\begin{array}{cccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}& 0\\ 0& 1& 0& 0\\ \sin \mathrm{\β}& 0& \cos \mathrm{\β}& 0\\ 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccc}\sqrt{6}/3& 0& -\sqrt{3}/3& 0\\ 0& 1& 0& 0\\ \sqrt{3}/3& 0& \sqrt{6}/3& 0\\ 0& 0& 0& 1\end{array}\right]---\left(1.9\right)$

变换矩阵T_Z，Ty进行组合变换，组合变换矩阵T，如式(1.10)所示：

$T{=T}_Z\·T_Y\left[\begin{array}{cccc}\sqrt{3}/3& -\sqrt{2}/2& -\sqrt{6}/6& 0\\ \sqrt{3}/3& \sqrt{2}/2& -\sqrt{6}/6& 0\\ \sqrt{3}/3& 0& \sqrt{6}/3& 0\\ 0& 0& 0& 1\end{array}\right]---\left(1.10\right)$

从组合变换矩阵可得出式(1.11)：

$\left\{\begin{array}{c}I=\sqrt{3}/3R+\sqrt{3}/3G+\sqrt{3}/3B\\ U=-\sqrt{2}/2R+\sqrt{2}/2G\\ V=-\sqrt{6}/6R-\sqrt{6}/6G+\sqrt{6}/3\end{array}\right.---\left(1.11\right)$

IUV色彩空间和RGB色彩空间只是同一物理量的不同表示法，因而它们之间存在着转换关系。RGB色彩空间转和IUV色彩空间转公式为下式(1.12)所示：

[I U V 1]＝[R G B 1]T                           (1.12)

当R=G=B=1时

U＝V＝0,此时I的大小和向量的模相等，可以得知此颜色模型是正确的可用，IUV颜色空间模型如图2（b）所示。

二、改进活动轮廓法

为了避免初始化和减少计算量的目标，对原有的活动轮廓进行改进，改进模型的能量泛函的水平集形式定义如下式(1.13)所示：

$E(c_1,c_2,\mathrm{\φ})={\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{|{\mathrm{\μ}}_0(x,y)-c_1|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|{\mathrm{\μ}}_0(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}$

$\left(1.13\right)$

$+\mathrm{\μ}{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}+\frac{1}{2}{\∫}_{\mathrm{\Ω}}{\left(\right|\▿\mathrm{\φ}(x,y)|-1)}^2\mathrm{dxdy}$

其中，H(φ)和δ(φ)是Heaviside函数和Dirac函数的规则，为了便于计算，将采用这两个函数的近似表示式(1.3)。设φ不变，根据最小化能量泛函式(1.13)，可以得到c₁和c₂的表达式如下式(1.14)所示：

$\left\{\begin{array}{c}{c}_1\left(\mathrm{\φ}\right)=\frac{{\∫}_{\mathrm{\Ω}}{\mathrm{\μ}}_0(x,y)H_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}{H}_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y0)\right)\mathrm{dxdy}}\\ c_2\left(\mathrm{\φ}\right)=\frac{{\∫}_{\mathrm{\Ω}}{\mathrm{\μ}}_0(x,y)(1-H_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}(1-H_{\mathrm{\ϵ}}\left(\mathrm{\φ}(x,y0)\right))\mathrm{dxdy}}\end{array}\right.---\left(1.14\right)$

设c₁和c₂不变，通过欧拉—拉格朗日方程和梯度下降流法，使改进模型的能量泛函式(1.14)最小化，可以得到水平集函数φ的演化方程如式(1.15)所示：

$\frac{\∂\mathrm{\φ}}{\∂t}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)[\mathrm{\μ}\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)-{\mathrm{\λ}}_1{({\mathrm{\μ}}_0-c_1)}^2+{\mathrm{\λ}}_2{({\mathrm{\μ}}_0-c_2)}^2]+\mathrm{\α}({\▿}^2\mathrm{\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))---\left(1.15\right)$

水平集函数中这一项

是为了避免反复初始化水平集函数，加快演化的速度。本改进算法采用有限差分法离散化式(1.15)，设Δt为时间步长，h是空间步长，则

可以离散化为如下式(1.16)所示：

$K_{i,j}^n=\mathrm{div}(\▿\mathrm{\φ}/|\▿\mathrm{\φ}\left|\right)=D_{i,j}^{0x}\left(\frac{D_{i,j}^{0x}{\mathrm{\φ}}_{i,j}^n}{\sqrt{{\left(D_{i,j}^{0x}{\mathrm{\φ}}_{i,j}^n\right)}^2+{\left(D_{i,j}^{0y}{\mathrm{\φ}}_{i,j}^n\right)}^2}}\right)+D_{i,j}^{0y}\left(\frac{D_{i,j}^{0y}{\mathrm{\φ}}_{i,j}^n}{\sqrt{{\left(D_{i,j}^{0x}{\mathrm{\φ}}_{i,j}^n\right)}^2+{\left(D_{i,j}^{0y}{\mathrm{\φ}}_{i,j}^n\right)}^2}}\right)---\left(1.16\right)$

其中，

和

是中心差分算子，计算如下式(1.17)所示：

$D_{i,j}^{0x}=\frac{{{\mathrm{\φ}}_{i+1,j}}^n-{{\mathrm{\φ}}_{i-1,j}}^n}{2h},$ $D_{i,j}^{0y}=\frac{{{\mathrm{\φ}}_{i,j+1}}^n-{{\mathrm{\φ}}_{i,j-1}}^n}{2h}---\left(1.17\right)$

则水平集函数φ的演化方程(1.15)可以离散为式(1.18)所示：

$\frac{{\mathrm{\φ}}_{i,j}^{n+1}-{\mathrm{\φ}}_{i,j}^n}{\mathrm{\Δt}}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left({\mathrm{\φ}}_{i,j}^n\right)[\mathrm{\μ}{K}_{i,j}^n-{\mathrm{\λ}}_1{({\mathrm{\μ}}_{0,i,j}-c_1\left({\mathrm{\φ}}^n\right))}^2+{\mathrm{\λ}}_2{({\mathrm{\μ}}_{0,i,j}-c_2\left({\mathrm{\φ}}^n\right))}^2]+---\left(1.18\right)$

$\mathrm{\μ}({\mathrm{\φ}}_{i+1,j}^n+{\mathrm{\φ}}_{i-1,j}^n+{\mathrm{\φ}}_{i,j+1}^n+{\mathrm{\φ}}_{i,j-1}^n-4{\mathrm{\φ}}_{i,j}^n-K_{i,j}^n)$

使用改进的活动模型检测图像边缘步骤如下:

第1步用公式(1.19)

$\mathrm{\φ}(x,t=0)=\left\{\begin{array}{c}-\mathrm{\ρ}(x\∈{\mathrm{\Ω}}_1)\\ 0(x\∈\∂{\mathrm{\Ω}}_1)\\ \mathrm{\ρ}(x\∈{\mathrm{\Ω}}_2)\end{array}\right.---\left(1.19\right)$

水平集函数初始化，其中ρ＞0是常数，轮廓线的内部区域为Ω₁，轮廓线的外部区域为Ω₂，

为区域的边界;

第2步根据迭代方程式(1.4)更新c₁和c₂;

第3步根据

可以离散化式(1.18)演化水平集函数;

第4步判断水平集函数的演化是否达到稳定状态，如果没有，则转向第2步，进行下一次迭代，否则算法结束。

三、几个特征提取

1、表征颗粒大小的参数包括颗粒周长、颗粒面积、等效面积直径、周长等效直径和等效粒径。

1)颗粒周长

根据前面介绍的改进的活动轮廓法提取颗粒的边缘，边缘是由一系列的象素点构成，象素点用P_i(i=1，2…n)表示，周长用链码表示，各象素点之间按照八连通方式排列，则微粒的周长L_P可表示为式(1.20)所示：

$L_p=N_E+\sqrt{2}{N}_o---\left(1.20\right)$

其中，N_E和N_o分别是边界链码(8方向)中走偶数步与走奇数步的数目。

2)颗粒面积

颗粒的面积就是颗粒所占的象素点数目，即颗粒的边界及边界内包含的象素点数。设颗粒μ₀(x,y)的大小为M×N，设目标颗粒μ₀(x,y)=1，背景为0，则目标颗粒的面积如下式(1.21)所示：

$A=\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{\mathrm{\μ}}_0(x,y)---\left(1.21\right)$

若有图像中n个颗粒，则上式表示所有目标颗粒的总面积。

3)面积等效直径如下式(1.22)所示：

$\mathrm{ds}=\sqrt{\frac{4A}{\mathrm{\π}}}---\left(1.22\right)$

其中，A表示颗粒截面积。

4)周长等效直径如下式(1.23)所示：

$\mathrm{dc}=\frac{L_p}{\mathrm{\π}}---\left(1.23\right)$

其中，L_p表示颗粒截面周长。

5)等效粒径

在实际的测量中，等效粒径一般是通过某种的等效方法来测量颗粒的粒径。实际生产过程中，颗粒的形状有些是不规则的，粒径只能采取统计的等效方式，这样统计出来的结果对于实际生产才有意义。对于球形的颗粒来说，粒径就是指颗粒的直径，等效粒径的表达式如下式(1.24)所示:

$d=\frac{4A}{L_P}---\left(1.24\right)$

其中，A代表颗粒的面积，L_p代表颗粒的周长。

2.表征颗粒的形状的参数有球型度、长宽比、庞大率和表面率。

1)球型度

球型度表征颗粒接近球形的程度,用Q表示，如公式(1.25)所示，可以看出球体的球型度等于1，其他颗粒球形度小于1，

2)长宽比

长宽比表征颗粒接近柱形的程度，用L表示，如公式(1.26)所示，一般颗粒长宽比恒大于1，

3)庞大率

庞大率表征颗粒填充最小外接矩形空间的程度，用P表示，如公式(1.27)所示，一般颗粒庞大率恒小于1，

4)表面率

表面率表征颗粒表面积的相对大小，用R表示，如公式(1.28)所示，球体的表面率为1，非球体的表面率大于1，

显微颗粒经过处理之后，选择“静态分析”功能菜单下的“每个颗粒数据”，然后点击一个颗粒，计算并显示所选中的单个颗粒的几何尺寸：面积等效直径、周长等效直径、X-切线、Y-切线、球形度和长宽比。

四、颜色识别

本发明使用最接近点方法解决曲面分割问题。最接近点方法是采用最接近的点来代表曲面，因此并不需要参数化设计，也不需要闭合曲面，甚至不需要知道曲面方向。该方法的一个显着特点是，它在一个小型的三维相邻曲面上完全将标准方法与插值点法相结合。这在实践中特别有吸引力，因为它能帮助我们利用现有的有限元差分软件将3维区域分割成2维曲面形状段，所有的分割都使用无边活动轮廓。曲面上的水平集方程梯度下降法。这些方程由最接近点方法演变而来，是一种用于解决曲面上的偏微分方程(PDE)的最新技术。这种算法有一个特别简单的形式：它只是在三维邻近曲面用一个时间步长来代替普通C-V模型中的一个插值点。该方法可以用于处理一般的曲面，因为它使用一个最接近点函数来表示基础表面。

设给定的曲面图像μ₀：M→R被水平集函数φ(x):M→R将曲面μ₀(x)分为两部分，水平集函数轮廓线内部为正，外部为负。使用能量最小化来进行分割，能量方程如下式(1.29)所示：

E(μ₁,μ₂,φ)＝μ·φ(x)+λ₁∫₊|μ₀(x)-μ₁|²dM+λ₂∫|-μ₀(x)-μ₂|²dM    (1.29)

式中μ₁和μ₂为未知量，且μ≥0,λ₁，λ₂>0是权重参数，整个曲面引入H(φ)和δ(φ)是Heaviside函数和Dirac函数的规则来将能量视为一个整体。我们使用μ₁和μ₂来保持φ(x)的稳定和能量函数E(μ₁,μ₂,φ)的最小化，μ₁(φ)和μ₂(φ)的表达式如下式(1.30)所示:

$\left\{\begin{array}{c}{\mathrm{\μ}}_1\left(\mathrm{\φ}\right)=\frac{\∫{\mathrm{\μ}}_0\left(x\right)H\left(\mathrm{\φ}\left(x\right)\right)\mathrm{dxdy}}{\∫H\left(\mathrm{\φ}\left(x\right)\right)\mathrm{dxdy}}\\ {\mathrm{\μ}}_2\left(\mathrm{\φ}\right)=\frac{\∫{\mathrm{\μ}}_0\left(x\right)(1-H\left(\mathrm{\φ}\left(x\right)\right))\mathrm{dxdy}}{\∫(1-H\left(\mathrm{\φ}\left(x\right)\right))\mathrm{dxdy}}\end{array}\right.---\left(1.30\right)$

也就是说，μ₁和μ₂分别是φ(x)对μ₀(x)轮廓的内部和外部的灰度的平均值。

式中H(φ)和δ(φ)是Heaviside函数和Dirac函数的规则近似为公式(1.31)所示，为了保持μ₁和μ₂来保持φ(x)的稳定和能量函数E(μ₁,μ₂,φ)的最小化，通过欧拉-拉格朗日方程和梯度下降流法，推导出水平集函数的演化方程式(1.25)所示：

$\frac{\∂\mathrm{\φ}}{\∂t}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)[\mathrm{\μ}{\▿}_M\·\left(\frac{{\▿}_M\left(\mathrm{\φ}\right)}{\left|{\▿}_M\mathrm{\φ}\right|}\right)]-{\mathrm{\λ}}_1{({\mathrm{\μ}}_0\left(x\right)-{\mathrm{\μ}}_1)}^2+{\mathrm{\λ}}_2{({\mathrm{\μ}}_0\left(x\right)-{\mathrm{\μ}}_2)}^2---\left(1.31\right)$

其中t≥0直到水平集演化到稳定状态。

最接近点法主要分为嵌入偏微分方程和离散化两部分。

(1)嵌入式偏微分方程

由于曲面导数的存在，曲面偏微分方程的计算较为复杂。最接近点法在三维空间中应用标准的笛卡尔网格法，这是一种计算偏微分方程和其他曲面处理过程的新方法。这种方法由曲面上距离最近的点表示。在已知的曲面M上，最接近点函数采用点x∈R³，并返回距离x最近的一个或多个点，由cp(x)∈M表示，可以用公式(1.32)所示：

cp(x)＝argmin_q∈M||x-q||₂                             (1.32)

最接近点法的主要思想：曲面上的很多微分算子可以在高维的空间中由相应的笛卡尔微分算子来替代，用最接近点cp(x)来计算数量。例如，在光滑曲面M上的点x，

但由于φcp(x)在正常方向保持不变，所以函数只能沿曲面变化。与分割问题相关的一个类似特性是适用于非线性的扩散条件。曲面曲率项满足下式(1.33)所示：

$\▿M\·\left(\frac{\▿\mathrm{M\φ}}{|\▿\mathrm{M\φ}|}\right)=\▿\·\left(\frac{\▿\mathrm{\φ}\left(\mathrm{cp}\right)}{|\▿\mathrm{\φ}\left(\mathrm{cp}\right)|}\right)---\left(1.33\right)$

其中，x是光滑曲面M上的所有点。再将这个问题带入三维空间，用上式替换曲面的曲率，我们获得了嵌入式的偏微分方程如下式(1.34)所示：

$\frac{\∂\mathrm{\φ}}{\∂t}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\left(\mathrm{cp}\right)\right)[\mathrm{\μ}\▿\·\left(\frac{\▿\mathrm{\φ}\left(\mathrm{cp}\right)}{|\▿\mathrm{\φ}\left(\mathrm{cp}\right)|}\right)-{\mathrm{\λ}}_1{({\mathrm{\μ}}_0\left(\mathrm{cp}\right)-{\mathrm{\μ}}_1)}^2+{\mathrm{\λ}}_2{({\mathrm{\μ}}_0\left(\mathrm{cp}\right)-{\mathrm{\μ}}_2)}^2]---\left(1.34\right)$

这是使用在三维空间分割的偏微分方程，用函数计算出最接近点cp(x)，函数式μ₁和μ₂也可以由嵌入空间中cp(x)计算出来。

梯度下降方程(1.33)中出现的诺伊曼齐次边界条件说明，在一个开曲面中不需要明确指出接近点，可以看出出现这样的边界条件是扩张过程的一部分。此外，不需要在计算域的边缘时引入人工介入的边界条件。

(2)离散化

为了离散化偏微分方程，首先,我们选择一个区域Ω_C中的离散网格点:这是没有必要去解决对整个领域,但只是包含曲面中的某一窄带。

对于区域Ω_C中的所有离散点而言，我们需要在曲面M上构造一个最接近点函数cp(x)，我们用它来代表曲面。对于简单的曲面，例如一个球体，圆环体或立方体，用最接近点函数来分析非常简单。为了解决嵌入式偏微分方程(1.34)，我们将以下两个步骤交替使用：

第一步是将计算域扩展至整个曲面，再使用最近点函数进行求解。然后，将计算域上的每个节点用φ(cp)代替φ。将最接近点扩展为插值点，并在三维空间Ω_C上用质心拉格朗日插值法将之展开。

在笛卡尔网格的一个时间步长上使用标准的有限元差分法计算嵌入式偏微分方程。这是C-V算法在三维空间区域分割时常用的一步。

第二步是差分法离散算法，我们可知下式(1.35)：

$\frac{{\mathrm{\φ}}_{\mathrm{ijk}}^{n+1}-{\mathrm{\φ}}_{\mathrm{ijk}}^n}{\mathrm{\Δt}}={\mathrm{\δ}}_h\left({\mathrm{\φ}}_{\mathrm{ijk}}^n\right)[\frac{\mathrm{\μ}}{h^2}{\mathrm{\Δ}}^x-\frac{{\mathrm{\Δ}}^x+{\mathrm{\φ}}_{\mathrm{ijk}}^{n+1}}{\sqrt{\frac{{({\mathrm{\Δ}}^x+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{h^2}+}\frac{{({\mathrm{\Δ}}_c^y+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{h^2}+\frac{{({\mathrm{\Δ}}_c^z+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{h^2}}$

$+\frac{\mathrm{\μ}}{h^2}{\mathrm{\Δ}}^y-\frac{{\mathrm{\Δ}}^y+{\mathrm{\φ}}_{\mathrm{ijk}}^{n+1}}{\sqrt{\frac{{({\mathrm{\Δ}}^y+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{h^2}+}\frac{{({\mathrm{\Δ}}_c^x+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{4h^2}+\frac{{({\mathrm{\Δ}}_c^z+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{4h^2}}---\left(1.35\right)$

$+\frac{\mathrm{\μ}}{h^2}{\mathrm{\Δ}}^z-\frac{\mathrm{\Δz}+{\mathrm{\φ}}_{\mathrm{ijk}}^{n+1}}{\sqrt{\frac{{({\mathrm{\Δ}}^z+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{h^2}+}\frac{{({\mathrm{\Δ}}_c^x+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{4h^2}+\frac{{({\mathrm{\Δ}}_c^y+{\mathrm{\φ}}_{\mathrm{ijk}}^n)}^2}{4h^2}}$

$-{\mathrm{\λ}}_1{({\mathrm{\μ}}_{0,\mathrm{ijk}}-{\mathrm{\μ}}_1\left({\mathrm{\φ}}^n\right))}^2+{\mathrm{\λ}}_2{({\mathrm{\μ}}_{0,\mathrm{ijk}}-{\mathrm{\μ}}_2\left({\mathrm{\φ}}^n\right))}^2$

此时，

代表方向，中心差分法和向后差分法，由上标α表示方向，h表示网格间距。在整个区域Ω_C的网格单元上用简单正交法计算μ₁和μ₂的平均值。

线性系统(1.35)使用高斯——赛德尔迭代方法来求出。因为我们只对稳定状态下式(1.34)的解决方案有兴趣，所以我们不必要保证每个时间段迭代的收敛性，只需要求能量下降的最小值。我们还要注意到，重新初始化有时可以避免水平集函数变得过于平坦或者过于陡峭^[37]。最终目的是为了梯度下降方程(1.28)的稳定状态。

在本文中，将C–V分割算法与最接近点法相结合，提出了一种定义在一般的曲面,甚至包括开曲面的分割方式。这种方法简单明了，在一个邻近的小型三维曲面中只需要在标准C-V算法计算时用时间步长代替插值。

五、颜色特征提取

1.颜色的量化

颜色量化是指将指定的颜色空间映射到一个给定大小的子集中，使其总体误差最小。颜色量化方法有均匀量化、非均匀量化和静态量化。常用的量化方法是非均匀量化法。本课题使用IUV颜色模型，下面以IUV为例进行非均匀量化：

1)将I空间分成3份，U空间分成8份，V空间分成3份。

2)根据色彩的不同范围和主观颜色感知进行量化，如下式(1.36)所示：

$I=\left\{\begin{array}{cc}0& I\∈\left[\mathrm{0,0.3}\right]\\ 1& I\∈\left[\mathrm{0.3,0.7}\right]\\ 2& I\∈\left[\mathrm{0.7,1}\right]\end{array}\right.$ $V=\left\{\begin{array}{cc}0& V\∈\left[\mathrm{0,0.3}\right]\\ 1& V\∈\left[\mathrm{0.3,0.7}\right]\\ 2& V\∈\left[\mathrm{0.7,1}\right]\end{array}\right.$

$U=\left\{\begin{array}{cccc}0& U\∈\left[\mathrm{316,20}\right]& 4& U\∈\left[\mathrm{156,190}\right]\\ 1& U\∈\left[\mathrm{21,40}\right]& 5& U\∈\left[\mathrm{191,270}\right]\\ 2& U\∈\left[\mathrm{41,75}\right]& 6& U\∈\left[\mathrm{271,295}\right]\\ 3& U\∈\left[\mathrm{76,155}\right]& 7& U\∈\left[\mathrm{276,315}\right]\end{array}\right.---\left(1.36\right)$

3)构造一维特征矢量。按照以上的量化级，把各颜色分量合成为一维特征矢量如下式(1.37)所示：

W＝UQ_IQ_V+IQ_V+V                              (1.37)

其中，O₁和O_V分别是I、V的量化级数，O₁=O_V=3。公式实际上为如下式(1.38)所示：

W＝9U+3I                                    (1.38)

这样，I,U,V三个分量在一维矢量上分布开来。

2.颜色特征的分类及表示

在图像颜色特征的提取中，主要包括全局颜色特征和空间颜色特征两个方面。全局颜色特征包括颜色直方图、累计直方图、颜色矩和颜色熵。颜色空间特征包括颜色聚合向量和颜色相关图。

1)全局颜色特征

(1)颜色直方图

在颜色特征的表达方法中最常用的颜色直方图。对于一幅彩色图像，已确定颜色空间的基础上，统计每一种颜色出现的象素数，把颜色作为横坐标，某种颜色出现的象素点总数作为纵坐标，以此统计画出的图形就是颜色直方图。图像的归一化颜色直方图C可以定义为如下式(1.39)所示：

$C=\left\{\left(c\right[n_1],c[n_2],....,c[n_k\left]\right)\right|\underset{k=1}{\overset{m}{\mathrm{\Σ}}}c\left[n_k\right]=\mathrm{1,0}\≤c\left[n_k\right]\≤1\}---\left(1.39\right)$

$c\left[n_k\right]=\frac{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}\underset{j-1}{\overset{W-1}{\mathrm{\Σ}}}\left\{\begin{array}{cc}1& I[i,j]=c_k\\ 0& \end{array}\right.}{L\×W}$

其中，L和W代表图像的长和宽。

(2)累积直方图

累积直方图是以颜色值作为横坐标，以颜色累积出现的频数为纵坐标。因此累积直方图的定义如下式(1.40)所示：

$\left(k\right)=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{n}_i}{N}$ k=0,1……L-1                                  (1.40)

其中，k是图像的颜色特征取值，L是颜色特征可取值的个数，n_i是图像中颜色特征值为k的像素的个数，N是像素的总数。在累积直方图中，相邻颜色在频数上是相关的。

(3)颜色矩

颜色矩是一种简单有效的表示方法，图像中任何的颜色分布均可以用它的矩表示，无需对特征进行向量化，降低了颜色特征的维数,分为一矩表示均值，二矩表示方差，三矩表示斜度，公式如下式(1.41)所示：

${\mathrm{\μ}}_i=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{p}_{\mathrm{ij}},$ ${\mathrm{\σ}}_1=\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(p_{\mathrm{ij}}-{\mathrm{\μ}}_i)}^2},$ ${\mathrm{\ξ}}_1=\sqrt[3]{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(p_{\mathrm{ij}}-{\mathrm{\μ}}_i)}^3}---\left(1.41\right)$

其中，p_ij表示象素j的颜色值i的概率，N为图像中象素点的个数。

(4)颜色熵

John提出采用图像颜色的信息熵来表示图像的颜色特征，从而将图像的颜色直方图由多维降低到一维。假设一幅图像中有N中颜色值，第i个颜色值在图像中出现的概率记为P_i，用概率密度函数表示不同颜色值的象素出现的概率，颜色熵可以表达式如下式(1.42)所示：

$C=-\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_i{\log }_2{P}_i---\left(1.42\right)$

2)空间颜色特征

(1)颜色聚合矢量

空间颜色特征用颜色聚合矢量表示，其核心思想就是把每个区间的像素分为两部分，这些区间属于颜色直方图中。这样统计图像所包含的每种颜色的聚合像素和非聚合像素的比率称为该图像的颜色聚合矢量^[57]，它在某种程度上保留了图像颜色的空间信息。

如果定义α为第i个聚合像素，β为第j个非聚合像素，那么颜色聚合矢量定义如下式(1.43)所示：

<(α₁,β₁),(α₂,β₂),....(α_N,β_N)>                     (1.43)

<(α₁+β₁),(α₂+β₂),....(α_N+β_N)>为图像的颜色累计直方图。

(2)颜色相关图

颜色空间分布的另一种表达方式是利用图像中各象素间的颜色关系来描述，这种特征不仅可以表示某一种颜色的象素的百分比，还可以反映不同颜色值之间的空间相关性。

本发明是以颜色模型理论基础，以识别微粒曲面和微粒全息特性为目标，研究微粒图像全息识别的方法，该全息识别方法是在三维流形中进行的，识别过程中各个颜色通道保持独立，没有产生降维，其识别结果是微粒的曲面片，包含微粒的颜色特征信息和几何特征信息，因此，本发明不仅解决了彩色图像的识别问题，提高了微粒识别精度和速度，为微粒图像识别提供了坚实的理论基础。

Claims (10)

1.一种用于显微图像全息识别的方法，其特征是，包括以下过程：

用RGB颜色模型表示彩色显微图像；

建立IUV颜色模型过程；

检测图像边缘过程；

提取显微图像几何特征的过程；

颜色识别过程；

提取显微图像颜色特征的过程。

2.根据权利要求1所述的一种用于显微图像全息识别的方法，其特征是，所述建立IUV颜色模型过程为将彩色显微图像从RGB颜色模型转换为IUV颜色模型的过程。

3.根据权利要求2所述的一种用于显微图像全息识别的方法，其特征是，所述IUV颜色模型与RGB颜色模型的转换转换关系为：

[I U V 1]＝[R G B 1]T

其中， $T=\left[\begin{array}{cccc}\sqrt{3}/3& -\sqrt{2}/2& -\sqrt{6}/6& 0\\ \sqrt{3}/3& \sqrt{2}/2& -\sqrt{6}/6& 0\\ \sqrt{3}/3& 0& \sqrt{6}/3& 0\\ 0& 0& 0& 1\end{array}\right].$

4.根据权利要求2所述的一种用于显微图像全息识别的方法，其特征是，所述IUV颜色模型的IUV分量的分布为：

$\left\{\begin{array}{c}I=\sqrt{3}/3R+\sqrt{3}/3G+\sqrt{3}/3B\\ U=-\sqrt{2}/2R+\sqrt{2}/2G\\ V=-\sqrt{6}/6R-\sqrt{6}/6G+\sqrt{6}/3\end{array}\right.$

其中，IUV分别为显微图像的使用色调、饱和度和透明度。

5.根据权利要求1所述的一种用于显微图像全息识别的方法，其特征是，所述检测图像边缘过程为提取IUV颜色模型的I分量并使用改进活动轮廓法完成微粒边缘检测的过程。

6.根据权利要求5所述的一种用于显微图像全息识别的方法，其特征是，所述使用改进活动轮廓法完成微粒边缘检测的过程包括以下步骤：

第1步用水平集函数初始化，所述水平集函数为

$\mathrm{\&phi;}(x,t=0)=\left\{\begin{array}{c}-\mathrm{\&rho;}(x\&Element;{\mathrm{\&Omega;}}_1)\\ 0(x\&Element;\&PartialD;{\mathrm{\&Omega;}}_1)\\ \mathrm{\&rho;}(x\&Element;{\mathrm{\&Omega;}}_2)\end{array}\right.$

其中，ρ＞O是常数，轮廓线的内部区域为Ω₁，轮廓线的外部区域为Q₂，

为区域的边界；

第2步根据迭代方程式更新c₁和c₂，所述迭代方程式为

$\left\{\begin{array}{c}{\mathrm{\&psi;}}_t=\mathrm{sign}\left(\mathrm{\&phi;}\left(t\right)\right)(1-|\&dtri;\mathrm{\&psi;}\left|\right)\\ \mathrm{\&psi;}(0,\&CenterDot;)=\mathrm{\&phi;}(t,\&CenterDot;)\end{array}\right.;$

第3步根据

通过离散化式来演化水平集函数，所述离散化式为

$\frac{{\mathrm{\&phi;}}_{i,j}^{n+1}-{\mathrm{\&phi;}}_{i,j}^n}{\mathrm{\&Delta;t}}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left({\mathrm{\&phi;}}_{i,j}^n\right)[\mathrm{\&mu;}{K}_{i,j}^n-{\mathrm{\&lambda;}}_1{({\mathrm{\&mu;}}_{0,i,j}-c_1\left({\mathrm{\&phi;}}^n\right))}^2+{\mathrm{\&lambda;}}_2{({\mathrm{\&mu;}}_{0,i,j}-c_2\left({\mathrm{\&phi;}}^n\right))}^2]+$

$\mathrm{\&mu;}({\mathrm{\&phi;}}_{i+1,j}^n+{\mathrm{\&phi;}}_{i-1,j}^n+{\mathrm{\&phi;}}_{i,j+1}^n+{\mathrm{\&phi;}}_{i,j-1}^n-4{\mathrm{\&phi;}}_{i,j}^n-K_{i,j}^n);$

第4步判断水平集函数的演化是否达到稳定状态，如果没有，则转向第2步，进行下一次迭代，否则算法结束。

7.根据权利要求1所述的一种用于显微图像全息识别的方法，其特征是，所述几何特征包括：表征颗粒大小的参数包括颗粒周长、颗粒面积、等效面积直径、周长等效直径和等效粒径，以及表征颗粒的形状的参数有球型度、长宽比、庞大率和表面率。

8.根据权利要求1所述的一种用于显微图像全息识别的方法，其特征是，所述根颜色识别过程为根据由IUV颜色模型I分量创建的曲面以及UV分量的分布，使用最接近点法完成颜色特征识别的过程。

9.根据权利要求8所述的一种用于显微图像全息识别的方法，其特征是，所述最接近点法为曲面上的很多微分算子可以在高维的空间中由相应的笛卡尔微分算子来替代，用最接近点cp(x)来计算数量。

10.根据权利要求1所述的一种用于显微图像全息识别的方法，其特征是，所述颜色特征包括：全局颜色特征和空间颜色特征，所述全局颜色特征包括颜色直方图、累计直方图、颜色矩和颜色熵；所述颜色空间特征包括颜色聚合向量和颜色相关图。

Images