CN101877122B - 扩散程度可控的各向异性扩散图像去噪增强方法 - Google Patents

Abstract

本发明提供一种扩散程度可控的基于迹模型的各向异性扩散图像去噪增强方法，包括以下步骤：计算每个像元点的结构张量和Hessian矩阵；对每个像元点的结构张量进行特征值分解；构建每个像元点的扩散张量，使该扩散张量的特征向量是结构张量的特征向量，特征值是结构向量特征值的函数，即扩散率函数，该函数可通过5个参数的调节使扩散程度可控，其中一个参数控制图像增强的效果；对迹模型迭代求解。采用该方法，可以有效对图像去噪增强。

Description

扩散程度可控的各向异性扩散图像去噪增强方法

技术领域

本发明涉及图像处理技术，尤其是图像去噪技术，可用于处理有加性或乘性噪声的图像。

背景技术

图像去噪，尤其对于抑制斑点噪声的滤波算法，通常使用的有Lee滤波、Kuan滤波、Frost滤波等降噪方法，这类方法的共同之处是根据图像局部统计特征，选取合适窗口大小及调整滤波函数进行滤波，只是不同方法窗口选取准则和采取的滤波器不同。但共同的弱点是：这类方法在边缘方向上都是各向同性滤波，模糊了图像结构和细节信息。后来发展了一系列基于各向异性扩散的滤波方法，其中方法之一是Tschumperle的文章“Vector-Valued ImageRegularization with PDE’s：A Common Framework for Different Applications”(发表在国际会议IEEE International Conference on Computer Vision and Pattern Recognition(CVPR’03)，Vol.1，pp.651--659，Madison/USA，June 2003)所揭示的基于迹算子模型的各向异性滤波方法。该方法提出了基于局部几何结构的各向异性扩散模型-迹算子模型滤波方法，其中基于迹算子模型的偏微分方程是：

$\frac{\∂I}{\∂t}=\mathrm{trace}\left(\mathrm{TH}\right)---\left(1\right)$

其中，trace(·)为矩阵的求迹运算，H为Hessian矩阵，T是扩散张量。扩散张量是驱动滤波过程的张量，它控制着局部的扩散过程，也是决定不同的各向异性扩散方法的本质变量。

已知的利用迭代方式求解迹模型的各向异性扩散方法包括以下步骤：

(a)读入二维(n＝2)初始灰度图像，记为I(0)。符号I(m)表示图像序列，其中m表示迭代次数。

(b)计算图像I(m)上每个像元的一阶差分I_x，I_y；并计算对应每个像元点的2×2的对称矩阵，即结构张量S。

(c)计算I(m)每一像元点的二阶差分I_xx，I_xy，I_yy，得到每一像元点的2×2的Hessian矩阵H。

(d)对每一像元点的结构张量进行矩阵的特征值分解，使得 $S=\left(\begin{array}{cc}{\mathrm{\θ}}_1& {\mathrm{\θ}}_2\end{array}\right)\left(\begin{array}{cc}{u}_1& 0\\ 0& u_2\end{array}\right)\left(\begin{array}{c}{\mathrm{\θ}}_1\\ {\mathrm{\θ}}_2\end{array}\right),$ 其中θ₁、θ₂分别是对应于特征值u₁、u₂的特征向量。

(e)对每一像元点计算2×2的扩散张量D，D具有形式 $D=g_1(u_1,u_2){\mathrm{\θ}}_1{\mathrm{\θ}}_1^T+g_2(u_1,u_2){\mathrm{\θ}}_2{\mathrm{\θ}}_2^T,$ 其中D的特征向量等于S的特征向量，D的特征值分别为扩散率函数g₁(u₁，u₂)、g₂(u₁，u₂)，上标T表示向量或矩阵的转置。

(e)使用公式(1)的离散方程I(m+1)＝I(m)+c×trace(DH)，从I(m)计算I(m+1)中的每一点。其中，c是预定义常数，0≤c≤1。

(f)重复步骤(b)至(e)共M次，获得结果图像。

上述步骤中Tschumperle基于平均曲率定义的扩散率函数g₁(u₁，u₂)、g₂(u₁，u₂)如下式：

$\left\{\begin{array}{c}{g}_1(u_1,u_2)=\frac{1}{1+u_1+u_2}\\ g_2(u_1,u_2)=\frac{1}{{(1+u_1+u_2)}^{1/2}}\end{array}\right.---\left(2\right)$

这样的扩散率函数使得扩散滤波随着平均曲率的增加而减弱，即在局部结构强的地方，滤波减弱，而在局部结构弱的地方，滤波增强，且在梯度方向上扩散滤波减弱速度快于边缘方向。这是这类扩散率函数带来的扩散滤波的优势之处。尽管如此，由该扩散率函数生成的滤波器仍然具有不足之处：它缺少可控性和松弛性，即对于任何噪声图像，不管噪声水平如何，扩散率函数仅与平均曲率相关；而且不具备同时增强边缘等特征的能力。

发明内容

本发明针对现有技术中存在的缺陷和不足，在迹算子模型的偏微分方程(1)各向异性滤波的框架下提出一种新的用于扩散滤波保持边界的扩散率函数，该函数具有很好的可控性和松弛性，并具备能力在去噪同时增强图像。

本发明的技术方案如下：

一种基于图像局部几何结构的新的扩散率函数的各向异性扩散的图像降噪方法，其特征在于包含如下步骤：

(a)读入二维(n＝2)初始灰度图像，记为I(0)。符号I(m)表示图像序列，其中m表示迭代次数。

(b)计算图像I(m)上每个像元的一阶差分I_x，I_y；并计算对应每个像元点的2×2的对称矩阵，即结构张量S。

(c)计算I(m)每一像元点的二阶差分I_xx，I_xy，I_yy，得到每一像元点的2×2的Hessian矩阵H。

(d)对每一像元点的结构张量进行矩阵的特征值分解，使得 $S=\left(\begin{array}{cc}{\mathrm{\θ}}_1& {\mathrm{\θ}}_2\end{array}\right)\left(\begin{array}{cc}{u}_1& 0\\ 0& u_2\end{array}\right)\left(\begin{array}{c}{\mathrm{\θ}}_1\\ {\mathrm{\θ}}_2\end{array}\right),$

其中θ₁、θ₂分别是对应于特征值u₁、u₂的特征向量。θ₁、θ₂为图像灰度值变化的最大、最小方向；u₁、u₂则反映了特征方向上的灰度变化强弱程度。

(e)在每一像元点计算新的扩散率函数g₁′(u₁，u₂)、g₂′(u₁，u₂)，其中u₁，u₂是步骤(d)中结构张量S矩阵分解得到的特征值。

(f)对每一像元点计算新的2×2的扩散张量D，D具有形式 $D={g_1}^{\′}(u_1,u_2){\mathrm{\θ}}_1{\mathrm{\θ}}_1^T+{g_2}^{\′}(u_1,u_2){\mathrm{\θ}}_2{\mathrm{\θ}}_2^T,$ 其中D的特征向量等于S的特征向量，D的特征值分别为扩散率函数g₁′(u₁，u₂)、g₂′(u₁，u₂)。

(g)使用公式(1)的离散方程I(m+1)＝I(m)+c×trace(DH)，从I(m)计算I(m+1)中的每一点。其中，c是预定义常数，0≤c≤1。

(h)重复步骤(b)至(g)共M次，获得结果图像。

所述图像上每个像元点对应的结构张量是2阶对称矩阵S，S是梯度算子的外积与某个

低通核函数的卷积，即将卷积平滑后的矩阵作为结构张量矩阵。梯度算子的外积具有形式 $\left(\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right),$ 其中，I_x、I_y是一阶微分或差分算子。

所述图像上每个像元点对应的Hessian矩阵是2阶对称矩阵H，其中 $H=\left(\begin{array}{cc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}\\ I_{\mathrm{yx}}& I_{\mathrm{yy}}\end{array}\right).$

所述的矩阵分解是指图像上每个像元点对应的结构矩阵的特征值分解。

所述的扩散张量是一个2阶矩阵，它直接利用矩阵特征值分解的形式构建，构建的特征值与特征向量具有以下特点：

(1)扩散张量T的特征向量与局部结构张量的特征向量相同，即扩散张量的特征向量选取局部结构张量的特征向量θ₁，θ₂；

(2)扩散张量D的特征值分别为扩散率函数g₁′(u₁，u₂)、g₂′(u₁，u₂)，分别表征在梯度方向和轮廓方向上的扩散滤波程度。

所述的新的扩散率函数具有形式：

$g_1^{\′}(u_1,u_2)=\mathrm{\α}(1-\frac{{(u_1+u_2)}^{b_1}}{{(u_1+u_2)}^{b_1+a}+{\mathrm{\β}}^{b_1}}),$ $g_2^{\′}(u_1,u_2)=\mathrm{\α}(1-\frac{{(u_1+u_2)}^{b_2}}{{(u_1+u_2)}^{b_2}+{\mathrm{\β}}^{b_2}})---\left(3\right)$

其中α，β，a，b₁，b₂是五个控制参数，均是大于零的实数，可依据图像噪声特点选择。

所述的预定义常数c在区间[0，1]，预先给出。

所述的迭代次数M预先给出。

本发明与现有技术相比有以下特点：本发明的扩散率函数有多个自由参数，这些参数从不同方面控制扩散滤波过程，使得可依据图像局部几何结构和噪声水平

附图说明

图1为迭代实现的基于偏微分方程图像扩散滤波的基本流程图

图2为基于迹算子模型的各向异性扩散滤波的基本流程图

图3分(a)、(b)、(c)图分别是控制参数a、b、β分别变化时扩散率函数变化的示意图

图4本发明去噪结果示例：(a)原噪声图像(含加性白噪声)；(b)Lee滤波去噪结果；(c)本方法迭代20次去噪结果；(d)本方法迭代40次去噪结果

具体实施方式

下面结合附图通过实施例对本发明作进一步的详细说明。

图1是迭代实现的基于偏微分方程图像扩散滤波的基本流程图。迭代次数是M，是预先指定的。当扩散滤波次数小于指定次数M时，将滤波结果图像作为当前图像再次滤波，否则滤波结果图像作为最终滤波结果。其中滤波单元100是本发明基于偏微分方程扩散滤波的主要部分。

图2是本发明实现的基于迹算子模型的各向异性扩散滤波的流程图。图2中读入的图像经过下列6个单元处理后获得结果。

单元110对读入的图像计算对应每个像元的结构张量S。图像上每个像元点对应的结构张量是2阶对称矩阵，它是梯度算子的外积与某个低通核函数的卷积，梯度算子的外积具有形式 $\left(\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right),$ 其中，I_x、I_y是一阶微分或差分算子。卷积运算在降低结构张量噪声影响的同时揭示图像的局部几何结构。合适的低通核函数的例子是高斯核函数。2维的高斯核函数的形式是

$K_{\mathrm{\σ}}=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-\frac{x^2+y^2}{{\mathrm{\σ}}^2})$

因此，实际应用时，局部结构张量通常由下式表达：

$S=K_{\mathrm{\σ}}*(\▿I\⊗\▿I)=K_{\mathrm{\σ}}*(\▿I\▿I^T)(\mathrm{\σ}\≥0)$

这里，符号*代表卷积运算，

代表外积运算，

表示梯度或一阶差分算子，上标T表示矩阵或向量的转置。

单元120计算I(m))每一像元点的二阶差分 $I_{\mathrm{xx}}=\frac{{\∂}^{2}I}{{\∂}^{2}x},$ $I_{\mathrm{xy}}=\frac{{\∂}^{2}I}{\∂x\∂y},$ $I_{\mathrm{yy}}=\frac{{\∂}^{2}I}{{\∂}^{2}y},$ 得到每一像元点的Hessian矩阵 $H=\left(\begin{array}{cc}{I}_{\mathrm{xx}}& I_{\mathrm{xy}}\\ I_{\mathrm{yx}}& I_{\mathrm{yy}}\end{array}\right).$

单元130对每一像元点的结构张量S进行特征值分解，使得 $S=\left(\begin{array}{cc}{\mathrm{\θ}}_1& {\mathrm{\θ}}_2\end{array}\right)\left(\begin{array}{cc}{u}_1& 0\\ 0& u_2\end{array}\right)\left(\begin{array}{c}{\mathrm{\θ}}_1\\ {\mathrm{\θ}}_2\end{array}\right),$ 其中θ₁、θ₂分别是对应于特征值u₁、u₂的特征向量。θ₁、θ₂为图像灰度值变化的最大、最小方向；u₁、u₂则反映了特征方向上的灰度变化强弱程度。

单元140按照公式(3)计算新的扩散率函数g₁′(u₁，u₂)、g₂′(u₁，u₂)。新的扩散率函数中有五个调控参数α，β，a，b₁，b₂。这五个参数分别从不同方面控制扩散过程：α控制着平坦区域、边缘方向扩散滤波程度和梯度方向上最大的扩散滤波程度。在平坦区u₁≈u₂＝0；g₁≈g₂＝α这时梯度方向和边缘方向上的扩散滤波程度相同，相当于对平坦区域进行高斯滤波；对于完全结构区域u₁＜u₂，g₁≈0，g₂≈α，这时在梯度方向不滤波，只进行边缘方向的滤波，相当于平均曲率运动扩散；对于不完全结构区域，a，b，β对扩散过程的控制可以由附图2的函数y(x)＝1-x^b/(x^b+a+β^b)的图示直观看出。附图3(a)是参数a，β固定b变化时的曲线的变化情况；附图2(b)是参数b，β固定a变化时的曲线的变化情况；附图2(c)是参数a，b固定β变化时的曲线的变化情况。由图2可见，参数b控制着扩散率函数从哪开始弯曲下降，作为一个软域值控制噪声信号区域的转换，控制转换的速度，b越大函数下降越快，随平均曲率增加，越快转换，滤波程度越快减弱；参数β控制从何处扩散率函数接近0；参数a控制着该扩散率函数最小值，即扩散率是否会过冲，给出小于0的值，扩散率小于零可增加边界对比度，起到图像增强的效果。实验中，令b₁＞b₂使扩散率函数在梯度方向下降快于边缘方向，并选择参数控制不完全结构区域的滤波。

单元150对每一像元点重构新的扩散张量D，其中D的特征向量等于S的特征向量，D的特征值为公式(3)定义的新的扩散率函数g₁′(u₁，u₂ )、g₂′(u₁，u₂)。

单元160使用公式(3)的离散方程I(m+1)＝I(m)+c*trace(DH)，从I(m)计算I(m+1)中的每一点。其中，c是预定义常数，0≤c≤1。

附图4是本发明应用的一个例子。图3(a)是添加了噪声的棋盘图；图3(b)是经典的Lee滤波的结果，显然去噪的同时模糊了边缘；图3(c)是本发明迭代20次处理的结果；图3(d)是本发明迭代40次处理的结果。图3(c)和图3(d)扩散滤波时选择的相关参数为：a＝0.1，b₁＝3，b₂＝1，α＝0.8，β＝0.1，迭代步长c＝0.05。图3结果表明对比Lee滤波结果，本发明方法在有效去噪的同时，边缘和细节信息保持得更充分。

Claims (2)

1.一种扩散程度可控的基于迹模型的各向异性扩散图像去噪增强方法，其特征在于包含以下步骤：计算每个像元点的结构张量和Hessian矩阵；对每个像元点的结构张量进行特征值分解；构建每个像元点的扩散张量，使该扩散张量的特征向量是结构张量的特征向量，特征值是结构张量特征值的函数，即扩散率函数，该扩散率函数可通过参数调节使扩散程度可控，该函数具有形式

$g_1^{\′}(u_1,u_2)=\mathrm{\α}(1-\frac{{(u_1+u_2)}^{b_1}}{{(u_1+u_2)}^{b_1+a}+{\mathrm{\β}}^{b_1}}),g_2^{\′}(u_1,u_2)=\mathrm{\α}(1-\frac{{(u_1+u_2)}^{b_2}}{{(u_1+u_2)}^{b_2}+{\mathrm{\β}}^{b_2}})$

其中α，β，a，b₁，b₂是五个控制参数，均是大于零的实数，通过改变参数的大小改变函数的形状，从而控制扩散程度；参数a控制着该扩散率函数最小值，即扩散率是否会过冲从而给出小于0的值，扩散率小于零可增加边界对比度，起到图像增强的效果；最后对迹模型迭代求解，获得去噪增强的结果；其中所述结构张量是该点梯度算子外积和一个合适的低通滤波核函数的卷积。

2.根据权利要求1所述的一种扩散程度可控的基于迹模型的各向异性扩散图像去噪增强方法，其特征在于：所述合适的低通滤波核函数是高斯核函数。

Images