CN103440644B - 一种基于最小描述长度的多尺度图像弱边缘检测方法 - Google Patents

Abstract

本发明公开了一种基于最小描述长度（Minimum Description Length，MDL）原则的自适应多尺度图像弱边缘检测方法，首先，利用多尺度高斯平滑构造线性尺度空间；然后，计算图像的局部描述长度，利用最小描述长度原则确定最佳局部平滑尺度；最后，对局部平滑后的图像进行边缘检测，得到弱边缘图像的全部边缘。本发明能有效地滤除常见噪声和提取真实的弱边缘，避免边缘断裂、偏移和虚假边缘响应等现象，而且本发明中的算法为非迭代地滤波方式，不仅极大地提高运算速度，还使得算法的稳定性大为增强。

Description

一种基于最小描述长度的多尺度图像弱边缘检测方法

技术领域

本发明涉及一种图像弱边缘检测方法，具体来说，它涉及一种利用局部相干扩散和最小描述长度准则(minimum description length，MDL)保护图像弱边缘并估计局部平滑尺度的自适应平滑滤波方式。

背景技术

图像弱边缘缘于图像的运动模糊和散焦模糊，其中蕴含的重要空间信息在机器视觉、图像恢复和图像压缩等领域中得到日益广泛的应用。为了提取这种空间信息首先要对这些弱边缘进行检测，而传统的边缘检测方法容易发生弱边缘的误判和漏检现象。本发明提出了一种基于最小描述长度的多尺度图像弱边缘检测技术，能够提高弱边缘的提取效率。

最先提出多尺度边缘检测概念的是Witkin【1984】。作者使用多尺度边缘识别和边缘跟踪，并采用了由粗到细的跟踪策略用于特征提取。理论基础源于尺度空间理论中的尺度不变性原理，该方法的前提是真实边缘存在于任何尺度下，且能被检测到。于是只要找到真实边缘在不同尺度之间的位置变化的联系，就能准确定位到所有真实边缘。

这种方法的缺陷在于，在多尺度平滑过程中，相邻的真实的边缘发生融合，而在跨尺度跟踪时，像素点的匹配问题难以解决。

在此基础上，Bergholm【1987】提出了一种多尺度联合感知的边缘检测算法，使用了由粗到细的边缘跟踪方法，称为边缘聚焦。Bergholm的研究实质是Witkin方法的一种实用性扩展，使用的是二维信号而非一维信号作为研究对象，有助于解决离散尺度带来的边缘融合问题。

Mallat和Zhong【1999】，运用了小波多尺度特性，先将图像分解然后进行边缘检测。其依据是：第一、Canny检测方法类似于在小波变换域中求极值；第二、多尺度的边缘意味着对一幅图像进行唯一性的表征。

Lindeberg在1993首次提出尺度空间理论，并在1998年提出了一种多尺度边缘检测的自动阈值选取方法，在此基础上，Lindeberg对基于高斯平滑的边缘检测做了更为深入研究，分别于2011年，2012年提出了一系列尺度空间新理论，为多尺度边缘检测奠定了理论基础。

Rissanen【1999】，首次提出了最小长度描述原则(MDL)，因其广泛而有效的信息描述能力，迅速成为信息论和机器学习领域的重要概念，MDL在其他领域的应用也日益受到重视。

目前的多尺度边缘检测技术主要存在算法复杂度高，尺度确定不合理等问题。

发明内容

针对现有算法缺乏合理的尺度选取标准等问题，本发明确定了一种最佳局部尺度计算标准，鉴于此标准，提出一种基于最小描述长度原则的局部多尺度自适应图像弱边缘检测方法，它包括：

a)对图像的预处理。对图像进行灰度化与边缘增强，提高图像中的弱边缘的对比度。先将图像转化为灰度图像，然后采用相干增强算法对图像的弱边缘进行增强，其计算方法为：

a1)计算结构张量S，可由下式估计出它的局部方向

$S=\left(\begin{array}{cc}{s}_{11}& s_{12}\\ s_{12}& s_{22}\end{array}\right)=\left(\begin{array}{cc}\frac{{\∂}^{2}L}{\∂x^2}*G_{\mathrm{\σ}}& \frac{{\∂}^{2}L}{\∂x\∂y}*G_{\mathrm{\σ}}\\ \frac{{\∂}^{2}L}{\∂x\∂y}*G_{\mathrm{\σ}}& \frac{{\∂}^{2}L}{\∂y^2}*G_{\mathrm{\σ}}\end{array}\right)---\left(1\right)$

其中G_σ表示局部尺度为σ的高斯滤波器，结构张量的特征向量为局部梯度的方向，而两个特征值的差，代表着局部邻域的非均匀性。

a2)构造扩散张亮D，如下所示，

$D=R^T\left(\begin{array}{cc}{c}_1& 0\\ 0& c_2\end{array}\right)R---\left(2\right)$

其中R为旋转矩阵，它的列向量由结构张量的特征向量表示，其中c₁和c₂为沿梯度方向的引导相干系数。

扩散张量D的元素可采用以下方式计算：

$d_{11}=\frac{1}{2}(c_1+c_2+\frac{(c_2-c_1)(s_{11}-s_{22})}{\mathrm{\α}})---\left(3\right)$

$d_{12}=\frac{(c_2-c_1)s_{12}}{s_{12}}---\left(4\right)$

$d_{22}=\frac{1}{2}(c_1+c_2-\frac{(c_2-c_1)(s_{11}-s_{22})}{\mathrm{\α}})---\left(5\right)$

其中，

$\mathrm{\α}=\sqrt{{(s_{11}-s_{22})}^2+4s_{12}^2}---\left(6\right)$

结构张量的特征值为：

${\mathrm{\λ}}_{1,2}=\frac{1}{2}(s_{11}+s_{22}\±\mathrm{\α})---\left(7\right)$

在计算中，本发明利用两个特征向量的差(即邻域的非均匀程度)来控制扩散速度c₁和c₂：

$c_1=max(0.01,1-e^{-{({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)}^2/k^2})---\left(8\right)$

c₂＝0.01

b)构造线性尺度空间。具体包括：

b1)构造一个线性尺度空间，G_σ表示标准差为σ的高斯滤波器。对于一幅R×C灰度图像I，利用一系列尺度σ递增的高斯滤波器与图像I卷积，得到平滑后的图像序列，即I_σ＝I₀*G_σ，其中σ＝σ₁,...,σ_n。

b2)在尺度空间中，搜索每个位置(x,y)处的σ_i。

c)计算图像的局部描述长度；

d)采用最小描述长度原则确定局部区域的最佳平滑尺度。

d1)利用最小描述长度原则作为方差σ的选取标准。

将一个高斯平滑的过程，用下面的公式简单描述

其中ε_σ为残差，表示原图像与平滑后的图像之间差值。

d2)将最小描述长度的思想应用到图像的平滑，本发明可以理解为：平滑量最多而残差最小，即最大限度地滤除噪声，压缩源图像的冗余信息的同时，保持原图像和平滑后的图像之间的相似性。本发明可以用长度描述算子(dl)来表示上述平滑过程：

${\mathrm{dl}}_{I_0(x,y)}={\mathrm{dl}}_{I_{\mathrm{\σ}}(x,y)}+{\mathrm{dl}}_{{\mathrm{\ϵ}}_{\mathrm{\σ}}(x,y)}---\left(10\right)$

d3)计算I_σ的描述长度

对于任意的σ大于0，当σ→∞时，I_σ的信息量一定是最小的，原因是，在尺度空间中，高斯滤波器的方差σ控制着信息的数量。由此，以比特数来衡量信息的长度，I₀的描述长度肯定大于I_σ，为了能够快速地得到有关I_σ信息量的近似表达式，本发明以一个理想的低通滤波器为模型开始分析：

由傅里叶变换可知，振幅为a的信号在空间域和频域之间满足缩放原理，因此得到如下关系：

$F\left\{s\left(ax\right)\right\}=\frac{1}{a}S\left(\frac{f}{a}\right);a\≠0---\left(11\right)$

那么，对于高斯分布而言，空域和频域之间满足如下关系：

$e^{(-\frac{{\mathrm{\ω}}^2{\mathrm{\σ}}^2}{2})}=e^{(-\frac{x^2}{2{\mathrm{\σ}}^2})};{\mathrm{\σ}}^2\≠0---\left(12\right)$

从式可以看出，高斯分布的空间尺度与频域的尺度成反比，即

由抽样定理可知，若原始信号的最高频率为f，为了能够精确的重建原始信号，采样频率至少为2f。而本发明已经知道，或者说，f与高斯滤波器的带宽成正比，设该比例为常量α时，而高斯滤波器的带宽由其标准差确定。因此，本发明可以将采样率重新表示为：

$s=n\left({\mathrm{\α\σ}}_{\mathrm{\ω}}^2\right);n\≥2---\left(13\right)$

那么结合公式(12)，得到空间尺度与频域尺度的反比关系：

$s=n\left(\frac{\mathrm{\α}}{{\mathrm{\σ}}_x^2}\right);n\≥2---\left(14\right)$

本发明中的算法无需计算出对于每一个s需要多少个位(bit)来表示。已知表示s所需位数(bits)与其信息量成正比，即s∝bits，若给定常量β表示某一个精确值，那么以bit数表示I_σ的信息量，由此本发明可以将上述关系具体表示为：

${\mathrm{dl}}_{I_{\mathrm{\σ}}}=n\left(\frac{\mathrm{\α}\mathrm{\β}}{{\mathrm{\σ}}_x^2}\right);n\≥2---\left(15\right)$

d4计算ε_σ的描述长度

即便本发明在未知噪声的分布情况下，依然可以依据中心极限定理，它的概率密度分布会近似满足零均值的高斯分布，因此，噪声的概率分布函数可以用下式表示：

$P_r\left(\mathrm{\ϵ}\right)=e^{-\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}}---\left(16\right)$

其中表示噪声方差，而ε²表示原始图像I₀与平滑后的图像I_σ之间的二次残差，正如Shannon和Rissanen等所指出，信息的度量与概率密切相关。因此，可以用对数log₂来表示公式(16)，得到残差信息量的描述长度如下所示：

${\mathrm{dl}}_{{\mathrm{\ϵ}}_{\mathrm{\σ}}}=k\left(\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}\right)---\left(17\right)$

其中一旦定义了有限的和那么本发明就可以得到局部信息量的描述长度：

${\mathrm{dl}}_{I_0}=n\left(\frac{\mathrm{\α}\mathrm{\β}}{{\mathrm{\σ}}_x^2}\right)+k\left(\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}\right)---\left(18\right)$

由此，简化后得到：

${\mathrm{dl}}_{I_0}=\left(\frac{\mathrm{\λ}}{{\mathrm{\σ}}_x^2}\right)+{\mathrm{\ϵ}}^2---\left(19\right)$

其中，所有参数均大于0。

e)对图像中的每个点采用最佳局部平滑尺度滤除噪声保持边缘。

e1)在每个位置(x,y)利用σ_i∈(σ_min,σ_max)遍历计算出各个描述长度dl_σ；

e2)依据MDL原则，选择dl_σ的最小值。该最小值min(dl_σ)对应(x,y)处的最优平滑尺度σ^*(x,y)。

e3)利用最优平滑尺度σ^*(x,y)，和位置(x,y)，计算出平滑后的图像I_σ(x,y)。

e4)由于局部方差σ与位置(x,y)存在对应关系，本发明可以构造自适应的高斯滤波器。对于每个位置(x,y)采用对应的特定的σ^*(x,y)进行平滑。本发明中的自适应高斯平滑可以用下面的公式表示，

$I_{\mathrm{\σ}(x,y)}=I_0(x,y)*e^{(-\frac{x^2+y^2}{2{\mathrm{\σ}}^{*}(x,y)})}---\left(20\right)$

其中，I₀为输入图像，I_σ为输出图像，σ^*(x,y)为基于图像中的每个位置，依据MDL原则自动计算出的最优平滑尺度。

e5)为了确定算法中的唯一参数λ，通过大量实际实验确定，该参数在3000左右，能够取得最佳的滤波效果与边缘保护的平衡。

f)对平滑后的图像进行边缘检测，得到图像的全部弱边缘。

本发明的有益效果

本发明针对近来实际生产生活中大量出现的弱边缘的检测进行了细致的研究，提出了一种图像弱边缘的检测方法，除了主观对比外，还应用了一种最新的广泛适用的客观评价标准，峰值信噪比(PSNR)客观评价方法。弱边缘进检测中普遍存在关键的问题：平滑滤波与保护弱边缘之间存在矛盾。

本发明为了解决这种矛盾，提出了多尺度的图像弱边缘检测的观点，在深入研究多尺度技术的基础上，受最小描述长度思想的启发，最终提出了一种基于最小描述长度的多尺度图像弱边缘检测算法，取得良好的滤波平滑效果与弱边缘的自适应增强效果。本发明中，在平滑阶段提出了对弱边缘增强和保护的平滑方式：首先对弱边缘进行局部增强，提高图像中弱边缘的对比度；然后由局部区域的一致性决定平滑尺度，能够在局部区域较为平缓的区域，即噪声比较集中的区域内，利用较大平滑尺度有效抑制噪声；而临近局部区域一致性变化较快的地方，即弱边缘和图像中的边界时，迅速减小尺度，避免对弱边缘产生较强抑制。本发明的算法无需迭代，稳定性高。

附图说明

为了对本发明的算法方案，或传统技术中的算法方案更清晰地予以阐述，下面将对本发明实施例和现有技术结合附图作简单而明了的介绍。

图1为本发明基于最小描述长度的多尺度图像弱边缘检测方法的整体流程图。

图2为本发明基于最小描述长度的多尺度平滑方法流程图。

图3为弱边缘增强效果对比图。

图4为无噪声弱边缘图像leaf。

图5为无噪声情况下，普通高斯平滑与基于最小描述长度的多尺度平滑方法效果对比图。

图6为添加高斯白噪声的弱边缘图像bird。

图7为高斯白噪声情况下，普通高斯平滑与基于最小描述长度的多尺度平滑方法效果对比图。

具体实施方式

下面结合附图对本发明进行进一步阐述。

如图1所示，传统的Canny算子包含下面四个基本步骤：

(1)平滑(Smoothing)：尽可能抑制噪声，同时不破坏真正的边缘。

(2)微分(Gradient)：应用微分算子求平滑图像的梯度，将其作为边缘性的指标。

(3)检测(Detection)：判别哪些边缘像素应该作为噪声而剔除，哪些必须予以保留，通常采用阈值法作为判别标准。

(4)定位(Localization)：确定边缘的准确位置(在有些特殊场合，定位要达到亚像素分辨率，如卫星地图和视觉测量)。

(5)细化与链接(thinning and linking)：保持边缘的单像素特性和完整性。

本发明针对现有算法复杂度过高，而传统Canny算子在弱边缘检测时存在大量漏检和误判现象，在普通边缘检测中增加了对弱边缘的增强与保护环节，以增强弱边缘检测的效果。

本发明边缘检测处理过程中先对弱边缘进行增强，然后进行自适应平滑滤波。本发明首先确立了最小描述长度原则，(minimum description length,MDL)，依据MDL原则确定局部增强尺度与局部滤波尺度在弱边缘增强阶段，本发明利用局部相干扩散(coherencediffusion)法来增强弱边缘对比度。而在局部尺度自适应平滑阶段，利用最小描述长度准则进行局部平滑尺度的估计，最后进行普通边缘检测。

本发明中基于最小描述长度原则的图像弱边缘检测方法执行的流程图如图2所示，它主要包括：

(1)依据尺度空间理论，对原图像进行连续高斯平滑，构造线性尺度空间；

(2)利用有关公式计算图像的局部描述长度；

(3)利用最小描述长度(MDL)准则计算出最优局部尺度；

(4)利用最优局部尺度σ^*作为弱边缘所需的增强尺度σ；

(5)基于局部邻域结构进行自适应弱边缘增强；

(6)基于图像的每个点位置采用最优局部尺度进行滤波去除噪声。

(7)完成弱边缘的增强与滤除噪声。

下面对各步骤进行详细阐述：

1、依据尺度空间理论，对原图像进行连续高斯平滑，构造线性尺度空间；

(1)、G_σ表示标准差为σ的高斯滤波器。对于一幅R×C灰度图像I，利用一系列尺度σ递增的高斯滤波器与图像I卷积，得到平滑后的图像序列，即I_σ＝I₀*G_σ，其中σ＝σ₁,...,σ_n。

(2)在尺度空间中，寻找每个位置(x,y)处的σ_i。本发明利用最小描述长度原则(Minimal Description Length principle)作为σ的选取标准。

本发明中，用下面的公式描述一般的高斯平滑过程，

其中ε_σ为残差，表示原图像与平滑后的图像之间差值。

(3)、将最小描述长度的思想应用到图像的平滑，本发明用长度描述算子(dl)来表示上述平滑过程：

${\mathrm{dl}}_{I_0(x,y)}={\mathrm{dl}}_{I_{\mathrm{\σ}}(x,y)}+{\mathrm{dl}}_{{\mathrm{\ϵ}}_{\mathrm{\σ}}(x,y)}---\left(10\right)$

2、计算图像的局部描述长度；

(1)、I_σ的描述长度

对于任意的σ大于0，当σ→∞时，I_σ的信息量最小。在尺度空间中，高斯滤波器的方差σ控制着提取信息的数量。由此，以比特数来衡量信息的长度，I₀的描述长度远大于I_σ。有关I_σ信息量的表达式，可以通过理想的低通滤波器为基本模型做如下分析：

由傅里叶变换可知，振幅为a的信号在空间域和频域之间满足缩放原理，于是得到：

$F\left\{s\left(ax\right)\right\}=\frac{1}{a}S\left(\frac{f}{a}\right);a\≠0---\left(11\right)$

对于满足高斯分布的信号，空域和频域之间满足如下关系：

$e^{(-\frac{{\mathrm{\ω}}^2{\mathrm{\σ}}^2}{2})}=e^{(-\frac{x^2}{2{\mathrm{\σ}}^2})};{\mathrm{\σ}}^2\≠0---\left(12\right)$

从上式可以看出，高斯分布的空间尺度与频域的尺度成反比，即

由抽样定理可知，若原始信号的最高频率为f，为了能够精确的重建原始信号，采样频率至少为2f。已知3σ_ω为高斯滤波器的最高截止频率，因此f与高斯滤波器的带宽成正比，而高斯滤波器的带宽由其标准差决定，因此可将二者关系记为设α为比例系数，本发明可以将采样信息重新表示为：

$s=n\left({\mathrm{\α\σ}}_{\mathrm{\ω}}^2\right);n\≥2---\left(13\right)$

结合空域和频域之间的关系式，得到空间尺度与频域尺度的反比关系：

$s=n\left(\frac{\mathrm{\α}}{{\mathrm{\σ}}_x^2}\right);n\≥2---\left(14\right)$

在本发明中，无需计算出对于每一个s需要多少个位(bit)来表示，已知表示s所需位数(bits)与其信息量成正比，即s∝bits。若给定常量β表示某一个精确值，那么以bit数表示I_σ的信息量，由此本发明可以将上述关系具体表示为：

${\mathrm{dl}}_{I_{\mathrm{\σ}}}=n\left(\frac{\mathrm{\α}\mathrm{\β}}{{\mathrm{\σ}}_x^2}\right);n\≥2---\left(15\right)$

(2)、ε_σ的描述长度

即便本发明在未知噪声的分布情况下，依然可以依据中心极限定理，它的概率密度分布会近似满足零均值的高斯分布，因此，噪声的概率分布函数可以用下式表示：

$P_r\left(\mathrm{\ϵ}\right)=e^{-\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}}---\left(16\right)$

其中表示噪声方差，而ε²表示原始图像I₀与平滑后的图像I_σ之间的二次残差，用对数log₂来表示噪声的概率分布函数，得到残差的信息量的描述长度如下所示：

${\mathrm{dl}}_{{\mathrm{\ϵ}}_{\mathrm{\σ}}}=k\left(\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}\right)---\left(17\right)$

其中一旦定义了有限的和就可以由下式得到局部信息量的描述长度：

${\mathrm{dl}}_{I_0}=n\left(\frac{\mathrm{\α}\mathrm{\β}}{{\mathrm{\σ}}_x^2}\right)+k\left(\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}\right)---\left(18\right)$

由此，简化后得到：

${\mathrm{dl}}_{I_0}=\left(\frac{\mathrm{\λ}}{{\mathrm{\σ}}_x^2}\right)+{\mathrm{\ϵ}}^2---\left(19\right)$

其中，所有参数均大于0。

3、利用最小描述长度(MDL)准则计算出最优局部尺度；

在每个位置(x,y)利用σ_i∈(σ_min,σ_max)遍历计算出各个描述长度dl_σ；然后，依据MDL原则，选择dl_σ的最小值。此时该min(dl_σ)对应着(x,y)处的最优平滑尺度σ^*。该最小描述长度，意味着最大的平滑量，同时最小的残差。这一步中，本发明得到了最优局部方差(σ^*(x,y))

4、利用最优局部尺度σ^*作为弱边缘增强的局部尺度σ；

5、基于局部结构进行自适应弱边缘增强；

先将图像转化为灰度图像，然后采用相干扩散算法对图像的弱边缘进行增强，具体计算方法为：

a1)计算结构张量S，可由下式估计它的局部方向：

$S=\left(\begin{array}{cc}{s}_{11}& s_{12}\\ s_{12}& s_{22}\end{array}\right)=\left(\begin{array}{cc}\frac{{\∂}^{2}L}{\∂x^2}*G_{\mathrm{\σ}}& \frac{{\∂}^{2}L}{\∂x\∂y}*G_{\mathrm{\σ}}\\ \frac{{\∂}^{2}L}{\∂x\∂y}*G_{\mathrm{\σ}}& \frac{{\∂}^{2}L}{\∂y^2}*G_{\mathrm{\σ}}\end{array}\right)---\left(1\right)$

其中G_σ表示局部尺度为σ的高斯滤波器，结构张量的特征向量为局部梯度的方向，而两个特征值的差，代表着局部邻域的非均匀性。

a2)构造扩散张亮D，如下所示：

$D=R^T\left(\begin{array}{cc}{c}_1& 0\\ 0& c_2\end{array}\right)R---\left(2\right)$

其中R为旋转矩阵，它的列向量由结构张量的特征向量表示，其中c₁和c₂为沿梯度方向的引导相干系数。

扩散张量D的元素可采用下式计算：

$d_{11}=\frac{1}{2}(c_1+c_2+\frac{(c_2-c_1)(s_{11}-s_{22})}{\mathrm{\α}})---\left(3\right)$

$d_{12}=\frac{(c_2-c_1)s_{12}}{s_{12}}---\left(4\right)$

$d_{22}=\frac{1}{2}(c_1+c_2-\frac{(c_2-c_1)(s_{11}-s_{22})}{\mathrm{\α}})---\left(5\right)$

其中，

$\mathrm{\α}=\sqrt{{(s_{11}-s_{22})}^2+4s_{12}^2}---\left(6\right)$

结构张量的特征值为：

${\mathrm{\λ}}_{1,2}=\frac{1}{2}(s_{11}+s_{22}\±\mathrm{\α})---\left(7\right)$

在计算中，本发明利用两个特征向量的差(即邻域的非均匀程度)来控制扩散速度c₁和c₂：

$c_1=max(0.01,1-e^{-{({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)}^2/k^2})---\left(8\right)$

c₂＝0.01

6、基于图像的每个点位置采用最优局部尺度进行滤波去除噪声。

对于每个位置(x,y)采用对应的特定的σ^*(x,y)进行模糊。本发明中的自适应高斯平滑可用如下公式表示：

$I_{\mathrm{\σ}(x,y)}=I_0(x,y)*e^{(-\frac{x^2+y^2}{2{\mathrm{\σ}}^{*}(x,y)})}---\left(20\right)$

其中，I₀为输入图像，I_σ为输出图像，σ^*为基于图像中的每个位置，依据MDL原则自动计算出的最优平滑尺度。

算法中存在的唯一参数λ，默认值为3000，为了能够取得更佳的滤波效果与边缘保护的平衡，只需通过观察实际图像略微调整该参数。

7、完成弱边缘的增强与滤除噪声。

由图3的结果可以明显地看出，原图中已经退化的边缘的局部对比度得到提升，弱边缘达到了很好的增强效果。经过尺度的自适应变化，所有的弱边缘都得到增强，提高了抗平滑抑制能力。

图4为无噪声的弱边缘图像，分别采用传统的平滑方法和本发明方法，对其进行平滑的效果如图5所示，左图显示普通高斯平滑的效果，右图显示了本发明的平滑效果。可以看出，普通高斯平滑使得部分弱边缘对比度提高，而部分弱边缘对比度降低；而本发明的多尺度MDL平滑方法使原图像的整体对比度都得到提高，边缘数量更多，细节更丰富。

图6为含有高斯白噪声的弱边缘图像，分别采用传统的平滑方法和本发明方法，对其进行平滑的效果如图7所示，左图显示普通高斯平滑的效果，右图显示了本发明的平滑效果。可以看出，传统方法在滤除噪声的同时，只保留了原来的强对比度的边缘，对弱边缘进行了抑制，使得部分弱边缘消失，产生平滑过度现象。而本发明中的多尺度MDL平滑方法，在有效滤除噪声的同时，提升了弱边缘的对比度，因此图像的整体对比度得到提高。图像的细节更丰富，弱边缘更清晰。

表1为采用PSNR(峰值信噪比)图像质量客观评价标准的对比数据，其中sigma为高斯白噪声的标准差，PSNR的单位为dB，表示相差一个数量级。PSNR值分别反映了平滑图像的质量：PSNR值越高，图像信噪比越高，图像质量越好，对应的图像平滑效果越好。

表1 采用PSNR客观评价两种平滑方法的图像质量

从表1中看出，在无噪声和有噪声情况下，本发明中基于最小描述长度的多尺度平滑算法的PSNR量化指标，始终优于普通高斯平滑方法。说明本发明比普通平滑方法的抗噪声性能更强，得到的图像质量更佳。

综合以上对比，本发明达到了滤除噪声同时保护了弱边缘的效果，显著提高了弱边缘图像的质量，增强了弱边缘的提取效果。

Claims (4)

1.一种基于最小描述长度的多尺度图像弱边缘检测方法，其特征在于，该方法包含以下步骤：

a)对图像的预处理，包括对图像进行灰度化与边缘增强，提高图像中的弱边缘的对比度；

b)构造线性尺度空间；

c)计算图像的局部描述长度；

d)采用最小描述长度MDL原则确定局部区域的最佳平滑尺度；

e)对图像中的每个点采用局部区域的最佳平滑尺度滤除噪声保持边缘；

f)对平滑后的图像进行边缘检测，得到图像的全部边缘；

步骤d)包括：

d1)利用最小描述长度原则作为方差σ的选取标准；

对于一般高斯平滑过程，可用下面的公式描述：

其中ε_σ为残差，表示原图像与平滑后的图像之间差值；

d2)先将最小描述长度应用于图像平滑，然后定义长度描述算子(dl)，则上述平滑过程可重新表示为：

${\mathrm{dl}}_{I_0(x,y)}={\mathrm{dl}}_{I_{\mathrm{\σ}}(x,y)}+{\mathrm{dl}}_{{\mathrm{\ϵ}}_{\mathrm{\σ}}(x,y)}$

d3)计算I_σ的描述长度

对于任意的σ大于0，当σ→∞时，I_σ的信息量最小，以比特数来衡量信息的长度，I₀的描述长度远大于I_σ，以低通滤波器表示信息采样模型：

由傅里叶变换可知，振幅为a的信号在空间域和频域之间满足缩放原理，从而得到：

$F\left\{s\left(ax\right)\right\}=\frac{1}{a}S\left(\frac{f}{a}\right);a\≠0$

对于满足高斯分布的图像信号，空域和频域之间满足如下关系：

$e^{(-\frac{{\mathrm{\ω}}^2{\mathrm{\σ}}^2}{2})}=e^{(-\frac{x^2}{2{\mathrm{\σ}}^2})},{\mathrm{\σ}}^2\≠0$

由上式可见，高斯分布的空间尺度与频域的尺度成反比，即由抽样定理可知，若原始信号的最高频率为f，采样频率至少为2f，而高斯滤波器的σ_ω决定最高采样频率，因此，采样信息可重新表示为：

$s=n\left({\mathrm{\α\σ}}_{\mathrm{\ω}}^2\right);n\≥2$

其中α为比例常量，结合频域与空域之间的反比关系，得到的空间尺度与频域尺度的反比关系：

$s=n\left(\frac{\mathrm{\α}}{{\mathrm{\σ}}_x^2}\right),n\≥2$

由上式可知s∝bits，若给定常量β表示某一个精确值，以bit数表示I_σ的信息量，将上述关系具体表示为：

${\mathrm{dl}}_{I_{\mathrm{\σ}}}=n\left(\frac{\mathrm{\α}\mathrm{\β}}{{\mathrm{\σ}}_x^2}\right);n\≥2$

d4)ε_σ的描述长度

依据中心极限定理，它的概率密度分布会近似满足零均值的高斯分布，噪声的概率分布函数可以用下式表示：

$P_r\left(\mathrm{\ϵ}\right)=e^{-\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}}$

其中表示噪声方差，而ε²表示原始图像I₀与平滑后的图像I_σ之间的二次残差，用对数log₂来表示上式，得到残差的信息量的描述长度如下所示

${\mathrm{dl}}_{{\mathrm{\ϵ}}_{\mathrm{\σ}}}=k\left(\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}\right)$

其中定义信息长度和后，局部信息的描述长度可以表示为：

${\mathrm{dl}}_{I_0}=n\left(\frac{\mathrm{\α}\mathrm{\β}}{{\mathrm{\σ}}_x^2}\right)+k\left(\frac{{\mathrm{\ϵ}}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}\right)$

由此，简化后得到：

${\mathrm{dl}}_{I_0}=\left(\frac{\mathrm{\λ}}{{\mathrm{\σ}}_x^2}\right)+{\mathrm{\ϵ}}^2$

其中，所有参数均大于0。

2.根据权利要求1所述的基于最小描述长度的多尺度图像弱边缘检测方法，其特征在于，步骤a)包括：

先将图像转化为灰度图像，然后采用相干增强算法对图像的弱边缘进行增强，其计算方法为：

a1)计算结构张量S，可由下式计算局部方向：

$S=\left(\begin{array}{cc}{s}_{11}& s_{12}\\ s_{12}& s_{22}\end{array}\right)=\left(\begin{array}{cc}\frac{{\∂}^{2}L}{\∂x^2}*G_{\mathrm{\σ}}& \frac{{\∂}^{2}L}{\∂x\∂y}*G_{\mathrm{\σ}}\\ \frac{{\∂}^{2}L}{\∂x\∂y}*G_{\mathrm{\σ}}& \frac{{\∂}^{2}L}{\∂y^2}*G_{\mathrm{\σ}}\end{array}\right)$

其中G_σ表示局部尺度为σ的高斯滤波器，结构张量的特征向量为局部梯度的方向，而两个特征值的差，代表局部邻域的非均匀性；

a2)构造扩散张量D，如下所示，

$D=R^T\left(\begin{array}{cc}{c}_1& 0\\ 0& c_2\end{array}\right)R$

其中R为旋转矩阵，它的列向量由结构张量的特征向量表示，其中c₁和c₂为沿梯度方向的引导相干系数，即扩散速度；

扩散张量D的元素采用以下方式计算：

$d_{11}=\frac{1}{2}(c_1+c_2+\frac{(c_2-c_1)(s_{11}-s_{22})}{\mathrm{\α}})$

$d_{12}=\frac{(c_2-c_1)s_{12}}{s_{12}}$

$d_{22}=\frac{1}{2}(c_1+c_2-\frac{(c_2-c_1)(s_{11}-s_{22})}{\mathrm{\α}})$

其中，

$\mathrm{\α}=\sqrt{{(s_{11}-s_{22})}^2+4s_{12}^2}$

结构张量的特征值为：

${\mathrm{\λ}}_{1,2}=\frac{1}{2}(s_{11}+s_{22}\±\mathrm{\α})$

邻域的非均匀程度由两个特征值之差λ₁-λ₂表示，可用来控制扩散速度c₁和c₂：

$c_1=max(0.01,1-e^{-{({\mathrm{\λ}}_1-{\mathrm{\λ}}_2)}^2/k^2})$

c₂＝0.01。

3.根据权利要求1所述的基于最小描述长度的多尺度图像弱边缘检测方法，其特征在于，步骤b)包括：

b1)构造一个线性尺度空间；G_σ表示标准差为σ的高斯滤波器，对于一幅R×C灰度图像I，利用一系列线性递增的尺度σ作为高斯滤波器的标准差，与图像I卷积，得到平滑后的图像序列：即I_σ＝I₀*G_σ，其中σ＝σ₁,...,σ_n，I₀为图像I；

b2)在尺度空间中，搜索每个位置(x,y)处的σ_i。

4.根据权利要求1所述的基于最小描述长度的多尺度图像弱边缘检测方法，其特征在于，步骤e)包括：

e1)在每个位置(x,y)利用σ_i∈(σ_min,σ_max)遍历计算出各点处描述长度dl_σ；

e2)依据MDL原则，选择dl_σ的最小值，该min(dl_σ)对(x,y)的最优平滑尺度为σ^*(x,y)；

e3)利用最优平滑尺度σ^*(x,y)，计算局部平滑图像I_σ(x,y)；

e4)利用局部方差σ与位置(x,y)存在对应关系，构造自适应的高斯滤波器，对每个位置(x,y)采用对应的最优平滑尺度σ^*(x,y)进行局部平滑；对于整幅图像，自适应高斯平滑用下面的公式表示：

$I_{\mathrm{\σ}(x,y)}=I_0(x,y)*e^{(-\frac{x^2+y^2}{2{\mathrm{\σ}}^{*}(x,y)})}$

其中，I₀为输入图像，I_σ为输出图像，σ^*(x,y)为基于图像中的每个位置，依据MDL原则自动计算出的最优平滑尺度；

e5)实验观察确定唯一参数λ，默认值为3000。