一种水稻谷粒计数方法

摘要

本发明公开了一种水稻谷粒计数方法，包括以下步骤：读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像；对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图；对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线；根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图；根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数。本发明的一种水稻谷粒计数方法，不需要精细昂贵的设备，成本低廉，操作简单，易于实现，拍摄谷粒图像时，允许谷粒之间有一定程度的粘连，无需人工干预，节省人力成本。

说明书

技术领域

本发明涉及农业生产技术领域，尤其涉及一种水稻谷粒计数方法。

背景技术

在农业生产和研究中，经常需要对水稻谷粒数量进行统计，以进一步检测谷粒的各种形态参数。人工统计数量是较为传统的计数方式。但人工计数的缺点在于，当计数的时间过长时，人眼在持续性高强度工作后很容易产生视觉疲劳，不仅难以保证计数的准确率，而且计数的效率较低，无法达到快速准确计数的要求。因此，实现谷粒数目的自动统计，提高计数的效率和准确率，具有很强的现实意义。

近年来，随着数字图像处理和模式识别技术在工程、农业、医学等领域的广泛应用，利用数字图像处理技术对谷粒实现自动化提取成为可能，并且其对所检测目标的准确率和效率有了更高的保证。

基于图像处理技术统计谷粒数目，一般需要人工将谷粒以一定间隔摆放，因为如果谷粒间存在粘连，则会导致统计算法出错，最终使得谷粒数目不准确。而这种摆放方式存在2个问题，一是仍旧需要人工参与，二是由于以一定间隔摆放，所以所放置的谷粒数目不能太多，这对于后续的谷粒研究来讲约束过多。

另外，还有一些方法需要专门配置的机器，其中包括多个固定位置的CCD相机，使用专门的设备可以提供计数准确率，但带来的问题是设备价格昂贵，并需要占据较多场地，这对用户来说也带来了诸多不便。另外，这些设备结构复杂，有一定的操作流程，对操作人员的要求较高。同时，在使用前，也需要人工将谷粒间较为严重的粘连去除，这也将耗费一定的人力成本。

发明内容

有鉴于现有技术的上述缺陷，本发明所要解决的技术问题是现有的谷粒计数通常为人工计数，使用图像处理统计谷粒计数时，需要人工将谷粒间隔摆放，同样耗费大量的人力，并且当出现谷粒黏连时，计算准确率不高；计数准确的图像处理方法则需要多个价格昂贵的设备辅助，操作要求比较高。因此，本发明提供了一种水稻谷粒计数方法，不需要精细昂贵的设备，成本低廉，操作简单，易于实现，拍摄谷粒图像时，允许谷粒之间有一定程度的粘连，无需人工干预，节省人力成本。

为实现上述目的，本发明提供了一种水稻谷粒计数方法，包括以下步骤：

读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像；

对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图；

对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线；

根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图；

根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数。

进一步地，读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像，具体包括以下步骤：

采集谷粒图像，并将谷粒图像导入到计算机中；

读取采集到的谷粒图像，将谷粒图像中各像素的属性表示为各像素的红色值-蓝色值和绿色值-蓝色值，当差值超出阈值，则将差值调整为30；

对谷粒图像中所有像素点利用kmeans聚类算法分为2个簇，确定谷粒对应的簇，得到谷粒二值图像。

进一步地，对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图，具体包括以下步骤：

统计谷粒二值图像中各谷粒连通区域面积，删除面积小于80的连通区域，并对连通区域中的孔洞进行填补；

统计各连通区域的灰度直方图。

进一步地，对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线，具体包括以下步骤：

利用高斯滤波对各灰度直方图进行卷积运算，得到平滑后的灰度直方图；

对平滑后的灰度直方图进行五次曲线拟合，得到拟合曲线。

进一步地，根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图，具体包括以下步骤：

通过寻找拟合曲线的最大值和第一个谷点，确定阈值thr，删除连通区域中灰度值小于阈值thr的像素，得到更为精确的谷粒连通区域；

对各谷粒连通区域按面积值从小到大进行排序，得到面积分布曲线图。

进一步地，根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数，具体包括以下步骤：

从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将第一个面积变化相对平稳的点作为起始端点i

从面积分布曲线的起始端点i

计算起始端点i

将谷粒二值图像中各连通区域的面积除以单个谷粒的面积值s，当商小于1且大于0.7，则连通区域的谷粒数目等于1，否则连通区域的谷粒数目等于商的整数部分；

对各个连通区域所包含的谷粒数作累加，从而得到图像中所包含的谷粒总数。

进一步地，从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将第一个面积变化相对平稳的点作为起始端点i

从面积分布曲线的起点出发，将起点对应的面积记作s(i)，位于起点右边连续三个点的面积分别记作s(i+1)，s(i+2)和s(i+3)，判断该起点是否同时满足以下4个条件：

条件1：

条件2：

条件3：

条件4：

如果该起点不能同时满足上述4个条件，则进一步判断该起点右边的一个点是否同时满足上述4个条件；这样一直进行下去，直到找到第一个同时满足上述4个条件的点，这个同时满足上述4个条件的点的面积变化相对平稳，可以被作为起始端点i

进一步地，从面积分布曲线的起始端点i

如果起始端点不能满足上述条件，则进一步判断该起始端点右边的点是否满足上述条件；这样一直进行下去，直到找到第一个满足上述条件的点，这个满足上述条件的点的面积变化相对剧烈，可以被作为终点e

进一步地，计算起始端点i

将起始端点i

将s值作为单个谷粒的面积值。

技术效果

拍摄时仅需准备黑色塑料盆，不反光的黑色背景纸和手机，即可完成谷粒图像的拍摄，不需要搭建专门的摄像设备，对于照明环境无任何特殊要求，具有成本低廉，操作简便，易于实现的优点。

在拍摄谷粒图像时，允许谷粒间存在一定程度的粘连，不需要在拍摄前人工对谷粒的摆放方式进行干预，从而节省了人力成本。

与现有的基于图像的水稻谷粒计数方法相比，本发明不要求分割出所有单个谷粒，通过面积分布曲线统计得到单个谷粒的最小，最大和平均面积值，并利用单个谷粒面积的平均面积值去除各连通区域面积，得到总的谷粒数目，很好地解决了粘连谷粒的计数难题，进一步提高了谷粒计数的准确率，同时节省了大量的人力成本。

以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。

附图说明

图1是本发明的一个较佳实施例的一种水稻谷粒计数方法的流程示意图；

图2是本发明的一个较佳实施例采集的谷粒图像；

图3是本发明的一个较佳实施例像素点kmeans聚类后得到的二值图像；

图4是本发明的一个较佳实施例谷粒图像某个连通区域的灰度直方图的趋势图；

图5是本发明的一个较佳实施例某个连通区域的灰度直方图经过高斯滤波平滑后的趋势图；

图6是本发明的一个较佳实施例对高斯滤波后的灰度直方图进行五次曲线拟合后的趋势图；

图7是本发明的一个较佳实施例各谷粒连通区域按照从小到大顺序排序后得到的曲线图。

具体实施方式

为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

以下描述中，为了说明而不是为了限定，提出了诸如特定内部程序、技术之类的具体细节，以便透彻理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。

本发明提供了一种水稻谷粒计数方法，包括以下步骤：

读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像；

对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图；

对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线；

根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图；

根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数。

进一步地，读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像，具体包括以下步骤：

采集谷粒图像，并将谷粒图像导入到计算机中；

读取采集到的谷粒图像，将谷粒图像中各像素的属性表示为各像素的红色值-蓝色值和绿色值-蓝色值，当差值超出阈值，则将该差值调整为30；

对谷粒图像中所有像素点利用kmeans聚类算法分为2个簇，确定谷粒对应的簇，得到谷粒二值图像。

进一步地，对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图，具体包括以下步骤：

统计谷粒二值图像中各谷粒连通区域面积，删除面积小于80的连通区域，并对连通区域中的孔洞进行填补；

统计各连通区域的灰度直方图。

进一步地，对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线，具体包括以下步骤：

利用高斯滤波对各灰度直方图进行卷积运算，得到平滑后的灰度直方图；

对平滑后的灰度直方图进行五次曲线拟合，得到拟合曲线。

进一步地，根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图，具体包括以下步骤：

通过寻找拟合曲线的最大值和第一个谷点，确定阈值thr，删除连通区域中灰度值小于阈值thr的像素，得到更为精确的谷粒连通区域；

对各谷粒连通区域按面积值从小到大进行排序，得到面积分布曲线图。

进一步地，根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数，具体包括以下步骤：

从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将第一个面积变化相对平稳的点作为起始端点i

从面积分布曲线的起始端点i

计算起始端点i

将谷粒二值图像中各连通区域的面积除以单个谷粒的面积值s，当商小于1且大于0.7，则连通区域的谷粒数目等于1，否则连通区域的谷粒数目等于商的整数部分；

对各个连通区域所包含的谷粒数作累加，从而得到图像中所包含的谷粒总数。

本发明实施例以粘连谷粒为例说明本发明的发明。

如图1所示的实施例中，一种粘连谷粒的自动计数方法，包括如下步骤：

S1：采集谷粒彩色图像，并将其导入到计算机中；其包括如下步骤：

首先，将谷粒置于黑色塑料盆内；

接着，将黑色塑料盆置于不反光的黑色背景纸上；

然后，利用手机拍摄谷粒，得到谷粒彩色图像，并将其导入计算机中。

根据上述方法采集到图2所示的谷粒图像。

S2：读取谷粒彩色图像，将图像中各像素的属性表示为该像素的红色值-蓝色值和绿色值-蓝色值，如差值超出阈值，则将该差值调整为30；其包括如下步骤：

首先，读取谷粒彩色图像I，I是大小为M×N×3的三维矩阵；

其中，M和N分别表示该谷粒彩色图像包括M行和N列，I(i,j,1)，I(i,j,2)和I(i,j,3)分别表示位于图像第i行，第j列的像素的红色值，绿色值和蓝色值；

接着，将图像中各像素的属性表示为该像素的红色值-蓝色值和绿色值-蓝色值，得到I(i,j,1)-I(i,j,3)，I(i,j,2)-I(i,j,3)两个属性；

然后，确定阈值为60，当I(i,j,1)-I(i,j,3)大于60，则将该差值调整为30；

同样，当I(i,j,2)-I(i,j,3)大于60，则将该差值调整为30。

S3：对图像中所有像素点利用kmeans聚类算法分为2个簇，确定谷粒对应的簇，得到谷粒二值图像；其包括如下步骤：

首先，确定2个初始聚类中心(40,40)和(0,0)；

其次，根据步骤S2得到的所有像素的两个属性，利用kmeans聚类算法将所有的像素点分为2个簇；

接着，分别计算各像素灰度值I_gray(i,j)：

其中，I_gray(i,j)＝(I(i,j,1)+I(i,j,2)+I(i,j,3))/3；

然后，分别计算2个簇中所有像素灰度值的均值，灰度值均值较大的簇为谷粒对应的簇，得到谷粒二值图像IB，IB的大小为M×N；

其中，像素值1表示谷粒，像素值0表示非谷粒。

该步骤中kmeans聚类算法是对二维数据对象进行聚类，输入为初始聚类中心(40,40)和(0,0)，聚类个数2以及n个二维数据对象，这里的二维数据对象即为图像中的像素，n为图片中像素点的总个数，将n个数据对象划分为2个簇，要求同一簇中的数据对象相似度较高，而不同簇的数据对象相似度较小。聚类相似度是利用聚类中心来计算的，第一次进行聚类操作时，由于指定了两个初始聚类中心，根据每个二维数据对象到两个初始聚类中心的距离，将各个数据对象归到与其距离较小的聚类中心所在的簇，从而得到两个簇，接着分别计算每个簇中的各个数据对象同一维度值的平均值，从而得到两个新的聚类中心，继续进行上述聚类操作直到聚类中心不再改变，此时得到的两个簇即为kmeans聚类的最终结果。图3为一幅经过上述kmeans聚类算法得到的二值图像。

S4：统计谷粒二值图像中各谷粒连通区域面积，删除面积小于80的连通区域，并对连通区域中的孔洞进行填补；其包括如下步骤：

首先，对于步骤S3得到的谷粒二值图像IB，统计像素值等于1的各连通区域中所包含的像素数目，并将像素数目作为连通区域的面积；

其次，删除面积小于80的连通区域；

然后，对连通区域中的孔洞进行填补。

S5：统计各连通区域的灰度直方图；其包括如下步骤：

对于S4中得到的各连通区域，根据连通区域中各像素的灰度值I_gray(i,j)，分别统计灰度值落在0，1，…，255上面的数目，得到各连通区域的灰度直方图。

例如，编号为2的连通区域内各像素的灰度值用如下大小为22*13的矩阵表示为：

灰度值为45的像素数目为2，那么该连通区域的灰度直方图中横坐标为45对应的纵坐标高度为2，依次计算该连通区域中的各个灰度值对应的像素数目，从而得到整个连通区域的灰度直方图。计算得到的上述编号为2的连通区域的灰度直方图用大小为1*256的矩阵表示为：

图4为上述编号为2的连通区域的灰度直方图得到的趋势线，可以清晰地看到趋势线上存在许多的峰值点和谷点，需要做进一步的平滑处理。

S6：利用高斯滤波对各灰度直方图进行卷积运算，得到平滑后的灰度直方图；其包括如下步骤：

首先，定义长度为100，方差为7.5的高斯低通滤波器；

其次，将高斯低通滤波器与S5中得到的各连通区域的灰度直方图进行卷积运算，得到平滑后的灰度直方图。

高斯滤波使用的是高斯函数，即满足正态分布的概率密度函数：

本发明中生成的高斯模板就是从上述公式得来的，例如，要生成一个大小为3，标准差为1的模板，则只需要代入公式计算(此处均值μ为0)：f(-1)、f(0)、f(1)就可以得到模板的值。本发明中，高斯滤波器参数设置为大小为100，标准差为7.5，代入公式f(x)可以得到：

从而可以得到以下1*100的高斯滤波器：

将S5得到的编号为2的连通区域的直方图矩阵histgram与大小为1*100的高斯滤波器h做卷积运算，得到新的灰度直方图矩阵histgram：

其中，

经过卷积之后，得到

卷积运算的方法如下：

例如，有一个序列1,2,3,4,5使用滤波器模板[1,2,1]对序列中的3进行卷积操作，那么计算过程为(2×1+3×2+4×1)÷(1+2+1)，得到对3卷积运算后的结果3，如果是对开头的1和结尾的5进行卷积运算，当模板与序列对齐时超出了序列的范围，那么需要为序列不对齐的位置补0，再进行卷积运算。

图5为编号为2的连通区域的灰度直方图与高斯低通滤波器进行卷积运算后得到的灰度直方图的趋势线，可以看到其相较于图4更加的光滑，但仍有一些波折，需要做更进一步的拟合处理。

S7：对平滑后的灰度直方图进行五次曲线拟合，得到拟合曲线；其包括如下步骤：

首先，定义五次拟合曲线y＝ax

其次，利用平滑后的灰度直方图计算五次拟合曲线中的6个参数，包括a，b，c，d，e和f，得到拟合曲线。

在利用多项式函数拟合数据点时，多项式函数的形式如下：

y(x,W)＝w

令

y(x,W)＝XW

根据这个原理，将n个灰度直方图坐标点带入五次曲线方程，写为矩阵形式为：

利用最小二乘法计算五次曲线的6个系数a,b,c,d,e和f，从而拟合五次曲线。

图6是对图5曲线进行五次曲线拟合后的结果，可以看到其相较于图5更加平滑。

S8：通过寻找拟合曲线的最大值和第一个谷点，确定阈值thr，删除连通区域中灰度值小于阈值thr的像素，得到更为精确的谷粒连通区域；其包括如下步骤：

首先，根据步骤七得到的五次拟合曲线，在灰度值为150～200之间寻找最大值；

其次，从最大值出发，沿拟合曲线左下角方向，判断曲线上各点与其前、后相邻两点之间的灰度值之差，当遇到点A，其灰度值比其前、后相邻两点的灰度值都要小，则将A点作为第一个谷点；

然后，将第一个谷点对应的灰度值作为阈值thr，删除连通区域中灰度值小于阈值thr的像素，即，将谷粒二值图像IB中该像素的数值调整为0，从而得到更为精确的连通区域。

S9：对各谷粒连通区域按其面积值从小到大进行排序，得到面积分布曲线图；其包括如下步骤：

首先，对于步骤S8得到的更为精确的各连通区域，统计像素值等于1的各连通区域中所包含的像素数目，并将像素数目作为连通区域的面积；

其次，对各连通区域按其面积值从小到大进行排序，以横坐标为各连通区域的排序顺序，纵坐标为各连通区域的面积，得到面积分布曲线图。

图7所示为各谷粒连通区域按面积值从小到大进行排序后得到的面积分布曲线图。

S10：从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将该点作为起始端点s0；其包括如下步骤：

从S9得到的面积分布曲线的起点出发，将该点对应的面积记作s(i)，位于该点右边的三个点分别记作s(i+1)，s(i+2)和s(i+3)，一旦该点同时满足以下4个条件：

条件1：

条件2：

条件3：

条件4：

如果该起点不能同时满足上述4个条件，则进一步判断该起点右边的一个点s(i+1)是否同时满足上述4个条件；这样一直进行下去，直到找到第一个同时满足上述4个条件的点，这个同时满足上述4个条件的点的面积变化相对平稳，可以被作为起始端点i

S11：从面积分布曲线的起始端点i

从S10得到的起始端点i

如果起始端点不能满足上述条件，则进一步判断该起始端点右边的点是否满足上述条件；这样一直进行下去，直到找到第一个满足上述条件的点，这个满足上述条件的点的面积变化相对剧烈，可以被作为终点e

S12：计算起始端点i

将起始端点i

将s值作为单个谷粒的面积值。

S13：将二值图像中各连通区域的面积除以单个谷粒的面积值s，当商小于1且大于0.7，则该连通区域的谷粒数目等于1，否则该连通区域的谷粒数目等于商的整数部分。

S14：对各个连通区域所包含的谷粒数作累加，从而得到图像中所包含的谷粒总数。

表1为对20幅不同谷粒图像使用本发明中的方法和装置进行测量得出的结果，这20幅图像包括不同尺度，不同光照，不同粘连程度的谷粒图像。从表1可见，在20幅谷粒图像中，有15幅得到的准确率为100％，准确率最低的是96.91％，总之，本发明提供的方法得到了较高的准确率。

表1 20幅不同谷粒图像谷粒数测量结果

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。