一种小麦条锈病菌夏孢子自动统计计数方法与流程

本发明属于农业图像处理技术领域，特别涉及一种小麦条锈病菌夏孢子自动统计计数方法。

背景技术：

由条形柄锈菌(pucciniastriiformisf.sp.tritici)引起的小麦条锈病是我国西北、西南、华北和淮北等冬麦区和西北春麦区小麦的重要病害，一般流行年份可致小麦减产10％～20％，特大流行年份可致小麦减产50％～60％，对我国粮食安全生产具有重大威胁。小麦条锈病是一种气传病害，病原菌主要以夏孢子随气流传播到小麦上完成周年侵染循环，因此，田间空气中的夏孢子数量与小麦条锈病的流行程度紧密相关。对田间空气中条锈病菌夏孢子的定量监测和准确计数，并利用计数结果进行病害的早期预测，是及早制定正确的防控策略和采取防治措施的一个重要依据。

目前，小麦条锈病菌夏孢子数量多采用孢子捕捉器取样和监测，然后将黏附有孢子的载玻片拿回实验室在显微镜下人工计数或分子生物学方法计数。传统的人工计数方法是在光学显微镜下通过肉眼观测载玻片确定孢子个数，或者是将黏着的孢子洗刷下来，配成悬浮液，制成临时玻片进行镜检计数，均存在工作量大、效率低、且随工作时间准确性降低等缺点。分子生物学方法需要根据病原孢子dna的量实现计数，过程复杂，对技术、试验条件和仪器要求高、花费高，很难大面积推广应用。因此，急需一种简便、快捷、准确、高效的小麦条锈病菌夏孢子计数方法。

近年来，基于显微图像处理技术被逐步应用到气传植物病原真菌孢子检测的相关研究上，并逐步应用到小麦条锈病菌夏孢子的相关研究中。齐龙等提出了基于距离变换和高斯滤波的改进分水岭算法的稻瘟病菌显微图像孢子自动检测和计数方法，但当粘连孢子的接触线长度超过单个孢子的宽度时，容易造成局部极小值点间的距离小于结构元素的长度，引起漏分割现象。李小龙等基于k-means聚类和分水岭分割算法，对小麦条锈病菌夏孢子显微图像进行处理，实现了对夏孢子的自动计数。然而图像中的噪声、结构特征和孢子表面细微的灰度变化等因素，易导致出现极小值点过多、分割位置不准确，从而产生过度分割现象。张荣标等提出了图像处理和支持向量机的圆褐固氮菌浓度快速检测方法，但方法仅适用于圆褐固氮菌相互粘连较少的简单情况。此外，上述方法中，大多数方法对未粘连的孢子能很好的计数，但对经常发生的多孢子粘连分割计数问题尚未解决。

技术实现要素：

为了克服上述现有技术的缺点，本发明的目的在于提供一种小麦条锈病菌夏孢子自动统计计数方法，利用一系列的显微图像处理方法实现夏孢子的自动统计与计数，以期提高小麦条锈病菌夏孢子计数精度，解决多孢子粘连分割计数的难题，并为在线式小麦条锈病夏孢子监测装备的开发提供技术支持。

为了实现上述目的，本发明采用的技术方案是：

一种小麦条锈病菌夏孢子自动统计计数方法，其特征在于包括以下步骤：

步骤1，夏孢子显微图像采集

用tpbz3型孢子捕捉仪模拟捕捉田间空气中的小麦条锈病菌夏孢子。在野外小麦田间，抽出捕捉仪的载波器并将表面均匀涂抹一薄层凡士林的载玻片放其上(涂面朝上)，捕捉一定时间后取出载玻片。为了使玻片上黏附的孢子密度不同，按照上述方法制备凡士林玻片30片，获取捕捉时间分别为120、180、240min的条锈病菌夏孢子载玻片，重复10次，共获不同孢子密度的载玻片30片；

用bx52型倒置显微镜对孢子临时玻片进行观察、拍照，放大倍数为10×20。在显微镜下，分别对30片载玻片进行显微图像采集，每片随机选取5个视野拍照，共获得150幅供试夏孢子图像，图像尺寸为4140×3096pixel，图像分辨率为72dpi，存贮为bmp格式；

步骤2，显微图像选取

由于显微图像尺寸过大，为减少图像处理运算量，提高检测和计数速度，将所述步骤s1的原始夏孢子显微图像用双线性插值方法缩小为911×682pixel。然后随机选取30幅用于算法训练，其余120幅图像用于测试；

s3、k-means聚类分割

将所述步骤s2的夏孢子显微图像用k-means聚类算法实现孢子和背景的分割，即将其分为孢子目标和背景2大类，然后采用全局阈值算法进行二值化处理得到二值图像；

s4、形态学处理

二值图像中有小面积区域等噪声，且夏孢子目标边缘有凸刺或内部含有孔洞，故首先用填充操作对所述步骤s3中的二值图像进行孔洞填充，再用完整孢子中最小面积值为阈值，去除小面积区域噪声，此时当边界上孢子所占面积小于完整孢子最小面积的均从图像中移除，计数时作为其他视野中的孢子。最后，用圆盘结构元素对二值图像进行形态学开运算，以消除夏孢子边界小的凸刺；

s5、基于形状因子和面积的粘连孢子判别

所述步骤s4中的二值图像中既有单个孢子，也有2个或多个孢子粘连在一起的孢子群。通过观察可知，粘连孢子群的区域轮廓要比单个孢子区域的轮廓复杂，故选择描述目标边界复杂程度的形状因子作为粘连孢子判别的依据。形状因子公式为：

sf＝4πs/l²

式中sf为形状因子；s为一个连通区域的面积像素值；l为连通区域的周长像素值；

经粘连孢子判别后，对判定为单个孢子的区域直接采用最小二乘椭圆拟合算法进行拟合并记录椭圆个数num供后续统计计数。对粘连孢子区域将每个孢子一一分割出来并自动计数是本发明下面拟解决的问题；

s6、基于凹度的粘连孢子轮廓分割

对所述步骤s5二值图像中的粘连孢子区域通过移除内部像素点操作提取出边缘轮廓，然后遍历图像跟踪边缘轮廓点，并将轮廓点坐标保存在一个有序表中。pt(xt,yt)是边缘轮廓上任意一点，两向量ptpt-k和ptpt+k之间的夹角称为pt的凹度，pt-k和pt+k表示点pt的相邻轮廓点。凹度concavity的计算公式为：

其中凹点是满足如下两个条件的轮廓点：

(1)凹度concavity(pt)在角angle(δ1,δ2)范围之内；

(2)直线不包含在粘连孢子内；

粘连孢子边缘轮廓由轮廓上的所有凹点和由凹点分割成的多个轮廓段组成，如下式：

c＝cs1+...+csi+...+csn+cp1+...+cpj+...+cpk

式中c为粘连孢子轮廓；csi为第i个轮廓段；n为轮廓段总数；cpj为第j个凹点；k为凹点总数；

删除凹点后，将粘连孢子边缘轮廓分割成多个轮廓段，通过连通区域标记，将轮廓段保存在一个有序的cell结构轮廓段集中，用于后续孢子轮廓段融合处理；

s7、粘连孢子轮廓段融合

经所述步骤s6的轮廓分割后同一孢子的轮廓可能被分割成多个轮廓段，所以，识别出同一孢子的轮廓段并融合，以获得各个孢子尽可能完整的椭圆拟合数据是本发明的关键。故本发明以轮廓段距离测量方法为准则确定候选轮廓段，对候选轮廓段进行最小二乘椭圆拟合并用偏移误差方法进行评价，将符合条件的轮廓段融合成新的轮廓段，对新的轮廓段进行椭圆拟合即为最后的正确椭圆；

s8、最后统计所有椭圆的总数即为夏孢子的个数。

所述步骤s4中的圆盘结构元素，大小为3×3。

所述步骤s5中的基于形状因子和面积的粘连孢子判别，形状因子阈值sf0＝0.8，面积阈值smax＝560，依次提取单个连通区域i，若i满足公式：

sfi>sf0&si<smax

则判定区域i为未粘连的单个夏孢子，否则为粘连孢子。

所述步骤s6中的凹度的角δ1和δ2分别设定为50°和150°。

所述步骤s7的轮廓段距离方法，第i和第j个轮廓段之间的距离dm定义为公式：

式中d(pi1,pj1)、分别为轮廓段csi和csj的2端点位置的欧氏距离；为轮廓段csi和csj的中点位置的欧式距离。

所述步骤s7的最小二乘椭圆拟合算法的具体步骤为：

在二维平面坐标系中，一般的椭圆曲线方程可以用2个向量相乘的隐式方程来表示：

f(α,x)＝x·α＝ax²+bxy+cy²+dx+ey+f＝0

式中，α＝[abcdef]^t为椭圆方程的系数；x和y分别是曲线上点横、纵坐标，x＝[x²xyy²xy1]；

加约束条件：α^tcα＝1，以保证拟合结果为椭圆，

其中c为一个6×6的矩阵，

方程f(α,xi)为二值图像中边缘点(xi,yi)到给定椭圆方程的代数距离。根据最小二乘原理，椭圆拟合问题通过将代数距离平方之和

最小化来解决。方程改写为向量形式：e＝||dα||²

其中d为一个m×6的矩阵，

引入拉格朗日系数并微分得：

式中λ为特征值；w为一个6×6的散射矩阵，w＝d^td；

根据广义特征值求解方法可得到：

式中λi和ui分别为特征值和特征向量。由此解得至多6个实数解(λi,αi)；

推导出数据点到椭圆的代数距离平方和公式为：

e＝||dα||²＝α^td^tdα＝α^twα＝λα^tcα＝λ

由上式可知所需要的是最小的正特征值λi所对应的特征向量αi。得到解αi之后，即可实现对椭圆的拟合。

所述步骤s7的偏移误差公式为：

dem(cs^#,ce)＝e/m^#＝λm^#

式中e为由数据点到椭圆的代数距离平方和公式求出的给定轮廓段cs^#到椭圆ce的最小二乘代数距离；m^#为轮廓段cs^#上点的总数。

若偏移误差小于阈值，则该轮廓段属于同一孢子的轮廓段，否则不是。

所述步骤s7的轮廓段融合步骤具体为：

(1)从轮廓段集中选取轮廓段中最长轮廓段cs1；

(2)将孢子轮廓段中剩余的轮廓段csi分别与cs1作距离测量，设定距离测量阈值ωdm，若满足dm(cs1,csi)&lt;ωdm条件，则有

cs^*＝{csi|dm(cs1,csi)&lt;ωdm,i＝1,2,...,k}；

(3)将cs1依次与cs^*中每一个轮廓段csi作椭圆拟合，同时求出相应的候选椭圆和偏移误差dem；

(4)若偏移误差dem小于预设的阈值σdem，则融合成新的轮廓段csnew：

csnew＝{cs1+csi|dem(cs1,csi)&lt;σdem,i＝1,2,...,k}；

(5)对csnew作椭圆拟合，作为一个孢子的正确拟合椭圆，孢子数量计数num+1，并且删除轮廓段集中的csnew；

(6)最后，判断轮廓段集中的轮廓段是否全部删除。若是，则输出所有正确拟合椭圆并且结束程序，否则返回到步骤1。

所述步骤s7中的距离测量阈值ωdm和偏移误差阈值σdem由30幅夏孢子图像样本的训练确定，分别为40和95。

与现有技术相比，本发明通过k-means聚类分割及形态学处理，能够准确地将夏孢子从背景中分割出来，且保留了夏孢子原有形态。提出基于夏孢子形状因子和面积特征的粘连孢子群判别方法，可正确地判别出未粘连的单个孢子和粘连孢子。

同时，本发明小麦条锈病菌夏孢子自动统计计数方法，基于凹度的轮廓分割和轮廓段融合方法对粘连孢子进行处理，可准确识别出粘连孢子中同一孢子的轮廓段，通过最小二乘椭圆拟合算法对夏孢子进行拟合，可有效对粘连孢子进行分割并计数，提高了孢子计数结果的准确性。经过试验结果表明，本发明最低计数准确率为92.7％，最高计数准确率为100％，总平均计数准确率为98.6％，具有较高的孢子计数精度。为在线式小麦条锈病夏孢子监测装备的开发提供技术支持。

附图说明

图1为本发明的自动统计计数方法处理流程图。

图2为条锈病菌夏孢子原始图像。

图3为k-means聚类后的夏孢子图像。

图4为夏孢子二值图像。

图5为形态学处理后的夏孢子二值图像。

图6为单个孢子二值图像。

图7为粘连孢子二值图像。

图8为本发明粘连夏孢子椭圆拟合和计数过程实例：a.预处理后的二值图像；b.粘连孢子的轮廓；c.轮廓上十字状黑色凹点；d.轮廓段；e.候选拟合椭圆；f.孢子正确拟合椭圆；g.最终分割结果。

图9为本发明分割及计数效果图。

具体实施方式

下面结合附图并通过实施例对本发明作进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于此。

本发明的小麦条锈病菌夏孢子自动统计计数方法处理流程如图1所示，具体包含以下步骤：

1、供试材料与设备

供试小麦条锈病菌夏孢子在西北农林科技大学植物保护学院小麦抗病种质资源创制与利用研究中心的东南窑低温温室内繁育，其他材料和仪器有凡士林、载玻片、洗耳球、便携式定量风流孢子捕捉仪(型号：tpbz3，浙江托普云农公司)、研究级倒置显微镜(型号：bx52，olympus，japan)、数码成像系统(型号：dp72，olympus，japan)、cellsensstandard显微图像采集软件(olympus，japan)和台式计算机。

2、夏孢子显微图像采集

用tpbz3型孢子捕捉仪模拟捕捉田间空气中的小麦条锈病菌夏孢子。在野外小麦田间，抽出捕捉仪的载波器并将表面均匀涂抹一薄层凡士林的载玻片放其上(涂面朝上)，捕捉一定时间后取出载玻片。为了使玻片上黏附的孢子密度不同，按照上述方法制备凡士林玻片30片，获取捕捉时间分别为60、120、180、240min的条锈病菌夏孢子载玻片，重复10次，共获不同孢子密度的载玻片30片。

用bx52型倒置显微镜对孢子临时玻片进行观察、拍照，放大倍数为10×20。在显微镜下，分别对30片载玻片进行显微图像采集，每片随机选取5个视野拍照，共获得150幅供试夏孢子图像，图像尺寸为4140×3096pixel，图像分辨率为72dpi，存贮为bmp格式。

3、显微图像选取

上述二个步骤之后，由于显微图像尺寸过大，为减少图像处理运算量，提高检测和计数速度，将原始图像用双线性插值方法缩小为911×682pixel，然后随机选取30幅用于算法训练，其余120幅图像用于测试。

4k-means聚类分割

为了对夏孢子进行计数，首先需要将夏孢子从背景中分割出来。夏孢子目标分割可看成聚类问题，即图像中像素点类别未知的前提下，根据像素点的特征值，将图像划分为若干个区域。考虑到k-means聚类算法是一种准确、高效的目标分割算法，故本发明用k-means聚类算法实现对孢子图像的分割，即将其分为孢子目标和背景2大类，图2所示的夏孢子图像经k-means聚类后的分割结果如图3所示，可见，该算法能够较准确地将夏孢子分割出来。对图3采用全局阈值算法进行二值化处理后的结果如图4所示。

5形态学处理

由图4可知，二值图像中有小面积区域等噪声，且夏孢子目标边缘有凸刺或内部含有孔洞，故首先用填充操作对图像进行孔洞填充，再用完整孢子中最小面积值为阈值，去除小面积区域噪声，此时当边界上孢子所占面积小于完整孢子最小面积的均从图像中移除，计数时作为其他视野中的孢子。最后，为消除夏孢子边界小的凸刺，经试验用3×3的圆盘结构元素对二值图像进行形态学开运算，处理结果如图5所示，可见，经上述处理后，更好地保留了图像中孢子目标的原有形态，有利于进行后续处理。

6基于形状因子和面积的粘连孢子判别

从图5可看出，二值图像中既有单个孢子，也有2个或多个孢子粘连在一起的孢子。通过观察可知，粘连孢子的区域轮廓要比单个孢子区域的轮廓复杂，故选择描述目标边界复杂程度的形状因子和面积作为粘连孢子判别的依据。形状因子公式为：

sf＝4πs/l²(1)

式中sf为形状因子；s为一个连通区域的面积像素值；l为连通区域的周长像素值。

当多个孢子相互粘连时，粘连孢子边界会出现凹陷而变得复杂，形状因子会相应变小，故sf对孢子是否粘连有很好的区分度。对30幅夏孢子图像中单个孢子和粘连孢子sf和s的统计可知，单个孢子sf范围为0.9080～1.0912，s为200～523pixel；粘连孢子的sf在0.2625～0.7606，s为600～2301pixel。故设置形状因子阈值sf0＝0.8、面积阈值smax＝560，作为判定单个孢子和粘连孢子的依据。依次提取单个连通区域i，若i满足

sfi>sf0&si<smax(2)

则判定区域i为未粘连的单个夏孢子，否则为粘连孢子。对图5进行粘连孢子判别后，单个夏孢子的二值图像如图6所示，对判定为单个孢子的区域直接采用椭圆拟合算法进行拟合并记录椭圆个数num供后续统计计数。粘连夏孢子的二值图像如图7所示。可见，单个夏孢子和粘连孢子被正确判别出来。

7基于凹度的粘连孢子轮廓分割

由图7可知，图像中孢子目标粘连或部分重叠，故轮廓往往是多个目标轮廓的混合集，由于孢子的外形似椭圆，椭圆边界上所有点的凹度变化是连续的，当椭圆间出现粘连或重叠时凹度会变小，基于此思想，本发明依次提取粘连孢子边界轮廓，并基于凹度提取出轮廓上的凹点，由凹点将轮廓分割成多个轮廓段。

图8为本发明粘连夏孢子椭圆拟合和计数过程，图7中矩形框内的粘连孢子放大后如图8a所示。对二值图像中的粘连孢子通过移除内部像素点操作提取出边缘轮廓(图8b)，然后遍历图像跟踪边缘轮廓点，并将轮廓点坐标保存在一个有序表中。pt(xt,yt)是边缘轮廓上任意一点，两向量ptpt-k和ptpt+k之间的夹角称为pt的凹度，pt-k和pt+k表示点pt的相邻轮廓点，根据预试验本文设定k值为3。凹度concavity的计算公式为：

其中凹点是满足如下两个条件的轮廓点：

(1)凹度concavity(pt)在角angle(δ1,δ2)范围之内；

(2)直线不包含在粘连孢子内；

本文根据30幅夏孢子图像样本的训练统计得出凹度的角δ1和δ2分别设定为50°和150°。

粘连孢子边缘轮廓由轮廓上的所有凹点和由凹点分割成的多个轮廓段组成，如式(4)：

c＝cs1+...+csi+...+csn+cp1+...+cpj+...+cpk(4)

式中c为粘连孢子轮廓；csi为第i个轮廓段；n为轮廓段总数；cpj为第j个凹点；k为凹点总数。

基于凹度选取的轮廓凹点如图8c所示，删除凹点后，将粘连孢子边缘轮廓分割成多个轮廓段，通过连通区域标记，将轮廓段保存在一个有序的cell结构轮廓段集中，用于后续孢子轮廓段融合处理(图8d)。

8粘连孢子轮廓段融合

轮廓分割后同一孢子的轮廓可能被分割成多个轮廓段，所以，识别出同一孢子的轮廓段并融合，以获得各个孢子尽可能完整的椭圆拟合数据是本发明的关键。故本发明以轮廓段距离测量方法为准则确定候选轮廓段，对候选轮廓段进行最小二乘椭圆拟合并用偏移误差方法评价，将符合条件的轮廓融合成新的轮廓段，对新的轮廓段进行椭圆拟合即为最后的正确椭圆。

8.1轮廓段的距离测量方法

若2个轮廓段之间的距离很大，则属于同一孢子的可能性就很小，反之亦然。基于这一思想，本发明提出轮廓段的距离测量方法，判断轮廓段之间的位置关系，以排除属于同一孢子可能性小的轮廓段，确定候选轮廓段组合。第i和第j个轮廓段之间的距离dm定义为：

式中d(pi1,pj1)、分别为轮廓段csi和csj的2端点位置的欧氏距离；为轮廓段csi和csj的中点位置的欧式距离。

8.2最小二乘椭圆拟合算法

小麦条锈病菌夏孢子形状类似椭圆，故本发明用最小二乘椭圆拟合算法进行夏孢子轮廓拟合。

在二维平面坐标系中，一般的椭圆曲线方程可以用2个向量相乘的隐式方程来表示：

f(α,x)＝x·α＝ax²+bxy+cy²+dx+ey+f＝0(6)

式中α＝[abcdef]^t为椭圆方程的系数；x和y分别是曲线上点横、纵坐标，x＝[x²xyy²xy1]。

为保证拟合结果为椭圆，需要加约束条件：

α^tcα＝1(7)

其中c为一个6×6的矩阵，

方程f(α,xi)为二值图像中边缘点(xi,yi)到给定椭圆方程的代数距离。根据最小二乘原理，椭圆拟合问题通过将代数距离平方之和

最小化来解决。方程(9)改写为向量形式：

e＝||dα||²(10)

其中d为一个m×6的矩阵，

引入拉格朗日系数并微分可得：

式中λ为方程(12)的特征值；w为一个6×6的散射矩阵，

w＝d^td(13)

根据广义特征值求解方法可得到：

式中λi和ui分别为方程的特征值和特征向量。由此解得至多6个实数解(λi,αi)。

推导出数据点到椭圆的代数距离平方和公式：

e＝||dα||²＝α^td^tdα＝α^twα＝λα^tcα＝λ(15)

由式(15)可知所需要的是最小的正特征值λi所对应的特征向量αi。得到解αi之后，即可实现对椭圆的拟合。

8.3轮廓段到拟合椭圆的偏移误差方法

通过候选轮廓段进行椭圆拟合时，并未考虑轮廓段与椭圆的拟合程度，导致得到的候选椭圆与实际有较大偏差，故需要对候选轮廓段进行评价筛选，本文以偏移误差作为去除错误候选轮廓段的评价条件，偏移误差公式为：

dem(cs^#,ce)＝e/m^#＝λ/m^#(16)

式中e为由式(15)求出的给定轮廓段cs^#到椭圆ce的最小二乘代数距离；m^#为轮廓段cs^#上点的总数。

若偏移误差小于阈值，则该轮廓段属于同一孢子的轮廓段，否则不是。

8.4轮廓段融合步骤

以图8d为例对融合算法进一步说明，具体步骤为：

(1)首先选取图8d中最长的轮廓段cs8；

(2)依次计算cs8与剩余的csi的距离dm，通过计算，cs8与cs4、cs5、cs6和cs7的距离测量值小于预设的阈值ωdm，故一个新的轮廓段集cs*包括cs4、cs5、cs6和cs7；

(3)然后依次作轮廓段(cs8和cs4、cs8和cs5、cs8和cs6、cs8和cs7)的候选拟合椭圆(如图8e所示)，同时求出各轮廓段的偏移误差值dem；

(4)无轮廓段满足式dm(cs1,csi)&lt;ωdm，则只有cs8融合成一个新的轮廓段；

(5)用新的轮廓段拟合椭圆(如图8f所示)，同时孢子数量num+1，并且删除轮廓段集中的cs8；

(6)重复上述步骤，直到所有轮廓段都拟合出正确的椭圆(如图8g所示)。

所述步骤8.4中的距离测量阈值ωdm和偏移误差阈值σdem由30幅夏孢子图像样本的训练确定，分别为40和95。

9最后统计所有椭圆的总数即为夏孢子的个数。本发明对图2的分割计数结果如图9所示，从图9可以看出，粘连孢子被很好地分割出来。本发明的小麦条锈病菌夏孢子计数准确率结果如表1所示。应用本发明方法，最低计数准确率为92.7％，最高计数准确率为100％，总平均计数准确率为98.6％，具有较高的孢子计数精度。为在线式小麦条锈病夏孢子监测装备的开发提供技术支持。

表1本发明的小麦条锈病菌夏孢子计数准确率结果

本说明书中所描述的以上内容仅仅是对本发明所作的举例说明。本发明所述技术领域的技术人员可以对所描述的具体实施例做各种修改或补充或采用类似的方式替代，只要不偏离本发明说明说的内容或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。