一种基于Legendre矩的图像亚像素边缘检测方法

摘要

本发明公开了一种基于Legendre矩的图像亚像素边缘检测方法，包括以下步骤：S1:读取图像信息，将图像灰度化并对灰度图像进行去噪处理；S2:采用Sobel算子对去噪后的图像进行像素级边缘定位：利用像素点的各向邻点灰度加权值在边缘点达到最大值这一现象进行边缘检测；S3:采用Legendre矩对图像进行亚像素边缘检测，输出边缘图像。其中Sobel算子对噪声具有平滑作用，提供较为准确的边缘方向信息，利用Legendre矩进行亚像素边缘检测，减少了运算所需要的模板的数量，降低了计算的复杂度，同时在抗噪方面具有更好的鲁棒性。

权利要求

1.一种基于Legendre矩的图像亚像素边缘检测方法，其特征在于包括以下步骤：S1:读取图像信息，将图像灰度化并对灰度图像进行去噪处理；

S2:采用Sobel算子对去噪后的图像进行像素级边缘定位：利用像素点的各向邻点灰度加权值在边缘点达到最大值这一现象进行边缘检测；

S3:采用Legendre矩对图像进行亚像素边缘检测，输出边缘图像；

所述步骤S3中具体采用如下方式：遍历检测到的所有边缘点，对每个边缘点进行如下处理：以得到的边缘点为中心，在灰度图像中选取N×N的窗口,N为奇数，采用如下公式(25)将N×N灰度图像窗口中的值与Legendre正交矩的掩码CLM11对应位置的系数相乘得到N×N的矩阵，将该矩阵求和得到Legendre正交矩LM11，同样的方式利用式(26)再求得一个Legendre正交矩LM31，其中，f(m,n)是像素边缘检测的位置的灰度值；

采用如下公式(18)求出 值：

其中 为亚像素边缘点的角度，

利用亚像素边缘点的角度 和如下公式(21)和(22)计算出亚像素边缘点离中心的位置l的值：其中,

利用如下公式(27)得到图像的亚像素边缘位置：

其中，x,y是Sobel算子进行检测得到的边缘点的位置,N代表掩码的窗口大小。

2.根据权利要求1所述的一种基于Legendre矩的图像亚像素边缘检测方法，其特征还在于：S2中具体采用如下方式：遍历原始灰度图像中所有的像素点，计算得到每个像素点的梯度值G[f′(x,y)]，将所得的梯度值归一化到[0,255]区间，采用最大类间方差法计算得到归一化梯度值的阈值T，对每个像素点归一化的梯度值进行判断，即当G[f′(x,y)]>T时，对应的像素点设定为255，否则设定为0至此得到图像的像素级粗定位。

3.根据权利要求2所述的一种基于Legendre矩的图像亚像素边缘检测方法，其特征还在于：在得到像素级粗定位后遍历图像中所有的边缘点，进行判断：若该边缘点是孤立的边缘点即在以该点为中心的3×3的邻域中、除本点外的边缘点的个数小于等于1，则将该点除去，即该点不作为边缘点、判断为噪声。

说明书

一种基于Legendre矩的图像亚像素边缘检测方法

技术领域

[0001] 本发明涉及图像边缘检测领域，尤其涉及一种基于Legendre矩的图像亚像素边缘检测方法。

背景技术

[0002] 工程中常基于图像实现非接触的几何尺寸精密测量，该方式以其非接触、全视场、高精度的特点而获得广泛应用。其原理就是通过处理被测物体图像的边缘而获得图像的几何参数。由此可见，边缘检测是图像测量的基础和关键。传统边缘检测方法多是基于图像像素灰度的变化，如Sobel算子、Laplacian算子和canny算子等。这些方法形式简单、易于实现、但定位精度不高，通常只有整数像素级的精度，且微分算子对噪声非常敏感，常会产生一些伪边缘。随着人们对检测精度要求的不断提高，像素级检测精度已经不能满足实际测量的要求。为了解决这个问题，人们提出了亚像素边缘检测方法。这些方法可以突破摄像机物理分辨率的限制，使图像的边缘定位精度达到亚像素级，从而提高了图像测量系统的检测精度。当算法的精度为0.1个像素，则相当于检测系统的硬件分辨率提高了10倍。目前亚像素边缘检测方法，在数学上可以归纳为插值法、拟合法、矩方法三种类型。拟合法通过对给定的边缘模型将图像中的灰度值进行拟合，这种方法具有很高的精度但是耗时，插值法通过对实际图像的灰度分布进行插值来得到亚像素的位置，但是对噪声很敏感。矩方法使用了对噪声不敏感的积分算子，因此应用最广泛。

[0003] 在文献[1]“Subpixel edge location based on orthogonal Fourier–Mellin moments”中，Bin提出了基于OFMM矩的亚像素边缘检测，该技术利用傅里叶-梅林矩对图像进行亚像素级的定位，采用了5×5的模板，求得亚像素坐标的四个参数，φ、l、k、h，然后对h进行判断，若大于阈值T，则判断该点为边缘点。该技术通过使用OFMM矩进行亚像素边缘检测，虽然该技术对噪声不敏感，克服了噪声的影响，但是由于需要三个实数模板，两个复数模板，计算复杂度大，影响求解速度。

发明内容

[0004] 根据现有技术存在的问题，本发明公开了一种基于Legendre矩的图像亚像素边缘检测方法，包括以下步骤：

[0005] S1:读取图像信息，将图像灰度化并对灰度图像进行去噪处理；

[0006] S2:采用Sobel算子对去噪后的图像进行像素级边缘定位：利用像素点的各向邻点灰度加权值在边缘点达到最大值这一现象进行边缘检测；

[0007] S3:采用Legendre矩对图像进行亚像素边缘检测，输出边缘图像。

[0008] S2中具体采用如下方式：遍历原始灰度图像中所有的像素点，计算得到每个像素点的梯度值G[f′(x,y)]，将所得的梯度值归一化到[0,255]区间，采用最大类间方差法计算得到归一化梯度值的阈值T，对每个像素点归一化的梯度值进行判断，即当G[f′(x,y)]>T时，对应的像素点设定为255，否则设定为0至此得到图像的像素级粗定位。

[0009] 进一步的，在得到像素级粗定位后遍历图像中所有的边缘点，进行判断：若该边缘点是孤立的边缘点即在以该点为中心的3×3的邻域中、除本点外的边缘点的个数小于等于1，则将该点除去，即该点不作为边缘点、判断为噪声。

[0010] 进一步的，S3中具体采用如下方式：遍历检测到的所有边缘点，对每个边缘点进行如下处理：以得到的边缘点为中心，在灰度图像中选取N×N的窗口,N为奇数，采用如下公式(25)将N×N灰度图像窗口中的值与Legendre正交矩的掩码CLM11对应位置的系数相乘得到N×N的矩阵，将该矩阵求和得到Legendre正交矩LM11，同样的方式利用式(26)再求得一个Legendre正交矩LM31，

[0011]

[0012]

[0013] 其中，f(m,n)是像素边缘检测的位置的灰度值；

[0014] 采用如下公式(18)求出 值：

[0015]

[0016] 其中 为亚像素边缘点的角度，

[0017] 利用亚像素边缘点的角度 和如下公式(21)和(22)计算出亚像素边缘点离中心的位置l的值：

[0018]

[0019] 其中,

[0020] 利用如下公式(27)得到图像的亚像素边缘位置：

[0021]

[0022] 其中，x,y是Sobel算子进行检测得到的边缘点的位置,N代表掩码的窗口大小。

[0023] 由于采用了上述技术方案，本发明提供的基于Legendre矩的图像亚像素边缘检测方法，首先对输入图像进行灰度化并采用自适应中值滤波器对图像进行去噪处理，然后用Sobel算子进行像素级边缘粗定位，最后利用Legendre矩进行图像的亚像素级的边缘检测，其中Sobel算子对噪声具有平滑作用，提供较为准确的边缘方向信息，利用Legendre矩进行亚像素边缘检测，减少了运算所需要的模板的数量，降低了计算的复杂度，同时在抗噪方面具有更好的鲁棒性。

附图说明

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

[0025] 图1为本发明的图像亚像素边缘检测方法的流程图；

[0026] 图2为本发明对图像进行去噪处理时边界扩展示意图；

[0027] 图3为本发明中理想2D边缘模型的示意图；

[0028] 图4(a)旋转前的阶跃边缘模型的示意图；

[0029] 图4(b)旋转后的阶跃边缘模型的示意图；

[0030] 图5为本发明中模板系数计算模型的示意图。

具体实施方式

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

[0032] 如图1所示的基于Legendre矩的图像亚像素边缘检测方法，具体包括以下步骤：读取图像信息，将图像灰度化并对灰度图像进行去噪处理：首先将输入RGB图像转换为灰度图像。采用如下方式：

[0033] Gray＝(28×B+151×G+77×R)>>8   (1)

[0034] 其中，“>>”表示二进制右移。R、G、B代表红、绿、蓝三个通道的颜色。遍历输入图像的所有像素点，对每个像素点采用式(1)，Gray为得到的灰度图像对应像素点的灰度值。

[0035] 对图像进行去噪声处理：

[0036] 本发明采用自适应中值滤波器，自适应中值滤波器可以处理含有大概率的脉冲噪声图像，在平滑非脉冲噪声时能保留细节。自适应中值滤波器工作于矩形窗口区域Sxy内，与传统的滤波器不同的是自适应中值滤波器在进行滤波处理时会根据某些条件而改变Sxy的尺寸。

[0037] 具体操作步骤：本发明中Smax＝10，首先是边界扩展，在图像上下左右各增加Smax个像素，如图2所示，设原图像的大小为m×n，区域1,2,3,4均是扩展的区域。区域1和区域2的大小为Smax×n，区域3和区域4的大小为(m+2×Smax)×Smax，首先扩展区域1，将原图像最左边与区域1对应大小的区域的数据复制到区域1，同理，扩展区域2，将原图像最右边与区域2对应大小的区域的数据复制到区域2。接着，扩展区域3，将原图像和区域1和区域2最上边(m+2×Smax)×Smax对应的区域的数据复制到区域3，同理，扩展区域4.得到的扩展图像进行后续的滤波处理。

[0038] 初始滤波半径r＝1，对应的初始矩形窗口Sxy的大小为(2r+1)×(2r+1),(即为3×3)，算法由两个进程表示，分别为进程A和进程B，依次从原图像的像素点开始遍历(不包括填充的像素点)。

[0039] Zmin表示在矩形窗口Sxy中的最小灰度值，Zmax表示在矩形窗口Sxy中的最大灰度值，Zmed表示在矩形窗口Sxy中的灰度中值，Zxy表示矩形窗口Sxy中心位置像素的灰度值，Smax为Sxy矩形窗口的最大滤波半径，在本文算法中，Smax＝10。

[0040] 进程A：A1＝Zmed-Zmin

[0041] A2＝zmed-zmax

[0042] 如果A1>0且A2<0，则转到进程B，否则滤波半径r＝r+1

[0043] 如果矩形窗口的滤波半径<＝Smax,则重复进程A

[0044] 否则输出Zmed

[0045] 进程B：B1＝Zxy-Zmin

[0046] B2＝Zxy-Zmax

[0047] 如果B1>0且B2<0，则输出Zxy

[0048] 否则输出Zmed

[0049] 输出的Zmed即为该像素点经过自适应中值滤波器滤波后的像素值。

[0050] 采用Sobel算子对去噪后的图像进行像素级边缘定位：利用像素点上下，左右邻点的灰度加权算法，根据在边缘点处达到极值现象进行图像边缘的检测。

[0051] Sobel算子很容易在空间上实现，Sobel边缘检测器不但产生较好的边缘检测效果，而且受噪声影响也比较小。Sobel算子利用像素点上下，左右邻点的灰度加权算法，根据在边缘点处达到极值这一现象进行边缘的检测。Sobel算子对噪声具有平滑作用，提供较为准确的边缘方向信息。

[0052] fx′(x,y)＝f(x-1,y+1)+2f(x,y+1)+f(x+1,y+1)

[0053] -f(x-1,y-1)-2f(x,y-1)-f(x+1,y-1)   (2)

[0054] fy′(x,y)＝f(x-1,y-1)+2f(x-1,y)+f(x-1,y+1)

[0055] -f(x+1,y-1)-2f(x+1,y)-f(x+1,y+1)   (3)

[0056] G[f′(x,y)]＝|fx′(x,y)|+|fy′(x,y)|   (4)

[0057] 其中，fx’(x,y)、fy’(x,y)分别是x(水平)方向和y(垂直)方向的一阶微分，G[f′(x,y)]为Sobel算子的梯度总和，f(x,y)是输入图像在坐标(x,y)点处的灰度值。

[0058] 在阈值T的设定方面，采用最大类间方差法(也称大津法，Otsu法)，它的主要思想是按照灰度特性将图像划分为背景和目标2部分，划分依据为选取门限值，使得背景和目标之间的方差最大。其主要实现原理如下：

[0059] 1)建立图像灰度直方图(共有L个灰度级，每个出现概率为p，ni为灰度值为i的像素的个数)

[0060]

[0061]

[0062] 2)计算背景和目标的出现概率，计算方法如下：

[0063]

[0064]

[0065] 其中，t为所选定的阈值，A代表背景(灰度级为0-t)，PA为背景出现的概率，同理B代表目标(灰度级为t+1-L-1)，PB为目标出现的概率。

[0066] 3)计算A和B两个区域的类间方差如下：

[0067]

[0068]

[0069] σ2＝PA(ωA-ωo)2+PB(ωB-ωo)2   (11)

[0070] 式(9)分别计算A和B区域的平均灰度值,ωA代表A区域的平均灰度值，ωB代表B区域的平均灰度值；式(10)计算灰度图像全局的灰度平均值ωo；式(11)计算A、B两个区域的类间方差σ2。

[0071] 4)以上几个步骤计算出了单个灰度值上的类间方差，因此最佳分割门限值应该是图像中能够使得A与B的类间灰度方差最大的灰度值。在程序中t的取值从0到255，依次计算式(11)的值，最大的σ2所对应的t值即为阈值T。在得到像素级粗定位后，遍历图像中所有的边缘点，并进行判断，若该边缘点是孤立的边缘点(即在以该点为中心的3×3的邻域中(除本点外)边缘点的个数小于等于1，)，则将该点除去，即该点不作为边缘点，判断为噪声。

[0072] 采用Legendre矩对图像进行亚像素边缘检测，输出边缘图像。

[0073] 自从TABATABAI等在1984年提出了利用灰度矩进行亚像素级别的边缘检测，经过20多年的研究，其他方法例如空间矩，Zernike矩，OFMM等等被提出，这些方法假设理想的边缘是阶跃模型,通过将图像映射到单位圆内，求得亚像素的4个参数,l(亚像素离中心的位置), (亚像素的角度)，k(灰度的阶跃值)，h(背景灰度)。为此，提出了基于Legendre矩的亚像素边缘检测。

[0074] (1)Legendre矩

[0075]

[0076] 其中，

[0077] 在单位圆内，Legendre矩可以定义为：

[0078]

[0079] 其中 是归一化系数，f(r,θ)为原灰度图像在(x,y)点处的灰度值

[0080] f(x,y)在极坐标下的表示形式.其中

[0081]

[0082] 核函数Tnm＝Qn(r)exp(-jmθ)使LMnm具有旋转不变性。

[0083]

[0084] 其中，LMnm表示原始图像的Legendre矩，LM'nm表示对图像旋转 角后的Legendre矩。如图3所示。

[0085] 基于Legendre矩的边缘检测:如图4所示，在旋转 的角度后，边缘与y轴垂直，旋转后的图像函数的积分具有如下的关系：

[0086]

[0087] 基于

[0088]

[0089] 根据式(9)，可以得到LM1'1的虚部为0，所以

[0090]

[0091] Re[LM11]和Im[LM11]分别是LM11的实部和虚部。

[0092] 因此，

[0093] 积分核函数可以表示为：

[0094]

[0095]

[0096] 注：LM1'1,LM3'1的解法可见[2]

[0097]

[0098] 其中,

[0099] 计算LM11和LM31模板的系数，本发明使用5×5的模板

[0100]

[0101]

[0102] 如图5和式(23),(24)可知，计算的系数实部关于y轴奇对称，关于x轴偶对称，虚部关于x轴奇对称，关于y轴偶对称。因此，仅需要计算图5中，方格1,2,3,6,7,8,11,12这八个系数，剩下的系数可通过对称性得到。

[0103] 首先利用式(23)计算CLM11

[0104] 对于方格1：

[0105]

[0106]

[0107]

[0108]

[0109]

[0110]

[0111]

[0112]

[0113] 根据对称性，可得剩下的系数，详细见表1.

[0114] 同理，利用式(24)可计算CLM31模板的系数，见表2

[0115] 表1 CLM11模板系数

[0116]-0.0147+0.0147j -0.0469+0.0933j 0.125j 0.0469+0.0933j 0.0147+0.0147j
-0.0933+0.0469j -0.064+0.064j 0.064j 0.064+0.064j 0.0933+0.0469j
-0.1253 -0.064 0.0 0.064 0.1253
-0.0933-0.0469j -0.064-0.064j -0.064j 0.064-0.064j 0.0933-0.0469j
-0.0147-0.0147j -0.0469-0.0933j -0.125j 0.0469-0.0933j 0.0147-0.0147j[0117] 表2 CLM31模板系数

[0118]-0.01116-0.01116j -0.018301-0.03672j -0.02712j 0.018301-0.03672j 0.01116-0.01116j-0.03672-0.0183017j 0.036267+0.036267j 0.061866j -0.036267+0.036267j 0.03672-0.018301j-0.0271204 0.061866 0.0 -0.061866 0.0271204
-0.03672+0.0183017j 0.036267-0.036267j -0.061866j -0.036267-0.036267j 0.03672+0.018301j-0.01116+0.01116j -0.018301+0.03672j 0.02712j 0.018301+0.03672j 0.01116+0.01116j[0119] 然后利用以下的公式，求LM11和LM31

[0120]

[0121]

[0122] 其中，f(m,n)是像素边缘检测的位置的灰度值。

[0123] 真实的边缘位置为：

[0124]

[0125] 其中，x,y是Sobel算子进行检测得到的边缘点的位置,N代表掩码的窗口大小。

[0126] 本发明的有益效果说明：

[0127] 为了验证本发明，进行了计算机仿真实验。在实验中，实验参数为CPU英特尔Pentium(奔腾)双核E53002.6GHz，2GB内存，显卡是英特尔G33/G31Express Chipset Family,操作系统为Window XP专业版32位SP2，软件编程环境为Matlab2010b，本发明实验的图像是利用人工合成的图像，对于人工合成的图片的大小是256像素×256像素。

[0128] 从[3]中，Feipeng Da推导了SGM,ZOM和OFMM之间的关系，可以得出计算SGM,ZOM和OFMM所得的0是一样的，而ZOM和OFMM的l值是一样的，由SGM和ZOM所计算得到的l值的差值是

[0129]

[0130] 所以将发明方法与SGM和ZOM进行了仿真对比实验。测试的图片为添加高斯白噪声2
的不同半径的圆，圆心为(128,128)，通过将计算得到的亚像素点进行拟合，求出(X-A) +(Y-B)2＝R2中A，B，R的值，拟合采用[6]中所提到的方法。

[0131]

[0132]

[0133] k是亚像素边缘点的个数，xt，yt代表了第i个亚像素边缘点的坐标，半径定义为亚像素边缘点到实际圆心的平均距离.

[0134] 表3 不同方法的位置边缘误差(拟合圆心与实际圆心的误差，注，误差计算的是欧式距离)

[0135]半径 本发明方法 基于SGM的算法 基于ZOM的算法
70 0.0077 0.9839 0.9968
75 0.0228 0.2288 0.2456
80 0.0176 0.6112 0.6246
85 0.0062 0.0878 0.0747
90 0.1804 1.5499 1.5317
95 0.5487 0.9952 0.9768
100 0.5842 0.8027 0.8240
105 0.6539 1.0609 1.0614
110 0.3456 0.1479 0.1478

[0136] 表4 不同方法的位置边缘误差(拟合半径与实际半径的误差，注，误差计算的是两者的差值)

[0137]半径 本文采用的算法 基于SGM的算法 基于ZOM的算法
70 0.048 1.5045 1.5203
75 0.0282 0.8873 0.9027
80 0.0015 0.4416 0.4390
85 0.0246 0.38 0.3936
90 0.0276 0.6257 0.6241
95 0.0414 0.3122 0.3011
100 0.0048 0.2015 0.2222
105 0.0195 0.1851 0.1915
110 0.0347 0.1907 0.1988

[0138] 参考文献：

[0139] (如专利/论文/标准)

[0140] [1]Bin T J,Lei A,Jiwen C,et al.Subpixel edge location based on orthogonal Fourier–Mellin moments[J].Image and Vision Computing,2008,26(4):563-569.

[0141] [2]Cui J,Feng K,Tan J B.Further improvement of edge location accuracy of double fiber spherical coupling sensor using orthogonal Jacobi–Fourier moments[J].Optik-International Journal for Light and Electron Optics,2014,125(1):353-359.

[0142] [3]Da F,Zhang H.Sub-pixel edge detection based on an improved moment[J].Image and Vision Computing,2010,28(12):1645-1658.

[0143] [4]Lyvers E P,Mitchell O R,Akey M L,et al.Subpixel measurements using a moment-based edge operator[J].Pattern Analysis and Machine Intelligence,IEEE Transactions on,1989,11(12):1293-1309.

[0144] [5]Ghosal S,Mehrotra R.Orthogonal moment operators for subpixel edge detection[J].Pattern recognition,1993,26(2):295-306.

[0145] [6]Fabijanska A.Gaussian-based approach to subpixel detection of blurred and unsharp edges[C]//Computer Science and Information Systems(FedCSIS),2014Federated Conference on.IEEE,2014:641-650.

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