CN101968881B - 一种运动模糊和散焦复合模糊的图像复原方法 - Google Patents

Abstract

本发明提供一种运动模糊和散焦复合模糊的图像复原方法，该方法能够对运动模糊和散焦同时出现的图像进行参数估计和图像复原，它包含如下步骤：(1)建立高斯白噪声模板，将降质图像与白噪声模板卷积，达到去除噪声的目的；(2)通过图像能量谱估计图像的模糊主方向和副方向；(3)计算图像的主方向导数矩阵和副方向导数矩阵；(4)分别对主方向导数矩阵和副方向导数矩阵实施自相关运算和方向累加运算；(5)根据主方向导数自相关的累加曲线和副方向导数自相关的累加曲线，估计主方向模糊长度和副方向模糊长度；(6)根据获得的主方向模糊长度和副方向模糊长度，建立复合模糊模型；(7)利用维纳滤波对降质图像进行复原。

Description

一种运动模糊和散焦复合模糊的图像复原方法

技术领域

本发明涉及一种运动模糊和散焦复合模糊的图像复原方法，其目的在于当图像中运动模糊和散焦同时存在时，建立复合降质模型，估计模糊参数，并实施图像复原，属于图像处理领域。

背景技术

为了改善图像质量，目前常采用基于模型的图像恢复方法进行处理，这种方法的核心是能够准确获知模糊模型和模糊参数。

运动模糊和散焦是两种常见的图像模糊，运动模糊是由拍摄图像与摄像器材间的相对运动造成的，散焦则是由拍摄器材对焦不准造成的。针对这两种模糊，目前都是分别处理，分别估计其模糊参数并建立模糊模型，常用的模糊参数估计方法有：基于频域的分析方法、基于边缘的分析方法、L-曲线法、自相关法等，其中基于自相关的方法由于抗干扰能力较强、计算量较小因而成为模糊参数估计中的一个热点。

通过对现有技术文献的检索发现现有的模糊参数估计方法都是根据标准运动模糊和散焦模型分别建立的。当运动模糊和散焦同时出现时，没有确定复合模型可以参考，由于模糊相互叠加，上述参数估计方法无法直接应用。本发明将建立一种运动模糊和散焦同时出现的模糊模型，并提供一种对于该模糊参数估计方法，并最终实现图像复原。

发明内容

技术问题

本发明提供一种运动模糊和散焦复合模糊的图像复原方法，能够对运动模糊与散焦同时出现的图像进行复原。

技术方案

一种运动模糊和散焦复合模糊的图像复原方法，包括如下步骤：

步骤1降质图像去噪

利用低通滤波对降质图像g(x，y)进行去噪处理，其中x和y分别为行坐标和列坐标，且都为大于0的整数，得到去噪后的降质图像f′(x，y)，其算法如下：

均值为零的高斯分布p(x，y)的函数形式p_z(i1，i2)可以表示为：

$p_z(i1,i2)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{({i1}^2+{i2}^2)}{{2\mathrm{\σ}}^2}}$

其中下标z表示均值为零的高斯分布的函数形式，i1和i2分别为横坐标和纵坐标，且都为实数，σ为标准差，e为自然对数函数的底数，π为圆周率，设成像***的加性噪声n(x，y)服从均值为零的高斯分布p(x，y)，这里σ取0.5，建立一个3×3的高斯白噪声模板n′(x，y)：

$n^{\′}(x,y)=\left[\begin{array}{ccc}0.0113& 0.0838& 0.0113\\ 0.0838& 0.6193& 0.0838\\ 0.0113& 0.0838& 0.0113\end{array}\right]$

将降质图像g(x，y)与高斯白噪声模板n′(x，y)进行卷积，得到M₁×M₂大小的去噪后的降质图像f′(x，y)：

f′(x，y)＝g(x，y)*n′(x，y)

其中*为卷积操作，

步骤2模糊方向识别

利用公知的傅立叶变换来估计去噪后的降质图像f′(x，y)的主方向θ和副方向θ′，假设此时主方向θ的角度为α，-90°≤α＜90°，副方向θ′的角度为β，0°≤β＜180°，旋转图像至主方向θ的角度为0°，具体方法如下：设散焦与运动模糊同时发生时运动模糊的方向为主方向θ，与主模糊方向垂直的方向为副方向θ′，计算步骤1得到的去噪后的降质图像f′(x，y)的傅立叶变换F(ω₁，ω₂)：

$F({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)=\∫\∫f^{\′}(x,y)e^{-j{\mathrm{\ω}}_{1}x}{e}^{-j{\mathrm{\ω}}_{2}y}\mathrm{dxdy}$

ω₁和ω₂为频率变量，F(ω₁，ω₂)为复数，j为虚数单位，去噪后的降质图像的傅立叶变换F(ω₁，ω₂)可以表示为：

F(ω₁，ω₂)＝R(ω₁，ω₂)+jI(ω₁，ω₂)

其中R(ω₁，ω₂)为实部，I(ω₁，ω₂)为虚部，

计算待估计降质图像f′(x，y)的能量谱E(ω₁，ω₂)：

E(ω₁，ω₂)＝R(ω₁，ω₂)²+I(ω₁，ω₂)²

能量谱E(ω₁，ω₂)中各元素为实数，将能量谱E(ω₁，ω₂)显示出来，在图中心位置会出现一个类似于椭圆的白色亮斑，亮斑周围会有连续的亮环，在所有穿过亮斑的弦中，最长的弦l_max所处的方向，取0°到180°之间的角做为β的值，与此垂直的方向角为α，α＝β-90°，此时主方向θ的角度为α，副方向θ′的角度为β，以垂直于去噪后的降质图像f′(x，y)，且过去噪后的降质图像f′(x，y)中心的直线为旋转轴，如果0°≤α＜90°，以α角顺时针旋转去噪后的降质图像f′(x，y)，如果-90°≤α＜0°，以α+90°角顺时针旋转去噪后的降质图像f′(x，y)，降质图像f′(x，y)旋转后，主方向θ的角度为0°，副方向θ′的角度为90°，得到M₃×M₄大小的旋转后的去噪后的降质图像f″(x，y)；

步骤3图像导数计算

对步骤2中得到的旋转后的去噪后的降质图像f″(x，y)，计算主方向的导数D_a(x，y)和副方向的导数D_b(x，y)：

D_a(x，y)＝f″(x，y)*L_a(x，y)

D_b(x，y)＝f″(x，y)*L_b(x，y)

其中下标a和b分别表示主方向和副方向，L_a(x，y)和L_b(x，y)分别为主方向和副方向的导数卷积模板，此时主方向θ的角度为0°，副方向θ′的角度为90°，主方向卷积模板L_a(x，y)＝[-1，1]，副方向卷积模板

此时主方向的导数D_a(x，y)和副方向的导数D_b(x，y)大小都为M₃×M₄；

步骤4图像导数的自相关计算

步骤4.1利用公知的自相关运算，对步骤3中得到的主方向导数D_a(x，y)和副方向导数D_b(x，y)，分别计算主方向导数D_a(x，y)的自相关C_a(x，y)和副方向导数D_b(x，y)的自相关C_b(x，y)：

$C_a(x,y)=D_a(x,y)\⊗D_a(x,y)$

$C_b(x,y)=D_b(x,y)\⊗D_b(x,y)$

其中下标a和b分别表示主方向和副方向，

表示相关操作，此时主方向导数的自相关C_a(x，y)和副方向导数的自相关C_b(x，y)的大小都为(2×M₃-1)×(2×M₄-1)；

步骤4.2对主方向导数的自相关C_a(x，y)按照与主方向θ垂直的方向进行累加得AC_a(x，y)，对副方向导数的自相关C_b(x，y)按照与副方向θ′垂直的方向进行累加得AC_b(x，y)，累加方法如下：

${\mathrm{AC}}_a(m,n)=C_a(1,n)+C_a(2,n)+L+C_a(2\×M_3-1,n){|}_{m=\mathrm{1,1}\≤n\≤2\×M_4-1}$

${\mathrm{AC}}_b(m,n)=C_b(m,1)+C_b(m,2)+L+C_b(m,2\×M_4-1){|}_{1\≤m\≤2\×M_3-1,n=1}$

其中m和n分别为行坐标和列坐标，且都为大于0的整数，主方向导数自相关的累加AC_a(x，y)的大小为1×(2×M₄-1)，副方向导数的自相关的累加AC_b(x，y)的大小为(2×M₃-1)×1；

步骤5模糊参数识别

设主方向模糊长度为d，副方向模糊半径为r，对步骤4中得到的主方向导数自相关的累加AC_a(x，y)和副方向导数的自相关的累加AC_b(x，y)分别绘制曲线图，在主方向导数自相关的累加AC_a(x，y)曲线图中将出现一个正的最高峰k_ap，在正最高峰两侧，出现负的左侧最低峰k_ap1和负的右侧最低峰k_ap2，下标ap，ap1和ap2分别表示主方向导数自相关的累加AC_a(x，y)曲线图中正的最高峰，负的左侧最低峰和负的右侧最低峰，两负最低峰k_ap1和k_ap2之间的距离s₁的值为2d+4r，在AC_a(x，y)曲线图中将出现一个正的最高峰k_bp，在正最高峰两侧，出现负的左侧最低峰k_bp1和负的右侧最低峰k_bp2，下标bp，bp1和bp2分别表示副方向导数自相关的累加AC_b(x，y)曲线图中正的最高峰，负的左侧最低峰和负的右侧最低峰，两负最低峰k_bp1和k_bp2之间的距离s₂的值为4r，计算主方向模糊长度d和副方向模糊半径r：

$\left\{\begin{array}{c}{S}_1=2d+4r\\ S_2=4r\end{array}\right.\⇒\left\{\begin{array}{c}d=(S_1-S_2)/2\\ r=S_2/4\end{array}\right.$

步骤6建立复合模糊模型

根据步骤5得到的主方向模糊长度d和副方向模糊半径r，通过对标准运动模糊模型h_q(x，y)和标准散焦模型h_p(x，y)进行线性运算和均一化处理来建立复合模糊模型h(x，y)，下标q和p分别表示运动模糊和散焦，建立复合模糊模型h(x，y)的方法如下：

步骤6.1对标准运动模糊模型h_q(x，y)与标准散焦模型h_p(x，y)进行卷积操作，得到复合模糊模型的理论模型h_v(x，y)，下标v表示复合模糊模型的理论模型，具体方法如下：

标准散焦模型h_p(x，y)的函数形式h_pt(i1，i2)为：

其中i1，i2为横坐标和纵坐标，且都为实数，下标pt表示标准散焦模型的数学模型，r为副方向模糊半径且以r为散焦模糊半径，r为大于零的实数，标准散焦模型h_p(x，y)是标准散焦模型的函数形式h_pt(i1，i2)的矩阵形式，标准散焦模型h_p(x，y)中的任意元素为大于或等于零的实数，

假设运动模糊方向水平，标准运动模糊模型h_p(x，y)的函数形式h_pt(i1，i2)为：

其中i1，i2为横坐标和纵坐标，且都为实数，下标qt表示标准运动模糊模型的数学模型，d为主方向模糊长度且以d为运动模糊长度，d为大于零的实数，标准散焦模型h_q(x，y)是标准运动模糊模型的函数模型h_qt(i1，i2)的矩阵形式，标准运动模糊模型h_q(x，y)中的任意元素为大于或等于零的实数，

对标准散焦模型h_p(x，y)和标准运动模糊模型h_q(x，y)进行卷积，得到复合模糊模型的理论模型h_v(x，y)，将复合模糊模型的理论模型h_v(x，y)表示为函数形式h_vt(i1，i2)：

其中下标vt表示复合模糊模型的理论模型的函数形式，

步骤6.2复合模糊模型的理论模型h_v(x，y)确定后，在能量守恒的前提下，对复合模糊模型的理论模型h_v(x，y)进行均一化处理，得到最终的复合模糊模型h(x，y)，方法如下：

步骤6.2.1

计算复合模糊模型的理论模型h_v(x，y)中大于零的点的个数w，

步骤6.2.2

将复合模糊模型的理论模型h_v(x，y)中大于零的点所对应的值替换为

得到最终的复合模糊模型h(x，y)，将复合模糊模型h(x，y)表示成函数形式h_t(i1，i2)：

其中下标t表示复合模糊模型的函数形式，

步骤7图像复原

根据步骤6得到的复合模糊模型h(x，y)，利用公知的维纳滤波进行复原，得到得到复原图像

与现有技术相比，本发明具有如下优点：

1.抗噪声能力强。由于采用了自相关方法估计模糊参数，并且在估计参数前进行了去噪预处理，因此本方法本身具有一定的抗干扰能力，能够在一定的噪声容忍度下进行估计。

2.参数估计精度高。在运动模糊和散焦同时存在的情况下，根据二者在模糊主方向和副方向的影响范围不同，分别计算相应导数，对导数进行自相关运算并按相应垂直方向进行累加，在保证估计正确性的同时，提高了估计的准确性。

3.根据标准运动模糊和散焦模型推导出复合模型的数学模型，并根据光的实际散射特性对复合模型进行修正，使该复合模糊模型能够反映运动模糊和散焦同时存在的情况。在实际环境中，同时出现运动模糊和散焦模糊时，对运动模糊和散焦模糊单独进行建模和复原效果不理想的情况下，使用该复合模糊模型进行复原能取得较好的效果。

总之，本发明所述的运动模糊和散焦复合模糊的图像复原方法的抗噪声能力和参数估计精度是令人满意的。

附图说明

图1运动模糊和散焦复合模糊图像复原算法***框架图。

图2主方向θ＝α识别。

图3主方向导数自相关的累加AC_a(x，y)曲线。

图4副方向导数自相关的累加AC_b(x，y)曲线。

图5复合模糊模型的理论模型h_v(x，y)。

图6复合模糊模型h(x，y)。

图7原始清晰图像f(x，y)。

图8降质图像g(x，y)。

图9复原后的图像

具体实施方式

在具体的实施方式中，将结合附图，清楚、完整地描述运动模糊和散焦复合模糊图像复原算法的详细过程：

一种运动模糊和散焦复合模糊的图像复原方法，包括如下步骤：

步骤1降质图像去噪

取f(x，y)为M₁×M₂大小的原始清晰图像，g(x，y)为M₁×M₂大小的降质后图像，h(x，y)为N₁×N₂大小的运动模糊和散焦复合模糊模型，n(x，y)为M₁×M₂大小的加性噪声，其中其中x和y分别为行坐标和列坐标，且都为大于0的整数，M₁、N₁为行数，M₂、N₂为列数，M₁、M₂、N₁和N₂均为大于零的整数，图像降质过程为：

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

其中*为公知的卷积操作，利用低通滤波对降质图像g(x，y)进行去噪处理，得到去噪后的降质图像f′(x，y)，其算法如下：

均值为零的高斯分布p(x，y)的函数形式p_z(i1，i2)可以表示为：

$p_z(i1,i2)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{({i1}^2+{i2}^2)}{{2\mathrm{\σ}}^2}}$

其中下标z表示均值为零的高斯分布的函数形式，i1和i2分别为横坐标和纵坐标，且都为实数，σ为标准差，e为自然对数函数的底数，π为圆周率，设成像***的加性噪声n(x，y)服从均值为零的高斯分布p(x，y)，这里σ取0.5，建立一个3×3的高斯白噪声模板n′(x，y)：

$n^{\′}(x,y)=\left[\begin{array}{ccc}0.0113& 0.0838& 0.0113\\ 0.0838& 0.6193& 0.0838\\ 0.0113& 0.0838& 0.0113\end{array}\right]$

将降质图像g(x，y)与高斯白噪声模板n′(x，y)进行卷积，得到M₁×M₂大小的去噪后的降质图像f′(x，y)：

f′(x，y)＝g(x，y)*n′(x，y)

其中*为卷积操作，

步骤2模糊方向识别

利用公知的傅立叶变换来估计去噪后的降质图像f′(x，y)的主方向θ和副方向θ′，假设此时主方向θ的角度为α，-90°≤α＜90°，副方向θ′的角度为β，0°≤β＜180°，旋转图像至主方向θ的角度为0°，具体方法如下：设散焦与运动模糊同时发生时运动模糊的方向为主方向θ，与主模糊方向垂直的方向为副方向θ′，计算步骤1得到的去噪后的降质图像f′(x，y)的傅立叶变换F(ω₁，ω₂)：

$F({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)=\∫\∫f^{\′}(x,y)e^{-j{\mathrm{\ω}}_{1}x}{e}^{-j{\mathrm{\ω}}_{2}y}\mathrm{dxdy}$

ω₁和ω₂为频率变量，F(ω₁，ω₂)为复数，j为虚数单位，去噪后的降质图像的傅立叶变换F(ω₁，ω₂)可以表示为：

F(ω₁，ω₂)＝R(ω₁，ω₂)+jI(ω₁，ω₂)

其中R(ω₁，ω₂)为实部，I(ω₁，ω₂)为虚部，

计算待估计降质图像f′(x，y)的能量谱E(ω₁，ω₂)：

E(ω₁，ω₂)＝R(ω₁，ω₂)²+I(ω₁，ω₂)²

能量谱E(ω₁，ω₂)中各元素为实数，将能量谱E(ω₁，ω₂)显示出来，图2为模糊方向识别的能量谱示意图，在图中心位置会出现一个类似于椭圆的白色亮斑，亮斑周围会有连续的亮环，在所有穿过亮斑的弦中，红色虚线被两黄色线段所截之间的部分为最长的弦l_max，l_max所处的方向，取0°到180°之间的角做为β的值，与此垂直的方向角为α，α＝β-90°，此时主方向θ的角度为α，副方向θ′的角度为β，以垂直于去噪后的降质图像f′(x，y)，且过去噪后的降质图像f′(x，y)中心的直线为旋转轴，如果0°≤α＜90°，以α角顺时针旋转去噪后的降质图像f′(x，y)，如果-90°≤α＜0°，以α+90°角顺时针旋转去噪后的降质图像f′(x，y)，降质图像f′(x，y)旋转后，主方向θ的角度为0°，副方向θ′的角度为90°，得到M₃×M₄大小的旋转后的去噪后的降质图像f″(x，y)；

步骤3图像导数计算

对步骤2中得到的旋转后的去噪后的降质图像f″(x，y)，计算主方向的导数D_a(x，y)和副方向的导数D_b(x，y)：

D_a(x，y)＝f″(x，y)*L_a(x，y)

D_b(x，y)＝f″(x，y)*L_b(x，y)

其中下标a和b分别表示主方向和副方向，L_a(x，y)和L_b(x，y)分别为主方向和副方向的导数卷积模板，此时主方向θ的角度为0°，副方向θ′的角度为90°，主方向卷积模板L_a(x，y)＝[-1，1]，副方向卷积模板

此时主方向的导数D_a(x，y)和副方向的导数D_b(x，y)大小都为M₃×M₄，主方向导数D_a(x，y)的函数形式D_at(i1，i2)和副方向导数D_b(x，y)的函数形式D_bt(i1，i2)为：

其中下标at和bt分别表示主方向导数的函数形式和副方向导数的函数形式，i1和i2分别为横坐标和纵坐标，且都为实数；

步骤4图像导数的自相关计算

步骤4.1利用公知的自相关运算，对步骤3中得到的主方向导数D_a(x，y)和副方向导数D_b(x，y)，分别计算主方向导数D_a(x，y)的自相关C_a(x，y)和副方向导数D_b(x，y)的自相关C_b(x，y)：

$C_a(x,y)=D_a(x,y)\⊗D_a(x,y)$

$C_b(x,y)=D_b(x,y)\⊗D_b(x,y)$

其中下标a和b分别表示主方向和副方向，

表示相关操作，此时主方向导数的自相关C_a(x，y)和副方向导数的自相关C_b(x，y)的大小都为(2×M₃-1)×(2×M₄-1)；

步骤4.2对主方向导数的自相关C_a(x，y)按照与主方向θ垂直的方向进行累加得AC_a(x，y)，对副方向导数的自相关C_b(x，y)按照与副方向θ′垂直的方向进行累加得AC_b(x，y)，主方向导数的自相关C_a(x，y)表示为：

$C_a(x,y)=\left[\begin{array}{cccc}{C}_a\left(\mathrm{1,1}\right)& C_a\left(\mathrm{1,2}\right)& L& C_a(\mathrm{1,2}\×M_4-1)\\ C_a\left(\mathrm{2,1}\right)& C_a\left(\mathrm{2,2}\right)& L& C_a(\mathrm{2,2}\×M_4-1)\\ M& M& M& M\\ C_a(2\×M_3-\mathrm{1,1})& C_a(2\×M_4-\mathrm{1,1})& L& C_a(2\×M_3-\mathrm{1,2}\×M_4-1)\end{array}\right]$

主方向导数自相关的累加AC_a(x，y)的大小为1×(2×M₄-1)，且主方向导数自相关的累加AC_a(x，y)中的元素都为实数，每个元素表示为：

AC_a(1，1)＝C_a(1，1)+C_a(2，1)+L+C_a(2×M₃-1，1)

AC_a(1，2)＝C_a(1，2)+C_a(2，2)+L+C_a(2×M₃-1，2)

AC_a(1，2×M₄-1)＝C_a(1，2×M₄-1)+C_a(2，2×M₄-1)+L+C_a(2×M₃-1，2×M₄-1)主方向导数自相关的累加AC_a(x，y)中的元素通项AC_a(m，n)可以表示为：

${\mathrm{AC}}_a(m,n)=C_a(1,n)+C_a(2,n)+L+C_a(2\×M_3-1,n){|}_{m=\mathrm{1,1}\≤n\≤2\×M_4-1}$

其中m和n分别为行坐标和列坐标，且都为大于0的整数，

主方向导数自相关的累加AC_a(x，y)表示为：

AC_a(x，y)＝[AC_a(1，1) AC_a(1，2) L AC_a(1，2×M₄-1)]

副方向导数的自相关C_b(x，y)表示为：

$C_b(x,y)=\left[\begin{array}{cccc}{C}_b\left(\mathrm{1,1}\right)& C_b\left(\mathrm{1,2}\right)& L& C_b(\mathrm{1,2}\×M_4-1)\\ C_b\left(\mathrm{2,1}\right)& C_b\left(\mathrm{2,2}\right)& L& C_b(\mathrm{2,2}\×M_4-1)\\ M& M& M& M\\ C_b(2\×M_3-\mathrm{1,1})& C_b(2\×M_4-\mathrm{1,1})& L& C_b(2\×M_3-\mathrm{1,2}\×M_4-1)\end{array}\right]$

副方向导数自相关的累加AC_b(x，y)的大小为(2×M₃-1)×1，且副方向导数自相关的累加AC_b(x，y)中的元素都为实数，每个元素表示为：

AC_b(1，1)＝C_b(1，1)+C_b(1，2)+L+C_b(1，2×M₄-1)

AC_b(2，1)＝C_b(2，1)+C_b(2，2)+L+C_b(2，2×M₄-1)

AC_b(2×M₃-1，1)＝C_b(2×M₃-1，1)+C_b(2×M₃-1，2)+L+C_b(2×M₃-1，2×M₄-1)

副方向导数自相关的累加AC_b(x，y)中的元素通项AC_b(m，n)可以表示为：

${\mathrm{AC}}_b(m,n)=C_b(m,1)+C_b(m,2)+L+C_b(m,2\×M_4-1){|}_{1\≤m\≤2\×M_3-1,n=1}$

其中m和n分别为行坐标和列坐标，且都为大于0的整数，

主方向导数自相关的累加AC_b(x，y)表示为：

${\mathrm{AC}}_b(x,y)=\left[\begin{array}{c}{\mathrm{AC}}_b\left(\mathrm{1,1}\right)\\ {\mathrm{AC}}_b\left(\mathrm{2,1}\right)\\ M\\ {\mathrm{AC}}_b(2\×M_3-\mathrm{1,1})\end{array}\right]$

步骤5模糊参数识别

设主方向模糊长度为d，副方向模糊半径为r，对步骤4中得到的主方向导数自相关的累加AC_a(x，y)和副方向导数的自相关的累加AC_b(x，y)分别绘制曲线图，主方向导数自相关的累加AC_a(x，y)曲线图如图3所示，在主方向导数自相关的累加AC_a(x，y)曲线图中将出现一个正的最高峰k_ap，在正最最高峰两侧，出现负的左侧最低峰k_ap1和负的右侧最低峰k_ap2，下标ap，ap1和ap2分别表示主方向导数自相关的累加AC_a(x，y)曲线图中正的最高峰，负的左侧最低峰和负的右侧最低峰，两负最低峰k_ap1和k_ap2的横坐标距离s₁的值为2d+4r，副方向导数自相关的累加AC_b(x，y)曲线图如图4所示，在副方向导数自相关的累加AC_b(x，y)曲线图中将出现一个正的最高峰k_bp，在正最高峰两侧，出现负的左侧最低峰k_bp1和负的右侧最低峰k_bp2，下标bp，bp1和bp2分别表示副方向导数自相关的累加AC_b(x，y)曲线图中正的最高峰，负的左侧最低峰和负的右侧最低峰，两负最低峰k_bp1和k_bp2之间的横坐标距离s₂的值为4r，计算主方向模糊长度d和副方向模糊半径r：

$\left\{\begin{array}{c}{S}_1=2d+4r\\ S_2=4r\end{array}\right.\⇒\left\{\begin{array}{c}d=(S_1-S_2)/2\\ r=S_2/4\end{array}\right.$

步骤6建立复合模糊模型

根据步骤5得到的主方向模糊长度d和副方向模糊半径r，通过对标准运动模糊模型h_q(x，y)和标准散焦模型h_p(x，y)进行线性运算和均一化处理来建立复合模糊模型h(x，y)，下标q和p分别表示运动模糊和散焦，建立复合模糊模型h(x，y)的方法如下：

步骤6.1对标准运动模糊模型h_q(x，y)与标准散焦模型h_p(x，y)进行卷积操作，得到复合模糊模型的理论模型h_v(x，y)，下标v表示复合模糊模型的理论模型，具体方法如下：

标准散焦模型h_p(x，y)的函数形式h_pt(i1，i2)为：

其中i1，i2为横坐标和纵坐标，且都为实数，下标pt表示标准散焦模型的数学模型，r为副方向模糊半径且以r为散焦模糊半径，r为大于零的实数，标准散焦模型h_p(x，y)是标准散焦模型的函数形式h_pt(i1，i2)的矩阵形式，标准散焦模型h_p(x，y)中的任意元素为大于或等于零的实数，

假设运动模糊方向水平，标准运动模糊模型h_p(x，y)的函数形式h_pt(i1，i2)为：

其中i1，i2为横坐标和纵坐标，且都为实数，下标qt表示标准运动模糊模型的数学模型，d为主方向模糊长度且以d为运动模糊长度，d为大于零的实数，标准散焦模型h_q(x，y)是标准运动模糊模型的函数模型h_qt(i1，i2)的矩阵形式，标准运动模糊模型h_q(x，y)中的任意元素为大于或等于零的实数，

对标准散焦模型h_p(x，y)和标准运动模糊模型h_q(x，y)进行卷积，得到复合模糊模型的理论模型h_v(x，y)，将复合模糊模型的理论模型h_v(x，y)表示为函数形式h_vt(i1，i2)：

其中下标vt表示复合模糊模型的理论模型的函数形式，复合模糊模型的理论模型h_v(x，y)的空间分布如图5所示；

步骤6.2复合模糊模型的理论模型h_v(x，y)确定后，在能量守恒的前提下，对复合模糊模型的理论模型h_v(x，y)进行均一化处理，得到最终的复合模糊模型h(x，y)，方法如下：

步骤6.2.1

计算复合模糊模型的理论模型h_v(x，y)中大于零的点的个数w：

w＝0；

假设x＝1:N₁

假设y＝1:N₂

如果h_v(x，y)＞0

w＝w+1；

或者w＝w；

结束

结束

结束

步骤6.2.2

将复合模糊模型的理论模型h_v(x，y)中大于零的点所对应的值替换为

假设x＝1:N₁

假设y＝1:N₂

如果h_v(x，y)＞0

$h(x,y)=\frac{1}{w};$

或者h(x，y)＝0；

结束

结束

结束

得到最终的复合模糊模型h(x，y)，将复合模糊模型h(x，y)表示成函数形式h_t(i1，i2)：

复合模糊模型h(x，y)的空间分布如图6所示；

步骤7图像复原

根据步骤6得到的复合模糊模型h(x，y)，利用公知的维纳滤波进行复原，得到得到复原图像方法如下：

$\hat{F}({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)=\frac{1}{H({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)}*\frac{{\left|H({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)\right|}^2}{{\left|H({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)\right|}^2+\mathrm{\γ}}G({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)$

其中ω₁和ω₂为频率变量，G(ω₁，ω₂)为降质模糊图像g(x，y)的傅立叶变换，|H(ω₁，ω₂)|²＝H*(ω₁，ω₂)H(ω₁，ω₂)，H(ω₁，ω₂)为复合模糊模型的傅立叶变换，H*(ω₁，ω₂)是H(ω₁，ω₂)的复共轭，γ是噪信功率比，这里取γ＝0.05，

是复原图像的傅立叶变换，

最后，将复原图像的傅立叶变换

进行傅立叶反变换，方法如下：

$\hat{F}(x,y)=\frac{1}{2\mathrm{\π}}\∫\∫F({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)e^{j{\mathrm{\ω}}_{1}x}{e}^{j{\mathrm{\ω}}_{2}y}d{\mathrm{\ω}}_{1}d{\mathrm{\ω}}_2$

其中j为虚数单位，原始图像f(x，y)如图7所示，降质后的图像g(x，y)如图8所示，复原后的图像

如图9所示。

Claims (1)

1.一种运动模糊和散焦复合模糊的图像复原方法，其特征在于，包括如下步骤：

步骤1降质图像去噪

利用低通滤波对降质图像g(x，y)进行去噪处理，其中x和y分别为行坐标和列坐标，且都为大于0的整数，得到去噪后的降质图像f′(x，y)，其算法如下：

均值为零的高斯分布p(x，y)的函数形式p_z(i1，i2)可以表示为：

$p_z(i1,i2)=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-\frac{({i1}^2+{i2}^2)}{{2\mathrm{\σ}}^2}}$

其中下标z表示均值为零的高斯分布的函数形式，i1和i2分别为横坐标和纵坐标，且都为实数，σ为标准差，e为自然对数函数的底数，π为圆周率，设成像***的加性噪声n(x，y)服从均值为零的高斯分布p(x，y)，这里σ取0.5，建立一个3×3的高斯白噪声模板n′(x，y)：

$n^{\′}(x,y)=\left[\begin{array}{ccc}0.0113& 0.0838& 0.0113\\ 0.0838& 0.6193& 0.0838\\ 0.0113& 0.0838& 0.0113\end{array}\right]$

将降质图像g(x，y)与高斯白噪声模板n′(x，y)进行卷积，得到M₁×M₂大小的去噪后的降质图像f′(x，y)：

f′(x，y)＝g(x，y)*n′(x，y)

其中*为卷积操作，

步骤2模糊方向识别

利用公知的傅立叶变换来估计去噪后的降质图像f′(x，y)的主方向θ和副方向θ′，假设此时主方向θ的角度为α，-90°≤α＜90°，副方向θ′的角度为β，0°≤β＜180°，旋转图像至主方向θ的角度为0°，具体方法如下：设散焦与运动模糊同时发生时运动模糊的方向为主方向θ，与主模糊方向垂直的方向为副方向θ′，计算步骤1得到的去噪后的降质图像f′(x，y)的傅立叶变换F(ω₁，ω₂)：

$F({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)=\∫\∫f^{\′}(x,y)e^{-j{\mathrm{\ω}}_{1}x}{e}^{-j{\mathrm{\ω}}_{2}y}\mathrm{dxdy}$

ω₁和ω₂为频率变量，F(ω₁，ω₂)为复数，j为虚数单位，去噪后的降质图像的傅立叶变换F(ω₁，ω₂)可以表示为：

F(ω₁，ω₂)＝R(ω₁，ω₂)+jI(ω₁，ω₂)

其中R(ω₁，ω₂)为实部，I(ω₁，ω₂)为虚部，

计算待估计降质图像f′(x，y)的能量谱E(ω₁，ω₂)：

E(ω₁，ω₂)＝R(ω₁，ω₂)²+I(ω₁，ω₂)²

能量谱E(ω₁，ω₂)中各元素为实数，将能量谱E(ω₁，ω₂)显示出来，在图中心位置会出现一个类似于椭圆的白色亮斑，亮斑周围会有连续的亮环，在所有穿过亮斑的弦中，最长的弦l_max所处的方向，取0°到180°之间的角做为β的值，与此垂直的方向角为α，α＝β-90°，此时主方向θ的角度为α，副方向θ′的角度为β，以垂直于去噪后的降质图像f′(x，y)，且过去噪后的降质图像f′(x，y)中心的直线为旋转轴，如果0°≤α＜90°，以α角顺时针旋转去噪后的降质图像f′(x，y)，如果-90°≤α＜0°，以α+90°角顺时针旋转去噪后的降质图像f′(x，y)，降质图像f′(x，y)旋转后，主方向θ的角度为0°，副方向θ′的角度为90°，得到M₃×M₄大小的旋转后的去噪后的降质图像f″(x，y)；

步骤3图像导数计算

对步骤2中得到的旋转后的去噪后的降质图像f″(x，y)，计算主方向的导数D_a(x，y)和副方向的导数D_b(x，y)：

D_a(x，y)＝f″(x，y)*L_a(x，y)

D_b(x，y)＝f″(x，y)*L_b(x，y)

其中下标a和b分别表示主方向和副方向，L_a(x，y)和L_b(x，y)分别为主方向和副方向的导数卷积模板，此时主方向θ的角度为0°，副方向θ′的角度为90°，主方向卷积模板L_a(x，y)＝[-1，1]，副方向卷积模板

此时主方向的导数D_a(x，y)和副方向的导数D_b(x，y)大小都为M₃×M₄；

步骤4图像导数的自相关计算

步骤4.1利用公知的自相关运算，对步骤3中得到的主方向导数D_a(x，y)和副方向导数D_b(x，y)，分别计算主方向导数D_a(x，y)的自相关C_a(x，y)和副方向导数D_b(x，y)的自相关C_b(x，y)：

$C_a(x,y)=D_a(x,y)\⊗D_a(x,y)$

$C_b(x,y)=D_b(x,y)\⊗D_b(x,y)$

其中下标a和b分别表示主方向和副方向，

表示相关操作，此时主方向导数的自相关C_a(x，y)和副方向导数的自相关C_b(x，y)的大小都为(2×M₃-1)×(2×M₄-1)；

步骤4.2对主方向导数的自相关C_a(x，y)按照与主方向θ垂直的方向进行累加得AC_a(x，y)，对副方向导数的自相关C_b(x，y)按照与副方向θ′垂直的方向进行累加得AC_b(x，y)，累加方法如下：

${\mathrm{AC}}_a(m,n)=C_a(1,n)+C_a(2,n)+L+C_a(2\×M_3-1,n){|}_{m=\mathrm{1,1}\≤n\≤2\×M_4-1}$

${\mathrm{AC}}_b(m,n)=C_b(m,1)+C_b(m,2)+L+C_b(m,2\×M_4-1){|}_{1\≤m\≤2\×M_3-1,n=1}$

其中m和n分别为行坐标和列坐标，且都为大于0的整数，主方向导数自相关的累加AC_a(x，y)的大小为1×(2×M₄-1)，副方向导数的自相关的累加AC_b(x，y)的大小为(2×M₃-1)×1；

步骤5模糊参数识别

设主方向模糊长度为d，副方向模糊半径为r，对步骤4中得到的主方向导数自相关的累加AC_a(x，y)和副方向导数的自相关的累加AC_b(x，y)分别绘制曲线图，在主方向导数自相关的累加AC_a(x，y)曲线图中将出现一个正的最高峰k_ap，在正最高峰两侧，出现负的左侧最低峰k_ap1和负的右侧最低峰k_ap2，下标ap，ap1和ap2分别表示主方向导数自相关的累加AC_a(x，y)曲线图中正的最高峰，负的左侧最低峰和负的右侧最低峰，两负最低峰k_ap1和k_ap2之间的距离s₁的值为2d+4r，在AC_b(x，y)曲线图中将出现一个正的最高峰k_bp，在正最高峰两侧，出现负的左侧最低峰k_bp1和负的右侧最低峰k_bp2，下标bp，bp1和bp2分别表示副方向导数自相关的累加AC_b(x，y)曲线图中正的最高峰，负的左侧最低峰和负的右侧最低峰，两负最低峰k_bp1和k_bp2之间的距离s₂的值为4r，计算主方向模糊长度d和副方向模糊半径r：

$\left\{\begin{array}{c}{s}_1=2d+4r\\ s_2=4r\end{array}\right.\⇒\left\{\begin{array}{c}d=(s_1-s_2)/2\\ r=s_2/4\end{array}\right.$

步骤6建立复合模糊模型

根据步骤5得到的主方向模糊长度d和副方向模糊半径r，通过对标准运动模糊模型h_q(x，y)和标准散焦模型h_p(x，y)进行线性运算和均一化处理来建立复合模糊模型h(x，y)，下标q和p分别表示运动模糊和散焦，建立复合模糊模型h(x，y)的方法如下：

步骤6.1对标准运动模糊模型h_q(x，y)与标准散焦模型h_p(x，y)进行卷积操作，得到复合模糊模型的理论模型h_v(x，y)，下标v表示复合模糊模型的理论模型，具体方法如下：

标准散焦模型h_p(x，y)的函数形式h_pt(il，i2)为：

其中il，i2为横坐标和纵坐标，且都为实数，下标pt表示标准散焦模型的数学模型，r为副方向模糊半径且以r为散焦模糊半径，r为大于零的实数，标准散焦模型h_p(x，y)是标准散焦模型的函数形式h_pt(il，i2)的矩阵形式，标准散焦模型h_p(x，y)中的任意元素为大于或等于零的实数，

假设运动模糊方向水平，标准运动模糊模型h_q(x，y)的函数形式h_qt(il，i2)为：

其中il，i2为横坐标和纵坐标，且都为实数，下标qt表示标准运动模糊模型的数学模型，d为主方向模糊长度且以d为运动模糊长度，d为大于零的实数，标准散焦模型h_q(x，y)是标准运动模糊模型的函数模型h_qt(il，i2)的矩阵形式，标准运动模糊模型h_q(x，y)中的任意元素为大于或等于零的实数，

对标准散焦模型h_p(x，y)和标准运动模糊模型h_q(x，y)进行卷积，得到复合模糊模型的理论模型h_v(x，y)，将复合模糊模型的理论模型h_v(x，y)表示为函数形式h_vt(il，i2)：

其中下标vt表示复合模糊模型的理论模型的函数形式，

步骤6.2复合模糊模型的理论模型h_v(x，y)确定后，在能量守恒的前提下，对复合模糊模型的理论模型h_v(x，y)进行均一化处理，得到最终的复合模糊模型h(x，y)，方法如下：

步骤6.2.1

计算复合模糊模型的理论模型h_v(x，y)中大于零的点的个数w，

步骤6.2.2

将复合模糊模型的理论模型h_v(x，y)中大于零的点所对应的值替换为

得到最终的复合模糊模型h(x，y)，将复合模糊模型h(x，y)表示成函数形式h_t(il，i2)：

其中下标t表示复合模糊模型的函数形式，

步骤7图像复原

根据步骤6得到的复合模糊模型h(x，y)，利用公知的维纳滤波进行复原，得到复原图像

Images