CN102393955B - 用于图像复原的全信息非局部约束全变分方法 - Google Patents

Abstract

本发明公开了一种用于图像复原的全信息非局部约束全变分方法，主要解决现有技术在图像复原时，不能锐化边缘和恢复高频细节的问题。本发明的技术方案为：(1)用非局部均值滤波法对模糊图像抑噪；(2)用维纳滤波法初始化复原结果；(3)计算全信息非局部权重系数矩阵；(4)用阈值迭代公式更新复原结果；(5)用全变分去噪法对复原结果抑噪；(6)判断是否要更新全信息非局部权重系数矩阵，如果是，返回步骤(3)，否则，执行步骤(7)；(7)判断是否满足停止条件，如果是，得到最终结果；否则，返回步骤(4)，直到满足停止条件为止。本发明在复原时，能够锐化图像边缘，恢复高频细节，可用于对已知模糊类型的模糊图像进行复原。

Description

用于图像复原的全信息非局部约束全变分方法

技术领域

本发明属于图像处理技术领域，具体地说是一种对模糊图像进行复原的方法，该方法可用于对各种已知模糊类型的模糊图像进行复原。

背景技术

图像复原是指去除或减轻在获取数字图像过程中发生的图像质量下降的现象，它是图像处理中重要而又富有挑战性的研究内容。对于图像复原问题，研究者已经提出了很多方法。

传统的复原方法有逆滤波，维纳滤波，卡尔曼滤波和广义逆的奇异值分解法等，这些方法已经被广泛地应用于图像复原上，但是这些方法要求模糊图像具有较高的信噪比，如逆滤波的方法仅适用于高信噪比的图像，这一点限制了传统的复原方法在实际中的应用。这些方法的另一个缺点就是在复原时，图像边缘不能很好地恢复，同时又丢失了一些细节信息。

上述经典的复原方法不但效果差，而且在实际应用中不能很好的实现。因此，目前国际上提出了一些改进上述缺点的图像复原方法。如，I.Daubechies等人提出基于小波的阈值迭代法，参见文献《An iterative thresholding algorithm for linear inverseproblems with a sparsity constraint》，Commun.Pure Appl.Math.，2004，Vol.57，No.11，pp.1413-1457。这种方法将两次迭代所得的复原结果的差值作为下一次迭代结果的补偿，是一种有效的复原方法。但是，这种方法是在小波域进行噪声抑制，容易产生振铃效应，且不能锐化图像边缘。此后，J.Bioucas-Dias等人将阈值迭代法进行了改进，参见文献《A new TwIST：two-step iterative shrinkage/thresholding algorithms for imagerestoration》，IEEE Trans.Image Process.，2007，Vol.16，No.12，pp.2992-3004。该方法的收敛速度比一般的阈值迭代法有所提高，同时，J.Bioucas-Dias等人在他们的代码示例中，将噪声系数转换到全变分域中进行抑制，去除了振铃效应，但是这种方法在图像的平滑区域容易产生阶梯效应，且不能很好的恢复图像高频细节。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种用于图像复原的全信息非局部约束全变分方法，以在图像复原时，锐化图像边缘，避免产生阶梯效应，恢复图像高频细节信息，提高模糊图像的恢复质量。

实现本发明目的的技术方案是：将非局部均值滤波的权重系数的产生方法进行改进，并将其作为约束加入全变分图像复原模型中，以阈值迭代法作为求解方法，来进行图像复原。其具体步骤包括如下：

(1)用非局部均值滤波法对输入的模糊图像y进行噪声抑制，得到抑噪后的模糊图像x^(-1)；

(2)设定迭代误差ε＝1×10^-6，设定当前的迭代次数k＝0，用维纳滤波法对抑噪后的模糊图像x^(-1)进行滤波，得到初始复原结果图x⁽⁰⁾，用初始复原结果图x⁽⁰⁾初始化迭代复原结果图x^(k)；

(3)计算迭代复原结果图x^(k)的非局部权重系数矩阵

其中，W₁(i，j)为非局部权重系数矩阵W₁的第i行，第j列的元素值，i＝1，2，...，N，j＝1，2，...，N，N为迭代复原结果图x^(k)的像素总数，令W₁(i，j)的计算公式为：

其中，

表示迭代复原结果图x^(k)的第i个f×f像素的图像小块x_i和第j个f×f像素的图像小块x_j之间的权重系数，f为图像小块的边长，f＝5，

c_i为归一化因子，

h为调节权重大小的参数，h＝65，a_i是一个列向量，它包含x_i所对应的最大的10个权重系数

(4)计算迭代复原结果图x^(k)的贝叶斯非局部权重系数矩阵

其中，W₂(i，j)为贝叶斯非局部权重系数矩阵W₂的第i行，第j列的元素值，令W₂(i，j)的计算公式为：

其中，

表示迭代复原结果图x^(k)的第i个f×f像素的图像小块x_i和第j个f×f像素的图像小块x_j之间的权重系数， $b_i^j=\frac{1}{z_i}\exp (-\frac{1}{2}{(\frac{\left|\right|x_i-x_j\left|\right|}{\mathrm{\σ}}-\sqrt{2f-1})}^2),$ z_i为归一化因子， $Z_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp (-\frac{1}{2}{(\frac{\left|\right|x_i-x_j\left|\right|}{\mathrm{\σ}}-\sqrt{2f-1})}^2),$ σ为已知噪声标准差，b_i是一个列向量，它包含x_i所对应的最大的10个权重系数

(5)计算迭代复原结果图x^(k)的全信息非局部权重系数矩阵W＝rW₁+(1-r)W₂，其中，r为调节参数，r＝0.4；

(6)按照如下的全信息非局部约束阈值迭代公式，计算迭代复原中间结果图x^(k+1/2)：

x^(k+1/2)＝x^(k)+(H^Ty-Ux^(k)-Vx^(k))

其中，U为变形的模糊核矩阵，U＝H^TH，H为已知的模糊核函数矩阵，H^T为H的转置矩阵，V为变形的非局部系数矩阵，V＝γ²W^TW，W^T为权重系数矩阵W的转置矩阵，γ为调节参数，γ＝0.1；

(7)用基于全变分模型的去噪方法对迭代复原中间结果图x^(k+1/2)进行噪声抑制，得到抑噪后的迭代复原结果图x^(k+1)；

(8)判断迭代次数k是否等于更新代数q，q＝100，200，...，1000，如果k＝q，则返回步骤(3)；否则，执行步骤(9)；

(9)判断抑噪后的迭代复原结果图x^(k+1)与迭代复原结果图x^(k)的差值

是否小于迭代误差ε，如果

则x^(k+1)就是最终的复原结果图；否则，令迭代次数k＝k+1，返回步骤(6)，直到满足

为止。

本发明构造了全信息非局部权重系数矩阵的计算公式，并将它作为约束计算复原中间结果图，与传统的基于全变分的图像复原方法相比，不仅解决了恢复图像在平坦区域容易产生阶梯效应的问题，而且可以恢复出更多的高频细节信息。

附图说明

图1是本发明的流程图；

图2是本发明在仿真实验中使用的Cameraman清晰图像；

图3是本发明在仿真实验中使用的Cameraman模糊图像；

图4是本发明对图3进行复原仿真得到的Cameraman复原结果图；

图5是现有的基于小波的阈值迭代法对图3进行复原仿真得到的Cameraman复原结果图；

图6是现有的基于全变分的阈值迭代法对图3进行复原仿真得到的Cameraman复原结果图。

具体实施方式

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

步骤1，用现有的“非局部均值滤波法”对输入的模糊图像y进行噪声抑制，得到抑噪后的模糊图像x^(-1)，其中，“非局部均值滤波法”由A.Buades等人在《A non-localalgorithm for image denoising》，IEEE Int.Conf.on Computer Vision and PatternRecognition(CVPR 2005)，San Diego，CA，USA，June 20-25，2005中提出，具体按如下公式计算：

$x^{(-1)}\left(i\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}g(i,j)y\left(j\right)$

其中，x^(-1)(i)为抑噪后的模糊图像x^(-1)的第i个像素值，y(j)为模糊图像y的第j个像素值，g(i，j)为像素值y(i)和像素值y(j)之间的权重系数，y(i)为模糊图像y的第i个像素值，

y_b(i)和y_b(j)分别为模糊图像y按照从上到下，从左到右的顺序，划分成的第i个和第j个大小为f×f的图像小块，f为图像小块的边长，f＝5，d(i)为归一化因子，

p为调节参数，p＝30。

步骤2，设定迭代误差ε＝1×10^-6，设定当前的迭代次数k＝0，用现有的“维纳滤波法”对抑噪后的模糊图像x^(-1)进行滤波，得到初步复原结果图x⁽⁰⁾，用初步复原结果图x⁽⁰⁾初始化迭代复原结果图x^(k)。

这里采用的“维纳滤波法”由Helstrom C.W.等人在《Image restoration by themethod of least squares》，J.Opt.Soc.Amer.，1967，Vol.57，No.3，pp：297～303中提出，其实现步骤如下：

(2a)用如下公式求解初始复原结果图x⁽⁰⁾的傅立叶变换结果

${\hat{x}}^{\left(0\right)}=\frac{{\hat{H}}^{*}{\hat{x}}^{(-1)}}{{\left|\hat{H}\right|}^2+k}$

式中，

为抑噪后的模糊图像x^(-1)的傅立叶变换结果，

为模糊核函数矩阵H的傅立叶变换结果，

为傅立叶变换结果

的共轭转置矩阵，k为调节参数，k＝0.1；

(2b)对傅立叶变换结果

进行逆傅立叶变换，得到初始复原结果图x⁽⁰⁾：

$x^{\left(0\right)}=F^{-1}\left({\hat{x}}^{\left(0\right)}\right)$

式中，

表示将傅立叶变换结果

进行逆傅立叶变换。

步骤3，计算迭代复原结果图x^(k)的非局部权重系数矩阵

其中，W₁(i，j)为非局部权重系数矩阵W₁的第i行，第j列的元素值，i＝1，2，...，N，j＝1，2，...，N，N为迭代复原结果图x^(k)的像素总数，令W₁(i，j)的计算公式为：

式中，

表示迭代复原结果图x^(k)的第i个f×f像素的图像小块x_i和第j个f×f像素的图像小块x_j之间的权重系数，f为图像小块的边长，f＝5，

c_i为归一化因子，h为调节权重大小的参数，h＝65，a_i是一个列向量，它包含x_i所对应的最大的10个权重系数

步骤4，计算迭代复原结果图x^(k)的贝叶斯非局部权重系数矩阵

其中，W₂(i，j)为贝叶斯非局部权重系数矩阵W₂的第i行，第j列的元素值，令W₂(i，j)的计算公式为：

其中，表示迭代复原结果图x^(k)的第i个f×f像素的图像小块x_i和第j个f×f像素的图像小块x_j之间的权重系数， $b_i^j=\frac{1}{z_i}\exp (-\frac{1}{2}{(\frac{\left|\right|x_i-x_j\left|\right|}{\mathrm{\σ}}-\sqrt{2f-1})}^2),$ z_i为归一化因子， $Z_i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\exp (-\frac{1}{2}{(\frac{\left|\right|x_i-x_j\left|\right|}{\mathrm{\σ}}-\sqrt{2f-1})}^2),$ σ为已知噪声标准差，b_i是一个列向量，它包含x_i所对应的最大的10个权重系数

步骤5，计算迭代复原结果图x^(k)的全信息非局部权重系数矩阵W＝rW₁+(1-r)W₂，其中，r为调节参数，r＝0.4。

步骤6，构建全信息非局部约束的阈值迭代公式，并利用该式计算复原中间结果图x^(k+1/2)。

(6a)将全信息非局部权重系数矩阵W和迭代复原结果图x^(k)的乘积的二范数模值

作为全信息非局部正则项，加入到由J.Bioucas-Dias等在文献《Totalvariation image deconvolution：A majorization-minimization approach》，presented at theIEEE Int.Conf.Acoustics，Speech，and Signal Processing，Toulouse，France，2006中提出的基于全变分的图像复原模型中，得到全信息非局部约束的全变分图像复原模型：

$x^{(k+1/2)}=\arg \underset{x^{\left(k\right)}}{\min }\{{\left|\right|y-{\mathrm{Hx}}^{\left(k\right)}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|x^{\left(k\right)}\left|\right|}_{\mathrm{TV}}+\mathrm{\&gamma;}{\left|\right|W{x}^{\left(k\right)}\left|\right|}_2^2\}---<<1>>$

其中，

表示求解当大括号{}里的式子达到最小值时，所得到迭代复原结果图x^(k)，

为惩罚项，λ为调节参数，||x^(k)||_TV为TV正则项；

(6b)将公式《1》中的惩罚项和TV正则项进行合并，得到：

$x^{(k+1/2)}=\arg \underset{x^{\left(k\right)}}{\min }\{{\left|\right|\left[\begin{array}{c}y\\ 0\end{array}\right]-\left[\begin{array}{c}H\\ \mathrm{\&gamma;W}\end{array}\right]x^{\left(k\right)}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|x^{\left(k\right)}\left|\right|}_{\mathrm{TV}}\}---<<2>>$

令

将公式《2》改写成：

(6c)阈值迭代公式求解公式《3》，得到：

其中，

为

的转置矩阵，该阈值迭代公式由I.Daubechies等在文献《Aniterative thresholding algorithm for linear inverse problems with a sparsity constraint》，Commun.Pure Appl.Math.，2004，Vol.57，No.11，pp.1413-1457中提出；

(6d)将公式《4》展开，得到计算迭代复原中间结果图x^(k+1/2)的展开式：

《5》

令变形的模糊核矩阵U＝H^TH，令变形的全信息非局部系数矩阵V＝γ²W^TW，其中，调节参数γ＝0.1，得到全信息非局部约束的阈值迭代公式：

x^(k+1/2)＝x^(k)+(H^Ty-Ux^(k)-Vx^(k))                    《6》

(6e)利用公式《6》，计算得到迭代复原中间结果图x^(k+1/2)，即将迭代复原结果图x^(k)、模糊核转置矩阵H^T、变形的模糊核矩阵U和变形的非局部系数矩阵V代入到公式《6》，计算后得到迭代复原中间结果图x^(k+1/2)。

步骤7，利用现有的“基于全变分模型的去噪方法”对迭代复原中间结果图x^(k+1/2)进行噪声抑制，得到抑噪后的迭代复原结果图x^(k+1)，其中，基于全变分模型的去噪方法由A.Chambolle在文献《An algorithm for total variation minimization andapplication》，Journal of Mathematical Imaging and Vision，2004，Vol.20，No.1-2，pp.89-97中给出。

步骤8，判断迭代次数k是否等于更新标记数q，q＝100，200，...，1000，如果k＝q，则返回步骤(3)；否则，执行步骤(9)；

步骤9，判断抑噪后的迭代复原结果图x^(k+1)与迭代复原结果图x^(k)的差值

是否小于等于迭代误差ε，如果则x^(k+1)即为最终的输出复原结果图；否则，令k＝k+1，并返回步骤(6)，直到满足为止。

本发明是一种迭代的图像复原方法，每次迭代得到的抑噪后的迭代复原结果图x^(k+1)均发生变化，当抑噪后的迭代复原结果图x^(k+1)与迭代复原结果图x^(k)的差值

时，则抑噪后的迭代复原结果图x^(k+1)即为最终的复原结果图。

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

1.实验条件：实验所用计算机的CPU为Intel Core2 Duo 2.33GHz，内存为2GB，编程平台为Matlab R2009a。实验所用到的图像来源于标准图像库，分别为Cameraman，Peppers，Boas，图像大小均为256×256。

实验参数设置为：模糊核函数矩阵H＝[1 4 6 4 1]^T[1 4 6 4 1]/256，噪声为高斯白噪声，标准差为σ＝2。

2.实验内容

用本发明和现有的基于小波的阈值迭代法、现有的基于全变分的阈值迭代法分别对图3所示的Cameraman模糊图像进行复原仿真，其中用本发明得到的Cameraman复原结果如图4，用现有的基于小波的阈值迭代法得到的Cameraman复原结果如图5；用现有的基于全变分的阈值迭代法得到的Cameraman复原结果如图6。

仿真实验中，应用峰值信噪比PSNR评价指标来评价复原结果的优劣，其PSNR的定义为：

$\mathrm{PSNR}=10{\log }_{10}\left(\frac{{255}^2\&times;M\&times;N}{\mathrm{\&Sigma;}{\left|\right|x-f\left|\right|}^2}\right)$

其中，f为清晰图像，x为复原后的图像，M和N为清晰图像f的像素行数和像素列数。

用本发明和现有的基于小波的阈值迭代法、基于全变分的阈值迭代法，分别对图像Cameraman，Peppers和Boats进行模糊图像复原仿真。应用峰值信噪比PSNR对复原结果图进行评价，评价结果如表1所示，其中，Alg1是本发明的方法，Alg2是基于小波的阈值迭代法，Alg3是基于全变分的阈值迭代法。

表1.本发明和两种对比方法在仿真实验中得到的PSNR值(单位为dB)

3.实验结果分析

对比本发明得到的复原结果图4、基于小波的阈值迭代法得到的复原结果图5和基于全变分的阈值迭代法得到的复原结果图6可以看出，图4所示的本发明得到的Cameraman复原结果不但有效地去除了模糊，使图像边缘清晰，同时还保留了更多的图像细节，更接近图2所示的Cameraman清晰图像；图5所示的基于小波的阈值迭代法得到的复原结果中残留了许多噪声，并且产生了严重的锯齿和振铃效应；图6所示的基于全变分的阈值迭代法能够有效地抑制噪声，但是其复原结果过于平滑，丢失了图像的细节。

从表1中可以看出，本发明比其它两种对比方法具有更高的PSNR值，有更好的复原性能。

Claims (4)

1.一种用于图像复原的全信息非局部约束全变分方法，包括如下步骤： 

(1)用非局部均值滤波法对输入的模糊图像y进行噪声抑制，得到抑噪后的模糊图像x^(-1)； 

(2)设定迭代误差ε＝1×10^-6，设定当前的迭代次数k＝0，用维纳滤波法对抑噪后的模糊图像x^(-1)进行滤波，得到初始复原结果图x⁽⁰⁾，用初始复原结果图x⁽⁰⁾初始化迭代复原结果图x^(k)； 

(3)计算迭代复原结果图x^(k)的非局部权重系数矩阵

其中，W₁(i,j)为非局部权重系数矩阵W₁的第i行，第j列的元素值，i＝1,2,...,N，j＝1,2,...,N，N为迭代复原结果图x^(k)的像素总数，令W₁(i,j)的计算公式为： 

其中，

表示迭代复原结果图x^(k)的第i个f×f像素的图像小块x_i和第j个f×f像素的图像小块x_j的非局部相似权重系数，f为图像小块的边长，f＝5，

为迭代复原结果图像小块的非局部相似权重系数，

c_i为归一化因子， 

h为调节权重大小的参数，h＝65，a_i是一个列向量，它包含x_i所对应的最大的10个非局部相似权重系数

(4)计算迭代复原结果图x^(k)的贝叶斯非局部权重系数矩阵

其中，W₂(i,j)为贝叶斯非局部权重系数矩阵W₂的第i行，第j列的元素值，令W₂(i,j)的 计算公式为： 

其中，

表示迭代复原结果图x^(k)的第i个f×f像素的图像小块x_i和第j个f×f像素的图像小块x_j的贝叶斯非局部相似权重系数，f为图像小块的边长，

为迭代复原结果图像小块的贝叶斯非局部相似权重系数，

z_i为归一化因子，

σ为已知噪声标准差，b_i是一个列向量，它包含x_i所对应的最大的10个贝叶斯非局部相似权重系数

(5)计算迭代复原结果图x^(k)的全信息非局部权重系数矩阵W=rW₁+(1-r)W₂，其中，r为调节参数，r＝0.4； 

(6)按照如下的全信息非局部约束阈值迭代公式，计算迭代复原中间结果图x^(k+1/2)： 

x^(k+1/2)＝x^(k)+(H^Ty-Ux^(k)-Vx^(k)) 

其中，U为变形的模糊核矩阵，U＝H^TH，H为已知的模糊核函数矩阵，H^T为H的转置矩阵，V为变形的非局部系数矩阵，V＝γ²W^TW，W^T为权重系数矩阵W的转置矩阵，γ为调节参数，γ＝0.1； 

(7)用基于全变分模型的去噪方法对迭代复原中间结果图x^(k+1/2)进行噪声抑制，得到抑噪后的迭代复原结果图x^(k+1)； 

(8)判断迭代次数k是否等于更新代数q，q＝100,200,...,1000，如果k＝q，则返回步骤(3)；否则，执行步骤(9)； 

(9)判断抑噪后的迭代复原结果图x^(k+1)与迭代复原结果图x^(k)的差值

是否小于迭代误差ε，如果

则x^(k+1)就是最终的复原结果图；否则，令 迭代次数k＝k+1，返回步骤(6)，直到满足为止。 

2.根据权利要求1所述的全信息非局部约束全变分方法，其中步骤(1)所述的用非局部均值滤波法对输入的模糊图像y进行噪声抑制，按如下公式计算： 

    《1》 

其中，x^(-1)(i)为抑噪后的模糊图像x^(-1)的第i个像素值，y(j)为模糊图像y的第j个像素值，g(i,j)为像素值y(i)和像素值y(j)之间的权重系数，y(i)为模糊图像y的第i个像素，g(i,j)为模糊图像小块的相似权重系数， 

y_b(i)和y_b(j)分别为模糊图像y按照从上到下，从左到右的顺序，划分成的第i个和第j个大小为f×f的图像小块，f为图像小块的边长，d(i)为归一化因子，

p为调节参数，p＝30。 

3.根据权利要求1所述的全信息非局部约束全变分方法，其中步骤(2)所述的用维纳滤波法对抑噪后的模糊图像x^(-1)进行滤波，按如下步骤进行： 

(3a)用如下公式求解初始复原结果图x⁽⁰⁾的傅立叶变换结果

    《2》 

其中，

为抑噪后的模糊图像x^(-1)的傅立叶变换结果，

为模糊核函数矩阵H的傅立叶变换结果，

为傅立叶变换结果

的共轭转置矩阵，k为调节参数，k＝0.1； 

(3b)对傅立叶变换结果

进行逆傅立叶变换，得到初始复原结果图x⁽⁰⁾： 

    《3》 

其中，

表示将傅立叶变换结果进行逆傅立叶变换。 

4.根据权利要求1所述的全信息非局部约束全变分方法，其中步骤(6)所述的计算迭代复原中间结果图x^(k+1/2)的公式，按如下步骤构建： 

(4a)将全信息权重系数矩阵W和迭代复原结果图x^(k)的乘积的二范数模值

作为非局部正则项，加入到基于全变分的图像复原模型中，得到全信息非局部约束全变分图像复原模型： 

    《4》 

其中，

表示求解当大括号{}里的式子达到最小值时，所得到迭代复原结果图x^(k)，

为惩罚项，λ为调节参数，

为TV正则项； 

(4b)将公式《4》中的惩罚项

和TV正则项

进行合并，得到： 

    《5》 

令

将公式《5》改写成： 

    《6》 

(4c)按照阈值迭代公式求解公式《6》，得到： 

    《7》 

其中，

为

的转置矩阵； 

(4d)将公式《7》展开，得到计算迭代复原中间结果图x^(k+1/2)的展开式： 

           《8》 

令变形的模糊核矩阵U＝H^TH，变形的非局部系数矩阵V＝γ²W^TW，得到步骤(6)中所述的全信息非局部约束的阈值迭代公式： 

x^(k+1/2)＝x^(k)+(H^Ty-Ux^(k)-Vx^(k))    《9》 。

Images