基于经验模态分解的边缘检测方法

摘要

本发明公开了一种基于经验模态分解的边缘检测方法，主要解决现有技术在噪声环境下，不能够很好地检测到清晰完整的图像边缘和虚假边缘较多的问题。其技术特点是：(1)在经验模态分解中，利用求解两个偏微分方程得到图像的极大值包络和极小值包络；(2)由图像的极大值包络和极小值包络得到图像的平均包络和差值包络；(3)继续不断地对步骤(1)和(2)迭代，直到满足迭代停止条件时，得到图像的固有模态函数和剩余图像；(4)对得到的剩余图像用两个Prewitt算子计算其梯度和门限，进而得到图像的边缘。本发明与传统Prewitt算子和Canny算子的检测效果相比，能得到更清晰完整的图像边缘，同时减少了虚假边缘以及噪声对边缘检测的影响。

说明书

技术领域

本发明属于图像处理领域，涉及边缘检测，具体地说是将经验模态分解(EMD)用于图像分解，将分解得到的剩余图像用边缘检测算子来检测其边缘，该方法可用于目标识别。

背景技术

边缘是图像局部灰度变化的不连续部分，是图像中灰度的急剧变化，主要存在于目标与目标、目标与背景、区域与区域之间，是图像分割、纹理特征提取和形状特征提取等图像分析的重要基础。图像分析的第一步常常是边缘检测，所以对边缘检测技术的研究十分重要。

边缘检测方法归纳起来分为三大类，第一类是经典的边缘检测方法，如微分算子法中的Sobel算子、Robert算子、Prewitt算子和Laplacian算子、最优算子法中的LOG算子和Canny算子；第二类是以能量最小化为准则的全局提取方法，其特征是运用严格的数学方法对此问题进行分析，给出一维代价函数作为最优提取依据，从全局最优的观点提取边缘，如松驰法；第三类是以小波变换、数学形态学、模糊数学、分形理论等近年来发展起来的新技术为基础的图像边缘提取方法，尤其是基于多尺度特性的小波变换提取图像边缘的方法是目前研究较多的课题。

实际应用中，图像数据往往含有噪声。因此，边缘检测方法要求既能检测到边缘的准确位置，又可以抑制无关的细节和噪声，而近年来出现的经验模态分解很适合用于处理非平稳的信号，而且可以用于图像的去噪，将得到的图像再进行边缘检测可以得到好的效果。

在1998年，时频分析领域出现了一种适合用于分析非线性和非平稳信号的信号处理方法Hilbert-Huang变换，该方法由Norden Huang提出，在这之后该方法迅速得到广泛的应用，这种方法的关键在于：任何复杂的数据都可以分解为一系列有限而且是少量的固有模态函数IMF，其中IMF定义为具有相同数量的零点和极值点而且有对称的极大值和极小值包络，还可以对每一个IMF进行Hilbert变换，进行时频分析。由于这种分解是自适应的，基于数据的局部特征时间尺度的，所以可以应用在非线性和非平稳数据处理过程中。近几年来EMD在各行各业的应用越来越多，如一维EMD用于机械损伤分析，二维EMD用于图像压缩，图像融合和图像识别等。

对于二维图像的EMD，如何得到图像的极大值点包络和极小值包络是至关重要的问题，已有的方法主要有基于Delaunay三角剖分EMD(DEMD)，基于径向基函数插值的EMD，自适应快速EMD，方向EMD和限邻域EMD等。但是DEMD由于对图像的极值点进行三角剖分，经常不能覆盖到全部的图像，所以会产生不确定的值，而且很容易产生过亮和过黑的点或者黑斑；而基于径向基插值的EMD虽然分解出的图像略优于DEMD，但是需要耗费巨大的计算量。

发明内容

本发明的目的在于针对上述问题，提出一种基于经验模态分解的边缘检测方法，以消除经验模态分解中过亮点和过黑点的产生，减少经验模态分解的计算量，提高边缘检测的清晰度。

实现本发明的技术关键是：利用求解两个偏微分方程得到在经验模态分解中求图像的极大值包络和极小值包络，然后对经验模态分解得到的剩余图像用Prewitt算子来检测边缘，其具体步骤包括如下：

(1)设置固有模态函数数目n＝0，初始化剩余图像r_n(x，y)为原图像f(x，y)，该剩余图像是指在经验模态分解中，从原图像f(x，y)减去固有模态函数之和得到的图像，x＝1，...h，y＝1，...w，x和y是图像的横坐标和纵坐标，h和w是图像的高和宽；

(2)对于图像r_n(x，y)，按8邻域求取其局部极大值点和局部极小值点，得到相对应的极大值标志矩阵IMax(x，y)与极小值标志矩阵IMin(x，y)；

(3)求剩余图像r_n(x，y)的极大值包络fmax(x，y)，即通过对下式的迭代得到fmax(x，y)＝f_t+1，设初始迭代次数t＝0，f_t(x，y)＝r_n(x，y)(为了方便，下式f_t(x，y)简记为f_t)，则

$f_{t+1}=f_t-\mathrm{dt}\times \mathrm{IMax}(x,y)\times \mathrm{sign}(\frac{{\partial }^4{f}_t}{{\partial x}^4}+\frac{{\partial }^4{f}_t}{{\partial y}^4})\sqrt{{\left(\frac{{\partial }^4{f}_t}{{\partial x}^4}\right)}^2+{\left(\frac{{\partial }^4{f}_t}{{\partial y}^4}\right)}^2}---1)$

其中，f_t为第t次迭代的包络，dt是迭代步长，sign()是符号函数，为f_t对x的四阶偏导数，为f_t对y的四阶偏导数，f_t+1为第t+1次迭代的包络，令t＝t+1，当t达到最大迭代次数maxiter_pde时，得到图像r_n(x，y)的极大值包络：fmax(x，y)＝f_t+1；

(4)求剩余图像r_n(x，y)的极小值包络fmin(x，y)，即通过对下式的迭代得到fmin(x，y)＝f_t+1，初始迭代次数t＝0，f_t(x，y)＝r_n(x，y)(为了方便，下式f_t(x，y)简记为f_t)，则

$f_{t+1}=f_t-\mathrm{dt}\times \mathrm{IMin}(x,y)\times \mathrm{sign}(\frac{{\partial }^4{f}_t}{{\partial x}^4}+\frac{{\partial }^4{f}_t}{{\partial y}^4})\sqrt{{\left(\frac{{\partial }^4{f}_t}{{\partial x}^4}\right)}^2+{\left(\frac{{\partial }^4{f}_t}{{\partial y}^4}\right)}^2}---2)$

其中，f_t为第t次迭代的包络，dt是迭代步长，sign()是符号函数，为f_t对x的四阶偏导数，为f_t对y的四阶偏导数，f_t+1为第t+1次迭代的包络，令t＝t+1，当t达到最大迭代次数maxiter_pde时得到图像r_n(x，y)的极小值包络：fmin(x，y)＝f_t+1；

(5)根据之前步骤得到的极小值包络fmin(x，y)和极大值包络fmax(x，y)，求得均值包络frean(x，y)和差值包络h₁(x，y)；

(6)用差值包络h₁(x，y)代替步骤(2)中的图像r_n(x，y)，重复步骤(2)-(5)，依次得到差值包络h₂(x，y)，h₃(x，y)，...h_k(x，y)，直到k达到最大迭代次数max_iter，设置固有模态函数数目n＝n+1，得到固有模态函数imf_n(x，y)＝h_k(x，y)和剩余图像r_n(x，y)；

(7)重复步骤(2)-(6)，当n达到最大固有模态函数层数max imf时，依次得到固有模态函数imf₂(x，y)，imf₃(x，y)...imf_n(x，y)和剩余图像r_n(x，y)；

(8)将剩余图像r_n(x，y)与两个Prewitt算子做卷积，分别得到剩余图像垂直梯度bx(x，y)和水平梯度by(x，y)，其中两个Prewitt算子是：

$p1=\left(\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right),p2=\left(\begin{array}{ccc}1& 0& -1\\ 1& 0& -1\\ 1& 0& -1\end{array}\right)$

(9)根据上步骤得到的垂直梯度bx(x，y)和水平梯度by(x，y)，计算剩余图像梯度b(x，y)和门限thresh：

b(x，y)＝bx(x，y)×bx(x，y)+by(x，y)×by(x，y)    3)

$\mathrm{thresh}=\frac{\mathrm{scale}}{h\times w}\times \underset{x=1}{\overset{h}{\Sigma }}\underset{y=1}{\overset{w}{\Sigma }}b(x,y)---4)$

其中scale＝4，当梯度b(x，y)大于门限thresh则认为边缘存在，设置edge(x，y)＝1，否则edge(x，y)＝0，edge(x，y)为图像边缘。

本发明具有以下优点：

(1)本发明通过求解两个偏微分方程得到经验模态分解中的极大值包络和极小值包络，由于偏微分方程可以用差分迭代法求解，因而计算复杂度比现有基于径向基函数插值的经验模态分解方法和现有基于Delaunay三角剖分的经验模态分解方法低很多；

(2)由于图像的极小值包络和极大值包络是通过周围的点差分迭代得到，不是由过远的点插值得到的，所以本发明的经验模态分解得到的图像比其他算法在清晰度方面有较大的提高，不会产生很模糊的图像，不会出现过亮点和过黑点，而且不需要对图像的边界设置，在边界处也不会产生黑色的区域。

(3)由于本发明采用的经验模态分解能够用于去除噪声，所以再对剩余图像进行边缘检测可以有效地抑制噪声的影响，比传统边缘检测方法Prewitt算子和Canny算子，得到更清晰准确的边缘和更少的虚假边缘。

附图说明

图1是本发明的总流程图；

图2是本发明中经验模态分解中求极大值包络和极小值包络的子流程图；

图3是本发明用于边缘检测的图像；

图4是利用现有Delaunay三角剖分的经验模态分解的图像分解结果；

图5是用本发明中经验模态分解的图像分解结果；

图6是在高斯噪声下本发明方法与现有Prewitt算子、Canny算子的边缘检测结果比较图；

图7是在椒盐噪声下本发明方法与现有Prewitt算子、Canny算子的边缘检测结果比较图。

具体实施方式

参照图1，本发明的具体实现步骤包括如下：

步骤1.设置固有模态函数数目n＝0，初始化剩余图像r_n(x，y)为原图像f(x，y)，该剩余图像是指在经验模态分解中，从原图像减去固有模态函数之和得到的图像。

步骤2.计算剩余图像r_n(x，y)的极大值包络和极小值包络。

参照图2，本步骤的具体实现如下：

2.1)求剩余图像r_n(x，y)的极大值标志矩阵和极小值标志矩阵。

(2.1.1)初始化与剩余图像r_n(x，y)同样大小的极大值标志矩阵IMax(x，y)和极小值标志矩阵IMin(x，y)，矩阵值都设为0，x＝1，...h，y＝1，...w，x和y是图像的横坐标和纵坐标，h和w是图像的高和宽；

(2.1.2)比较r_n(x，y)与其八邻域r_n(i，j)(i＝x-1，x，x+1，j＝y-1，y，y+1)的大小，如果r_n(x，y)比任何一个都大，设置极大值标志矩阵IMax(x，y)＝1；如果比任何一个都小，设置极小值标志矩阵IMin(x，y)＝1；如果(x，y)在图像边界，则比较r_n(x，y)与其五邻域或者三邻域的点的大小；

2.2)求剩余图像r_n(x，y)的极大值包络fmax(x，y)。

(2.2.1)初始迭代次数t＝0，设置f_t＝r_n(x，y)；

(2.2.2)计算f_t(x，y)对x的二阶偏导数和f_t(x，y)对y的二阶偏导数$\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y):$

$\frac{{\partial }^2{f}_t}{{\partial x}^2}(x,y)=[f_t(x+1,y)+f_t(x-1,y)-2\times f_t(x,y)]---5)$

$\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y)=[f_t(x,y+1)+f_t(x,y-1)-2\times f_t(x,y)]---6)$

其中，f_t(x，y)是第t次迭代(x，y)点的包络；

(2.2.3)计算f_t(x，y)对x的四阶偏导数和f_t(x，y)对y的四阶偏导数$\frac{{\partial }^4{f}_t}{{\partial y}^4}(x,y):$

$\frac{{\partial }^4{f}_t}{{\partial x}^4}(x,y)=[\frac{{\partial }^2{f}_t}{{\partial x}^2}(x+1,y)+\frac{{\partial }^2{f}_t}{{\partial x}^2}(x-1,y)-2\times \frac{{\partial }^2{f}_t}{{\partial x}^2}(x,y)]---7)$

$\frac{{\partial }^4{f}_t}{{\partial y}^4}(x,y)=[\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y+1)+\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y-1)-2\times \frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y)]---8)$

(2.2.4)计算下一次迭代的极大值包络f_t+1：

$f_{t+1}=f_t-\mathrm{dt}\times \mathrm{IMax}(x,y)\times \mathrm{sign}(\frac{{\partial }^4{f}_t}{{\partial x}^4}+\frac{{\partial }^4{f}_t}{{\partial y}^4})\sqrt{{\left(\frac{{\partial }^4{f}_t}{{\partial x}^4}\right)}^2+{\left(\frac{{\partial }^4{f}_t}{{\partial y}^4}\right)}^2}---9)$

其中，f_t为第t次迭代的极大值包络，dt是迭代步长，sign()是符号函数(上式中的f_t实际为f_t(x，y)，为了记号方便简记为f_t)；

(2.2.5)令t＝t+1，按步骤(2.2.2)-步骤(2.2.4)迭代，当t达到最大迭代次数maxiter_pde时，得到剩余图像r_n(x，y)的极大值包络：fmax(x，y)＝f_t+1。

2.3)求剩余图像r_n(x，y)的极小值包络fmin(x，y)。

(2.3.1)初始迭代次数t＝0，f_t＝f_n(x，y)；

(2.3.2)计算f_t(x，y)对x的二阶偏导数和f_t(x，y)对y的二阶偏导数$\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y):$

$\frac{{\partial }^2{f}_t}{{\partial x}^2}(x,y)=[f_t(x+1,y)+f_t(x-1,y)-2\times f_t(x,y)]---10)$

$\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y)=[f_t(x,y+1)+f_t(x,y-1)-2\times f_t(x,y)]---11)$

其中，f_t(x，y)是第t次迭代(x，y)点的包络；

(2.3.3)计算f_t(x，y)对x的四阶偏导数和f_t(x，y)对y的四阶偏导数$\frac{{\partial }^4{f}_t}{{\partial y}^4}(x,y):$

$\frac{{\partial }^4{f}_t}{{\partial x}^4}(x,y)=[\frac{{\partial }^2{f}_t}{{\partial x}^2}(x+1,y)+\frac{{\partial }^2{f}_t}{{\partial x}^2}(x-1,y)-2\times \frac{{\partial }^2{f}_t}{{\partial x}^2}(x,y)]---12)$

$\frac{{\partial }^4{f}_t}{{\partial y}^4}(x,y)=[\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y+1)+\frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y-1)-2\times \frac{{\partial }^2{f}_t}{{\partial y}^2}(x,y)]---13)$

(2.3.4)计算下一次迭代的极小值包络f_t+1：

$f_{t+1}=f_t-\mathrm{dt}\times \mathrm{IMin}(x,y)\times \mathrm{sign}(\frac{{\partial }^4{f}_t}{{\partial x}^4}+\frac{{\partial }^4{f}_t}{{\partial y}^4})\sqrt{{\left(\frac{{\partial }^4{f}_t}{{\partial x}^4}\right)}^2+{\left(\frac{{\partial }^4{f}_t}{{\partial y}^4}\right)}^2}---14)$

其中，f_t为第t次迭代的极小值包络，dt是迭代步长，sign()是符号函数(上式中的f_t实际为f_t(x，y)，为了记号方便简记为f_t)；

(2.3.5)令t＝t+1，按步骤(2.3.2)-步骤(2.3.4)迭代，当t达到最大迭代次数maxiter_pde时，得到剩余图像r_n(x，y)的极小值包络：fmin(x，y)＝f_t+1。

步骤3.计算剩余图像的均值包络fmean(x，y)和差值包络h₁(x，y)：

fmean(x，y)＝(fmax(x，y)+fmin(x，y))/2    15)

h₁(x，y)＝f(x，y)-fmean(x，y)。           16)

步骤4.检查内层迭代是否满足迭代停止条件和得到固有模态函数，即用h₁(x，y)代替步骤2中的剩余图像r_n(x，y)，按步骤2-步骤3迭代，依次得到h₂(x，y)，h₃(x，y)，...h_k(x，y)，直到k达到最大迭代次数max_iter时停止迭代，设置固有模态函数数目n＝n+1，得到固有模态函数imf_n(x，y)和剩余图像r_n(x，y)：

imf_n(x，y)＝h_k(x，y)        17)

r_n(x，y)＝f(x，y)-imf_n(x，y)18)

步骤5.检查外层迭代是否满足迭代停止条件和更新剩余图像，即重复步骤2-步骤4，当固有模态函数数目n达到最大固有模态函数层数max_imf时停止迭代，依次得到固有模态函数imf₂(x，y)，imf₃(x，y)...imf_n(x，y)和剩余图像r_n(x，y)。

步骤6.计算剩余图像梯度和门限。

6.1)将剩余图像r_n(x，y)与两个Prewitt算子做卷积，分别得到剩余图像垂直梯度bx(x，y)和水平梯度by(x，y)：

$\mathrm{bx}(x,y)=\underset{i=x-1}{\overset{x+1}{\Sigma }}\underset{j=y-1}{\overset{y+1}{\Sigma }}{r}_n(i,j)\times p1(i-x+2,j-y+2)---19)$

$\mathrm{by}(x,y)=\underset{i=x-1}{\overset{x+1}{\Sigma }}\underset{j=y-1}{\overset{y+1}{\Sigma }}{r}_n(i,j)\times p2(i-x+2,j-y+2)---20)$

其中，两个Prewitt算子如下：

$p1=\left(\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right),p2=\left(\begin{array}{ccc}1& 0& -1\\ 1& 0& -1\\ 1& 0& -1\end{array}\right)$

6.2)根据剩余图像垂直梯度bx(x，y)和水平梯度by(x，y)，计算剩余图像梯度b(x，y)和门限thresh：

b(x，y)＝bx(x，y)×bx(x，y)+by(x，y)×by(x，y)    21)

$\mathrm{thresh}=\frac{\mathrm{scale}}{h\times w}\times \underset{x=1}{\overset{h}{\Sigma }}\underset{y=1}{\overset{w}{\Sigma }}b(x,y)---22)$

其中，scale＝4。

步骤7.确定检测图像边缘。

当剩余图像梯度b(x，y)大于门限thresh则认为边缘存在，设置edge(x，y)＝1，否则edge(x，y)＝0，edge(x，y)为被检测图像的边缘。

本发明的效果可通过以下实验仿真进一步说明：

1.仿真条件

本实验用cameraman图像作为边缘检测对象，如图3所示。首先比较本发明的经验模态分解和现有Delaunay三角剖分的经验模态分解的效果，然后分别向cameraman图像添加高斯噪声和椒盐噪声，比较本发明的边缘检测方法与Prewitt算子和Canny算子的边缘检测效果。在经验模态分解中设置最大固有模态函数层数max_imf＝3，偏微分方程迭代步长dt＝0.01，偏微分方程最大迭代次数maxiter_pde＝20，求固有模态函数最大迭代次数max_iter＝3，设置高斯噪声的幅度为0.02和椒盐噪声的幅度为0.09。

2.仿真内容和结果分析

2.1)分别对现有基于Delaunay三角剖分的经验模态分解和本发明的经验模态分解进行仿真。

基于Delaunay三角剖分的经验模态分解仿真结果如图4所示，其中图4(a)为固有模态函数1，图4(b)为固有模态函数2，图4(c)为固有模态函数3，图4(d)为剩余图像。

本发明的经验模态分解仿真结果如图5所示，其中图5(a)为固有模态函数1，图5(b)为固有模态函数2，图5(c)为固有模态函数3，图5(d)为剩余图像。

从图5可以看出，本发明的经验模态分解得到的图像相对比较清晰，从图4中看出，除图4(a)清晰的表示图像的边缘外，图4(b)和图4(c)都很模糊，而且有过黑和过亮的点，而这在本发明的图5(a)、图5(b)、图5(c)中是没有的，从图4(c)、图4(d)中看出，基于Delaunay三角剖分的经验模态分解的图像总会在边缘出现黑色的块，而这在本发明方法中是不会出现的。

2.2)对比分析在高斯噪声下，本发明方法和现有Prewitt算子、Canny算子的边缘检测仿真效果，仿真结果如图6所示，其中图6(a)为用本发明方法的检测结果，图6(b)为用Prewitt算子的检测结果，图6(c)为用Canny算子的检测结果。从图6(c)可以看出Canny算子检测到很多虚假边缘，图6(b)的Prewitt算子能够很好地检测到边缘，而本发明的图6(a)比用Prewitt算子得到的边缘更好，噪声也少很多。

2.3)对比分析在椒盐噪声下，本发明方法和现有Prewitt算子、Canny算子的边缘检测仿真效果，仿真结果如图7所示，其中图7(a)为本发明方法的检测结果，图7(b)为Prewitt算子的检测结果，图7(c)为Canny算子的检测结果。从图7(c)可以看出Canny算子仍然检测到很多虚假边缘，图7(b)的Prewitt算子能够很好地检测到边缘，但是噪声点相对较多，而本发明的图7(a)比用Prewitt算子得到的边缘更好，噪声和虚假边缘相对更少，但是比在高斯噪声下的本发明方法边缘检测效果稍差。