CN103308430B - 一种测量种子千粒重的方法及装置 - Google Patents

Abstract

本发明公开了一种测量种子千粒重的方法，包括：1）将待测试的多粒种子平铺，并采集种子的图像和重量；2）将所述的图像转换为灰度图像，并对灰度图像进行滤波处理；3）对滤波处理后的灰度图像进行二值化处理得到二值化图像；4）提取二值化图像中种子的轮廓线，并依据轮廓线按连通区域划分，选取每个连通区域内的轮廓线上的特征点，匹配特征点对，并在两特征点对间建立分割线，所述的分割线将连通区域划分为多个子区域；5）统计所有连通区域内的子区域数量，得到所述种子的数量；6）根据种子的粒数和重量，计算种子千粒重。本发明还公开了一种测量种子千粒重的装置。本发明装置结构简单、操作方便、测量精确高。

Description

一种测量种子千粒重的方法及装置

技术领域

本发明涉及粮食、油料等作物种子的千粒重测量，尤其涉及一种测量种子千粒重的方法及装置。

背景技术

千粒重就是一千粒（或折合成一千粒）谷物的重量，它是体现种子大小与饱满程度的一项指标，种子千粒重是衡量种子质量和丰产性能的重要指标，是评价品种和经营技术措施的重要标准，也是研究数量遗传的必要数据。一直以来，粮食和油料种子千粒重的测定都采用人工的方法，费工费时，且因烦琐的重复劳动极易使人疲劳而产生过误差，不仅效率低且不准确。

随着图像处理、计算机视觉和模式识别技术的迅速发展，国内外对基于计算机视觉技术的类圆物体的检测方法也进行了一些有益的探索，例如:基于圆周模板匹配方法[宋强，徐科，徐金梧等提出的基于图像处理的棒材自动记数技术[J].钢铁,2004,39(5):34-37.]、腐蚀法[魏瑾,韩斌和张其亮发表的粘连蚕卵分离算法研究[J].江苏科技大学学报:自然科学版,2006,20(1):46-50.]、图像距离法[尤育塞,于慧敏和刘圆圆提出的基于粒度测量的重叠圆形颗粒图像分离方法[J].浙江大学学报:工学版,2005,39(7):962-966.]等。

江苏大学的陈静等[陈静,朱伟兴,苏威；基于机器视觉的种子千粒重测定系统的设计，传感器与微系统，2008.27(9):97-99.]就在此基础上对这方面进行了有益的尝试。该研究主要采用图像处理法，提出一种结合了最大方差二值化和欧氏距离变换的方法来解决千粒重测量中出现的问题。其主要步骤是先由摄像机摄取种子图像，进行信号转换后采集到电脑中并通过分析显示颗粒数。天平对该图像中的种子进行称重，结合所得计数结果得出谷物千粒重。该方法速度快，清洁易携带。但对技术要求较高，价格不菲，在计算千粒重时往往需要不断换算，在一定程度上降低了精度。

重力式种子精选机又叫比重式精选机，是国内种子精选分级的主要设备。其中，5XZW-1.5[李法德,张晓辉,康景峰等.1997.5XZW-1.5型重力精选机-主机部分.农机与食品机械.(2):13～24]是用处较为广泛的机型。该机采用结构简单的正压式主机结构，主要组成部分有提升上料机构、长方形台面以及调节机构等。工作原理是由提升机均匀、连续地向工作台输送种子，在重力、振动以及自上而下的气流作用下，使种子不断下落。山东大学的李光提等（李光提,李汝莘,侯存良,迟淑筠,樊桂菊.基于虚拟仪器技术的重力式种子精选机测控系统的研制.农业工程学报.2004,20(4):99-103）提出将VI（虚拟仪器技术）跟5XZW-1.5结合起来，加入光电传感器与荷重传感器，实现系统自我调整，能实现6千粒重的精确测量。该方法测量精确，结构紧凑，但机型较大，需要很大安装空间，不太利于和相关机械如干燥机的结合。

何有根[何有根.基于图像识别的千粒重仪的研究.西北农林科技大学，硕士学位论文.2007.]进行了基于图像识别的千粒重仪的研究，电磁振动排钟装置对玉米进行排序，然后通过图像处理的方法进行单个玉米种子的分类识别，采用红外光电耦合型传感器和高灵敏度重量传感器对被测对象进行监测与称重，处理速度较慢，有较多需要改进的地方。

王刚进行了基于机器视觉的玉米千粒重快速检测仪的研制，由于该仪器设计时的工作对象为玉米，所以应用面相对较窄。使用摄像头不能连续采样，测定速度慢。机械平铺结构的工作情况还有待改进，单次采集样品量较小，存在籽粒重叠现象，影响千粒重测定准确度。

发明内容

本发明提供了一种测量种子千粒重的方法及装置，实现了种子千粒重的快速测量，提高了种子千粒重的测量精度和自动化水平。

一种测量种子千粒重的方法，包括以下几个步骤：

1）将待测试的多粒种子平铺，并采集种子的图像和重量；

2）将所述的图像转换为灰度图像，并对灰度图像进行滤波处理；

图像在形成、传输、接收和处理的过程中，不可避免地存在着外部干扰和内部干扰，图像的常见噪声主要有加性噪声、乘性噪声、量化噪声和椒盐噪声等。

一般图像噪声的来源有三种：

a）在光电、电磁转换的过程中引入的人为噪声；

b）大气层（磁）暴、闪电、电压等引起的强脉冲性冲击噪声的干扰；

c）自然起伏性噪声，由物理量的不连续或粒子性引起的。

反应在图像上，噪声使原本均匀和连续变化的颜色突然变大或减小，形成一些虚假的物体边缘或轮廓。抑制或消除这些噪声的平滑方法有邻域平均法、低通滤波法、中值滤波法、自适应滤波法等。其中邻域平均法的算法简单，但在降低噪声的同时会使图像模糊；中值滤波是一较好的方法，滤除噪声的同时也可以保护物体边缘。自适应滤波也可保留一定的边缘信息。如果平滑不当，就会使图像本身的细节如边界轮廓、线条等变的模糊不清。为了既平滑掉噪声又尽量保持图像细节，采用中值滤波法对灰度图像进行滤波处理，达到消除噪声的目的。

3）对滤波处理后的灰度图像进行二值化处理得到二值化图像；

所述的步骤3）中，在进行二值化处理之前，采用Ostu算法得到滤波处理后灰度图像的阈值。

需要确定阈值才能进行图像二值化处理，Ostu方法计算简单，稳定有效，是实际应用中经常采用的确定阈值的方法，Ostu算法也称最大类间差法，有时也称之为大津算法，被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。它是按图像的灰度特性，将图像分成背景和前景两部分。背景和前景之间的类间方差越大，说明构成图像的两部分的差别越大，当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此，使类间方差最大的分割意味着错分概率最小。

在步骤3）中，还需对所述的二值化图像进行腐蚀算法处理。

腐蚀是一种消除边界点，使边界向内部收缩的过程。可以用来消除二值图像中小且无意义的噪声像素。

4）提取二值化图像中种子的轮廓线，并依据轮廓线按连通区域划分，选取每个连通区域内的轮廓线上的特征点，匹配特征点对，并在两特征点对间建立分割线，所述的分割线将连通区域划分为多个子区域；

在步骤4）中，对所述的轮廓线平滑处理，计算平滑处理后的轮廓线上各点的曲率，并对由各点的曲率构建的曲率曲线一阶求导，选取求导结果为零的点，并结合求导结果为零的点所对应的曲率中心，得到所述轮廓线上的特征点。

原始轮廓线上的噪音太多，会影响特征点的计算，所有要先对原始轮廓线进行平滑，通过平滑得到了平滑的轮廓线，计算平滑的轮廓线上各点的曲率，将这些点的曲率值依次连接作为一条曲线，计算该曲线的一阶导数，一阶导数为零的点就是特征点。该特征点只是平滑的轮廓线上的特征点，而不是原始轮廓线的特征点，所以要结合特征点对应的曲率中心，两者结合来实现定位，从而确定原始轮廓线上的特征点。

所述的轮廓线包括外轮廓线和内轮廓线，针对轮廓线上的特征点A，与特征点A匹配的特征点B应满足以下几个条件：

a）建立通过特征点A的法线，特征点B与所述法线之间的距离最小；

b）特征点B与特征点A之间的距离小于种子的长度；

c）特征点B与特征点A处于不同的内轮廓线上；

d）特征点A与特征点B之间的距离最短。

5）统计所有连通区域内的子区域数量，得到所述种子的粒数；

在步骤5）中，将所述分割线的像素值设置为零，选取子区域中像素值不为零的像素a，运用八邻域算子对该像素a周围八个像素进行搜索，选择周围八个像素与像素a的像素值相同的像素进行标记，遍历所有子区域中的每个像素，完成对所有子区域的标记，并计算标记的子区域的个数。

针对子区域的标记，陈柏生提出了一种二值图像连通区域标记的新方法[陈柏生.一种二值图像连通区域标记的新方法[J].计算机工程与应用,2006,42(25):46-47.]。

6）根据种子的粒数和重量，计算种子千粒重。

种子千粒重计算公式：

Q=w/n*1000

式中：Q为千粒重，单位为克，w为种子的重量，n为种子的粒数，

本发明还提供了一种测量种子千粒重的装置，包括工作台，该工作台上放置有称重装置和位于所述称重装置上方的图像采集装置，以及分别与所述称重装置和图像采集装置连接的计算机。

所述工作台上竖有支撑杆，所述图像采集装置位于该支撑杆的顶部，所述支撑杆的中部还设有照明灯。图像采集装置对平铺在称重装置上的种子正面拍摄，防止各种子的影像重叠；照明灯选用环形灯，方便图像采集装置从种子的正上方拍摄，且还能缩小种子在图像中的阴影。

本发明具有下述有益效果：

（1）装置结构简单、操作方便、成本较低；

（2）测量精确、速度快，检测结果不受被测种子大小和形状的限制，且被测种子形状最好为圆形、椭圆或类似圆形。

附图说明

图1为本发明测量种子千粒重的装置的结构示意图。

图2为本发明测量种子千粒重的方法的流程图。

图3为本发明中灰度图像。

图4为本发明中的二值化图像。

图5为提取轮廓线上的曲率极值点的示意图。

图6为确定曲率极值点的流程图。

图7为选取匹配点的原理示意图。

图8为基点与匹配点的匹配图。

图9为匹配过程流程图。

图10为长型籽粒的分割图。

图11为圆型籽粒的分割图。

具体实施方式

如图1所示，一种测量种子千粒重的装置，包括工作台3、称重装置5、照明灯2、图像采集器1、支撑杆4、传输线6、传输线8和平板电脑7。

称重装置5用于测量种子的重量，其包括放置种子的托盘和位于托盘下面的重量传感器。支撑杆4竖直安装在工作台上，图像采集器1和照明灯2沿支撑杆4上下布置。图像采集器1位于称重装置5的上方，用于从上方采集种子图像，其优选用USB接口数码摄像头，无需配备图像采集卡，使用方便，降低装置开发成本。照明灯2为种子提供光照，一般为环形灯，既能减小种子的图像中留下的阴影，又能满足摄像头从正上方采集图像的要求。图像采集器1通过传输线6接入平板电脑7，称重装置5通过传输线8连接至平板电脑7，平板电脑7接受图像采集器1和称重装置5传输的信息，对种子图像进行处理，并计算种子数量，再结合种子重量计算千粒重。

如图2所示，利用上述装置计算种子千粒重的步骤如下：

第一步：种子图像预处理

将待测种子平铺在托盘内，摄像头拍摄JPEG格式的图像并传输至平板电脑，平板电脑将接受到的RGB图像转换为灰度图像。转换常用的方法为平均值法，即取红绿蓝三色的平均值为灰度，公式如下：

Gray=（R+G+B）/3

式中，Gray表示灰度值，R、G、B分别表示红绿蓝三色，得到的灰度图像见图3。

图像在形成、传输、接收和处理的过程中，不可避免地存在着外部干扰和内部干扰，图像的常见噪声主要有加性噪声、乘性噪声、量化噪声和椒盐噪声等，反应在图像上，噪声使原本均匀和连续变化的颜色突然变大或减小，形成一些虚假的物体边缘或轮廓，了既平滑掉噪声又尽量保持图像细节，采用中值滤波法消除灰度图像中的噪音。

中值滤波法采用统计滤波器，是一种非线性的空间滤波器，它的响应基于图像滤波器包围的图像区域中像素的排序，然后由统计排序结果决定的值代替中心像素的值。统计滤波器中最常见的例子就是中值滤波器，正如其名，它就是将像素（在中值计算中包括的原像素值）邻域内灰度的中值代替该像素的值。中值滤波器的使用非常普遍，这是因为对于一定类型的随机噪声，它提供了一种优秀的去噪能力，比小尺寸的线性平滑滤波器的模糊程度明显要低。中值滤波器对处理脉冲噪声（也称椒盐噪声）非常有效，因为这种噪声是以黑白点叠加在图像上的。一个数值集合的中值e是这样的数值，即数值集合中，有一半小于或者等于e，还有一半大于等于e。为了对灰度图像上的某个点做中值滤波处理，必须先将掩模内欲求的像素及其邻域的像素值排序，确定其中值，并将中值赋予该像素点。例如，对一个3×3的邻域，其中值是第5个值，而对于5×5的邻域来说，中值就是第13个值，等等。当一个邻域中的一些像素值相同时，它们中的任何一个都可以作为中值。例如，在一个3×3邻域内有一系列像素值（10，20，20，20，15，20，20，25，100），对这些值排序后为（10，15，20，20，20，20，20，25，100），那么其中值就是20。中值滤波器的主要功能是使拥有不同灰度的点看起来更接近与它的邻近值，消除灰度图像中虚假的物体边缘或轮廓。

第二步：将灰度图像转化为二值化图像

采用Ostu算法确定灰度图像的阈值，设灰度图像灰度级是L，则灰度范围为[0,L-1]，利用Ostu算法计算灰度图像的最佳阈值为：

t=Max[w₀(t)×(u₀(t)-u)²+w₁(t)×(u₁(t)-u)²)]

其中的变量说明：当分割的阈值为t时，w₀为背景比例，u₀为背景均值，w₁为前景比例，u₁为前景均值，u为整幅图像的均值。使以上表达式值最大的t，即为分割图像的最佳阈值。

用阈值t分割图像对原图像f(x,y)使用上面的准则找到阈值t，将图像分割为2个部分，也就是我们通常所说的二值化。如图4所示，二值化的图像为：

$g(x,y)=\left\{\begin{array}{cc}0& f(x,y)<T\\ 1& f(x,y)\&GreaterEqual;T\end{array}\right.$

式中，x,y表示图像中各像素坐标，T表示阈值。

对得到的二值图像运用腐蚀算法进行处理，减小干扰。腐蚀算法：用5*5的结构元素，扫描图像的每一个像素用结构元素与其覆盖的二值图像做“与”操作，如果都为1，结果图像的该像素为1。否则为0。定义：

$G^{\&prime;}=G\&CircleTimes;B$

其中B是结构元素，即全是1的五阶矩阵；G为二值图像，为腐蚀运算，G’为腐蚀后的二值图像。

第三步：图像分割

提取二值图像中种子的轮廓线，并依据轮廓线按连通区域划分，选取每个连通区域内的轮廓线上的特征点，匹配特征点对，并在两特征点间建立分割线，分割线将连通区域划分为多个子区域。

连通区域存在一个或多个，且每个连通区域的轮廓线包括一个外轮廓线和多个闭合的内轮廓线。原始的轮廓线上可能存在噪音，为了进一步去除噪音并平滑原始轮廓线，每一条原始轮廓线都用一维高斯核函数进行卷积，一维高斯核函数为：

$g(t,\mathrm{\&sigma;})e^{-\left(\frac{t^2}{2{\mathrm{\&sigma;}}^2}\right)}$

式中：σ为卷积核宽度，t为一维高斯核函数的自变量。卷积后得到的轮廓线由Φ表示：

$\mathrm{\&Phi;}=z\left(t\right)\&CircleTimes;g(t,\mathrm{\&sigma;})$

$=x\left(t\right)\&CircleTimes;g(t,\mathrm{\&sigma;})+\mathrm{jy}\left(t\right)\&CircleTimes;g(t,\mathrm{\&sigma;})$

$=X\left(t\right)+\mathrm{jY}\left(t\right)$

$X\left(t\right)=x\left(t\right)\&CircleTimes;g(t,\mathrm{\&sigma;})$

$Y\left(t\right)=y\left(t\right)\&CircleTimes;g(t,\mathrm{\&sigma;})$

其中，t为一维高斯核函数的自变量，z(t)为原始轮廓线，表示卷积，x(t)表示x方向上的变量值，y(t)表示y方向上的变量值，j表示虚函数，X(t)表示x方向上的卷积值，Y(t)表示y方向上的卷积值。

得到平滑后的轮廓线Φ后，轮廓线的特征点可以由曲率表征出来。计算平滑处理后的轮廓线上各点的曲率，并对由各点的曲率构建的曲率曲线一阶求导，选取求导结果为零的点，并结合求导结果为零的点所对应的曲率中心，得到所述轮廓线上的特征点。

曲率的绝对值大小表示曲线的弯曲程度，曲率绝对值小的表示曲线弯曲程度低。曲率的计算按照以下的公式得到：

$\mathrm{\&kappa;}\left(t\right)=\frac{X^{\&prime;}\left(t\right)Y^{\&prime;\&prime;}\left(t\right)-X^{\&prime;\&prime;}\left(t\right)Y^{\&prime;}\left(t\right)}{{(X^{\&prime;2}\left(t\right)+Y^{\&prime;\&prime;2}\left(t\right))}^{3/2}}$

其中，X'(t)、X''(t)和Y'(t)、Y''(t)分别是X、Y的一阶和二阶导。曲率的正负分别表征曲线的凸和凹。

曲率中心可以按照以下公式计算：

$\left\{\begin{array}{c}\mathrm{Xc}\left(t\right)=X\left(t\right)-\frac{Y^{\&prime;}\left(t\right)(X^{\&prime;2}\left(t\right)+Y^{\&prime;2}\left(t\right))}{Y^{\&prime;\&prime;}\left(t\right)X^{\&prime;}\left(t\right)-X^{\&prime;\&prime;}\left(t\right)Y^{\&prime;}\left(t\right)}\\ \mathrm{Yc}\left(t\right)=Y\left(t\right)-\frac{X^{\&prime;}\left(t\right)(X^{\&prime;2}\left(t\right)+Y^{\&prime;2}\left(t\right))}{Y^{\&prime;\&prime;}\left(t\right)X^{\&prime;}\left(t\right)-X^{\&prime;\&prime;}\left(t\right)Y^{\&prime;}\left(t\right)}\end{array}\right.$

其中，X'(t)、X''(t)和Y'(t)、Y''(t)分别是X、Y的一阶和二阶导；(Xc(t)，Yc(t))是曲率中心的位置，(X(t)，Y(t))是轮廓线上的点。

对曲率曲线κ(t)一阶求导得到：

t_MCCP＝arg{κ′(t)＝0:|κ(t)|＞λ}

式中，t_MCCP表示局部极值点，κ′(t)为κ(t)的一阶导数，λ∈R⁺是设定的阈值用于去除曲率曲线上不明显的极值点和噪声点，λ常设置为30。

根据平滑后轮廓曲线上的极值点位置和相应的曲率中心位置，原始轮廓线上的特征点可以按如下方法找到。

首先建立一个新的笛卡尔坐标系：

O'-x'y'是O-xy经过平移和旋转得到的笛卡尔坐标系。(x₀,y₀)是O'点在原始坐标系中的坐标。θ∈[0,2π]是新坐标系沿着顺时针方向的旋转角。平滑轮廓曲线上的每个极值点都将作为转换后新坐标系的原点。y'轴的正方形是由MCCP指向其曲率中心点D，对应的x'轴的正方向是由y'轴的正方向沿着顺时针方向旋转90度得到，如图5所示。带黑点的折线为原始轮廓线，不带点的平滑曲线为经过高斯核卷积滤波得到的平滑轮廓线。为了确定对应原始轮廓线上的特征点B，先设定一个搜索范围。原始轮廓线上的A点和C点分别分布在x'≥0和x'<0轴的正负半轴上，可以分别用一下方程表示：

$\left\{\begin{array}{c}{t}_A=\underset{t\&Element;\mathrm{\&Delta;}}{\arg \min }\{{\left|\right|(x^{\&prime;}\left(t\right),y^{\&prime;}\left(t\right))-\left(\mathrm{0,0}\right)\left|\right|}_2:|y^{\&prime;}\left(t\right)|<\mathrm{\&lambda;},x^{\&prime;}\left(t\right)\&GreaterEqual;0\}\\ t_C=\underset{t\&Element;\mathrm{\&Delta;}}{\arg \min }\{{\left|\right|(x^{\&prime;}\left(t\right),y^{\&prime;}\left(t\right))-\left(\mathrm{0,0}\right)\left|\right|}_2:|y^{\&prime;}\left(t\right)|<\mathrm{\&lambda;},x^{\&prime;}\left(t\right)<0\}\end{array}\right.$

式中，t_A表示A点的最小值，t_C表示C点的最小值，‖·‖表示2范数，λ是沿着x'轴方向搜索的范围，△代表点集。

如图5所示，原始轮廓线上的特征点B在A点和C点之间与x'轴有最远的距离d。该点可以表示为

t_MCCP＝argmax{|y′(t)|,t∈[min(t_A,t_C),Max(t_A,t_C)]}

式中，t_MCCP表示局部极大值，|y′(t)|表示y方向上的变量y(t)的一阶导数，t_A,t_C表示A、C点的极值。

如图6所示，在原始轮廓线上确定特征点的过程如下：

对原始轮廓线平滑处理，计算平滑处理后的轮廓线上的曲率得到曲率曲线，对曲率曲线一阶求导得到轮廓线段上的曲率极值点，并结合该曲率极值点对应的曲率中心，根据上述的方法可得到原始轮廓线上的特征点。

为了将连通区域内每颗种子的所在的区域标记出来，需要确定每颗种子的边缘，根据种子间边缘的交点（即上面得到的特征点），对特征点进行配对，并在配对成功的两个特征点间建立分割线，即为每颗种子在连通区域内的边界线。

选定其中一个特征点为基点BP（basicpoint），与其对应的形成分割线的特征点成为匹配点MP（matchingpoint）。为了正确的找到与基点对应的匹配点，该匹配点应同时满足以下4个要求：

条件1：

如图7所示，设O点为BP，与其相对应的MP为C点。从图中可以看出这两点的曲率法线夹角很小，因此我们可以根据BP点的法线（l₁）位置设定一个辐射角范围使得与其正确对应的MP点落在此范围内。法线l₁是O点与该点曲率中心点A连接线的延长线。将l₁线沿着O点左右旋转θ角，得到的辐射区域角2θ，如图中阴影部分所示就是C点所处的范围，我们称其为置信区间。θ角称为置信角，常设置为45度。法线l₁与连接线l₂的夹角可以通过以下方程式计算：

$\mathrm{\&alpha;}=\arctan ({(-1)}^{\mathrm{Heaviside}}\&CenterDot;\frac{k_2-k_1}{1+k_1{k}_2})+\mathrm{\&pi;}\&CenterDot;\mathrm{Heaviside}$

$\mathrm{Heaviside}=\left\{\begin{array}{c}0,x\cos \mathrm{\&phi;}+y\sin \mathrm{\&phi;}-(x_0\cos \mathrm{\&phi;}+y_0\sin \mathrm{\&phi;})\&GreaterEqual;0\\ 1,x\cos \mathrm{\&phi;}+y\sin \mathrm{\&phi;}-(x_0\cos \mathrm{\&phi;}+y_0\sin \mathrm{\&phi;})<0\end{array}\right.$

式中，k₁和k₂分别是l₁和l₂的斜率，α∈[0,π]是l₁和l₂的夹角，(x₀,y₀)和(x,y)分别是基点O和匹配点C的坐标，Φ是x'与x轴的夹角。夹角α小于置信角θ，说明C点有可能是O点的正确匹配点。而β角大于置信角θ，说明B点不是O点的正确匹配点。

条件2：

如图8所示，同一个连通区域中，该连通区域内同一闭合的内轮廓线上的特征点，彼此之间不能相互匹配，只能与其他内轮廓线或外轮廓线上的特征点进行匹配。如果外轮廓线里面包含内轮廓线，则优先选择外轮廓线上的特征点与内轮廓线上的特征点配对。如果外轮廓线的闭合区间里没有包含内轮廓线，则外轮廓线上的特征点相互配对。

条件3：

满足条件1和条件2的匹配点存在多个，利用所有的匹配点进行逐一验证，选取与基点距离最小的匹配点作为正确的匹配点。验证公式如下：

$\{t_i,t_j\}=\underset{\mathrm{\&zeta;},\mathrm{\&xi;}\&Element;\mathrm{\&Delta;}}{\arg \min }\{{\left|\right|(x\left(\mathrm{\&zeta;}\right),y\left(\mathrm{\&zeta;}\right))-(x\left(\mathrm{\&xi;}\right),y\left(\mathrm{\&xi;}\right))\left|\right|}_2,\mathrm{\&zeta;}\&NotEqual;\mathrm{\&xi;}\}$

式中，‖·‖表示2范数，△表示点集，ξ和ζ分别表示BP和MP，点（x(t_i),y(t_i)）和点（x(t_j),y(t_j)）分别为基点和距离基点最近的匹配点。

条件4：

连接线的长度要小于种子的纵径的长度，这可以在一定程度上防止误匹配。

$d<C\&CenterDot;\underset{i\&Element;\&lsqb;l,n\&rsqb;}{\max }\left\{L\left(i\right)\right\}$

式中，d为基点和匹配点之间的距离，L(i)是第i个种子的最大长度，是由该种子轮廓线上距离最大的两个像素点之间的距离决定，n是种子的个数，C∈R⁺是一个常数。

如图9所示，根据上面四个条件，对每个特征点进行逐一匹配，直至所有的特征点都找到相应的匹配特征点。

所有的匹配特征点对找出后，在相匹配的两个特征点之间建立分割线，完成对每个连通区域内单粒种子的分割，每个连通区域被轮廓线和分割线划分为单颗种子所处的子区域，图10表示长型种子籽粒分割后的示意图，图11表示圆型种子籽粒分割后的示意图。

第四步：对分割得到的子区域进行标记和计数

对粘连的种子组成的连通区域分割后，每颗种子在图像上被轮廓线和分割线划分为单据的子区域，将分割线的像素值设置为0，从图像的左上角进行检测，第一个不为零的像素值极为第一个标记子区域中的像素点，针对该像素点运用八邻域算子进行八个方向的搜索，若其周围八个方向上的像素值为1，则对为1像素进行标记，不为1的像素不用标记，再以新标记为1的像素出发，用八邻域算子进行搜索，直到第一个子区域标记完毕。按行进行搜索，以同样的方法进行第二个子区域的标记，直到所有子区域标记完毕，结束标记，得到标记的子区域个数为n，n即是种子的粒数。

第五步：根据种子的重量和个数，计算种子千粒重。

利用平板电脑读取重量传感器的数据，种子的重量为w，单位为克，通过上述的步骤计算得到种子的籽粒数为n。根据种子千粒重计算公式：

Q=w/n*1000

式中：Q为千粒重，单位为克，w为种子的重量，n为种子的籽粒数，并将算得的千粒重显示在平板电脑的显示器上。

Claims (5)

1.一种测量种子千粒重的方法，其特征在于，包括以下几个步骤：

1)将待测试的多粒种子平铺，并采集种子的图像和重量；

2)将所述的图像转换为灰度图像，并对灰度图像进行滤波处理；

3)对滤波处理后的灰度图像进行二值化处理得到二值化图像；

4)提取二值化图像中种子的轮廓线，并依据轮廓线按连通区域划分，选取每个连通区域内的轮廓线上的特征点，匹配特征点对，并在两特征点对间建立分割线，所述的分割线将连通区域划分为多个子区域；

所述的轮廓线包括外轮廓线和内轮廓线，针对轮廓线上的特征点A，与特征点A匹配的特征点B应满足以下几个条件：

a)建立通过特征点A的法线，特征点B与所述法线之间的距离最小；

b)特征点B与特征点A之间的距离小于种子的长度；

c)特征点B与特征点A处于不同的内轮廓线上；

d)特征点A与特征点B之间的距离最短；

在步骤4)中，对所述的轮廓线平滑处理，计算平滑处理后的轮廓线上各点的曲率，并对由各点的曲率构建的曲率曲线一阶求导，选取求导结果为零的点，并结合求导结果为零的点所对应的曲率中心，得到所述轮廓线上的特征点；

求得轮廓线上的特征点的过程如下：

首先建立一个新的笛卡尔坐标系，O'-x'y'是O-xy经过平移和旋转得到的笛卡尔坐标系，(x₀,y₀)是O'点在原始坐标系中的坐标，θ∈[0,2π]是新坐标系沿着顺时针方向的旋转角，平滑轮廓曲线上的每个极值点都将作为转换后新坐标系的原点，y'轴的正方向是由MCCP指向其曲率中心点D，对应的x'轴的正方向是由y'轴的正方向沿着顺时针方向旋转90度得到，为了确定对应原始轮廓线上的特征点B，先设定一个搜索范围，原始轮廓线上的A点和C点分别分布在x'≥0和x'<0轴的正负半轴上，可以分别用以下方程表示：

$\left\{\begin{array}{c}{t}_A=\underset{t\&Element;\mathrm{\&Delta;}}{\arg \min }\left\{\left|\right|\left(x^{\&prime;}\left(t\right),y^{\&prime;}\left(t\right)\right)-\left(0,0\right)|{|}_2:\left|y^{\&prime;}\left(t\right)\right|<\mathrm{\&lambda;},x^{\&prime;}\left(t\right)\&GreaterEqual;0\right\}\\ t_C=\underset{t\&Element;\mathrm{\&Delta;}}{\arg \min }\left\{\left|\right|\left(x^{\&prime;}\left(t\right),y^{\&prime;}\left(t\right)\right)-\left(0,0\right)|{|}_2:\left|y^{\&prime;}\left(t\right)\right|<\mathrm{\&lambda;},x^{\&prime;}\left(t\right)<0\right\}\end{array}\right.$

式中，t_A表示A点的最小值，t_C表示C点的最小值，‖·‖表示2范数，λ是沿着x'轴方向搜索的范围，△代表点集；原始轮廓线上的特征点B在A点和C点之间与x'轴有最远的距离d，该特征点B可以表示为

t_MCCP＝argmax{|y′(t)|,t∈[min(t_A,t_C),Max(t_A,t_C)]}

式中，t_MCCP表示局部极大值，|y′(t)|表示y方向上的变量y(t)的一阶导数，t_A,t_C表示A、C点的极值；

5)统计所有连通区域内的子区域数量，得到所述种子的粒数；

6)根据种子的粒数和重量，计算种子千粒重。

2.如权利要求1所述的测量种子千粒重的方法，其特征在于，所述的步骤2)中，采用中值滤波法对灰度图像进行滤波处理。

3.如权利要求1所述的测量种子千粒重的方法，其特征在于，所述的步骤3)中，在进行二值化处理之前，采用Ostu算法得到滤波处理后灰度图像的阈值。

4.如权利要求3所述的测量种子千粒重的方法，其特征在于，在步骤3)中，对所述的二值化图像进行腐蚀算法处理。

5.如权利要求1所述的测量种子千粒重的方法，其特征在于，在步骤5)中，将所述分割线的像素值设置为零，选取子区域中像素值不为零的像素a，运用八邻域算子对该像素a周围八个像素进行搜索，选择周围八个像素与像素a的像素值相同的像素进行标记，遍历所有子区域中的每个像素，完成对所有子区域的标记，并计算标记的子区域的个数。