CN113139934B - 一种水稻谷粒计数方法 - Google Patents

一种水稻谷粒计数方法 Download PDF

Info

Publication number
CN113139934B
CN113139934B CN202110325343.5A CN202110325343A CN113139934B CN 113139934 B CN113139934 B CN 113139934B CN 202110325343 A CN202110325343 A CN 202110325343A CN 113139934 B CN113139934 B CN 113139934B
Authority
CN
China
Prior art keywords
grain
area
point
value
connected region
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110325343.5A
Other languages
English (en)
Other versions
CN113139934A (zh
Inventor
马燕
熊斌斌
朱旻
雍开成
魏鑫
黄学辉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Normal University
Original Assignee
Shanghai Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Normal University filed Critical Shanghai Normal University
Priority to CN202110325343.5A priority Critical patent/CN113139934B/zh
Publication of CN113139934A publication Critical patent/CN113139934A/zh
Application granted granted Critical
Publication of CN113139934B publication Critical patent/CN113139934B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/40—Analysis of texture
    • G06T7/41—Analysis of texture based on statistical description of texture
    • G06T7/44—Analysis of texture based on statistical description of texture using image operators, e.g. filters, edge density metrics or local histograms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10004—Still image; Photographic image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20224—Image subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation
    • G06T2207/30188—Vegetation; Agriculture
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30242—Counting objects in image

Abstract

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

Description

一种水稻谷粒计数方法
技术领域
本发明涉及农业生产技术领域，尤其涉及一种水稻谷粒计数方法。
背景技术
在农业生产和研究中，经常需要对水稻谷粒数量进行统计，以进一步检测谷粒的各种形态参数。人工统计数量是较为传统的计数方式。但人工计数的缺点在于，当计数的时间过长时，人眼在持续性高强度工作后很容易产生视觉疲劳，不仅难以保证计数的准确率，而且计数的效率较低，无法达到快速准确计数的要求。因此，实现谷粒数目的自动统计，提高计数的效率和准确率，具有很强的现实意义。
近年来，随着数字图像处理和模式识别技术在工程、农业、医学等领域的广泛应用，利用数字图像处理技术对谷粒实现自动化提取成为可能，并且其对所检测目标的准确率和效率有了更高的保证。
基于图像处理技术统计谷粒数目，一般需要人工将谷粒以一定间隔摆放，因为如果谷粒间存在粘连，则会导致统计算法出错，最终使得谷粒数目不准确。而这种摆放方式存在2个问题，一是仍旧需要人工参与，二是由于以一定间隔摆放，所以所放置的谷粒数目不能太多，这对于后续的谷粒研究来讲约束过多。
另外，还有一些方法需要专门配置的机器，其中包括多个固定位置的CCD相机，使用专门的设备可以提供计数准确率，但带来的问题是设备价格昂贵，并需要占据较多场地，这对用户来说也带来了诸多不便。另外，这些设备结构复杂，有一定的操作流程，对操作人员的要求较高。同时，在使用前，也需要人工将谷粒间较为严重的粘连去除，这也将耗费一定的人力成本。
发明内容
有鉴于现有技术的上述缺陷，本发明所要解决的技术问题是现有的谷粒计数通常为人工计数，使用图像处理统计谷粒计数时，需要人工将谷粒间隔摆放，同样耗费大量的人力，并且当出现谷粒黏连时，计算准确率不高；计数准确的图像处理方法则需要多个价格昂贵的设备辅助，操作要求比较高。因此，本发明提供了一种水稻谷粒计数方法，不需要精细昂贵的设备，成本低廉，操作简单，易于实现，拍摄谷粒图像时，允许谷粒之间有一定程度的粘连，无需人工干预，节省人力成本。
为实现上述目的，本发明提供了一种水稻谷粒计数方法，包括以下步骤：
读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像；
对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图；
对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线；
根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图；
根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数。
进一步地，读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像，具体包括以下步骤：
采集谷粒图像，并将谷粒图像导入到计算机中；
读取采集到的谷粒图像，将谷粒图像中各像素的属性表示为各像素的红色值-蓝色值和绿色值-蓝色值，当差值超出阈值，则将差值调整为30；
对谷粒图像中所有像素点利用kmeans聚类算法分为2个簇，确定谷粒对应的簇，得到谷粒二值图像。
进一步地，对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图，具体包括以下步骤：
统计谷粒二值图像中各谷粒连通区域面积，删除面积小于80的连通区域，并对连通区域中的孔洞进行填补；
统计各连通区域的灰度直方图。
进一步地，对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线，具体包括以下步骤：
利用高斯滤波对各灰度直方图进行卷积运算，得到平滑后的灰度直方图；
对平滑后的灰度直方图进行五次曲线拟合，得到拟合曲线。
进一步地，根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图，具体包括以下步骤：
通过寻找拟合曲线的最大值和第一个谷点，确定阈值thr，删除连通区域中灰度值小于阈值thr的像素，得到更为精确的谷粒连通区域；
对各谷粒连通区域按面积值从小到大进行排序，得到面积分布曲线图。
进一步地，根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数，具体包括以下步骤：
从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将第一个面积变化相对平稳的点作为起始端点i0；
从面积分布曲线的起始端点i0出发，找到第一个面积变化相对剧烈的点，则将第一个面积变化相对剧烈的点作为终点e0；
计算起始端点i0和终点e0之间面积的平均值s，并将平均值作为单个谷粒的面积值；
将谷粒二值图像中各连通区域的面积除以单个谷粒的面积值s，当商小于1且大于0.7，则连通区域的谷粒数目等于1，否则连通区域的谷粒数目等于商的整数部分；
对各个连通区域所包含的谷粒数作累加，从而得到图像中所包含的谷粒总数。
进一步地，从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将第一个面积变化相对平稳的点作为起始端点i0，具体包括以下步骤：
从面积分布曲线的起点出发，将起点对应的面积记作s(i)，位于起点右边连续三个点的面积分别记作s(i+1)，s(i+2)和s(i+3)，判断该起点是否同时满足以下4个条件：
条件1：
条件2：
条件3：
条件4：
如果该起点不能同时满足上述4个条件，则进一步判断该起点右边的一个点是否同时满足上述4个条件；这样一直进行下去，直到找到第一个同时满足上述4个条件的点，这个同时满足上述4个条件的点的面积变化相对平稳，可以被作为起始端点i0。
进一步地，从面积分布曲线的起始端点i0出发，找到第一个面积变化相对剧烈的点，则将第一个面积变化相对剧烈的点作为终点e0，其特征在于，从起始端点i0出发，将起始端点的面积记作s(i)，另将起始端点右边第3个点的面积记作s(i+3)，判断起始端点是否满足以下条件：
如果起始端点不能满足上述条件，则进一步判断该起始端点右边的点是否满足上述条件；这样一直进行下去，直到找到第一个满足上述条件的点，这个满足上述条件的点的面积变化相对剧烈，可以被作为终点e0。
进一步地，计算起始端点i0和终点e0之间面积的平均值s，并将平均值作为单个谷粒的面积值，具体包括以下步骤：
将起始端点i0和终点e0之间n个点的面积分别记作s(1),s(2),…,s(n)，计算面积的平均值s：
将s值作为单个谷粒的面积值。
技术效果
拍摄时仅需准备黑色塑料盆，不反光的黑色背景纸和手机，即可完成谷粒图像的拍摄，不需要搭建专门的摄像设备，对于照明环境无任何特殊要求，具有成本低廉，操作简便，易于实现的优点。
在拍摄谷粒图像时，允许谷粒间存在一定程度的粘连，不需要在拍摄前人工对谷粒的摆放方式进行干预，从而节省了人力成本。
与现有的基于图像的水稻谷粒计数方法相比，本发明不要求分割出所有单个谷粒，通过面积分布曲线统计得到单个谷粒的最小，最大和平均面积值，并利用单个谷粒面积的平均面积值去除各连通区域面积，得到总的谷粒数目，很好地解决了粘连谷粒的计数难题，进一步提高了谷粒计数的准确率，同时节省了大量的人力成本。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
图1是本发明的一个较佳实施例的一种水稻谷粒计数方法的流程示意图；
图2是本发明的一个较佳实施例采集的谷粒图像；
图3是本发明的一个较佳实施例像素点kmeans聚类后得到的二值图像；
图4是本发明的一个较佳实施例谷粒图像某个连通区域的灰度直方图的趋势图；
图5是本发明的一个较佳实施例某个连通区域的灰度直方图经过高斯滤波平滑后的趋势图；
图6是本发明的一个较佳实施例对高斯滤波后的灰度直方图进行五次曲线拟合后的趋势图；
图7是本发明的一个较佳实施例各谷粒连通区域按照从小到大顺序排序后得到的曲线图。
具体实施方式
为了使本发明所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
以下描述中，为了说明而不是为了限定，提出了诸如特定内部程序、技术之类的具体细节，以便透彻理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。
本发明提供了一种水稻谷粒计数方法，包括以下步骤：
读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像；
对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图；
对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线；
根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图；
根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数。
进一步地，读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像，具体包括以下步骤：
采集谷粒图像，并将谷粒图像导入到计算机中；
读取采集到的谷粒图像，将谷粒图像中各像素的属性表示为各像素的红色值-蓝色值和绿色值-蓝色值，当差值超出阈值，则将该差值调整为30；
对谷粒图像中所有像素点利用kmeans聚类算法分为2个簇，确定谷粒对应的簇，得到谷粒二值图像。
进一步地，对谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计谷粒连通区域的灰度直方图，具体包括以下步骤：
统计谷粒二值图像中各谷粒连通区域面积，删除面积小于80的连通区域，并对连通区域中的孔洞进行填补；
统计各连通区域的灰度直方图。
进一步地，对灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线，具体包括以下步骤：
利用高斯滤波对各灰度直方图进行卷积运算，得到平滑后的灰度直方图；
对平滑后的灰度直方图进行五次曲线拟合，得到拟合曲线。
进一步地，根据获得的拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图，具体包括以下步骤：
通过寻找拟合曲线的最大值和第一个谷点，确定阈值thr，删除连通区域中灰度值小于阈值thr的像素，得到更为精确的谷粒连通区域；
对各谷粒连通区域按面积值从小到大进行排序，得到面积分布曲线图。
进一步地，根据面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数，具体包括以下步骤：
从面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将第一个面积变化相对平稳的点作为起始端点i0；
从面积分布曲线的起始端点i0出发，找到第一个面积变化相对剧烈的点，则将第一个面积变化相对剧烈的点作为终点e0；
计算起始端点i0和终点e0之间面积的平均值s，并将平均值作为单个谷粒的面积值；
将谷粒二值图像中各连通区域的面积除以单个谷粒的面积值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＝ax5+bx4+cx3+dx2+ex+f；
其次，利用平滑后的灰度直方图计算五次拟合曲线中的6个参数，包括a，b，c，d，e和f，得到拟合曲线。
在利用多项式函数拟合数据点时，多项式函数的形式如下：
y(x,W)＝w0+w1x+…+wmxm
令那么多项式函数可以化为线性代数的形式：
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个条件的点的面积变化相对平稳，可以被作为起始端点i0。
S11：从面积分布曲线的起始端点i0出发，找到第一个面积变化相对剧烈的点，则将该点作为终点e0；其包括如下步骤：
从S10得到的起始端点i0出发，将该点的面积记作s(i0)，另将该点右边第3个点的面积记作s(i0+3)，一旦该点满足以下条件：
如果起始端点不能满足上述条件，则进一步判断该起始端点右边的点是否满足上述条件；这样一直进行下去，直到找到第一个满足上述条件的点，这个满足上述条件的点的面积变化相对剧烈，可以被作为终点e0。
S12：计算起始端点i0和终点e0之间面积的平均值s，并将该值作为单个谷粒的面积值；其包括如下步骤：
将起始端点i0和终点e0之间n个点的面积分别记作s(1),s(2),…,s(n)，计算面积的平均值s：
将s值作为单个谷粒的面积值。
S13：将二值图像中各连通区域的面积除以单个谷粒的面积值s，当商小于1且大于0.7，则该连通区域的谷粒数目等于1，否则该连通区域的谷粒数目等于商的整数部分。
S14：对各个连通区域所包含的谷粒数作累加，从而得到图像中所包含的谷粒总数。
表1为对20幅不同谷粒图像使用本发明中的方法和装置进行测量得出的结果，这20幅图像包括不同尺度，不同光照，不同粘连程度的谷粒图像。从表1可见，在20幅谷粒图像中，有15幅得到的准确率为100％，准确率最低的是96.91％，总之，本发明提供的方法得到了较高的准确率。
表1 20幅不同谷粒图像谷粒数测量结果
/>
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (1)

1.一种水稻谷粒计数方法，其特征在于，包括以下步骤：
S1、读取谷粒图像，对谷粒图像进行预处理，得到谷粒二值图像；具体包括：
S11、采集谷粒图像，并将所述谷粒图像导入到计算机中；
S12、读取采集到的所述谷粒图像，将所述谷粒图像中各像素的属性表示为所述各像素的红色值-蓝色值和绿色值-蓝色值，当差值超出阈值，则将所述差值调整为30；
S13、对所述谷粒图像中所有像素点，根据步骤S12得到的所有像素的两个属性，利用kmeans聚类算法分为2个簇，确定谷粒对应的簇，得到谷粒二值图像；
S2、对所述谷粒二值图像统计谷粒连通区域，并对谷粒连通区域进行删除和填补，统计所述谷粒连通区域的灰度直方图；具体包括：
S21、统计所述谷粒二值图像中各谷粒连通区域面积，删除面积小于80的连通区域，并对连通区域中的孔洞进行填补；
S22、统计各连通区域的灰度直方图，其特征在于，依次计算各连通区域中的各个灰度值对应的像素数目，从而得到整个连通区域的灰度直方图；
S3、对所述灰度直方图进行卷积运算，然后进行曲线拟合，得到拟合曲线；具体包括：
S31、利用高斯滤波对各灰度直方图进行卷积运算，得到平滑后的灰度直方图；其包括如下步骤：
首先，定义长度为100，方差为7.5的高斯低通滤波器；
其次，将高斯低通滤波器与步骤S2中得到的各连通区域的灰度直方图进行卷积运算，得到平滑后的灰度直方图；
S32、对平滑后的灰度直方图进行五次曲线拟合，得到拟合曲线；其包括如下步骤：
首先，定义五次拟合曲线y＝ax5+bx4+cx3+dx2+ex+f；
其次，利用平滑后的灰度直方图计算五次拟合曲线中的6个参数，包括a，b，c，d，e和f，得到拟合曲线；
S4、根据获得的所述拟合曲线，确定阈值，计算谷粒连通区域的面积，并得到面积分布曲线图；具体包括：
S41、通过寻找所述拟合曲线的最大值和第一个谷点，确定阈值thr，删除连通区域中灰度值小于阈值thr的像素，得到更为精确的谷粒连通区域；其包括如下步骤：
首先，根据步骤S3得到的五次拟合曲线，在灰度值为150～200之间寻找最大值；
其次，从最大值出发，沿拟合曲线左下角方向，判断曲线上各点与其前、后相邻两点之间的灰度值之差，当遇到点A，其灰度值比其前、后相邻两点的灰度值都要小，则将A点作为第一个谷点；
然后，将第一个谷点对应的灰度值作为阈值thr，删除连通区域中灰度值小于阈值thr的像素，即，将谷粒二值图像中该像素的数值调整为0，从而得到更为精确的连通区域；
S42、对各谷粒连通区域按面积值从小到大进行排序，得到面积分布曲线图；其包括如下步骤：
首先，对于所述更为精确的各连通区域，统计像素值等于1的各连通区域中所包含的像素数目，并将像素数目作为连通区域的面积；
其次，对各连通区域按其面积值从小到大进行排序，以横坐标为各连通区域的排序顺序，纵坐标为各连通区域的面积，得到面积分布曲线图；
S5、根据所述面积分布曲线图，计算单个谷粒的面积值，通过谷粒连通区域的面积除以单个谷粒的面积，得到图像中所包含的谷粒总数；具体包括：
S51、从所述面积分布曲线的起点出发，找到第一个面积变化相对平稳的点，将所述第一个面积变化相对平稳的点作为起始端点i0；具体包括以下步骤：
从所述面积分布曲线的起点出发，将所述起点对应的面积记作s(i)，位于所述起点右边连续三个点的面积分别记作s(i+1)，s(i+2)和s(i+3)，判断所述起点是否同时满足以下4个条件：
条件1：
条件2：
条件3：
条件4：
如果所述起点不能同时满足上述4个条件，则进一步判断所述起点右边的点是否同时满足上述4个条件；直到找到第一个同时满足上述4个条件的点，所述同时满足上述4个条件的点的面积变化相对平稳，可以被作为起始端点i0；
S52、从所述面积分布曲线的所述起始端点i0出发，找到第一个面积变化相对剧烈的点，则将所述第一个面积变化相对剧烈的点作为终点e0；其特征在于，从所述起始端点i0出发，将所述起始端点的面积记作s(i0)，另将所述起始端点右边第3个点的面积记作s(i0+3)，判断起始端点是否满足以下条件：
如果起始端点不能满足上述条件，则进一步判断该起始端点右边的点是否满足上述条件；这样一直进行下去，直到找到第一个满足上述条件的点，这个满足上述条件的点的面积变化相对剧烈，可以被作为终点e0；
S53、计算所述起始端点i0和所述终点e0之间面积的平均值s，并将所述平均值作为单个谷粒的面积值，具体包括以下步骤：
将所述起始端点i0和所述终点e0之间n个点的面积分别记作s(1),s(2),…,s(n)，计算面积的平均值s：
将s值作为单个谷粒的面积值；
S54、将所述谷粒二值图像中各连通区域的面积除以单个谷粒的面积值s，当商小于1且大于0.7，则所述连通区域的谷粒数目等于1，否则所述连通区域的谷粒数目等于商的整数部分；
S55、对各个连通区域所包含的谷粒数作累加，从而得到图像中所包含的谷粒总数。
CN202110325343.5A 2021-03-26 2021-03-26 一种水稻谷粒计数方法 Active CN113139934B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110325343.5A CN113139934B (zh) 2021-03-26 2021-03-26 一种水稻谷粒计数方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110325343.5A CN113139934B (zh) 2021-03-26 2021-03-26 一种水稻谷粒计数方法

Publications (2)

Publication Number Publication Date
CN113139934A CN113139934A (zh) 2021-07-20
CN113139934B true CN113139934B (zh) 2024-04-30

Family

ID=76810530

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110325343.5A Active CN113139934B (zh) 2021-03-26 2021-03-26 一种水稻谷粒计数方法

Country Status (1)

Country Link
CN (1) CN113139934B (zh)

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08145871A (ja) * 1994-09-19 1996-06-07 Hitachi Ltd 粒子画像の領域分割方法及び装置
JP2005215712A (ja) * 2002-01-21 2005-08-11 Nisca Corp 個体数カウント用の撮像装置及び制御手段
CN101281112A (zh) * 2008-04-30 2008-10-08 浙江理工大学 一种对网状粘连稻米的图像式自动分析方法
CN101339118A (zh) * 2008-08-08 2009-01-07 华中科技大学 谷粒参数自动测量装置及方法
CN101458204A (zh) * 2008-12-24 2009-06-17 华中科技大学 谷物实粒数的自动测量装置及方法
CN102750584A (zh) * 2012-04-18 2012-10-24 中国农业大学 玉米粒在穗计数方法
CN104966085A (zh) * 2015-06-16 2015-10-07 北京师范大学 一种基于多显著特征融合的遥感图像感兴趣区域检测方法
CN108388853A (zh) * 2018-02-09 2018-08-10 重庆东渝中能实业有限公司 针对白细胞与血小板共存全息图的分步重建与计数方法
CN108765448A (zh) * 2018-05-28 2018-11-06 青岛大学 一种基于改进tv-l1模型的虾苗计数分析方法
CN108961208A (zh) * 2018-05-21 2018-12-07 江苏康尚生物医疗科技有限公司 一种聚集白细胞分割计数系统及方法
CN109636824A (zh) * 2018-12-20 2019-04-16 巢湖学院 一种基于图像识别技术的多目标计数方法
CN109843034A (zh) * 2016-10-19 2019-06-04 巴斯夫农化商标有限公司 用于麦田的产量预测

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08145871A (ja) * 1994-09-19 1996-06-07 Hitachi Ltd 粒子画像の領域分割方法及び装置
JP2005215712A (ja) * 2002-01-21 2005-08-11 Nisca Corp 個体数カウント用の撮像装置及び制御手段
CN101281112A (zh) * 2008-04-30 2008-10-08 浙江理工大学 一种对网状粘连稻米的图像式自动分析方法
CN101339118A (zh) * 2008-08-08 2009-01-07 华中科技大学 谷粒参数自动测量装置及方法
CN101458204A (zh) * 2008-12-24 2009-06-17 华中科技大学 谷物实粒数的自动测量装置及方法
CN102750584A (zh) * 2012-04-18 2012-10-24 中国农业大学 玉米粒在穗计数方法
CN104966085A (zh) * 2015-06-16 2015-10-07 北京师范大学 一种基于多显著特征融合的遥感图像感兴趣区域检测方法
CN109843034A (zh) * 2016-10-19 2019-06-04 巴斯夫农化商标有限公司 用于麦田的产量预测
CN108388853A (zh) * 2018-02-09 2018-08-10 重庆东渝中能实业有限公司 针对白细胞与血小板共存全息图的分步重建与计数方法
CN108961208A (zh) * 2018-05-21 2018-12-07 江苏康尚生物医疗科技有限公司 一种聚集白细胞分割计数系统及方法
CN108765448A (zh) * 2018-05-28 2018-11-06 青岛大学 一种基于改进tv-l1模型的虾苗计数分析方法
CN109636824A (zh) * 2018-12-20 2019-04-16 巢湖学院 一种基于图像识别技术的多目标计数方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Matlab图像处理在水稻谷粒计数中的应用;焦雁翔;唐玉琴;黄成志;黄仁军;;南方农业(第04期);全文 *
Ravi Kiran Boggavarapu ; Pushpendra Kumar Pateriya."Efficient method for counting and tracking of moving objects through region segmentation".《2016 IEEE International Conference on Computational Intelligence and Computing Research (ICCIC)》.2016,全文. *
一种基于Matlab的谷物颗粒计数方法;赵欣;贾晓剑;;河南科技学院学报(自然科学版)(第05期);全文 *
基于细胞面积估计的蓝藻细胞计数;胡洋洋;王鑫;;江南大学学报(自然科学版)(第05期);全文 *

Also Published As

Publication number Publication date
CN113139934A (zh) 2021-07-20

Similar Documents

Publication Publication Date Title
CN109636784B (zh) 基于最大邻域和超像素分割的图像显著性目标检测方法
Guijarro et al. Discrete wavelets transform for improving greenness image segmentation in agricultural images
CN110120042B (zh) 一种基于slic超像素和自动阈值分割的农作物图像病虫害区域提取方法
CN109447945B (zh) 基于机器视觉和图形处理的小麦基本苗快速计数方法
CN110415208B (zh) 一种自适应目标检测方法及其装置、设备、存储介质
CN106650580B (zh) 基于图像处理的货架快速清点方法
CN109685045A (zh) 一种运动目标视频跟踪方法及系统
CN112712518B (zh) 鱼类计数方法、装置、电子设备及存储介质
CN111695373B (zh) 斑马线的定位方法、系统、介质及设备
CN108829711B (zh) 一种基于多特征融合的图像检索方法
CN111382658B (zh) 一种基于图像灰度梯度一致性的自然环境下道路交通标志检测方法
CN116664559A (zh) 基于机器视觉的内存条损伤快速检测方法
CN112233076B (zh) 基于红色圆标靶图像处理的结构振动位移测量方法及装置
CN111709305B (zh) 一种基于局部图像块的人脸年龄识别方法
CN111476804A (zh) 托辊图像高效分割方法、装置、设备及存储介质
CN110287847A (zh) 基于Alexnet-CLbpSurf多特征融合的车辆分级检索方法
Paschos et al. A color texture based visual monitoring system for automated surveillance
CN114862902A (zh) 一种基于四叉树的光照自适应orb特征提取和匹配方法
CN106529441A (zh) 基于模糊边界分片的深度动作图人体行为识别方法
CN112053371A (zh) 一种遥感图像中的水体提取方法和装置
CN113139934B (zh) 一种水稻谷粒计数方法
CN111369497B (zh) 一种行走式树上果实连续计数方法及装置
CN111126303B (zh) 一种面向智能停车的多车位检测方法
CN109544614B (zh) 一种基于图像低频信息相似度的匹配图像对识别的方法
CN111105394B (zh) 一种检测发光小球特征信息的方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant