CN101359398A - 运动模糊图像的盲恢复方法 - Google Patents

Abstract

运动模糊图像的盲恢复方法，步骤如下：(1)采用倒谱法对图像进行变换，求出模糊图像的模糊尺度和模糊方向；(2)利用步骤(1)求出的模糊图像的模糊尺度和模糊方向，采用基于TV总变分方法对模糊图像进行恢复。基于TV总变分方法为：采用定点迭代法作为外循环，共轭梯度法作为内循环，进行循环迭代，得到恢复的图像。本发明具有强自适应，强抗噪能力和强鲁棒性等优点，仿真图片和实拍图片的处理效果表明本方法具有有效性和实用性。

Description

运动模糊图像的盲恢复方法

技术领域

本发明涉及一种运动模糊图像的盲恢复方法，属于图像处理领域。

背景技术

随着信息技术和计算机的迅猛发展，数字图像处理得到了广泛的应用和发展，但由于目标与成像系统间的相对运动，造成图像降晰，特别是当相对运动速度过快时，会严重影响成像的质量，造成图像的模糊拖尾，因而使图像识别，目标跟踪和目标检测等增加了难度和不确定性。

目前，尽管很多文献资料上也在研究和讨论解决运动模糊图像的方法，如去卷积(逆滤波)方法、维纳滤波、投影法、最大熵恢复等，这些方法对噪声比较敏感，鲁棒性较差，其中致命的弱点是都要知道点扩散函数PSF(Point Spread Function)才能进行复原，而在实际过程中点扩散函数PSF一般都是未知的，这就限制了这些方法的适用范围。

针对此种情况，图像处理的研究工作者提出了一系列的图像盲复原方法，如基于最大似然函数的Richardson-Lucy算法的图像盲复原方法，参见Richardson，W.H，“Bayesian-based Iterative method of image restoration”，Journal of opticalsociety of America，1972，Vol.62，P55-59和L.B.Lucy，“An iterative techniquefor rectification of observed distributions”，The Astronomical journal，1974，Vol.79，P745-754；基于小波变换的图像盲复原方法，参见汪雪林、韩华、彭思龙，“基于小波域局部高斯模型的图像复原”，软件学报，2004，Vol.15，No.3，443-450。其中基于最大似然函数的图像盲复原方法对初始的PSF的估计要求较高，当初始的估计值与真实的PSF相差较远时，图像复原效果很差。基于小波变换的图像盲复原是假设图像的小波变换系数是服从局部高斯模型的，但对于航空航天图像和医学图像而言这个假设一般不成立，并且该方法计算量大，耗时长。基于最大最小错误率EMM的图像盲复原方法得到的解是局部极小值，不一定是全局最优，并且解和初始化时的初值有很大的关系，参见李红阳、卓晴、王文渊，“EMM盲运动模糊图像的恢复”，计算机工程与应用，2003，P111-115。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种运动模糊图像的盲恢复方法，该方法不但能使复原图像的清晰度提高，同时具有良好的自适应强能力，抗噪能力，和鲁棒性，而且实际证明本方法具有强的实用性。

本发明的技术解决方案：运动模糊图像的盲恢复方法，其步骤如下：

(1)采用倒谱法对图像进行变换，求出模糊图像的模糊尺度和模糊方向；

(2)采用基于TV总变分方法进行对步骤(1)求出的模糊图像的模糊尺度和模糊方向进行模糊图像恢复。

所述步骤(1)中采用倒谱法对图像进行变换，求出模糊图像的模糊尺度和模糊方向的步骤如下：

(a)取G(u，v)的幅度谱|G(u，v)|，为了保证|G(u，v)|＝0时有意义，按下式取对数得

$\hat{G}(u,v)=\ln (1+|G(u,v)\left|\right)$

(b)对

进行反傅里叶得模糊图像g(x，y)的倒谱

(c)将

分成大小相等的四部分，并按对角线方向进行相互交换，使原点移到中心位置。

(d)确定模糊尺度和模糊方向。

取交换后的右半部分数据用来确定倒谱最小值位置(i，j)即负最大值的位置，则模糊尺度L：

$L=\mathrm{integer}(\sqrt{i^2+j^2}-1+0.5)$

模糊方向θ：

$\mathrm{\θ}=\mathrm{integer}(\arctan \left(\frac{j}{i}\right)\×\left(\frac{180}{\mathrm{\π}}\right)+0.5)$

所述的步骤(2)基于TV总变分方法为：采用定点迭代法作为外循环，共轭梯度法作为内循环，进行循环迭代，得到恢复的图像，具体实现步骤为：

(a)利用得到的模糊尺度和模糊方向，采用基于运动模糊的模型计算出点扩散函数PSF；

(b)令f＝g，确定定点迭代次数m、共轭梯度次数j、正则化参数α、常数β，其中：f为复原图像，g为原始图像；

(c)定点迭代，判断定点迭代次数是否等于0，是则转步骤(e)，否则令；

其中： $\▿J\left(f^m\right)=h^{\⊗}*(h*f^m-g)+\mathrm{\αL}\left(f^m\right)f^m,$ h为点扩散函数，

表示h的伴随算子，f^m是第m次定点迭代的估计图像， $L\left(f^m\right)f^m=-\▿\·\left[\frac{\▿f^m}{\sqrt{{\left|{\▿f}^m\right|}^2+{\mathrm{\β}}^2}}\right],$ ▽·v表示向量v的散度， $\▿f^m=(f_x^m,f_y^m),$ $\tilde{H}\left(f^m\right)=h^{\⊗}*h+\mathrm{\αL}\left(f^m\right).$

(d)判定共轭梯度次数j是否等于0，不等于0，则j＝j-1，同时继续循环迭代修改，否则令

并转步骤(c)；

(e)输出复原图像f。

本发明与现有技术相比的优点在于：

(1)本发明不同于一般的盲复原方法，由于采用了基于倒谱法的运动参数估计和基于总变分的图像复原方法，克服了由于运动模糊造成的模糊拖尾问题，使得图像的盲复原的清晰度大大提高。

(2)本发明采用的倒谱估计方法具有快速的衰减特性和抑制噪声干扰的特点决定了该方法具有估计准确，误差少，实用性强的优点。而现有的基于频率的估计方法由于对噪声比较敏感，使得在实际过程中估计的准确率低，误差大，导致了实用性差。而基于小波的方法由于它是假设图像的小波系数是遵循局部高斯模型的，但实际中大部分图像不遵循这个模型，例如医学图像、卫星图像等，对这些图像处理会导致恢复效果较差，因此具有很大的局限性，不具有实用性。

(3)本发明采用了总变分的方法。因为总变分方法不但能保留图像细节和纹理，而且能抑制噪声的影响，而一般的复原方法如逆滤波、维纳滤波等受噪声的影响大，同时这些方法会使复原后的图像变得平滑，不能保留图像的细节和纹理。

(4)本发明对运动模糊图像而言不要求任何先验知识，具有很强的适用性。

(5)本发明能保留图像细节和纹理特征，更适合于医学图像、卫星图像等运动图像的模糊处理。

总之，本发明具有自适应强，抗噪声能力好和强鲁棒性等优点，仿真图片和实拍图片的处理效果表明本发明具有有效性和实用性。

附图说明

图1为本发明方法的流程图；

图2为图像模糊模型；

图3为H(u)的频谱L＝30，N＝256；

图4为h(x，y)的倒谱L＝30，N＝256；

图5为本发明中的数据交换示意图；

图6为原始图像；

图7为本发明的仿真图像及处理效果图，其中图7a是模糊尺度为30，方向为135°仿真图像，图7b是其处理后的图像，图7c是模糊尺度为30，方向为135°，所加的高斯白噪声的方差是58.65的仿真图像，图7d是其处理后的图像；

图8为本发明的运动速度为50cm/s和60cm/s，曝光时间为0.8s和1s时的实拍和复原效果图，其中8a、8b、8c、8d、8e、8f、8g、8h为实拍图，8a’、8b’、8c’、8d’、8e’、8f’、8g’、8h’为复原图。

具体实施方式

下面结合具体实施方式及附图对本发明进一步详细说明。

如图1所示，本发明包括两个部分：采用倒谱法对图像进行变换，求出模糊图像的模糊尺度和模糊方向；采用基于TV总变分方法进行模糊图像恢复。

1.运动图像模糊模型

本发明采用的图像模糊模型如图2所示，其公式如1所示：

            g(x，y)＝h(x，y)*f(x，y)+n(x，y)                       (1)

其中g(x，y)表示模糊图像，h(x，y)表示点扩散函数，f(x，y)表示原始图像，n(x，y)表示噪声，*表示卷积运算，以下同。

对(1)式进行傅里叶变换可得其频率模型如下：

            G(u，v)＝H(u，v)F(u，v)+N(u，v)                        (2)

在实际拍摄过程中，由于相机的曝光时间相当短，可以认为运动物体做匀速直线运动，则运动方向为θ，模糊长度为L的点扩散函数可以表示为

$h(x,y)=\left\{\begin{array}{cc}\frac{1}{L}& \mathrm{if}\sqrt{x^2+y^2}\≤\mathrm{L\; and}\tan \mathrm{\θ}=\frac{y}{x}\\ 0& \mathrm{otherwise}\end{array}\right.---\left(3\right)$

其中L表示模糊长25度，θ为与x轴正方向成的角度。

2.基于倒谱法的点扩散函数辨识机理

2.1倒谱的定义和性质

根据应用背景的不同，倒谱的定义也是不同的，本文采用的倒谱定义如下：

                                                            (4)

其中|C(u，v)| $\hat{c}(x,y)=F^{-1}\left(\ln \right|\left(C(u,v)\right|)$ 表示c(x，y)的傅里叶变换的幅度谱，

表示与c(x，y)相对应的倒谱，为了表述方便将其记为：

倒谱具有很多重要的性质，受篇幅的限制，只讨论跟本文有关的二维倒谱的性质并假定其二维倒谱是存在的。

性质1：将二维卷积运算变成加法运算

若

g(x，y)＝h(x，y)*f(x，y)                                 (6)

则

$\hat{g}(x,y)=\hat{h}(x,y)+\hat{f}(x,y)---\left(7\right)$

性质2：旋转性

若在极坐标系中，有

则

该性质表明函数c(r，θ)旋转角度α后所得的倒谱等于将其倒谱

朝同一个方向旋转同样的角度α。

性质3：原点对称性

如果c(x，y)是实函数，则

具有原点对称性。即：

$\hat{c}(x,y)=\hat{c}(-x,-y)---\left(9\right)$

性质4：周期性

即 $\hat{c}(x+M,y+N)=\hat{c}(x,y)---\left(10\right)$

其中M，N分别为图像的行数和列数。

良好的性质使得倒谱方法得到了广泛的应用，本文所提的对称性和周期性，不仅可以简化计算，同时也为本文算法的可靠性提供了理论了论据。

2.2基于倒谱法的点扩散函数辨识原理

为分析方便，假设运动模糊图像的模糊长度为L，角度θ为0，则

对式(11)进行傅里叶变换得

$H_c\left(u\right)=\frac{\sin \left(\frac{\mathrm{uL}}{2}\right)}{\frac{\mathrm{uL}}{2}}---\left(12\right)$

其离散的傅里叶变换为

$H\left(u\right)=\frac{e^{-j\frac{\mathrm{\πLu}}{N}}\sin \left(\frac{\mathrm{Lu\π}}{N}\right)}{e^{-j\frac{\mathrm{\πu}}{N}}L\sin \left(\frac{\mathrm{u\π}}{N}\right)},0\≤u\≤N-1---\left(13\right)$

其频谱图如图3所示，从图中可以看出频谱的零点近似周期性地出现。

据式12可知，它的零点出现以

为周期，因此h(x，y)的倒谱会出现较大的负最大值，负最大值会在L与N-L处出现，据此可以确定其模糊尺度。图4中显示的是L＝30，N＝256时h(x，y)的倒谱图，从图4中可看出x＝30和x＝226处出现两个最大的负巅峰，并关于中心对称。

令式(1)中的n(x，y)＝0，由倒谱的性质1可知：

$\hat{g}(x,y)=\hat{h}(x,y)+\hat{f}(x,y)---\left(14\right)$

一般原始清晰图像的傅里叶变换不会周期地出现零点，因此它的倒谱不会在L+1处和N-L+1处产生较大的负最大值，原始图像对

的负最大值的影响可以忽略不计，模糊图像倒谱的负最大值表征了点扩散的特征，所以通过确定

的负最大值位置就可以确定运动模糊图像的模糊尺度。

对于模糊方向为θ的运动模糊，根据以上分析，由性质2可知，在模糊长度为L，方向为θ的地方将出现负最大值，据此可以确定模糊方向和模糊尺度。

3.基于倒谱的系统辨识算法

基于上面的分析，本发明提出如下的系统辨识算法：

步骤1：对模糊图像g(x，y)进行二维傅里叶变换fft2得到G(u，v)。

步骤2：取G(u，v)的幅度谱|G(u，v)|，为了保证|G(u，v)|＝0时有意义，按式(15)取对数得

$\hat{G}(u,v)=\ln (1+|G(u,v)\left|\right)---\left(15\right)$

步骤3：通过反傅里叶得模糊图像g(x，y)的倒谱

步骤4：将

从左到右从上到下分成大小相等的四部分，分别记为A、B、C、D，将A与D、B与C分别进行数据交换，见图5，使原点移到中心位置。

步骤5：基于性质2，确定模糊尺度和模糊方向。

根据性质3和性质4，取图5所示交换后的右半部分数据，认为是平面坐标系中的1、4象限，用这部分数据确定倒谱最小值位置(i，j)即负最大值的位置，则模糊尺度L：

$L=\mathrm{integer}(\sqrt{i^2+j^2}-1+0.5)---\left(16\right)$

模糊方向θ：

$\mathrm{\θ}=\mathrm{integer}(\arctan \left(\frac{j}{i}\right)\×\left(\frac{180}{\mathrm{\π}}\right)+0.5)---\left(17\right)$

4.基于TV总变分方法的图像复原原理

总变分最小化方法作为一种非线性规整化方法具有良好的棱边保持特性，因而得到广泛关注，被认为是一种对噪声抑制、图像重建和复原有前景的方法。基于总变分的图像复原方法可以归结为条件最小化问题：

min J_T(f)并满足‖h*f-g‖²＝‖n‖²                         (18)

其中

$J_T\left(f\right)={\∫}_{D_f}\sqrt{{\left|\right|\▿f\left|\right|}^2+{\mathrm{\β}}^2}={\∫}_{D_f}\sqrt{f_x^2+f_y^2+{\mathrm{\β}}^2}\mathrm{dxdy}---\left(19\right)$

$f_x=\frac{\∂f}{\∂x},$ $f_y=\frac{\∂f}{\∂y},$ D_f是f的支持域，β是一个为了避免J_T(f)在f_x＝f_y＝0处不可微的常数，利用拉格朗日乘子法，最小化式(18)等价干最小化下面的式子：

$\underset{f}{\min }J\left(f\right)=\underset{f}{\min }\frac{1}{2}{\left|\right|h*f-g\left|\right|}^2+\mathrm{\α}{J}_T\left(f\right)---\left(20\right)$

其中α作为正则化参数，它的取值应满足约束条件‖h*f-g‖²＝‖n‖²。

最小化(20)式导致一个Eular-Lagrange方程，连同Van Neumann边界条件可以表示为

$\left\{\begin{array}{ccc}\▿J\left(f\right)=h^{\⊗}*(h*f-g)+\mathrm{\αL}\left(f\right)f=0,& {\mathrm{\χ}\∈D}_f& \left(21a\right)\\ \frac{\∂f}{\∂n}=0,& \mathrm{\χ}\∈\∂D_f& \left(21b\right)\end{array}\right.$

式中

表示h的伴随算子，x＝[x，y]^T，L(f)是个微分算子，定义为：

$L\left(f\right)w=-\▿\·\left[\frac{\▿w}{\sqrt{{|\▿f|}^2+{\mathrm{\β}}^2}}\right]---\left(22\right)$

其中w为一个函数，▽·v表示向量v的散度，▽w、▽u分别表示其梯度。

式(21)能表示为一个非线性的一阶系统：

$\left\{\begin{array}{c}{h}^{\⊗}*h*f-\mathrm{\α}\▿\·\stackrel{\→}{v}=h^{\⊗}*g\\ -\▿f+\sqrt{{|\▿f|}^2+{\mathrm{\β}}^2}\stackrel{\→}{v}=\stackrel{\→}{0}\end{array}\right.---\left(23\right)$

其中

$\stackrel{\→}{v}=\frac{\▿f}{\sqrt{{|\▿f|}^2+{\mathrm{\β}}^2}}---\left(24\right)$

采用Vogel提出了定点迭代法解式(23)得：

$h^{\⊗}*h*f^{m+1}+\mathrm{\αL}\left(f^m\right)f^{m+1}=h^{\⊗}*g,m=\mathrm{0,1},\·\·\·---\left(25\right)$

上式以半牛顿方式表示为：

$f^{m+1}=f^m-{\left[\tilde{H}\left(f^m\right)\right]}^{-1}\▿J\left(f^m\right)---\left(26\right)$

其中

$\tilde{H}\left(f\right)=h^{\⊗}*h+\mathrm{\αL}\left(f\right)---\left(27\right)$

令

$-{\left[\tilde{H}\left(f^m\right)\right]}^{-1}\▿J\left(f^m\right)=\mathrm{\Δ}{f}^m---\left(28\right)$

则

$\▿J\left(f^m\right)=-\tilde{H}\left(f^m\right)\mathrm{\Δf}---\left(29\right)$

将式(21a)代入式(29)得：

$-\tilde{H}\left(f^m\right)\mathrm{\Δf}=h^{\⊗}*(h*f^m-g)+\mathrm{\αL}\left(f^m\right)f^m---\left(30\right)$

固定f^m对上式采用共轭梯度法求出最优的Δf，则

f^m+1＝f^m+Δf                                       (31)

(31)式表明通过采用定点迭代方法可以恢复出原始图像。

5.基于总变分法的图像复原算法

根据以上理论，本发明提出了基于定点迭代法的恢复方法，具体步骤如下：

步骤1：利用得到的模糊尺度和模糊方向，采用基于运动模糊的模型式(3)计算出点扩散函数PSF；

步骤2：令f＝g，确定定点迭代次数m、共轭梯度次数j、正则化参数α、常数β，其中：f为复原图像，g为原始图像；

步骤3：定点迭代，判断定点迭代次数是否等于0，是则转步骤(5)，否则令

步骤4：利用共轭梯度法求解式(30)中的Δf

Δf＝0

判定共轭梯度次数j是否等于0，不等于0，则j＝j-1，同时继续循环迭代修改

Δf f＝f+Δf，m＝m-1

否则令，并转步骤(3)；

步骤5：输出复原图像f，结束。

6.实验结果及结论

为了验证本发明所提算法的有效性和可靠性，进行了大量的实验。原始图像如图6所示，图7a是模糊尺度为30，方向为135°仿真图像，图7b是其处理后的图像，对比图6可以看出图7b中的天线细节都能得到很好的恢复，说明本文的方法能够很好地保留图像细节和纹理。图7c是模糊尺度为30，方向为135°，所加的高斯白噪声的方差是58.65的仿真图像，图7d是其处理后的图像，从图中可以看出在有较强噪声的情况下，图像细节和纹理仍然能够很好地恢复，说明本发明的抗噪能力强，鲁棒性好。

为了进一步验证本算法的可靠性和适用性，将相机固定在平动或转动的平台上，其相应的曝光时间、平动或转动速度、所拍的图像及处理结果如图8，运动速度为50cm/s到60cm/s，曝光时间为0.8s到1s时的实拍和复原效果图，其中a、b、c、d、e、f、g、h为实拍图，a’、b’、c’、d’、e’、f’、g’、h’为复原图。从图8中可以看出本发明所提的模糊图像盲复原方法效果显著，图像细节和纹理特征恢复得较好，其相应的复原效果与之相对应。

根据以上处理结果可知：本发明所提的方法在图像细节和纹理能恢复有很强的优势，抗噪能力强，鲁棒性好，同时在恢复过程中，不易产生振铃和鬼影效应，之所以具有这么优点，一是采用了新的基于倒谱理论的点扩散函数辨识方法，保证了点扩散函数的准确估计；二是在图像恢复过程中采用了基于总变分的恢复方法，这就保证了图像细节和纹理得到有效地恢复。因此本发明所提的方法具有可靠性和广泛的适用性，势必产生很好的经济效益。

Claims (5)

1、运动模糊图像的盲恢复方法，其特征在于步骤如下：

(1)采用倒谱法对图像进行变换，求出模糊图像的模糊尺度和模糊方向；

(2)采用基于TV总变分方法进行对步骤(1)求出的模糊图像的模糊尺度和模糊方向进行模糊图像恢复。

2、根据权利要求1所述的运动模糊图像的盲恢复方法，其特征在于：所述步骤(1)中采用倒谱法对图像进行变换，求出模糊图像的模糊尺度和模糊方向的步骤如下：

(1)对模糊图像g(x，y)进行二维傅里叶变换fft2得到G(u，v)；

(2)取G(u，v)的幅度谱|G(u，v)|，为了保证|G(u，v)|＝0时有意义，按下式取对数得

$\hat{G}(u,v)=\ln (1+|G(u,v)\left|\right)$

(3)对

进行反傅里叶得模糊图像g(x，y)的倒谱

(4)将

分成大小相等的四部分，并按对角线方向进行相互交换，使原点移到中心位置；

(5)确定模糊尺度和模糊方向

取交换后的右半部分或左半部分数据用来确定倒谱最小值位置(i，j)即负巅峰的位置，对于取右半部分数据的情形，模糊尺度L：

$L=\mathrm{integer}(\sqrt{i^2+j^2}-1+0.5)$

模糊方向θ：

。 $\mathrm{\θ}=\mathrm{integer}(\arctan \left(\frac{j}{i}\right)\×\left(\frac{180}{\mathrm{\π}}\right)+0.5)$

对于取左半部分数据与上类似。

3、根据权利要求1所述的运动模糊图像的盲恢复方法，其特征在于：所述的步骤(2)基于TV总变分方法为：采用定点迭代法作为外循环，共轭梯度法作为内循环，进行循环迭代，得到恢复图像。

4、根据权利要求3所述的运动模糊图像的盲恢复方法，其特征在于：所述的基于TV总变分方法的具体实现步骤为：

(1)利用得到的模糊尺度和模糊方向，采用基于运动模糊的模型计算出点扩散函数PSF；

(2)令f＝g，确定定点迭代次数m、共轭梯度次数j、正则化参数α、常数β，其中：f为复原图像，g为原始图像；

(3)定点迭代，判断定点迭代次数是否等于0，是则转步骤(5)，否则令；

(4)Δf＝0利用共轭梯度法求解下式中的Δf

$\mathrm{\Δf}=-{\left[\tilde{H}\left(f^m\right)\right]}^{-1}\▿J\left(f^m\right)$

其中： $\▿J\left(f^m\right)=h^{\⊗*}(h^{*}{f}^m-g)+\mathrm{\αL}\left(f^m\right)f^m,$ *表示卷积运算，h为点扩散函数，表示h的伴随算子，f^m是第m次定点迭代的估计图像， $L\left(f^m\right)f^m=-\▿\·\left[\frac{{\▿f}^m}{{|\▿f^m|}^2+{\mathrm{\β}}^2}\right],$ 表示向量v的散度， $\▿f^m=(f_x^m,f_y^m),$ $\tilde{H}\left(f^m\right)=h^{\⊗}*h+\mathrm{\αL}\left(f^m\right)$

判定共轭梯度次数j是否等于0，若不等于0，则j＝j-1，同时继续循环迭代修改；否则令Δff＝f+Δf，m＝m-1，并转步骤(3)；

(5)输出复原图像f。

5、根据权利要求4所述的运动模糊图像的盲恢复方法，其特征在于：所述的步骤(1)利用得到的模糊尺度和模糊方向，采用基于运动模糊的模型计算出点扩散函数PSF为：

$h(x,y)=\left\{\begin{array}{cc}\frac{1}{L}& \mathrm{if}\sqrt{x^2+y^2}\≤\mathrm{Land}\tan \mathrm{\θ}=\frac{y}{x}\\ 0& \mathrm{otherwise}\end{array}.\right.$

Images