CN101777122B - 一种角毛藻显微图像细胞目标提取方法 - Google Patents

Abstract

本发明提供一种角毛藻显微图像细胞目标提取方法，包括以下步骤：a.将角毛藻原始图像转换为灰度图像；b.构建灰度图像方向角模型，分别构成两幅灰度图像；c.对两幅灰度图像分别采用增强大津法自动阈值二值化后再进行逻辑与运算；d.对所得到的二值图像做双向滤波和中值滤波以去掉二值图像中小的噪声点；e.进行形态学膨胀操作以连接目标细胞断裂的边缘；f.提取二值图像的最大轮廓并填充后再进行形态学腐蚀操作；g.用腐蚀后的二值图像与原始图像做逻辑与运算，得到细胞目标图像。本发明方法能够从带有泥沙、残骸等干扰物的角毛藻显微图像中较为精确的提取目标细胞，提高效率和识别准确率。

Description

一种角毛藻显微图像细胞目标提取方法

技术领域

本发明涉及浮游植物细胞目标提取方法技术领域，特别涉及一种角毛藻显微图像细胞目标提取方法。

背景技术

角毛藻属是硅藻门中心纲盒形藻目角毛藻科下的属，是硅藻门下较大的属之一，其属下包含了大量种类的中国海常见浮游植物，且部分种类属于中国海常见有害赤潮藻种。

角毛藻细胞在壳面顶轴两端各生一细长的角毛，其断面大多为圆形。亦有四角形或多角形等。角毛的表面常有横裂点纹，及螺旋排列的实心小刺。少数种细胞单独生活，大部分种相邻细胞的相近角毛交接，联成直的、弯的或螺旋状的链状群体。链两端细胞的端角毛生出的位置不同，有的种链两端角毛伸出的位置也不同，它们向链外伸出的方向不同。

一直以来，生物外形是生物学家对生物进行分类的最重要，也是最基本的手段，而对角毛藻而言，角毛的形态、生长位置、伸出方向，色素体的形态、数量及其在细胞内的分布等都是角毛藻鉴定的重要依据，因此，从角毛藻显微图像中精确地提取目标细胞(尤其是角毛特征)就成为藻类鉴定和研究工作中重要的一步，而且还可以大幅度降低对角毛藻自动分类识别工作的技术要求，提高效率和识别准确率。

发明内容

本发明所要解决的技术问题在于，提供一种角毛藻显微图像细胞目标提取方法，以此为角毛藻的自动化鉴定提供一种手段，为后续的特征提取、分类和计数等工作打好基础，可以比较精确的提取角毛藻的角毛特征，为藻类研究者提供有效地指导和帮助。

为解决上述技术问题，本发明提供了一种角毛藻显微图像细胞目标提取方法，所述方法包含以下步骤：

a.将角毛藻原始图像转换为灰度图像；

b.构建灰度图像方向角模型，计算灰度曲面的法线方向矢量与三维坐标系中三个坐标轴X、Y、Z方向的矢量方向角，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像；

c.对两幅灰度图像分别采用增强大津法自动阈值二值化后再进行逻辑与运算；

d.对所得到的二值图像做双向滤波和中值滤波以去掉二值图像中小的噪声点；

e.进行形态学膨胀操作以连接目标细胞断裂的边缘；

f.提取二值图像的最大轮廓并填充后再进行形态学腐蚀操作；

g.用腐蚀后的二值图像与原始图像做逻辑与运算，得到细胞目标图像。

在步骤b中，灰度图像方向角模型优选是将图像中任意4个像素点的坐标(x，y)及灰度值f(x，y)在三维坐标系中确定的4个点(x，y，f(x，y))构成灰度曲面，计算灰度曲面的法线方向矢量与三维坐标系中三个坐标轴X、Y、Z方向的矢量方向角，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像。

在步骤c中，增强大津法优选是在运用大津法之前做一个线性拉伸以解决图像主体部分与背景灰度差较小的问题；选取使类间方差最大和类内方差最小的图像灰度作为最佳阈值。

在步骤d中，对步骤c所得到的二值图像优选先进行双向滤波再进行中值滤波以去掉二值图像中小的噪声点。

在步骤f中，最大轮廓优选是指寻找到的轮廓包含的像素数最多，填充可以是将寻找到的最大轮廓内部填充为白色。

在步骤g中，优选将腐蚀后的图像与原始图像做逻辑与运算，提取细胞图像的边缘并保留内部纹理等细胞特征。

所述步骤b优选为：

对于一幅图像I，由任意4个像素点的坐标及灰度值在三维坐标系中确定的4个点构成灰度曲面，建立空间直角坐标系，令相邻两像素点间的距离为1，并将每个像素点的灰度值除以255，利用空间点坐标的概念确定每个像素点灰度在该坐标系下的坐标；

设4个相邻像素(i，j)，(i+1，j)，(i，j+1)，(i+1，j+1)在Z轴上的灰度值分别为I(i，j)，I(i+1，j)，I(i，j+1)，I(i+1，j+1)，在空间坐标系中A′B′C′、A′C′D′、A′B′D′和B′C′D′为四个互不重叠的三角形平面；

在A′C′D′平面中，D′C′对应的矢量为f_D′C′＝[0，1，I(i，j+1)-I(i，j)]，A′D′对应的矢量为f_A′D′＝[1，0，I(i，j)-I(i+1，j)]，由矢量积的概念，该三角形面的法线方向就是此两矢量的矢量积，f_A′C′D′＝f_A′D′×f_D′C′；

同样求出其它三角形平面的法线方向f_B′C′D′、f_A′B′D′和f_A′B′C′，则灰度曲面的法线方向能够用四个三角平面的法线方向的均值近似得到，为：

f_{A′B′C′D′}＝(f_B′C′D′+f_A′B′C′+f_A′B′D′+f_A′C′D′)/4

这就是灰度曲面A′B′C′D′的法线方向矢量；

再利用矢量方向角及方向余弦的概念，得到该法线矢量与X、Y、Z三个坐标轴的夹角分别如下：

θ_x(i，j)＝360°arccos(f_{A′B′C′D′}(1)/abs(f_{A′B′C′D′}))/2π

θ_y(i，j)＝360°arccos(f_{A′B′C′D′}(2)/abs(f_{A′B′C′D′}))/2π

θ_z(i，j)＝360°arccos(f_{A′B′C′D′}(3)/abs(f_{A′B′C′D′}))/2π

即为图像灰度分别在X、Y、Z三个方向上的矢量方向角；

对于每个方向上得到的夹角值分别进行灰度映射处理：

${\mathrm{map}}_x(i,j)=\left[\frac{255\×({\mathrm{\θ}}_x(x,y)-\min ({\mathrm{\θ}}_x(x,y))}{\max \left({\mathrm{\θ}}_x(x,y)\right)-\min \left({\mathrm{\θ}}_x(x,y)\right)}\right]$

${\mathrm{map}}_y(i,j)=\left[\frac{255\×({\mathrm{\θ}}_y(x,y)-\min ({\mathrm{\θ}}_y(x,y))}{\max \left({\mathrm{\θ}}_y(x,y)\right)-\min \left({\mathrm{\θ}}_y(x,y)\right)}\right]$

${\mathrm{map}}_z(i,j)=\left[\frac{255\×({\mathrm{\θ}}_z(x,y)-\min ({\mathrm{\θ}}_z(x,y))}{\max \left({\mathrm{\θ}}_z(x,y)\right)-\min \left({\mathrm{\θ}}_z(x,y)\right)}\right]$

最后，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像：

${\mathrm{map}}_{\mathrm{xz}}(i,j)=\sqrt{{\mathrm{map}}_x^2+{\mathrm{map}}_z^2}$

${\mathrm{map}}_{\mathrm{yz}}(i,j)=\sqrt{{\mathrm{map}}_y^2+{\mathrm{map}}_z^2}.$

所述步骤c优选为：

令原图像的灰度级为X，期望处理后的图像灰度级为Y，原始图像和期望图像的灰度级的分布范围极值分别为X_max，X_min和Y_max，Y_min，期望变换前后的图像对比度保持线性关系，即满足下式：

$\frac{Y-Y_{\min }}{Y_{\max }-Y_{\min }}=\frac{X-X_{\min }}{X_{\max }-X_{\min }}$

经整理线性拉伸的简单数学表达式为：Y＝aX+b

其中 $a=\frac{Y_{\max }-Y_{\min }}{X_{\max }-X_{\min }},$ $b=\frac{{X_{\max }{Y}_{\min }-Y}_{\max }-X_{\min }}{X_{\max }-X_{\min }}$

当a＞1时，Y_max-Y_min＞X_max-X_min，则变换后图像的对比度增大，视觉感应为亮处更亮，暗处更暗；而通过设置b值，使整幅图像灰度级的均值改变相应的b值；本发明中令Y_max＝255，Y_min＝0；

原始灰度图像的灰度级为L灰度级为i的像素点数为n_i，全部像素为N，归一化直方图：

阈值t将灰度级划分为两类：C₀＝(0，1，...，t)和C₁＝(t+1，t+2，...L-1)

这两类的概率分别为：

${\mathrm{\ω}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{P}_i=\mathrm{\ω}\left(t\right),$ ${\mathrm{\ω}}_1=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}{P}_i=1-\mathrm{\ω}\left(t\right)$

均值为：

${\mathrm{\μ}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}i\frac{P_i}{{\mathrm{\ω}}_0}=\frac{\mathrm{\μ}\left(t\right)}{\mathrm{\ω}\left(t\right)},$ ${\mathrm{\μ}}_1=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}i\frac{P_i}{{\mathrm{\ω}}_1}=\frac{{\mathrm{\μ}}_T\left(t\right)-\mathrm{\μ}\left(t\right)}{1-\mathrm{\ω}\left(t\right)}$

上式中 $\mathrm{\μ}\left(t\right)=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}i{P}_i,$ ${\mathrm{\μ}}_T\left(t\right)=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}i{P}_i$

方差：

${\mathrm{\σ}}_0^2=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_0)}^2\frac{P_i}{{\mathrm{\ω}}_0},$ ${\mathrm{\σ}}_1^2=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_1)}^2\frac{P_i}{{\mathrm{\ω}}_1}$

类内方差为： ${\mathrm{\σ}}_{\mathrm{\ω}}^2={\mathrm{\ω}}_0{\mathrm{\σ}}_0^2+{\mathrm{\ω}}_1{\mathrm{\σ}}_1^2$

类间方差为：

${\mathrm{\σ}}_B^2={\mathrm{\ω}}_0{({\mathrm{\μ}}_0-{\mathrm{\μ}}_T)}^2+{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2={\mathrm{\ω}}_0{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_0)}^2$

总体方差为： ${\mathrm{\σ}}_T^2={\mathrm{\σ}}_B^2+{\mathrm{\σ}}_{\mathrm{\ω}}^2$

改变t的值，使类间方差取得最大值时的t值即为最佳阈值；

用求得的最佳阈值对灰度图像进行二值化。

本发明角毛藻显微图像细胞目标提取方法能够从带有泥沙、残骸等干扰物的角毛藻显微图像中较为精确的提取目标细胞，从而可以大幅度降低对角毛藻自动分类识别工作的技术要求，提高效率和识别准确率；而且也可以比较精确的提取角毛藻的角毛特征，为藻类研究者提供有效地指导和帮助。

附图说明

图1是本发明角毛藻显微图像细胞目标提取方法的流程图；

图2为灰度图像方向角模型示意图；

图3为具体实施例中选取的原始图像；

图4为具体实施例中经过转换后的灰度图像；

图5为具体实施例中水平方向角XZ矢量距灰度图像；

图6为具体实施例中垂直方向角YZ矢量距灰度图像；

图7为图5经过二值化后的二值图像；

图8为图6经过二值化后的二值图像；

图9为图7和图8逻辑与运算后的二值图像；

图10为具体实施例中对二值图像做双向滤波和中值滤波后的结果；

图11为具体实施例中对提取到的最大轮廓并填充后再进行形态学腐蚀的图像；

图12为具体实施例中腐蚀后的二值图像与原始图像做逻辑与运算后得到的细胞目标图像。

具体实施方式

本发明角毛藻显微图像细胞目标提取方法，所述方法利用计算机图像处理技术，首先对角毛藻显微图像进行灰度变换，然后构建灰度图像方向角模型，得到原始图像中的目标细胞部分。

下面结合附图和具体实施例对本发明作进一步的说明：

如图1所示，为本发明无角毛类赤潮藻显微图像细胞目标提取方法，所述方法包含以下步骤：(1)将角毛藻原始图像转换为灰度图像；(2)构建灰度图像方向角模型，计算灰度曲面的法线方向矢量与三维坐标系中三个坐标轴X、Y、Z方向的矢量方向角，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像；(3)对两幅灰度图像分别采用增强大津法自动阈值二值化后再进行逻辑与运算；(4)对所得到的二值图像做双向滤波和中值滤波以去掉二值图像中小的噪声点；(5)进行形态学膨胀操作以连接目标细胞断裂的边缘；(6)提取二值图像的最大轮廓并填充后再进行形态学腐蚀操作；(7)用腐蚀后的二值图像与原始图像做逻辑与运算，得到细胞目标图像。

在得到原始图像后，进入步骤(1)，将原始图像转换为灰度图像，因为微生物在颜色上差别不大，同时灰度图像可以大大减少图像处理的运算量并且更有利于对图像进行特征的提取等。图3为原始图像，图4为经过转换后的灰度图像。

然后进入步骤(2)，构建灰度图像方向角模型，计算灰度曲面的法线方向矢量与三维坐标系中三个坐标轴X、Y、Z方向的矢量方向角，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像。具体的原理如下所述：

对于一幅图像I，由任意4个像素点的坐标及灰度值在三维坐标系中确定的4个点构成灰度曲面，如图2建立空间直角坐标系，令相邻两像素点间的距离为1，并将每个像素点的灰度值除以255，利用空间点坐标的概念确定每个像素点灰度在该坐标系下的坐标。

设4个相邻像素(i，j)，(i+1，j)，(i，j+1)，(i+1，j+1)在Z轴上的坐标(灰度值)分别为I(i，j)，I(i+1，j)，I(i，j+1)，I(i+1，j+1)，因此在空间坐标系中A′B′C′、A′C′D′、A′B′D′和B′C′D′四个互不重叠的三角形平面。

以A′C′D′平面为例，D′C′对应的矢量为f_D′C′＝[0，1，I(i，j+1)-I(i，j)]，A′D′对应的矢量为f_A′D′＝[1，0，I(i，j)-I(i+1，j)]，由矢量积的概念，该三角形面的法线方向就是此两矢量的矢量积，f_A′C′D′＝f_A′D′×f_D′C′。

同理，可求出其它三角形平面的法线方向f_B′C′D′、f_A′B′D′和f_A′B′C′，则灰度曲面的法线方向可以用四个三角平面的法线方向的均值近似得到，即为

f_{A′B′C′D′}＝(f_B′C′D′+f_A′B′C′+f_A′B′D′+f_A′C′D′)/4

这也就是灰度曲面A′B′C′D′的法线方向矢量。

再利用矢量方向角及方向余弦的概念，得到该法线矢量与X、Y、Z三个坐标轴的夹角分别如下：

θ_x(i，j)＝360°arccos(f_{A′B′C′D′}(1)/abs(f_{A′B′C′D′}))/2π

θ_y(i，j)＝360°arccos(f_{A′B′C′D′}(2)/abs(f_{A′B′C′D′}))/2π

θ_z(i，j)＝360°arccos(f_{A′B′C′D′}(3)/abs(f_{A′B′C′D′}))/2π

即为图像灰度分别在X、Y、Z三个方向上的矢量方向角。对于每个方向上得到的夹角值分别进行灰度映射处理：

${\mathrm{map}}_x(i,j)=\left[\frac{255\×({\mathrm{\θ}}_x(x,y)-\min ({\mathrm{\θ}}_x(x,y))}{\max \left({\mathrm{\θ}}_x(x,y)\right)-\min \left({\mathrm{\θ}}_x(x,y)\right)}\right]$

${\mathrm{map}}_y(i,j)=\left[\frac{255\×({\mathrm{\θ}}_y(x,y)-\min ({\mathrm{\θ}}_y(x,y))}{\max \left({\mathrm{\θ}}_y(x,y)\right)-\min \left({\mathrm{\θ}}_y(x,y)\right)}\right]$

${\mathrm{map}}_z(i,j)=\left[\frac{255\×({\mathrm{\θ}}_z(x,y)-\min ({\mathrm{\θ}}_z(x,y))}{\max \left({\mathrm{\θ}}_z(x,y)\right)-\min \left({\mathrm{\θ}}_z(x,y)\right)}\right]$

最后，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像：

${\mathrm{map}}_{\mathrm{xz}}(i,j)=\sqrt{{\mathrm{map}}_x^2+{\mathrm{map}}_z^2}$

${\mathrm{map}}_{\mathrm{yz}}(i,j)=\sqrt{{\mathrm{map}}_y^2+{\mathrm{map}}_z^2}$

图5为水平方向角XZ矢量距灰度图像，图6为垂直方向角YZ矢量距灰度图像。

接着进入步骤(3)，对两幅灰度图像分别采用增强大津法自动阈值二值化后再进行逻辑与运算。增强大津法是在运用大津法之前做一个线性拉伸以解决图像主体部分与背景灰度差较小的问题，具体的原理如下所述：

若原图像的灰度级为X，期望处理后的图像灰度级为Y，原始图像和期望图像的灰度级的分布范围极值分别为X_max，X_min和Y_max，Y_min，我们期望变换前后的图像对比度保持线性关系，即满足下式：

$\frac{Y-Y_{\min }}{Y_{\max }-Y_{\min }}=\frac{X-X_{\min }}{X_{\max }-X_{\min }}$

经整理线性拉伸的简单数学表达式为：Y＝aX+b

其中 $a=\frac{Y_{\max }-Y_{\min }}{X_{\max }-X_{\min }},$ $b=\frac{{X_{\max }{Y}_{\min }-Y}_{\max }-X_{\min }}{X_{\max }-X_{\min }}$

当a＞1时，Y_max-Y_min＞X_max-X_min，则变换后图像的对比度增大，视觉感应为亮处更亮，暗处更暗；而通过设置b值，可以使整幅图像灰度级的均值改变相应的b值。本发明中令Y_max＝255，Y_min＝0。

大津法依据的原理是利用类别方差作为判据，选取使类间方差最大和类内方差最小的图像灰度值作为最佳阈值。大津法可作如下理解：因方差是灰度分布均匀性的一种度量，方差值越大，说明构成图像的两部分差别越大，当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小，因此使类间方差最大的分割意味着错分概率最小。

原始灰度图像的灰度级为L灰度级为i的像素点数为n_i，全部像素为N，归一化直方图：

阈值t将灰度级划分为两类：C₀＝(0，1，...，t)和C₁＝(t+1，t+2，...L-1)

这两类的概率分别为：

${\mathrm{\ω}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{P}_i=\mathrm{\ω}\left(t\right),$ ${\mathrm{\ω}}_1=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}{P}_i=1-\mathrm{\ω}\left(t\right)$

均值为：

${\mathrm{\μ}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}i\frac{P_i}{{\mathrm{\ω}}_0}=\frac{\mathrm{\μ}\left(t\right)}{\mathrm{\ω}\left(t\right)},$ ${\mathrm{\μ}}_1=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}i\frac{P_i}{{\mathrm{\ω}}_1}=\frac{{\mathrm{\μ}}_T\left(t\right)-\mathrm{\μ}\left(t\right)}{1-\mathrm{\ω}\left(t\right)}$

上式中 $\mathrm{\μ}\left(t\right)=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}i{P}_i,$ ${\mathrm{\μ}}_T\left(t\right)=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}i{P}_i$

方差：

${\mathrm{\σ}}_0^2=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_0)}^2\frac{P_i}{{\mathrm{\ω}}_0},$ ${\mathrm{\σ}}_1^2=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_1)}^2\frac{P_i}{{\mathrm{\ω}}_1}$

类内方差为： ${\mathrm{\σ}}_{\mathrm{\ω}}^2={\mathrm{\ω}}_0{\mathrm{\σ}}_0^2+{\mathrm{\ω}}_1{\mathrm{\σ}}_1^2$

类间方差为：

${\mathrm{\σ}}_B^2={\mathrm{\ω}}_0{({\mathrm{\μ}}_0-{\mathrm{\μ}}_T)}^2+{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2={\mathrm{\ω}}_0{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_0)}^2$

总体方差为： ${\mathrm{\σ}}_T^2={\mathrm{\σ}}_B^2+{\mathrm{\σ}}_{\mathrm{\ω}}^2$

改变t的值，使类间方差取得最大值时的t值即为最佳阈值。用求得的最佳阈值对灰度图像进行二值化。

图7、图8是利用步骤(3)中的增强大津法方法选取的阈值对灰度图像进行二值化的结果。

图9为图7和图8逻辑与运算后的二值图像。

在得到二值图像之后，进入步骤(4)，对所得到的二值图像分别采用3×3的掩模做双向滤波和中值滤波以去掉二值图像中小的噪声点。

图10为对二值图像做双向滤波和中值滤波后的结果。

之后进入步骤(5)，进行形态学膨胀以连接目标细胞断裂的边缘。

然后进入步骤(6)，提取二值图像的最大轮廓并填充后再进行形态学腐蚀操作。本步骤中说的最大轮廓是指所包含的像素点数最多的轮廓，在找到最大轮廓后将最大轮廓内部填充为白色。图11为对提取到的最大轮廓并填充后再进行形态学腐蚀的图像。

最后进入步骤(7)，将步骤(6)得到的二值图像与原始图像做逻辑与运算，提取到角毛藻显微图像细胞的边缘并保留了内部纹理等细胞特征，得到细胞目标图像。

图12为最终得到的细胞目标图像。

所有上述为这一知识产权的首要实施系统，并没有设定限制以其他形式实施这种新系统。本领域技术人员将利用这一重要信息，对上述内容修改，以实现类似的执行情况。但是，所有基于本发明的修改或改造新方法，属于保留的权利。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims (5)

1.一种角毛藻显微图像细胞目标提取方法，其特征在于，包括以下步骤：

a.将角毛藻原始图像转换为灰度图像；

b.对于一幅图像I，由任意4个像素点的坐标及灰度值在三维坐标系中确定的4个点构成灰度曲面，建立空间直角坐标系，令相邻两像素点间的距离为1，并将每个像素点的灰度值除以255，利用空间点坐标的概念确定每个像素点灰度在该坐标系下的坐标；

设4个相邻像素(i，j)，(i+1，j)，(i，j+1)，(i+1，j+1)在Z轴上的灰度值分别为I(i，j)，I(i+1，j)，I(i，j+1)，I(i+1，j+1)，在空间坐标系中A′B′C′、A′C′D′、A′B′D′和B′C′D′为四个互不重叠的三角形平面；

在A′C′D′平面中，D′C′对应的矢量为f_D′C′＝[0，1，I(i，j+1)-I(i，j)]，A′D′对应的矢量为f_A′D′＝[1，0，I(i，j)-I(i+1，j)]，由矢量积的概念，该三角形面的法线方向就是此两矢量的矢量积，f_A′C′D′＝f_A′D′×f_D′C′；

同样求出其它三角形平面的法线方向f_B′C′D′、f_A′B′D′和f_A′B′C′，则灰度曲面的法线方向能够用四个三角平面的法线方向的均值近似得到，为：

f_{A′B′C′D′}＝(f_B′C′D′+f_A′B′C′+f_A′B′D′+f_A′C′D′)/4

这就是灰度曲面A′B′C′D′的法线方向矢量；

再利用矢量方向角及方向余弦的概念，得到该法线矢量与X、Y、Z三个坐标轴的夹角分别如下：

θ_x(i，j)＝360°arccos(f_{A′B′C ′D′}(1)/abs(f_{A′B′C′D′}))/2π

θ_y(i，j)＝360°arccos(f_{A′B′C ′D′}(2)/abs(f_{A′B′C′D′}))/2π

θ_z(i，j)＝360°arccos(f_{A′B′C ′D′}(3)/abs(f_{A′B′C′D′}))/2π

即为图像灰度分别在X、Y、Z三个方向上的矢量方向角；

对于每个方向上得到的夹角值分别进行灰度映射处理：

${\mathrm{map}}_x(i,j)=\left[\frac{255\×({\mathrm{\θ}}_x(x,y)-\min ({\mathrm{\θ}}_x(x,y))}{\max \left({\mathrm{\θ}}_x(x,y)\right)-\min \left({\mathrm{\θ}}_x(x,y)\right)}\right]$

${\mathrm{map}}_y(i,j)=\left[\frac{255\×({\mathrm{\θ}}_y(x,y)-\min \left({\mathrm{\θ}}_y(x,y\right))}{\max \left({\mathrm{\θ}}_y(x,y)\right)-\min \left({\mathrm{\θ}}_y(x,y)\right)}\right]$

${\mathrm{map}}_z(i,j)=\left[\frac{255\×({\mathrm{\θ}}_z(x,y)-\min \left({\mathrm{\θ}}_z(x,y\right))}{\max \left({\mathrm{\θ}}_z(x,y)\right)\text{-min}\left({\mathrm{\θ}}_z(x,y)\right)}\right]$

最后，再利用水平(X)方向角和垂直(Y)方向角在XZ和YZ平面的矢量距作为灰度值分别构成两幅灰度图像：

${\mathrm{map}}_{\mathrm{xz}}(i,j)=\sqrt{{\mathrm{map}}_x^2+{\mathrm{map}}_z^2}$

${\mathrm{map}}_{\mathrm{yz}}(i,j)=\sqrt{{\mathrm{map}}_y^2+{\mathrm{map}}_z^2};$

c.令原图像的灰度级为X，期望处理后的图像灰度级为Y，原始图像和期望图像的灰度级的分布范围极值分别为X_max，X_min和Y_max，Y_min、期望变换前后的图像对比度保持线性关系，即满足下式：

$\frac{Y-Y_{\min }}{Y_{\max }-Y_{\min }}=\frac{X-X_{\min }}{X_{\max }-X_{\min }}$

经整理线性拉伸的简单数学表达式为：Y＝aX+b

其中 $a=\frac{Y_{\max }-Y_{\min }}{X_{\max }-X_{\min }},$ $b=\frac{X_{\max }{Y}_{\min }-Y_{\max }{X}_{\min }}{X_{\max }-X_{\min }}$

当a＞1时，Y_max-Y_min＞X_max-X_min，则变换后图像的对比度增大，视觉感应为亮处更亮，暗处更暗；而通过设置b值，使整幅图像灰度级的均值改变相应的b值；令Y_max＝255，Y_min＝0；

原始灰度图像的灰度级为L灰度级为i的像素点数为n_i，全部像素为N，归一化直方图： $P_i=\frac{n_i}{N},$ $\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{P}_i=1$

阈值t将灰度级划分为两类：C₀＝(0，1，...，t)和C₁＝(t+1，t+2，...L-1)这两类的概率分别为：

${\mathrm{\ω}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{P}_i=\mathrm{\ω}\left(t\right),$ ${\mathrm{\ω}}_1=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}{P}_i=1-\mathrm{\ω}\left(t\right)$

均值为：

${\mathrm{\μ}}_0=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}i\frac{P_i}{{\mathrm{\ω}}_0}=\frac{\mathrm{\μ}\left(t\right)}{\mathrm{\ω}\left(t\right)},$ ${\mathrm{\μ}}_1=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}i\frac{P_i}{{\mathrm{\ω}}_1}=\frac{{\mathrm{\μ}}_T\left(t\right)-\mathrm{\μ}\left(t\right)}{1-\mathrm{\ω}\left(t\right)}$

上式中 $\mathrm{\μ}\left(t\right)=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{\mathrm{iP}}_i,$ ${\mathrm{\μ}}_T\left(t\right)=\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{iP}}_i$

方差：

${\mathrm{\σ}}_0^2=\underset{i=0}{\overset{t}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_0)}^2\frac{P_i}{{\mathrm{\ω}}_0},$ ${\mathrm{\σ}}_1^2=\underset{i=t+1}{\overset{L-1}{\mathrm{\Σ}}}{(i-{\mathrm{\μ}}_1)}^2\frac{P_i}{{\mathrm{\ω}}_1}$

类内方差为： ${\mathrm{\σ}}_{\mathrm{\ω}}^2={\mathrm{\ω}}_0{\mathrm{\σ}}_0^2+{\mathrm{\ω}}_1{\mathrm{\σ}}_1^2$

类间方差为：

${\mathrm{\σ}}_B^2={\mathrm{\ω}}_0{({\mathrm{\μ}}_0-{\mathrm{\μ}}_T)}^2+{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2={\mathrm{\ω}}_0{\mathrm{\ω}}_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_0)}^2$

总体方差为： ${\mathrm{\σ}}_T^2={\mathrm{\σ}}_B^2+{\mathrm{\σ}}_{\mathrm{\ω}}^2$

改变t的值，使类间方差取得最大值时的t值即为最佳阈值；

用求得的最佳阈值对灰度图像进行二值化；

d.对所得到的二值图像做双向滤波和中值滤波以去掉二值图像中小的噪声点；

e.进行形态学膨胀操作以连接目标细胞断裂的边缘；

f.提取二值图像的最大轮廓并填充后再进行形态学腐蚀操作；

g.用腐蚀后的二值图像与原始图像做逻辑与运算，得到细胞目标图像。

2.根据权利要求1所述角毛藻显微图像细胞目标提取方法，其特征在于，在步骤c中选取使类间方差最大和类内方差最小的图像灰度作为最佳阈值。

3.根据权利要求1所述角毛藻显微图像细胞目标提取方法，其特征在于，在步骤d中，对步骤c所得到的二值图像先进行双向滤波，再进行中值滤波，用于去掉二值图像中小的噪声点。

4.根据权利要求1所述角毛藻显微图像细胞目标提取方法，其特征在于，在步骤f中最大轮廓是指：寻找到的轮廓包含的像素数最多；填充是指：将寻找到的最大轮廓内部填充为白色。

5.根据权利要求1所述角毛藻显微图像细胞目标提取方法，其特征在于，在步骤g中，将腐蚀后的图像与原始图像做逻辑与运算，提取细胞图像的边缘并保留内部纹理。

Images