CN102663406A - 一种基于显微图像的角毛藻和非角毛藻自动分类方法 - Google Patents

Abstract

本发明提供一种基于显微图像的角毛藻和非角毛藻自动分类方法，其中，包括以下步骤：(1)从藻类的显微图像中提取出藻类目标细胞二值化图像；(2)将二值化图像的高度即行数固定为300像素，列数则按原比例进行缩放；(3)采用半径为5像素宽度的圆形结构元素对图像进行闭运算(4)骨架提取：对二值图像采用距离变换细化算法进行骨架提取；(5)提取骨架的节点数和端点数：提取骨架的节点集合和端点集合，从而获取节点和端点的数量；(6)类判决：根据是否满足(节点数＞5且端点数＞5)条件作为是否为角毛藻的判据。本发明自动分类方法准确率高，显像清晰，为进一步的藻种分类识别打下基础。提高进一步藻种自动分类的识别率；也可以为广大第一线藻类监测人员、藻种研究者提供有效地指导和帮助，经济和社会效益都将十分显著。

Description

一种基于显微图像的角毛藻和非角毛藻自动分类方法

技术领域

本发明涉及一种浮游植物分类方法，特别涉及一种基于显微图像的角毛藻和非角毛藻自动分类方法。

背景技术

浮游植物指悬浮于水中的微小藻类植物，是浮游生物中的自养生物部分。浮游植物广泛存在于河流、湖泊和海洋中，多分布于水域的上层，个体极小，需要用显微镜才能观察到，繁殖极快。在淡水中主要是蓝藻、绿藻、硅藻等，在海水中主要是硅藻、甲藻。浮游植物为水生动物(主要是鱼类)的食料；地球上约一半的光合作用是由浮游植物进行的。传统的海洋浮游植物调查，一般是由海洋调查船进行现场采样，将其带回实验室，依赖经验丰富的藻类专家，依据种的生物形态学特征在显微镜下通过人工目视进行分类。该方法不仅费时、费力，且需要研究者具有丰富的海洋浮游植物专业知识和分类经验，难以满足快速分析和现场监测的需求。

为了适应海洋科学和环境检测中快速检测的需要，解决目前显微镜人工鉴定中存在的专业水平高、分类人员断层、耗时等问题，至今已有不少学者对各种浮游植物自动分类方法进行了研究，包括基于藻细胞形态的与计算机技术相结合的显微图像法；基于藻类色素组成的吸收光谱法、荧光光谱法和高压液相色谱法；基于藻细胞大小、色素组成、DNA等的流式细胞仪法；基于细胞基因序列的分子探针法等。其中，日新月异的计算机信息处理技术为基于图像信息并结合生物形态学分类特征的自动化藻类鉴定提供了基础，图像信号经过计算机图像分析与模式识别处理，是传统分类方法实现自动化的有效手段。

角毛藻细胞在壳面顶轴两端各生一细长的角毛，其断面大多为圆形。亦有四角形或多角形等。角毛的表面常有横裂点纹，及螺旋排列的实心小刺。少数种细胞单独生活，大部分种相邻细胞的相近角毛交接，联成直的、弯的或螺旋状的链状群体。链两端细胞的端角毛生出的位置不同，有的种链两端角毛伸出的位置也不同，它们向链外伸出的方向不同。一直以来，生物外形是生物学家对生物进行分类的最重要，也是最基本的手段，而对角毛藻而言，角毛的形态、生长位置、伸出方向，色素体的形态、数量及其在细胞内的分布等都是角毛藻鉴定的重要依据，因此，从角毛藻显微图像中精确地提取目标细胞(尤其是角毛特征)就成为藻类鉴定和研究工作中重要的一步，而且还可以大幅度降低对角毛藻自动分类识别工作的技术要求，提高效率和识别准确率。

中国专利号201010115606.1公开了一种角毛藻显微图像细胞目标提取方法。其主要为图样二值化的方法。这种方法包括如下步骤：

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

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

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

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

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

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

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

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

角毛藻属是硅藻门中心纲盒形藻目角毛藻科下的属，是硅藻门下较大的属之一，其属下包含了大量种类的中国海常见浮游植物，且部分种类属于中国海常见有害赤潮藻种。因此，基于藻种显微图像进行角毛藻和非角毛藻的自动分类，可以提高进一步藻种分类的识别准确率。本发明基于这种细胞目标提取方法，进一步改进计算方法，提高角毛藻和非角毛藻自动分类的识别准确率。本发明的方法准确率高，显像清晰，为进一步的藻种分类识别打下基础。为广大第一线藻类监测人员、藻种研究者提供有效地指导和帮助，经济和社会效益都将十分显著。

发明内容

本发明所要解决的技术问题在于，提供一种角毛藻与非角毛藻自动分类的识别方法，该方法准确率高，显像清晰，为进一步的藻种分类识别打下基础。本发明能够通过藻种显微图像较为准确的区分角毛藻和非角毛藻，从而可以提高进一步藻种自动分类的识别率；而且也可以为广大第一线藻类监测人员、藻种研究者提供有效地指导和帮助，经济和社会效益都将十分显著。

为解决上述技术问题，本发明提供了一种基于显微图像的角毛藻和非角毛藻自动分类方法，其中，包括以下步骤：

(1)从藻类的显微图像中提取出藻类目标细胞二值化图像；

(2)将二值化图像调整为如下尺寸：将二值化图像的高度即行数固定为300像素，列数则按原比例进行缩放；

(3)采用结构元素对图像进行闭运算：所述结构元素是指半径为5像素宽度的圆形结构元素；

(4)骨架提取：对二值图像采用距离变换细化算法进行骨架提取；

(5)提取骨架的节点数和端点数：提取骨架的节点集合和端点集合，从而获取节点和端点的数量；

(6)类判决：根据是否满足(节点数＞5且端点数＞5)条件作为是否为角毛藻的判据：如果满足则为角毛藻，否则为非角毛藻。

步骤(1)从藻类的显微图像中提取出藻类目标细胞二值化图像；具体为：

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

b.构建灰度图像方向角模型，分别构成两幅灰度图像；

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

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

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

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

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

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

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

在步骤c中增强大津法是在运用大津法之前做一个线性拉伸，用于解决图像主体部分与背景灰度差较小的问题；

在步骤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 \left({\mathrm{\θ}}_x(x,y)\right)}{\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)-\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{t}{\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值即为最佳阈值；

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

步骤(2)，用于证后续操作得以快速、有效的进行，防止因原图像尺寸过大引起的计算复杂度过高，或原图像尺寸过小引起的后续操作失真；

步骤(3)，通过闭运算操作连接邻近物体、填充物体内细小空洞、平滑其边界：若结构元素过小则图像改善效果不明显，结构元素过大则会丢失大量图像细节：

图像经二值化操作后会出现断裂、细小空洞、边界粗糙等问题，利用形态学闭运算可以连接邻近物体、填充物体内细小空洞、平滑其边界且不明显改变其面积；用结构元素B对集合A进行闭运算，表示为A·B，定义为：

$A\·B=(A\⊕B)\mathrm{\ΘB}$

其中，符号

表膨胀，Θ表腐蚀；这个公式表明，使用结构元素B对集合A的闭运算就是用B对A进行膨胀，然后用B对结果进行腐蚀，相当于在集合A边界的外部转动“转球”B，A·B的边界为B中的点所能到达的A的边界的最近点集合；

为避免结构元素过小造成的图像改善效果不明显，以及结构元素过大造成的大量图像细节丢失，优选为，半径为5像素宽度的圆形结构元素，用该结构元素对二值化图像进行闭运算。

步骤(4)中基于距离变换细化算法进行骨架提取是指，首先利用距离变换得到藻类目标双像素宽度的骨架，然后采用基于击中-击不中变换的形态学细化法进行进一步的细化处理，最终得到单像素宽度的骨架；

对二值化图像做距离变换，距离变换方法如下：

假设待处理二值化图像I的大小为m×n，i为图像I上的任意像素点，图像上的目标内容用F表示，

为背景内容，B为目标的任意边界点p构成的集合。则目标内容F的距离变换定义为：

$\mathrm{DT}\left(F\right)=\left\{\begin{array}{cc}\underset{p\∈B}{\min }\left\{\mathrm{dist}(i,p)\right\}& i\∈F\\ 0& i\∈\stackrel{\‾}{F}\end{array}\right.$

式中，dist(i，p)表示i、p两点的欧式距离；

为得到目标内容的骨架需要先确定起始骨架点，起始骨架点定义为距离变换值最大的像素点sk，其他骨架点就是在起始骨架点的基础上生长得到的。接着获得当前骨架点及其八邻域对应的边界点，对这9个边界点进行聚类，取聚类中心作为最大圆与边界的切点；其中最大圆是包含在目标内容F中，且不是任何其他包含于F中的圆的真子集的圆；

由当前骨架点与切点的连线将八邻域划分为若干连通区域，理论上每个区域内有且仅有一个邻域骨架点，为了有效消除细微毛刺的影响，实际上采用了多尺度控制的方法来判断存在骨架点的区域：沿最大圆圆周计算相邻切点之间的弦长D，设定弦长的阈值DT，认为当弦长D小于阈值DT的情况下其对应的区域不存在骨架点，在实际操作中，设定弦长的阈值DT＝35；

根据当前骨架点判断邻域骨架点：首先要由当前骨架点与切点的连线将八邻域分成若干区域，在由多尺度控制方法判断存在骨架点的区域内，若已有了骨架点则直接考虑下一个区域；若只有一个像素则该元素即为骨架点；若含有≥2个像素则将与至少2个边界点距离相等的元素作为骨架点；将找到的骨架点作为当前骨架点，继续进行上述操作，直到找出所有的骨架点；

通过距离变换获得的骨架一般为双像素宽度，接着采用基于击中-击不中变换的形态学细化法进行进一步的细化处理，细化处理为：

击中-击不中变换的基本思想是抽取图像中与给定邻域结构相匹配的所有像素，假设二值图像集合用A表示，邻域结构用结构元素B表示，对于二值图像而言，B包含两个结构元素，其中一个结构元素用B¹表示，定义为与前景目标匹配的像素集合；另一个以B²表示，定义为与背景相匹配的像素集合，即B＝{B¹，B²}；

利用结构元素B对集合A作击中-击不中变换(HMT)的结果如下式所示：

HMT_B(A)＝(AΘB¹)∩(A^cΘB²)

其中，AΘB表集合A与结构元素B的差，A^c表集合A的补集；

细化就是在原始图像中减去由击中-击不中变换得到的匹配集合，来寻找图像的骨架。利用结构元素B对图像A作细化表示为

或THIN_B(A)，定义为图像A与HMT_B(A)的算术差：

$A\⊗B={\mathrm{THIN}}_B\left(A\right)=A-{\mathrm{HMT}}_B\left(A\right)$

需要指出的是，仅用一个结构元素对来对任意形状进行细化是不可能的，实际上对于图像A更有效的细化是以结构元素序列为基础的：

{B}＝{B¹，B²，B³，...，Bⁿ}

这里Bⁱ是B^i-1旋转后的形式，因此二值图像A用结构元素序列进行细化可表示如下：

$A\⊗\left\{B\right\}=((\·\·\·((A\⊗B^1)\⊗B^2)\·\·\·)\⊗B^n)$

随着迭代的进行，得到的集合也不断细化，迭代过程不断在这些结构元素中循环，所有结构元素遍历一遍称为一个完整的循环；一个完整的循环结束后，如果所得结果还没有完全细化，继续进行下一个完整循环，直到所有结果不再变化则迭代过程结束，最终得到单像素宽度的骨架；

由于通过距离变换已经得到了双像素宽度的骨架，因此细化操作仅需较少的迭代就可完成，通过距离变换细化算法提取骨架，保证了骨架的连贯性、精确性，与单纯利用细化操作提取骨架相比，大大提高了运行效率。

步骤(5)中利用骨架点邻域矩阵参数打分的方法获取骨架的节点和端点集合，并记录节点及端点的数量；

对提取的骨架利用骨架点邻域矩阵参数打分法获取骨架的节点集合junction_point和端点集合end_point，具体原理如下：

定义一个3×3的掩模SE：

$\mathrm{SE}=\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$

在骨架图像中逐点地移动掩模，在每个像素点(x，y)处，其邻域矩阵参数打分的结果由掩模系数与掩模扫过区域的相应像素值的乘积之和给出；由于掩模系数均为1，因此在像素点(x，y)处，用该掩模得到的邻域矩阵参数打分结果R为：

R＝f(x-1，y-1)+f(x-1，y)+f(x-1，y+1)+

f(x，y-1)+f(x，y)+f(x，y+1)+

f(x+1，y-1)+f(x+1，y)+f(x+1，y+1)

对于经过二值化处理的图像来说，其骨架点的像素值为1，背景点的像素值为0；若某骨架点(x，y)为端点，即其像素值为1，则在其8邻域范围内必定还有一个像素点的值为1，按照上述骨架点邻域矩阵参数打分法得到的结果R为2，若某骨架点(x，y)为节点，即其像素值为1，则在其8邻域范围内必定还至少有3个像素点的值为1，按照上述骨架点邻域矩阵参数打分法得到的结果R至少为4；

判断某骨架点是否为节点或端点的依据为：

$\left\{\begin{array}{ccc}(x,y)\∈\mathrm{end}\_\mathrm{point}& \mathrm{if}& R=2\\ (x,y)\∈\mathrm{junction}\_\mathrm{point}& \mathrm{if}& R\≥4\end{array}\right.$

由节点集合junction_point和端点集合end_point即可得到端点数量e_num及节点数量j_num。

在步骤(6)中角毛藻与非角毛藻的节点数和端点数具有明显差异，将(节点数＞5且端点数＞5)条件作为角毛藻的判据，可以有效的区分角毛藻与非角毛藻：

首先将大量已知角毛藻和非角毛藻图像作为训练样本集，按照步骤(1)～(5)分别提取每幅图像的节点数与端点数，利用统计方法，分析角毛藻和非角毛藻图像节点数与端点数的统计差异，在此基础上设定节点数与端点数的判决门限j_T、e_T，使其对训练样本集达到最高的正确分类率，

最终设定j_T＝5，e_T＝5；

对测试样本集，按照步骤(1)～(5)分别提取每幅图像的节点数j_num与端点数e_num，通过分析j_num、e_num与判决门限j_T、e_T的关系，做出该藻类为角毛藻或非角毛藻的判决，具体的判决依据如下：

本发明为角毛藻与非角毛藻自动分类方法，该方法准确率高，显像清晰，为进一步的藻种分类识别打下基础。本发明能够通过藻种显微图像较为准确的区分角毛藻和非角毛藻，从而可以提高进一步藻种自动分类的识别率；而且也可以为广大第一线藻类监测人员、藻种研究者提供有效地指导和帮助，经济和社会效益都将十分显著。

附图说明

图1为本发明中角毛藻和非角毛藻自动分类方法的流程图；

图2为具体实施中选取的角毛藻(上)与非角毛藻(下)原始图像；

图3为角毛藻(上)与非角毛藻(下)二值化图像；

图4为经尺寸调整的角毛藻(上)与非角毛藻(下)二值化图像；

图5为经闭运算的角毛藻(上)与非角毛藻(下)二值化图像；

图6为对二值图像提取的角毛藻(上)与非角毛藻(下)的骨架；

图7为端点邻域矩阵参数打分结果；

图8为R＝4时的节点邻域矩阵参数打分结果；

图9为角毛藻(上)与非角毛藻(下)节点和端点的提取结果。

具体实施方式

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

参考图1，本发明中角毛藻和非角毛藻自动分类方法的流程图。本发明提供一种基于显微图像的角毛藻和非角毛藻自动分类方法，所述方法包含以下步骤：(1)将藻类显微图像转换为二值化图像；(2)将二值图像调整为适合处理的尺寸；(3)采用合适的结构元素对图像进行闭运算；(4)对二值图像采用距离变换细化算法进行骨架提取；(5)提取骨架的节点数和端点数；(6)根据是否满足(节点数＞5且端点数＞5)条件作为是否为角毛藻的分类判据。

在得到原始图像后，(如图2所示，图2为具体实施中选取的角毛藻(上)与非角毛藻(下)原始图像；)进入步骤(1)，从藻类显微图像中提取藻类目标细胞二值化图像。具体为：

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

b.构建灰度图像方向角模型，分别构成两幅灰度图像；

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

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

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

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

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

如图3所示，图3为角毛藻(上)与非角毛藻(下)的二值化图像。

然后进入步骤(2)，为保证后续操作得以快速、有效的进行，将二值图像自动调整为适合处理的尺寸：图像的高度即行数固定为300像素左右，列数则按原比例进行缩放。图4为经尺寸调整的角毛藻(上)与非角毛藻(下)二值化图像。

之后进入步骤(3)，图像经二值化操作后会出现断裂、细小空洞、边界粗糙等问题，利用形态学闭运算可以连接邻近物体、填充物体内细小空洞、平滑其边界且不明显改变其面积。用结构元素B对集合A进行闭运算，表示为A·B，定义为：

$A\·B=(A\⊕B)\mathrm{\ΘB}$

其中，符号

表膨胀，Θ表腐蚀。这个公式表明，使用结构元素B对集合A的闭运算就是用B对A进行膨胀，然后用B对结果进行腐蚀，相当于在集合A边界的外部转动“转球”B，A·B的边界为B中的点所能到达的A的边界的最近点集合。

为避免结构元素过小造成的图像改善效果不明显，以及结构元素过大造成的大量图像细节丢失，在实际操作中选取半径为5像素宽度的圆形结构元素，用该结构元素对二值图像进行闭运算。图5为经闭运算的角毛藻(上)与非角毛藻(下)二值化图像。

继而进入步骤(4)，首先对步骤(3)得到的二值化图像做距离变换，距离变换的原理如下：

假设待处理二值化图像I的大小为m×n，i为图像I上的任意像素点，图像上的目标内容用F表示，

为背景内容，B为目标的任意边界点p构成的集合。则目标内容F的距离变换定义为：

$\mathrm{DT}\left(F\right)=\left\{\begin{array}{cc}\underset{p\∈B}{\min }\left\{\mathrm{dist}(i,p)\right\}& i\∈F\\ 0& i\∈\stackrel{\‾}{F}\end{array}\right.$

式中，dist(i，p)表示i、p两点的欧式距离。

为得到目标内容的骨架需要先确定起始骨架点，起始骨架点定义为距离变换值最大的像素点sk，其他骨架点就是在起始骨架点的基础上生长得到的。接着获得当前骨架点及其八邻域对应的边界点，对这9个边界点进行聚类，取聚类中心作为最大圆与边界的切点。其中最大圆是包含在目标内容F中，且不是任何其他包含于F中的圆的真子集的圆。

由当前骨架点与切点的连线将八邻域划分为若干连通区域，理论上每个区域内有且仅有一个邻域骨架点。为了有效消除细微毛刺的影响，实际上采用了多尺度控制的方法来判断存在骨架点的区域：沿最大圆圆周计算相邻切点之间的弦长D，设定弦长的阈值DT，认为当弦长D小于阈值DT的情况下其对应的区域不存在骨架点。在实际操作中，设定弦长的阈值DT＝35。

根据当前骨架点判断邻域骨架点：首先要由当前骨架点与切点的连线将八邻域分成若干区域，在由多尺度控制方法判断存在骨架点的区域内，若已有了骨架点则直接考虑下一个区域；若只有一个像素则该元素即为骨架点；若含有≥2个像素则将与至少2个边界点距离相等的元素作为骨架点。将找到的骨架点作为当前骨架点，继续进行上述操作，直到找出所有的骨架点。

显然，由这种方法获得的骨架一般为双像素宽度，在实验中采用基于击中-击不中变换的形态学细化法进行进一步的细化处理。

击中-击不中变换的基本思想是抽取图像中与给定邻域结构相匹配的所有像素。假设二值图像集合用A表示，邻域结构用结构元素B表示。对于二值图像而言，B包含两个结构元素，其中一个结构元素用B¹表示，定义为与前景目标匹配的像素集合；另一个以B²表示，定义为与背景相匹配的像素集合，即B＝{B¹，B²}。

利用结构元素B对集合A作击中-击不中变换(HMT)的结果如下式所示：

HMT_B(A)＝(AΘB¹)∩(A^cΘB²)

其中，AΘB表集合A与结构元素B的差，A^c表集合A的补集。

细化就是在原始图像中减去由击中-击不中变换得到的匹配集合，来寻找图像的骨架。利用结构元素B对图像A作细化表示为或THIN_B(A)，定义为图像A与HMT_B(A)的算术差：

$A\⊗B={\mathrm{THIN}}_B\left(A\right)=A-{\mathrm{HMT}}_B\left(A\right)$

需要指出的是，仅用一个结构元素对来对任意形状进行细化是不可能的。实际上对于图像A更有效的细化是以结构元素序列为基础的：

{B}＝{B¹，B²，B³，...，Bⁿ}

这里Bⁱ是B^i-1旋转后的形式。因此二值图像A用结构元素序列进行细化可表示如下：

$A\⊗\left\{B\right\}=((\·\·\·((A\⊗B^1)\⊗B^2)\·\·\·)\⊗B^n)$

随着迭代的进行，得到的集合也不断细化，迭代过程不断在这些结构元素中循环，所有结构元素遍历一遍称为一个完整的循环。一个完整的循环结束后，如果所得结果还没有完全细化，继续进行下一个完整循环，直到所有结果不再变化则迭代过程结束，最终得到单像素宽度的骨架。

由于通过距离变换已经得到了双像素宽度的骨架，因此细化操作仅需较少的迭代就可完成。通过距离变换细化算法提取骨架，保证了骨架的连贯性、精确性，与单纯利用细化操作提取骨架相比，大大提高了运行效率。图6为对二值图像提取的角毛藻(上)与非角毛藻(下)的骨架。

随后进入步骤(5)，对提取的骨架利用骨架点邻域矩阵参数打分法获取骨架的节点集合junction_point和端点集合end_point，具体原理如下：

定义一个3×3的掩模SE：

$\mathrm{SE}=\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$

在骨架图像中逐点地移动掩模，在每个像素点(x，y)处，其邻域矩阵参数打分的结果由掩模系数与掩模扫过区域的相应像素值的乘积之和给出。由于掩模系数均为1，因此在像素点(x，y)处，用该掩模得到的邻域矩阵参数打分结果R为：

R＝f(x-1，y-1)+f(x-1，y)+f(x-1，y+1)+

f(x，y-1)+f(x，y)+f(x，y+1)+

f(x+1，y-1)+f(x+1，y)+f(x+1，y+1)

对于经过二值化处理的图像来说，其骨架点的像素值为1，背景点的像素值为0。若某骨架点(x，y)为端点，即其像素值为1，则在其8邻域范围内必定还有一个像素点的值为1，按照上述骨架点邻域矩阵参数打分法得到的结果R为2，如图7所示；若某骨架点(x，y)为节点，即其像素值为1，则在其8邻域范围内必定还至少有3个像素点的值为1，按照上述骨架点邻域矩阵参数打分法得到的结果R至少为4，如图8所示。即判断某骨架点是否为节点或端点的依据为：

$\left\{\begin{array}{ccc}(x,y)\∈\mathrm{end}\_\mathrm{point}& \mathrm{if}& R=2\\ (x,y)\∈\mathrm{junction}\_\mathrm{point}& \mathrm{if}& R\≥4\end{array}\right.$

图9为角毛藻(上)与非角毛藻(下)节点和端点的提取结果。由节点集合junction_point和端点集合end_point即可得到端点数量e_num及节点数量j_num。

最后进入步骤(6)，首先将大量已知角毛藻和非角毛藻图像作为训练样本集，按照步骤(1)～(5)分别提取每幅图像的节点数与端点数。利用统计方法，分析角毛藻和非角毛藻图像节点数与端点数的统计差异，在此基础上设定节点数与端点数的判决门限j_T、e_T，使其对训练样本集达到最高的正确分类率。

经过对训练样本集的大量实验，最终设定j_T＝5，e_T＝5。当然，需要指出的是：由于显微镜参数的不同、拍摄环境和水平的差异等因素的影响，会造成节点数、端点数判决门限的改变，因此j_T＝5，e_T＝5只是在一定实验条件下得出的判决门限，在具体实施中需实验者严格按照步骤(1)～(5)得出精确的判决门限。

对测试样本集，仍按照步骤(1)～(5)分别提取每幅图像的节点数j_num与端点数e_num，通过分析j_num、e_num与判决门限j_T、e_T的关系，做出该藻类为角毛藻或非角毛藻的判决。具体的判决依据如下：

经大量测试样本集验证，本专利提供了一种对角毛藻和非角毛藻进行分类的有效方法。所有上述为这一知识产权的首要实施系统，并没有设定限制以其他形式实施这种新系统。本领域技术人员将利用这一重要信息，对上述内容修改，以实现类似的执行情况。但是，所有基于本发明的修改或改造新方法，属于保留的权利。

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

Claims (6)

1.一种基于显微图像的角毛藻和非角毛藻自动分类方法，其特征在于，包括以下步骤：

(1)从藻类的显微图像中提取出藻类目标细胞二值化图像；

(2)将二值化图像调整为如下尺寸：将二值化图像的高度即行数固定为300像素，列数则按原比例进行缩放；

(3)采用结构元素对图像进行闭运算：所述结构元素是指半径为5像素宽度的圆形结构元素；

(4)骨架提取：对二值图像采用距离变换细化算法进行骨架提取；

(5)提取骨架的节点数和端点数：提取骨架的节点集合和端点集合，从而获取节点和端点的数量；

(6)类判决：根据是否满足(节点数＞5且端点数＞5)条件作为是否为角毛藻的判据：如果满足则为角毛藻，否则为非角毛藻。

2.根据权利要求1所述基于显微图像的角毛藻和非角毛藻自动分类方法，其特征在于，步骤(1)从藻类的显微图像中提取出藻类目标细胞二值化图像；具体为：

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

b.构建灰度图像方向角模型，分别构成两幅灰度图像；

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

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

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

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

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

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

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

在步骤c中增强大津法是在运用大津法之前做一个线性拉伸，用于解决图像主体部分与背景灰度差较小的问题；

在步骤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 \left({\mathrm{\θ}}_x(x,y)\right)}{\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)-\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将灰度级划分为两类：C0＝(0，1，...，t)和C1＝(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{t}{\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值即为最佳阈值；

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

3.根据权利要求1所述基于显微图像的角毛藻和非角毛藻自动分类方法，其特征在于，步骤(2)，用于以保证后续操作得以快速、有效的进行，防止因原图像尺寸过大引起的计算复杂度过高，或原图像尺寸过小引起的后续操作失真；

步骤(3)，通过闭运算操作连接邻近物体、填充物体内细小空洞、平滑其边界：若结构元素过小则图像改善效果不明显，结构元素过大则会丢失大量图像细节：

图像经二值化操作后会出现断裂、细小空洞、边界粗糙等问题，利用形态学闭运算可以连接邻近物体、填充物体内细小空洞、平滑其边界且不明显改变其面积；用结构元素B对集合A进行闭运算，表示为A·B，定义为：

$A\·B=(A\⊕B)\mathrm{\ΘB}$

其中，符号表膨胀，Θ表腐蚀；这个公式表明，使用结构元素B对集合A的闭运算就是用B对A进行膨胀，然后用B对结果进行腐蚀，相当于在集合A边界的外部转动“转球”B，A·B的边界为B中的点所能到达的A的边界的最近点集合；

为避免结构元素过小造成的图像改善效果不明显，以及结构元素过大造成的大量图像细节丢失，优选为，半径为5像素宽度的圆形结构元素，用该结构元素对二值化图像进行闭运算。

4.根据权利要求1所述基于显微图像的角毛藻和非角毛藻自动分类方法，其特征在于，步骤(4)中基于距离变换细化算法进行骨架提取是指，首先利用距离变换得到藻类目标双像素宽度的骨架，然后采用基于击中-击不中变换的形态学细化法进行进一步的细化处理，最终得到单像素宽度的骨架；

对二值化图像做距离变换，距离变换方法如下：

假设待处理二值化图像I的大小为m×n，i为图像I上的任意像素点，图像上的目标内容用F表示，为背景内容，B为目标的任意边界点p构成的集合。则目标内容F的距离变换定义为：

$\mathrm{DT}\left(F\right)=\left\{\begin{array}{cc}\underset{p\∈B}{\min }\left\{\mathrm{dist}(i,p)\right\}& i\∈F\\ 0& i\∈\stackrel{\‾}{F}\end{array}\right.$

式中，dist(i，p)表示i、p两点的欧式距离；

为得到目标内容的骨架需要先确定起始骨架点，起始骨架点定义为距离变换值最大的像素点sk，其他骨架点就是在起始骨架点的基础上生长得到的。接着获得当前骨架点及其八邻域对应的边界点，对这9个边界点进行聚类，取聚类中心作为最大圆与边界的切点；其中最大圆是包含在目标内容F中，且不是任何其他包含于F中的圆的真子集的圆；

由当前骨架点与切点的连线将八邻域划分为若干连通区域，理论上每个区域内有且仅有一个邻域骨架点，为了有效消除细微毛刺的影响，实际上采用了多尺度控制的方法来判断存在骨架点的区域：沿最大圆圆周计算相邻切点之间的弦长D，设定弦长的阈值DT，认为当弦长D小于阈值DT的情况下其对应的区域不存在骨架点，在实际操作中，设定弦长的阈值DT＝35；

根据当前骨架点判断邻域骨架点：首先要由当前骨架点与切点的连线将八邻域分成若干区域，在由多尺度控制方法判断存在骨架点的区域内，若已有了骨架点则直接考虑下一个区域；若只有一个像素则该元素即为骨架点；若含有≥2个像素则将与至少2个边界点距离相等的元素作为骨架点；将找到的骨架点作为当前骨架点，继续进行上述操作，直到找出所有的骨架点；

获得的骨架一般为双像素宽度，接着采用基于击中-击不中变换的形态学细化法进行进一步的细化处理，细化处理为：

击中-击不中变换的基本思想是抽取图像中与给定邻域结构相匹配的所有像素，假设二值图像集合用A表示，邻域结构用结构元素B表示，对于二值图像而言，B包含两个结构元素，其中一个结构元素用B¹表示，定义为与前景目标匹配的像素集合；另一个以B²表示，定义为与背景相匹配的像素集合，即B＝{B¹，B²}；

利用结构元素B对集合A作击中-击不中变换(HMT)的结果如下式所示：

HMT_B(A)＝(AΘB¹)∩(A^cΘB²)

其中，AΘB表集合A与结构元素B的差，A^c表集合A的补集；

细化就是在原始图像中减去由击中-击不中变换得到的匹配集合，来寻找图像的骨架。利用结构元素B对图像A作细化表示为

或THIN_B(A)，定义为图像A与HMT_B(A)的算术差：

$A\⊗B={\mathrm{THIN}}_B\left(A\right)=A-{\mathrm{HMT}}_B\left(A\right)$

需要指出的是，仅用一个结构元素对来对任意形状进行细化是不可能的，实际上对于图像A更有效的细化是以结构元素序列为基础的：

{B}＝{B¹，B²，B³，...，Bⁿ}

这里Bⁱ是B^i-1旋转后的形式，因此二值图像A用结构元素序列进行细化可表示如下：

$A\⊗\left\{B\right\}=((\·\·\·((A\⊗B^1)\⊗B^2)\·\·\·)\⊗B^n)$

随着迭代的进行，得到的集合也不断细化，迭代过程不断在这些结构元素中循环，所有结构元素遍历一遍称为一个完整的循环；一个完整的循环结束后，如果所得结果还没有完全细化，继续进行下一个完整循环，直到所有结果不再变化则迭代过程结束，最终得到单像素宽度的骨架；

由于通过距离变换已经得到了双像素宽度的骨架，因此细化操作仅需较少的迭代就可完成，通过距离变换细化算法提取骨架，保证了骨架的连贯性、精确性，与单纯利用细化操作提取骨架相比，大大提高了运行效率。

5.根据权利要求1所述基于显微图像的角毛藻和非角毛藻自动分类方法，其特征在于，步骤(5)中利用骨架点邻域矩阵参数打分的方法获取骨架的节点和端点集合，并记录节点及端点的数量；

对提取的骨架利用骨架点邻域矩阵参数打分法获取骨架的节点集合junction_point和端点集合end_point，具体原理如下：

定义一个3×3的掩模SE：

$\mathrm{SE}=\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]$

在骨架图像中逐点地移动掩模，在每个像素点(x，y)处，其邻域矩阵参数打分的结果由掩模系数与掩模扫过区域的相应像素值的乘积之和给出；由于掩模系数均为1，因此在像素点(x，y)处，用该掩模得到的邻域矩阵参数打分结果R为：

R＝f(x-1，y-1)+f(x-1，y)+f(x-1，y+1)+

f(x，y-1)+f(x，y)+f(x，y+1)+

f(x+1，y-1)+f(x+1，y)+f(x+1，y+1)

对于经过二值化处理的图像来说，其骨架点的像素值为1，背景点的像素值为0；若某骨架点(x，y)为端点，即其像素值为1，则在其8邻域范围内必定还有一个像素点的值为1，按照上述骨架点邻域矩阵参数打分法得到的结果R为2，若某骨架点(x，y)为节点，即其像素值为1，则在其8邻域范围内必定还至少有3个像素点的值为1，按照上述骨架点邻域矩阵参数打分法得到的结果R至少为4；

判断某骨架点是否为节点或端点的依据为：

$\left\{\begin{array}{ccc}(x,y)\∈\mathrm{end}\_\mathrm{point}& \mathrm{if}& R=2\\ (x,y)\∈\mathrm{junction}\_\mathrm{point}& \mathrm{if}& R\≥4\end{array}\right.$

由节点集合junction_point和端点集合end_point即可得到端点数量e_num及节点数量j_num。

6.根据权利要求1所述基于显微图像的角毛藻和非角毛藻自动分类方法，其特征在于，在步骤(6)中角毛藻与非角毛藻的节点数和端点数具有明显差异，将(节点数＞5且端点数＞5)条件作为角毛藻的判据，可以有效的区分角毛藻与非角毛藻：

首先将大量已知角毛藻和非角毛藻图像作为训练样本集，按照步骤(1)～(5)分别提取每幅图像的节点数与端点数，利用统计方法，分析角毛藻和非角毛藻图像节点数与端点数的统计差异，在此基础上设定节点数与端点数的判决门限j_T、e_T，使其对训练样本集达到最高的正确分类率，

最终设定j_T＝5，e_T＝5；

对测试样本集，按照步骤(1)～(5)分别提取每幅图像的节点数j_num与端点数e_num，通过分析j_num、e_num与判决门限j_T、e_T的关系，做出该藻类为角毛藻或非角毛藻的判决，具体的判决依据如下：

Images