CN104715491A - 一种基于一维灰度矩的亚像素边缘检测方法 - Google Patents

Abstract

本发明公开了一种基于一维灰度矩的亚像素边缘检测方法，包括以下步骤：S1：采用中值滤波方式对待处理图像进行去噪处理；S2：采用Canny边缘检测算子方法对待处理图像进行像素级边缘检测；S3：采用一维灰度矩方法对图像的像素进行边缘检测，完成该图像的亚像素边缘检测。本发明首先利用中值滤波器进行去噪的处理，其次利用Canny算子进行像素级边缘检测，然后在空间域的笛卡尔坐标下利用一维灰度矩进行像素边缘检测。

Description

一种基于一维灰度矩的亚像素边缘检测方法

技术领域

本发明涉及图形处理技术领域，尤其涉及一种基于一维灰度矩的亚像素边缘检测方法。

背景技术

图像边缘包含了大量的目标对象的有用信息，准确地提取图像的边缘对后续图像的处理，如物体配准、物体的尺寸测量、物体的检测与识别等也有着决定性的影响，因此，边缘检测技术在基于计算机视觉的检测系统中有着十分重要的作用。基于计算机视觉的边缘检测系统具有非接触、精度高、成本低、适用范围广、检测速度快、自动化程度高等优点。在工业生产与制造中，如螺母、轴承等零件的精确识别和配准，都将应用到基于计算机视觉的边缘检测技术。

随着工业的发展对检测精度的不断提高，传统像素级的边缘检测已经不能满足工业生产的实际需求，人们据此提出了更高精度的亚像素边缘检测技术，该技术将图像中的目标对象定位到亚像素级上，等效地提高了系统的分辨率，例如当检测的精度是0.1像素时，系统的分辨率就提高了10倍。

现有技术中公开号CN101477685B的中国专利中，申请人公开了一种具有景深零件加工质量的亚像素级图像检测方法。该发明首先对机器视觉系统进行分层标定，其次，对原始图像进行插值计算，通过粗精两步法实现零件边缘的精确定位，最后利用建立的图像各层面与零件各层面映射关系计算具有景深零件的形状，关键尺寸参数，通过比较分析获得质量数据。但是该技术是利用插值的方法来进行亚像素边缘检测的，由于插值技术本身很容易受噪声的影响，那么该技术也很容易受噪声的干扰，这样会造成边缘检测的精度降低，进而给后面的操作造成不好的影响。

2014年，Chen等^[3]在文献中“Sub-pixel dimensional measurement withLogistic edge model”中利用基于sigmod内核建立一个模糊边缘模型，然后利用建立的边缘模型所得到的灰度值与实际的灰度相减，然后使其均方误差最小的参数就是边缘所在的位置。该方案利用sigmod函数内核拟合的方法来进行亚像素边缘检测，虽然这种方法对噪声不敏感，但是精度不够高，并且模糊边缘模型需要进行积分，特别是要进行迭代运算来得到边缘的参数，这都将增加计算的复杂度。

发明内容

根据现有技术存在的问题，本发明公开了一种基于一维灰度矩的亚像素边缘检测方法，具体采用如下步骤：

S1：采用中值滤波方式对待处理图像进行去噪处理；

S2：采用Canny边缘检测算子方法对待处理图像进行像素级边缘检测：对待处理的图像进行Gauss滤波处理，计算图像梯度的幅度和方向，采用非极大值抑制算法对图像的边缘进行细化，通过双阈值递归方法进行图像边缘点的检测，按照上述方式对待处理图像的所有像素进行像素级边缘检测；

S3：采用一维灰度矩方法对图像的像素进行边缘检测：建立一维阶跃边缘模型和二次曲线边缘模型，以待处理像素为中心将其周围的像素分成n行三列构成的三个区域：待处理像素所在的列为第二区域，第二区域的左右两列分别是第一区域和第三区域，计算每个区域所对应像素的一阶灰度矩和二阶灰度矩，利用每个区域所对应的一阶灰度矩和二阶灰度矩的结果来计算该像素的二次曲线边缘参数，按照上述方式对待处理图像的所有像素进行处理，完成该图像的亚像素边缘检测。

S1中对输入图像进行去除噪处理采用如下方式：

S11：以待处理图像中待处理像素为中心把一个滑动窗口内的多个像素灰度排序；

S12：用已排好序的中间值代替窗口中心像素原来的灰度值，对待处理图像的所有像素进行上述操作完成去噪处理。

S3中采用一维灰度矩方法对图像进行亚像素边缘检测时采用如下方法：

设二次曲线边缘模型为：

y＝ax²+bx+c    (9)

其中，a、b、c是二次曲线的参数，也是一维阶跃边缘模型边缘的参数，所述二次曲线边缘模型将一维阶跃边缘模型分成两个部分；

在三个区域分别利用一阶灰度矩和二阶灰度矩进行边缘参数的计算，具体的计算过程如下：

S31：对于第一区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(1\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(1\right)}A+p_2^{\left(1\right)}B=B+p_1^{\left(1\right)}(A_1-B_1)---\left(10\right)$

$M_2^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(1\right)}{A}^2+p_2^{\left(1\right)}{B}^2=B^2+p_1^{\left(1\right)}({A_1}^2-{B_1}^2)---\left(11\right)$

$p_1^{\left(1\right)}=\left[{\∫}_{s\∈S\left(1\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{-1.5}{\overset{-0.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c-b+\frac{13}{12}a+\frac{n}{2})/n---\left(12\right)$

其中，分别是第一区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第一区域面积的百分比，其中A₁为一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₁是二次曲线边缘模型下方区域的灰度值，(a,b,c)是边缘的参数，S(1)代表第一区域的面积，灰度值A₁＝A,B₁＝B。

S32：对于第二区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(2\right)}A+p_2^{\left(2\right)}B=B+p_1^{\left(2\right)}(A_2-B_2)---\left(13\right)$

$M_2^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(2\right)}{A}^2+p_2^{\left(2\right)}{B}^2=B^2+p_1^{\left(2\right)}({A_2}^2-{B_2}^2)---\left(14\right)$

$p_1^{\left(2\right)}=\left[{\∫}_{s\∈S\left(2\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{-0.5}{\overset{-0.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c+\frac{1}{12}a+\frac{n}{2})/n---\left(15\right)$

其中，分别是第二区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第二区域面积的百分比，其中A₂为一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₂是二次曲线边缘模型下方区域的灰度值，(a,b,c)是边缘的参数，S(2)代表第二区域的面积，灰度值A₂＝A,B₂＝B。

S33：对于第三区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(3\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(3\right)}A+p_2^{\left(3\right)}B=B+p_1^{\left(3\right)}(A_3-B_3)---\left(16\right)$

$M_2^{\left(3\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(3\right)}{A}^2+p_2^{\left(3\right)}{B}^2=B^2+p_1^{\left(3\right)}({A_3}^2-{B_3}^2)---\left(17\right)$

$p_1^{\left(3\right)}=\left[{\∫}_{s\∈S\left(3\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{0.5}{\overset{1.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c+b+\frac{13}{12}a+\frac{n}{2})/n---\left(18\right)$

其中，分别是第三区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第三区域面积的百分比，其中A₃一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₃是二次曲线边缘模型下方区域的灰度值(a,b,c)是边缘的参数，S(3)代表第三区域的面积，灰度值A₃＝A,B₃＝B。

S34：利用上述每个区域所对应的一阶灰度矩和二阶灰度矩计算结果和公式(10)～(18)得到像素的二次曲线边缘模型的参数(a,b,c),具体的算法公式详见式(27)～(31)：

$a=\frac{1}{4}(\frac{M_2^{\left(1\right)}+M_2^{\left(3\right)}-2M_2^2}{(A^2-B^2)}+\frac{M_1^{\left(1\right)}+M_1^{\left(3\right)}-2M_1^2}{(A-B)})---\left(27\right)$

$b=\frac{1}{4}(\frac{M_2^{\left(3\right)}-M_2^{\left(1\right)}}{(A^2-B^2)}+\frac{M_1^{\left(3\right)}-M_1^{\left(1\right)}}{2(A-B)})---\left(28\right)$

$c=\frac{1}{2}(\frac{2M_2^2-n(A^2+B^2)}{2(A^2-B^2)}+\frac{2M_1^2-n(A+B)}{2(A-B)})-\frac{1}{12}a---\left(29\right)$

$A=\frac{T_1+\sqrt{{T_1}^2-4T_2}}{2}---\left(30\right)$

$B=\frac{T_1-\sqrt{{T_1}^2-4T_2}}{2}---\left(31\right)$

其中分别是第一区域，第二区域，第三区域的一阶灰度矩，分别是第一区域，第二区域，第三区域的二阶阶灰度矩，中间参数值T₁，T₂为 $T_1=\frac{M_2^{\left(3\right)}-M_2^{\left(1\right)}}{M_1^{\left(3\right)}-M_1^{\left(1\right)}},T_2=\frac{T_1{M}_1^{\left(2\right)}-M_2^{\left(2\right)}}{n},$ 区域长度n＝9。

由于采用了上述技术方案，本发明提供的一维灰度矩的亚像素边缘检测方法，首先利用中值滤波器进行去噪的处理，其次利用Canny算子进行像素级边缘检测，然后在空间域的笛卡尔坐标下利用一维灰度矩进行像素边缘检测。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明中基于一维灰度矩的亚像素边缘检测方法的流程图；

图2为本发明中一维阶跃边缘模型的示意图；

图3为本发明中一维阶跃边缘模型与二次曲线边缘模型相结合的示意图；

图4(a)图4(b)为本发明采用的方法对不同角度的直线的检测结果示意；

图5为本发明中不同角度的合成图像的亚像素边缘检测结果示意图；

图6为本发明中实际图像的亚像素边缘检测效果示意图。

具体实施方式

为使本发明的技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：

如图1所示的一种基于一维灰度矩的亚像素边缘检测方法，具体包括以下步骤：

S1：采用中值滤波方式对待处理图像进行去噪处理：

输入图像若含有噪声，将影响像素级边缘检测结果，所以在采集到图像后，先要对图像进行中值滤波器去噪的预处理。中值滤波器^[4]是局部平均平滑滤波器，也是最常用的非线性平滑滤波器。中值滤波的原理是把一个滑动窗口内的诸像素灰度排序，用中间值代替窗口中心像素原来的灰度：

$g_{\mathrm{median}}(x,y)=\underset{(s,t)\∈N(x,y)}{\mathrm{median}}\left[f(s,t)\right]---\left(1\right)$

其中，median代表取中值操作。

S1中对输入图像进行去除噪处理采用如下方式：

S11：以待处理图像中待处理像素为中心把一个滑动窗口内的多个像素灰度排序；

S12：用已排好序的中间值代替窗口中心像素原来的灰度值，对待处理图像的所有像素进行上述操作完成去噪处理。

S2：采用Canny边缘检测算子方法对待处理图像进行像素级边缘检测：对待处理的图像进行Gauss滤波处理，计算图像梯度的幅度和方向，采用非极大值抑制算法对图像的边缘进行细化，通过双阈值递归方法进行图像边缘点的检测，按照上述方式对待处理图像的所有像素进行像素级边缘检测。

Canny边缘检测算子是一类最优边缘检测算子^[5]。该算子的检测基本步骤如下所示：

2.1、对要处理的图像进行先进行Gauss滤波，这样做的目的是滤除高斯等有害的噪声；

2.2、按照式(2)、式(3)计算图像的梯度的幅度和方向

$f(i,j)=\sqrt{M_x^2(i,j)+M_y^2(i,j)}---\left(2\right)$

$H(i,j)=\arctan \frac{M_y(i,j)}{M_x(i,j)}---\left(3\right)$

其中，M(i,j)是图像的灰度值，M_x(i,j)和M_y(i,j)分别是纵轴和横轴邻域差，f(i,j)是梯度的幅度，H(i,j)是梯度的方向；

2.3、得到梯度的幅度值以后，不能直接地将幅度和阈值进行比较，应该采用为“非极大值抑制(Non-maxima Suppression)”的技术来使边缘细化；

2.4、虽然进过了极大值抑制，但是仍然存在一定虚假边缘，那么就需要利用阈值技术进行处理：通过双阈值递归方法来进行边缘点的检测，从而实现了真边缘的提取。

S3：采用一维灰度矩方法对图像的像素进行边缘检测：建立一维阶跃边缘模型和二次曲线边缘模型，以待处理像素为中心将其周围的像素分成n行三列构成的三个区域：待处理像素所在的列为第二区域，第二区域的左右两列分别是第一区域和第三区域，计算每个区域所对应像素的一阶灰度矩和二阶灰度矩，利用每个区域所对应的一阶灰度矩和二阶灰度矩的结果来计算该像素的二次曲线边缘参数，按照上述方式对待处理图像的所有像素进行处理，完成该图像的亚像素边缘检测。

1984年，Tabatabai和Mitchell首先建立阶跃模型，然后利用二维灰度矩变量来求得边缘变量的四个参数(位置，方向，背景的灰度值，目标的灰度值)。经过三十多年的理论和实际的研究，人们在其基础上进了大量的实验，并且取得了较好的检测精度，但是这些研究成果基本上都是在映射圆上并且是在极坐标上进行计算的。而本发明是在空间域笛卡尔坐标(直角坐标)下利用一维灰度矩进行计算的。

(1)建立一维阶跃边缘模型，其表达式为

$f\left(x\right)=\left\{\begin{array}{cc}{h}_1,& x\≤k\\ h_2,& x>k\end{array}\right.---\left(4\right)$

其中，f(x)是在x处的灰度值，h₁、h₂是边缘两侧的灰度值，k是边缘所在的位置。对于连续的图像，一维连续q阶灰度矩的定义详见下式，

M_q＝∫f^q(x)dx    (5)

因为本发明处理的图像数据是离散数据，所以本发明将详细介绍一下离散灰度矩，其一维阶跃边缘模型详见图2。

设p₁和p₂分别是灰度值为h₁和h₂的像素点所占的比例，且p₁+p₂＝1，则

$p_1=\frac{k}{n}---\left(6\right)$

$p_2=\frac{n-k}{n}=1-p_1---\left(7\right)$

利用上面的阶跃模型和灰度矩的定义可得，

$m_q=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{p}_j{h}_j^q=\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{g}_j^q---\left(8\right)$

其中，g_j(j＝1,2,…,n)是图像中实际的灰度值，n是总的像素点数。

(2)基于一维灰度矩的边缘参数的求解本发明所利用的边缘模型详见图3，★代表第三区域；其中图中A代表着一维阶跃边缘模型中二次曲线边缘模型上方的灰度值；B代表着一维阶跃边缘模型中二次曲线边缘模型下方的灰度值；n代表着区域的长度；■代表着其正下面的区域即第一区域；①代表着边缘二次曲线下面■所对应的区域；●代表着其正下面的区域即第二区域；②代表着边缘二次曲线下面●所对应的区域；★代表着其正下面的区域即第三区域；③代表着边缘二次曲线下面★所对应的区域。

如图3所示，本发明利用的二次曲线边缘模型为，其定义如下所示，

y＝ax²+bx+c    (9)

其中，其中，a、b、c是二次曲线的参数，也是一维阶跃边缘模型边缘的参数，所述二次曲线边缘模型将一维阶跃边缘模型分成两个部分。

下面将详解地介绍边缘的参数的求解。

S31：对于第一区域如图中的■区域，其一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(1\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(1\right)}A+p_2^{\left(1\right)}B=B+p_1^{\left(1\right)}(A_1-B_1)---\left(10\right)$

$M_2^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(1\right)}{A}^2+p_2^{\left(1\right)}{B}^2=B^2+p_1^{\left(1\right)}({A_1}^2-{B_1}^2)---\left(11\right)$

$p_1^{\left(1\right)}=\left[{\∫}_{s\∈S\left(1\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{-1.5}{\overset{-0.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c-b+\frac{13}{12}a+\frac{n}{2})/n---\left(12\right)$

其中，分别是第一区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第一区域面积的百分比，其中A₁为一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₁是二次曲线边缘模型下方区域的灰度值，(a,b,c)是边缘的参数，S(1)代表第一区域的面积，灰度值A₁＝A,B₁＝B。

S32：对于第二区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(2\right)}A+p_2^{\left(2\right)}B=B+p_1^{\left(2\right)}(A_2-B_2)---\left(13\right)$

$M_2^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(2\right)}{A}^2+p_2^{\left(2\right)}{B}^2=B^2+p_1^{\left(2\right)}({A_2}^2-{B_2}^2)---\left(14\right)$

$p_1^{\left(2\right)}=\left[{\∫}_{s\∈S\left(2\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{-0.5}{\overset{-0.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c+\frac{1}{12}a+\frac{n}{2})/n---\left(15\right)$

其中，分别是第二区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第二区域面积的百分比，其中A₂为一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₂是二次曲线边缘模型下方区域的灰度值，(a,b,c)是边缘的参数，S(2)代表第二区域的面积，灰度值A₂＝A,B₂＝B。

S33：对于第三区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(3\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(3\right)}A+p_2^{\left(3\right)}B=B+p_1^{\left(3\right)}(A_3-B_3)---\left(16\right)$

$M_2^{\left(3\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(3\right)}{A}^2+p_2^{\left(3\right)}{B}^2=B^2+p_1^{\left(3\right)}({A_3}^2-{B_3}^2)---\left(17\right)$

$p_1^{\left(3\right)}=\left[{\∫}_{s\∈S\left(3\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{0.5}{\overset{1.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c+b+\frac{13}{12}a+\frac{n}{2})/n---\left(18\right)$

其中，分别是第三区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第三区域面积的百分比，其中A₃一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₃是二次曲线边缘模型下方区域的灰度值(a,b,c)是边缘的参数，S(3)代表第三区域的面积，灰度值A₃＝A,B₃＝B。

S34：利用上述每个区域所对应的一阶灰度矩和二阶灰度矩计算结果和公式(10)～(18)得到像素的二次曲线边缘模型的参数(a,b,c),具体的算法公式详见式(27)～(31)：

$a=\frac{1}{4}(\frac{M_2^{\left(1\right)}+M_2^{\left(3\right)}-2M_2^2}{(A^2-B^2)}+\frac{M_1^{\left(1\right)}+M_1^{\left(3\right)}-2M_1^2}{(A-B)})---\left(27\right)$

$b=\frac{1}{4}(\frac{M_2^{\left(3\right)}-M_2^{\left(1\right)}}{(A^2-B^2)}+\frac{M_1^{\left(3\right)}-M_1^{\left(1\right)}}{2(A-B)})---\left(28\right)$

$c=\frac{1}{2}(\frac{2M_2^2-n(A^2+B^2)}{2(A^2-B^2)}+\frac{2M_1^2-n(A+B)}{2(A-B)})-\frac{1}{12}a---\left(29\right)$

$A=\frac{T_1+\sqrt{{T_1}^2-4T_2}}{2}---\left(30\right)$

$B=\frac{T_1-\sqrt{{T_1}^2-4T_2}}{2}---\left(31\right)$

其中分别是第一区域，第二区域，第三区域的一阶灰度矩，分别是第一区域，第二区域，第三区域的二阶阶灰度矩，中间参数值T₁，T₂为 $T_1=\frac{M_2^{\left(3\right)}-M_2^{\left(1\right)}}{M_1^{\left(3\right)}-M_1^{\left(1\right)}},T_2=\frac{T_1{M}_1^{\left(2\right)}-M_2^{\left(2\right)}}{n},$ 区域长度n＝9。

这样通过上面的式(27)～式(31)，二次曲线边缘的参数(a,b,c,A,B)就得到了。

实施例：

为了验证本发明的有效性，进行了计算机仿真实验。在实验中，实验参数为CPU IntelR CoreTM i32.4GHz，2G内存，显卡是ATI Mobility Radeon HD5470，系统为Window7家庭版，软件编程环境为Matlab2010b。本发明实验的图像是利用人工合成的图像和实际图像，对于人工合成的图片的大小是400像素×400像素，而对于实际图像大小是512像素×512像素。

噪声比的定义：

$\mathrm{NR}=\frac{{\mathrm{\σ}}_n}{k}\×100\%---\left(32\right)$

其中，σ_n是加入的高斯白噪声的标准方差，k是目标对象灰度值与背景灰度的差值。

边缘参数(位置，方向)的误差定义：

${\mathrm{\ϵ}}_i={\stackrel{\‾}{r}}_i-r_i,i=\mathrm{1,2,3}...M---\left(33\right)$

其中，是利用基于伪Zernike矩计算出来的第i个边缘的估计值，r_i是第i个边缘的实际值，M是检测出来的总的边缘点总数。

边缘参数(位置，方向)的最大误差，平均误差，误差的标准方差的定义：

ε_max＝max{|ε_i|}  i＝1,2,3...M    (34)

${\mathrm{\μ}}_{\mathrm{\ϵ}}=\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|{\mathrm{\ϵ}}_i\right|---\left(35\right)$

${\mathrm{\σ}}_{\mathrm{\ϵ}}=\sqrt{\frac{M\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i^2-{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i\right)}^2}{M(M-1)}}---\left(36\right)$

其中，ε_i是边缘参数(位置，方向)的误差，M是检测的边缘点总数。

(2)不同方法的仿真结果

将本发明方法与文献[2，3，7]中技术进行了对比实验，具体的仿真结果，详见图4、图5和图6所示，以及表1—表3所示。

表1不同方法的边缘位置的误差

表2不同方法的边缘方向的误差

表3不同方法的边缘检测时间

从表1可见，在不同角度直线的实验情况下，文献[3]的位置检测的精度比文献[2]的检测精度高，而文献[7]比文献[3]的检测精度高，而本文所提出的方法比文献[7]更高，也就是说本文的方法在这几种方法中是较好的；从表2可以看出，对不同角度直线的检测，不同的方法有着不同的结果，文献[3]的方向检测的精度比文献[2]和文献[7]的检测精度高，这是因为文献[3]利用了迭代的方法，这在一定程度上改善了方向检测精度，而文献[7]比文献[3]的检测精度高，而本文的方法比文献[3]高。从表3可以看出本文所提的方法的计算复杂度最小。

如图4所示，本文的方法对不同角度的直线都有着很好的检测结果的示意图，并且从图4(a)可以看出，本文所提的方法最大误差在0.45像素以下，平均误差在0.1像素以下，标准误差在0.07像素左右，那么本文所提的方法其位置检测精度到达了0.08像素左右的精度；从图4(b)可以看出，本文所提的方法最大误差在1.4度以下，平均误差在0.4°以下，标准误差在0.2°左右，那么本文所提的方法其方向检测精度到达了0.3°左右的精度。

图4随角度变化下一维灰度矩的检测结果：(a)在不同角度下位置的检测结果；(b)在不同角度下，角度的检测结果。

图5为不同角度的合成图像的亚像素边缘检测结果：(a)～(f)没有加白噪声的边缘检测结果；(a)角度为零的检测结果；(b)角度为零的检测结果的局部放大图；(c)角度为-25度的检测结果；(d)角度为-25度的检测结果的局部放大图；(e)角度为65度的检测结果；(f)角度为65度的检测结果的局部放大图；(g)～(l)加入噪声的边缘检测结果；(g)角度为零的检测结果；(h)角度为零的检测结果的局部放大图；(i)角度为-25度的检测结果；(j)角度为-25度的检测结果的局部放大图；(k)角度为65度的检测结果；(l)角度为65度的检测结果的局部放大图。

图6为实际图像的亚像素边缘检测效果：(a)～(c)具有丰富边缘的图像检测结果；(a)b图的局部放大图；(b)边缘的检测结果；(c)b图的局部放大图；(d)～(f)简单圆环工件的检测图像检测结果；(d)e图的局部放大图；(e)圆环的检测结果；(f)e图的局部放大图。

从图5可以看出，本发明中基于一维灰度矩的方法对合成的直线的灰度图像有着很好的检测精度，并且有一定的抗噪能力。而从图6可以看出，本发明方法在提取实际的物体边缘也有着较好的效果，图6(a)显示对于有着较多边缘的玩具虎，其提取的边缘精度较高，从图6(b)显示而对于简单的环形工件，同样也有着很好的检测结果。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

参考文献(如专利/论文/标准)

[1]曹建福，史博，汪霖，张家良.具有景深零件加工质量的亚像素级图像检测方法.中国，公开号：CN101477685B.2011.06.01.

[2]Sun Q,Hou Y,Tan Q,et al.A robust edge detection method with sub-pixelaccuracy[J].Optik International Journal for Light and Electron Optics,2014,125(14):3449-3453.

[3]Chen P,Chen F,Han Y,et al.Sub-pixel dimensional measurement with Logisticedge model[J].Optik-International Journal for Light and Electron Optics,2014,125(9):2076-2080.

[4]章毓晋.图像处理和分析教程[M].北京：人民邮电出版社，2009

[5]Canny J.A computational approach to edge detection[J].IEEE Transactions onPattern Analysis and Machine Intelligence.1986(6):679-698.[5]Sun Q,Hou Y,Tan Q,et al.A robust edge detection method with sub-pixel accuracy[J].OptikInternational Journal for Light and Electron Optics,2014,125(14):3449-3453.

[6]Ye J,Fu G,Poudel U P.High-accuracy edge detection with blurred edge model[J].Image and Vision Computing,2005,23(5):453-467.

[7]Kaur A,Singh C.Sub-pixel edge detection using pseudo Zernike moment[J].Int.J.Signal Process.Image Process.Pattern Recognit.2011,4(2):107-118.

Claims (3)

1.一种基于一维灰度矩的亚像素边缘检测方法，其特征在于：包括以下步骤：

S1：采用中值滤波方式对待处理图像进行去噪处理；

S2：采用Canny边缘检测算子方法对待处理图像进行像素级边缘检测：对待处理的图像进行Gauss滤波处理，计算图像梯度的幅度和方向，采用非极大值抑制算法对图像的边缘进行细化，通过双阈值递归方法进行图像边缘点的检测，按照上述方式对待处理图像的所有像素进行像素级边缘检测；

S3：采用一维灰度矩方法对图像的像素进行边缘检测：建立一维阶跃边缘模型和二次曲线边缘模型，以待处理像素为中心将其周围的像素分成n行三列构成的三个区域：待处理像素所在的列为第二区域，第二区域的左右两列分别是第一区域和第三区域，计算每个区域所对应像素的一阶灰度矩和二阶灰度矩，利用每个区域所对应的一阶灰度矩和二阶灰度矩的结果来计算该像素的二次曲线边缘参数，按照上述方式对待处理图像的所有像素进行处理，完成该图像的亚像素边缘检测。

2.根据权利要求1一种基于一维灰度矩的亚像素边缘检测方法，其特征还在于：S1中对输入图像进行去除噪处理采用如下方式：

S11：以待处理图像中待处理像素为中心把一个滑动窗口内的多个像素灰度排序；

S12：用已排好序的中间值代替窗口中心像素原来的灰度值，对待处理图像的所有像素进行上述操作完成去噪处理。

3.根据权利要求1所述的一种基于一维灰度矩的亚像素边缘检测方法，其特征还在于：S3中采用一维灰度矩方法对图像进行亚像素边缘检测时采用如下方法：

设二次曲线边缘模型为：

y＝ax²+bx+c    (9)

其中，a、b、c是二次曲线的参数，也是一维阶跃边缘模型边缘的参数，所述二次曲线边缘模型将一维阶跃边缘模型分成两个部分；

在三个区域分别利用一阶灰度矩和二阶灰度矩进行边缘参数的计算，具体的计算过程如下：

S31：对于第一区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(1\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(1\right)}A+p_2^{\left(1\right)}B=B+p_1^{\left(1\right)}(A_1-B_1)---\left(10\right)$

$M_2^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(1\right)}{A}^2+p_2^{\left(1\right)}{B}^2=B^2+p_1^{\left(1\right)}(A_1^2-B_1^2)---\left(11\right)$

$p_1^{\left(1\right)}=\left[{\∫}_{s\∈S\left(1\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{-1.5}{\overset{-0.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c-b+\frac{13}{12}a+\frac{n}{2})/n---\left(12\right)$

其中，分别是第一区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第一区域面积的百分比，其中A₁为一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₁是二次曲线边缘模型下方区域的灰度值，(a,b,c)是边缘的参数，S(1)代表第一区域的面积，灰度值A₁＝A,B₁＝B。

S32：对于第二区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(2\right)}A+p_2^{\left(2\right)}B=B+p_1^{\left(2\right)}(A_2-B_2)---\left(13\right)$

$M_2^{\left(2\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(2\right)}{A}^2+p_2^{\left(2\right)}{B}^2=B^2+p_1^{\left(2\right)}(A_2^2-B_2^2)---\left(14\right)$

$p_1^{\left(2\right)}=\left[{\∫}_{s\∈S\left(2\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{-0.5}{\overset{-0.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c+b+\frac{1}{12}a+\frac{n}{2})/n---\left(15\right)$

其中，分别是第二区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第二区域面积的百分比，其中A₂为一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₂是二次曲线边缘模型下方区域的灰度值，(a,b,c)是边缘的参数，S(2)代表第二区域的面积，灰度值A₂＝A,B₂＝B；

S33：对于第三区域的一阶灰度矩，二阶灰度矩分别为，

$M_1^{\left(3\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f(i,j)=p_1^{\left(3\right)}A+p_2^{\left(3\right)}B=B+p_1^{\left(3\right)}(A_3-B_3)---\left(16\right)$

$M_2^{\left(3\right)}=\frac{1}{n}\underset{i}{\mathrm{\Σ}}f{(i,j)}^2=p_1^{\left(3\right)}{A}^2+p_2^{\left(3\right)}{B}^2=B^2+p_1^{\left(3\right)}(A_3^2-B_3^2)---\left(17\right)$

$p_1^{\left(3\right)}=\left[{\∫}_{s\∈S\left(3\right)}({\mathrm{ax}}^2+\mathrm{bx}+c)\mathrm{dx}\right]/n=\left[\underset{0.5}{\overset{1.5}{\∫}}({\mathrm{ax}}^2+\mathrm{bx}+c+\frac{n}{2})\mathrm{dx}\right]/n=(c+b+\frac{13}{12}a+\frac{n}{2})/n---\left(18\right)$

其中，分别是第三区域中二次曲线边缘模型下方所对应的面积和二次曲线边缘模型上方所对应的面积占第三区域面积的百分比，其中A₃一维阶跃边缘模型中二次曲线边缘模型上方区域的灰度值，B₃是二次曲线边缘模型下方区域的灰度值(a,b,c)是边缘的参数，S(3)代表第三区域的面积，灰度值A₃＝A,B₃＝B；

S34：利用上述每个区域所对应的一阶灰度矩和二阶灰度矩计算结果和公式(10)～(18)得到像素的二次曲线边缘模型的参数(a,b,c),具体的算法公式详见式(27)～(31)：

$a=\frac{1}{4}(\frac{M_2^{\left(1\right)}+M_2^{\left(3\right)}-{2M}_2^2}{(A^2-B^2)}+\frac{M_1^{\left(1\right)}+M_1^{\left(3\right)}{-2M}_1^2}{(A-B)})---\left(27\right)$

$b=\frac{1}{4}(\frac{M_2^{\left(3\right)}-M_2^{\left(1\right)}}{(A^2-B^2)}+\frac{M_1^{\left(3\right)}-M_1^{\left(1\right)}}{2(A-B)})---\left(28\right)$

$c=\frac{1}{2}(\frac{{2M}_2^2-n(A^2+B^2)}{2(A^2-B^2)}+\frac{{2M}_1^2-n(A+B)}{2(A-B)})-\frac{1}{12}a---\left(29\right)$

$A=\frac{T_1+\sqrt{{T_1}^2-{4T}_2}}{2}---\left(30\right)$

$B=\frac{T_1-\sqrt{{T_1}^2-{4T}_2}}{2}---\left(31\right)$

其中分别是第一区域，第二区域，第三区域的一阶灰度矩，分别是第一区域，第二区域，第三区域的二阶阶灰度矩，中间参数值T₁，T₂为 $T_1=\frac{M_2^{\left(3\right)}-M_2^{\left(1\right)}}{M_1^{\left(3\right)}-M_1^{\left(1\right)}},T_2=\frac{T_1{M}_1^{\left(2\right)}-M_2^{\left(2\right)}}{n}.$ 区域长度n＝9。