CN105139348A - 一种运动模糊图像的双迭代混合盲复原方法 - Google Patents

Abstract

本发明提供一种运动模糊图像的双迭代混合盲复原方法，包括计算待复原的运动模糊图像的三种频谱，包括常规频谱、扩展对数频谱和去主线的扩展对数频谱；用通用迭代图像盲复原法对模糊图像进行初步复原，并用强空域约束霍夫变换和弱空域约束霍夫变换对运动模糊点扩散函数进行粗估计；用双迭代法进行多轮粗估计；用对称性与功率联合度量搜索的方法，对运动方向进行精估计；用在去主线的扩展对数频谱中寻找全局转折处的方法，对运动距离进行精估计；创建光学传输函数，并用常规图像复原方法对模糊图像进行复原。本发明对运动模糊点扩散函数的估计准确，复原效果好，并且兼有抗噪性强、稳定性好、无收敛性问题、参数设置简单等优点，可以广泛应用于各种场合下运动模糊图像的复原。

Description

一种运动模糊图像的双迭代混合盲复原方法

技术领域

本发明属于图像处理技术领域，具体涉及一种运动模糊图像的双迭代混合盲复原方法。

背景技术

图像复原是近年来的一个热点技术，已经广泛地用于天文、医学、军事、影视等领域。它的目的是从降质的观测图像中恢复出原来的清晰场景，其数学模型可以描述为清晰图像与点扩散函数(pointspreadfunction，PSF)的卷积再与噪声相加，如下式：

$g(x,y)=f(x,y)\⊗h(x,y)+n(x,y)---\left(1\right)$

其中，x,y是图像的空间坐标，g(x,y)是模糊图像的灰度，f(x,y)清晰图像的灰度，h(x,y)是PSF的灰度，n(x,y)是加性噪声，代表卷积运算。通过傅立叶变换，该模型可在频域表示为：

G(u,v)＝F(u,v)·H(u,v)+N(u,v)(2)

其中，u,v表示离散频率坐标，G(u,v)、F(u,v)、H(u,v)和N(u,v)分别是g(x,y)、f(x,y)、h(x,y)和n(x,y)的离散傅立叶变换，H(u,v)又称为光学传输函数(opticaltransferfunction，OTF)。图像复原也就是从g(x,y)中把f(x,y)恢复出来。然而在实际情况中，h(x,y)和n(x,y)往往是未知的，以致恢复f(x,y)非常困难，这种情况被称为“盲复原”。

运动模糊是一种常见的模糊类型，无论是在日常的摄影、摄像还是在交通、监控等领域中，都极其普遍地存在。当目标运动较快或镜头本身发生抖动时，在相机的曝光时间内目标与相机之间容易发生相对位移，这时就会出现一种“拖尾”现象，也就是运动模糊。根据运动性质，运动模糊可分为很多种类，如匀速运动、变速运动、直线运动、曲线运动、旋转运动等等。但一般匀速直线运动最有代表性，其他类型的运动模糊可以被分解成分段的匀速直线运动，因此本发明主要考虑匀速直线运动模糊。

设模糊图像的行数和列数分别为M和N，运动方向与x轴正方向的夹角为运动距离为则运动模糊的OTF模型可表示为：

$H(u,v)=\sin c\[\hat{L}(ucos\widehat{\mathrm{\θ}}/M+vsin\widehat{\mathrm{\θ}}/N)\]\·e^{-j\mathrm{\π}\hat{L}(ucos\widehat{\mathrm{\θ}}/M+vsin\widehat{\mathrm{\θ}}/N)}---\left(3\right)$

式中，j表示虚数单位。可见，只要知道了运动的方向和距离，就可以唯一确定运动模糊的OTF或PSF。

由于运动模糊的OTF是一个sinc函数，故有许多特殊的性质。若考虑|H(u,v)|，当任意指定一个离散频率坐标u₀时有：

|H(u₀,v)|＝sinc[L(u₀cosθ/M+vsinθ/N)](4)

从中可知，在|H(u₀,v)|中每隔N/(Lsinθ)将会出现一个零点。同理，|H(u,v₀)|每隔M/(Lcosθ)将会出现一个零点。所有的零点会组成一系列等间隔的平行直线，在|H(u,v)|中则显现为一系列平行的暗条纹，这些暗条纹与u轴的夹角φ为：

$\widehat{\mathrm{\φ}}=-arctan\left(\frac{N}{Mtan\widehat{\mathrm{\θ}}}\right)---\left(5\right)$

在噪声较小的运动模糊图像中，常常能在频谱中观察到这些平行的暗条纹。许多盲复原算法就是通过辨识暗条纹的方向和间隔来估计运动参数。然而，由于实际情况复杂，噪声普遍较大，这些暗条纹常常被淹没在噪声之中，很难直接观测，导致参数难以准确估计，复原效果不佳。

目前，运动模糊图像的盲复原方法大致分为两类。第一类是通用的迭代图像盲复原法(简称“迭代法”)，比如最大似然估计算法、迭代盲反卷积算法、约束递归逆滤波算法、R-L算法等。这类算法并非专门针对运动模糊图像，或者说运动模糊是这类算法的一个特例。这类算法的优点是无需PSF或OTF模型，抗噪性和稳定性较好；缺点是有收敛性问题，复原结果不唯一，不能确保有较好的复原效果，并且计算量大。大量实验证明，这类算法对运动模糊图像的复原效果并不令人满意，原因在于未使用运动模糊的先验信息，是一个笼统的复原结果，并且可能未达到最佳收敛。虽然有专门针对运动模糊的改进方法，但是仍然没有解决迭代过程中不可避免的收敛性问题。

第二类是先辨识运动模糊参数再用经典复原法进行复原的方法(简称“参数法”)。这类方法的优点在于更具针对性，在参数准确的条件下可以得到较好的复原效果，无收敛性问题，并且计算量小；而缺点在于很难辨识出准确的参数，稳定性和适用性较差。迄今为止，人们提出了多种方法来估计运动的方向和运动距离这两个参数，但适用范围总是有限。比如频谱分析法，通过计算频谱中暗条纹的方向和间隔来对运动的方向和距离进行估计。但在信噪比较低的真实图像中，暗条纹往往被噪声完全淹没，从而很难准确找出暗条纹的位置。再如倒谱分析法，利用图像的倒谱来辨识参数。但该方法抗噪性能很低，极少的噪声就能使倒谱失去运动模糊的特征。又如空间分析法，利用图像本身的轨迹来计算参数。该方法计算复杂，具有一定抗噪性，但对低信噪比的图像或细节不丰富的图像仍然束手无策。另外，还有的方法通过对复原后的图像进行评价的方式来搜索参数，但图像的客观评价与主观认知还有很大差距，从而导致该类方法只能适用于某些场合的运动模糊图像。

综上，虽然运动模糊图像的盲复原方法层出不穷，但大多方法的抗噪性、稳定性和适用性不能兼得，无法满足实际应用的需求，故亟待一种更加有效、更加可靠、更加实用的新方法。

发明内容

本发明要解决的技术问题为：克服现有运动模糊图像盲复原方法的抗噪性、稳定性和适用性不足的问题，提供一种运动模糊图像的双迭代混合盲复原方法。该方法将迭代法和参数法进行结合，将二者取长补短，既利用了迭代法的抗噪性强和参数法复原效果好的优点，又避开了迭代法的收敛性问题和参数法稳定性差的缺点，从而可以广泛应用于各种场合下运动模糊图像的复原。

本发明的技术方案：一种运动模糊图像的双迭代混合盲复原方法，其实现步骤如下：

步骤一、计算待复原的运动模糊图像(以下简称“模糊图像”)的三种频谱，包括常规频谱、扩展对数频谱和去主线的扩展对数频谱。具体实现过程为：

(1)用常规的快速傅里叶变换计算常规频谱。

(2)将模糊图像转换为灰度图像，再进行补0扩展，然后进行快速傅里叶变换，最后进行中心化和对数化得到扩展对数频谱。

(3)计算扩展对数频谱高频区域的竖直和水平投影曲线，再在该曲线中用适当的阈值找出频谱主线的对应点，最后将竖直和水平的主线合并成需要去除的区域，从而得到去主线的扩展对数频谱。

步骤二、用通用迭代图像盲复原法对模糊图像进行初步复原，得到全局点扩散函数(PSF)，并用强空域约束霍夫变换和弱空域约束霍夫变换对运动模糊PSF进行粗估计。具体实现过程为：

(1)用最大似然估计算法对模糊图像进行初步复原，得到全局PSF。

(2)用最大类间方差法对全局PSF进行二值化。

(3)以全局PSF的最亮点为强空域约束点，计算全局PSF的强空域约束霍夫变换。

(4)以全局PSF中按灰度排序靠前的若干点为弱空域约束点，计算全局PSF的弱空域约束霍夫变换。

(5)用上述强空域约束霍夫变换和弱空域约束霍夫变换对全局PSF提取直线，作为运动模糊PSF的粗估计。

步骤三、用双迭代法进行多轮粗估计，对运动模糊PSF的支持域和粗估计结果进行优化。具体实现过程为：

(1)用符合模糊图像半宽的正方形作为步骤二中全局PSF的初始支持域。

(2)用步骤二粗估计得到的直线尺寸为新的支持域，再次重复步骤二，如此循环。

(3)当两轮粗估计的直线长度差小于指定阈值时终止循环。

步骤四、用对称性与功率联合度量搜索的方法，对运动方向进行精估计。具体实现过程为：

(1)以步骤三最后得到直线的角度附近范围作为待搜索范围。

(2)对于待搜索范围内的任一角度，计算其经过强空域约束点的直线两侧附近的全局PSF对称性。

(3)对于待搜索范围内的任一角度，计算其经过强空域约束点的直线的功率。

(4)将上述对称性和功率进行归一化并加权和，得到待搜索角度的评价函数，在待搜索范围中使该评价函数最大的角度就是运动方向的精估计。

步骤五、用在去主线的扩展对数频谱中寻找全局转折处的方法，对运动距离进行精估计。具体实现过程为：

(1)用运动方向的精估计值计算暗条纹的角度估计值。

(2)根据暗条纹的角度估计值，将去主线的扩展对数频谱顺时针旋转，使可能存在的暗条纹呈竖直方向，并沿竖直方向求均值，得到其竖直均值曲线。

(3)计算上述竖直均值曲线的全局转折点，得到第一条暗条纹位置。

(4)利用第一条暗条纹位置与运动距离的关系，计算运动距离的精估计值。

步骤六、根据运动方向和运动距离的精估计结果，创建光学传输函数，并用常规图像复原方法对模糊图像进行复原。

本发明与现有技术相比的有益效果在于：

(1)本发明对运动模糊PSF的估计准确，复原效果好。

(2)抗噪性强，本发明对信噪比极低的运动模糊图像也能较好的估计其PSF。

(3)稳定性好，本发明估计运动模糊PSF时可以克服各种干扰。

(4)无收敛性问题，本发明迭代过程的收敛条件很容易达成，能确保收敛。

(5)参数设置简单，本发明的参数设置对结果影响很小，故无需频繁尝试。

(6)适用范围广，本发明可以广泛应用于各种场合下运动模糊图像的复原。

发明内容

图1是实施例中待复原的运动模糊图像；

图2是实施例步骤一中扩展对数频谱(左)和去主线的扩展对数频谱(右)示意图；

图3是实施例步骤二中初步得到的全局PSF示意图；

图4是实施例步骤三中每轮的全局PSF以及提取的直线示意图；

图5是实施例步骤四中待搜索角度的对称性、功率和评价函数曲线示意图；

图6是实施例步骤五中去主线的扩展对数频谱的定向均值曲线示意图；

图7是实施例步骤六中对模糊图像的复原结果。

具体实施方式

以下是本发明的具体实施办法。但以下的实施例仅限于解释本发明，本发明的保护范围应包括权利要求的全部内容，而且通过以下实施例对该领域的技术人员即可以实现本发明权利要求的全部内容。

本发明包含六个步骤，其具体实现方法如下：

步骤一、计算待复原的运动模糊图像的三种频谱，包括常规频谱、扩展对数频谱和去主线的扩展对数频谱。

本实施例中待复原的运动模糊图像如图1所示，可见无法辨认其中的文字，并且信噪比极低。在本发明的实施过程中，一共将用到三种频谱，详细如下：

一是常规频谱G(u,v)，用于最后复原图像。该频谱用常规的快速傅里叶变换(FFT)即可算得，此处不再赘述。如果模糊图像是彩色图像，则要对每个色彩面分别进行变换。

二是扩展对数频谱，用G_e(u,v)表示，其计算公式为：

G_e(u,v)＝ln|FFT[(-1)^x+yg_e(x,y)]|(6)

其中，FFT[·]表示快速傅里叶变换，(-1)^x+y因子的作用是将频谱原点移到中心(中心化)，ln表示取自然对数(对数化)，g_e(x,y)是模糊图像经过如下转换而来的扩展图像：

1)如果模糊图像是彩色图像，则转换为灰度图像；

2)如果模糊图像尺寸较小(在本发明中设定行或列以128为界限)，为了获得较好的频谱分辨率，则对转换后的灰度图像进行补0扩展，使行数和列数都不小于该界限。后面M和N表示扩展图像的行数和列数。

由于本实施例的模糊图像是行列都大于128的灰度图像，故无需进行转换，只需中心化和对数化即可，其结果如图2中左图所示。

三是去主线的扩展对数频谱，用表示。对于自然图像，通常在G_e(u,v)中原点附近的横纵方向上有若干条贯穿频谱、亮度极高的直线，即频谱的主线，会对参数的估计产生不良影响，必须予以去除。频谱主线的特点是在高频区能量极大，故可以利用该特性来找出其位置。先考虑竖直方向的主线，选取如下的竖直高频区域Ω_vp：

Ω_vp＝{(u,v)|u∈[-M/2,-M/4]∪[M/4,M/2],v∈[-N/4,N/4]}(7)

再在该区域中沿竖直方向求和，得到G_e(u,v)的竖直投影曲线C_vp(v)：

$C_{vp}\left(v\right)=\underset{u}{\Σ}{G}_e(u,v),(u,v)\∈{\mathrm{\Ω}}_{vp}---\left(8\right)$

在该曲线中可以看到原点附近有若干明显凸出的点，每个点也就对应一条竖直主线。现在需要用一个适当的阈值(用A_vp表示)把这些点区分开来，本发明提供一种自动阈值法，其计算公式为：

$A_{vp}=\underset{\left|v\right|\≤{\mathrm{\ω}}_0}{max}\[{\tilde{C}}_{vp}\left(v\right)\]+E\[|{\tilde{C}}_{vp}\left(v\right)-C_{vp}\left(v\right)|\]---\left(9\right)$

式中，ω₀表示主线偏离原点的最大距离，max(·)表示取最大值，表示对C_vp(v)进行低通滤波得到的平滑曲线，E[·]表示取均值。本发明选用11点移动平均法对C_vp(v)进行低通滤波，并取ω₀为5。从而得到需去除的竖直主线区域为：

同理，对于水平方向的主线，首先选取如下的水平高频区域Ω_up：

Ω_up＝{(u,v)|v∈[-N/2,-N/4]∪[N/4,N/2],u∈[-M/4,M/4]}(11)

再在该区域中沿水平方向求和，得到频谱的水平投影曲线C_up(u)：

$C_{up}\left(u\right)=\underset{v}{\Σ}{G}_e(u,v),(u,v)\∈{\mathrm{\Ω}}_{up}---\left(12\right)$

然后用相同的方法将C_up(u)进行低通滤波得到平滑曲线以及相应阈值：

$A_{up}=\underset{\left|u\right|\≤{\mathrm{\ω}}_0}{\max }\[{\tilde{C}}_{up}\left(u\right)\]+E\[|{\tilde{C}}_{up}\left(u\right)-C_{up}\left(u\right)|\]---\left(13\right)$

从而得到需去除的水平主线区域为：

频谱的主线对参数识别的影响主要体现在高频区。而对于低频区来说，去除主线可能反而不利于参数的识别，因此要有所取舍。首先原点不能去除，其次水平和竖直主线的交集部分可予以保留，从而整个频谱需要去除的区域Ω°为：

于是可以表示为：

本实施例中去主线的扩展对数频谱如图2中右图所示。

步骤二、用通用迭代图像盲复原法对模糊图像进行初步复原，得到全局PSF，并用强空域约束霍夫变换和弱空域约束霍夫变换对运动模糊PSF进行粗估计。

一般的通用迭代图像盲复原法(比如最大似然估计算法、迭代盲反卷积算法、约束递归逆滤波算法、R-L算法等)，可以估计出模糊图像的PSF，但该PSF是包含各种模糊因素在内且有噪声影响的“全局PSF”，因此复原效果往往不够理想。对于含有运动模糊的图像，全局PSF中必会出现在某一方向上明显较亮的特征。如果能从其中提取运动模糊的信息，则可以更加有针对性地去除运动模糊。所以该步骤的目的就是从全局PSF中提取出基本符合运动模糊PSF的直线，即对运动模糊PSF进行粗估计。

在该步骤中，首先需要选用一种通用迭代图像盲复原法对模糊图像进行初步复原，得到全局PSF。在该实施例中选用的是“最大似然估计算法”，且固定迭代10次。详细的算法公式可参考相关文献，此处不再赘述。由于只是为了得到全局PSF，故不需要收敛，从而避免了迭代法的收敛性问题。下面用表示上述迭代法得到的全局PSF。本实施例初步得到的全局PSF如图3所示，可见全局模糊方向掩盖了运动模糊方向。

然后，需要对全局PSF进行二值化。这里采用一般的阙值分割即可，本实施例采用的是最大类间方差法。下面用B_i(x,y)表示将分割后的二值图像。

再后，就是计算全局PSF的强空域约束霍夫变换和弱空域约束霍夫变换。霍夫变换是一种常用的直线提取方法。因为直线可以表示为如下的极坐标形式：

ρ＝xcosα+ysinα(17)

式中，ρ是原点到直线的距离，α是直线的法线与x轴正方向的夹角，所以霍夫变换的原理就是把(x,y)空间的所有点变换到(ρ,α)空间，(x,y)空间的共线点则会在(ρ,α)空间中形成一个峰值点，检测出该峰值点就可以反变换为直线。

然而，由于有各种模糊因素和噪声的干扰，故直接用霍夫变换来提取全局PSF中的直线往往是非常困难的，比如二值化误差、车灯效应、全局模糊方向掩盖运动模糊方向等等。为了克服这些干扰，本发明提出了“强空域约束霍夫变换”和“弱空域约束霍夫变换”来实现运动模糊PSF的粗估计。

强空域约束霍夫变换，就是使霍夫变换提取出的直线必须经过空间中的某些指定点(下称“强空域约束点”)的方法。设的支持域为Ω，其常规霍夫变换为强空域约束点的集合为Ω_s，点数为n_s。首先建立一幅与B_i(x,y)尺寸相同的强空域约束二值图像B_s(x,y)：

$B_s(x,y)=\left\{\begin{array}{cc}1,& (x,y)\∈{\mathrm{\Ω}}_s\\ 0,& (x,y)\∈{\stackrel{\‾}{\mathrm{\Ω}}}_s\∩\mathrm{\Ω}\end{array}\right.---\left(18\right)$

然后，对B_s(x,y)进行霍夫变换，即得到强空域约束霍夫变换并定义强空域约束有效域∑_s为：

那么所有位于∑_s中的就对应所有经过Ω_s的直线。因为要提取的直线必然要经过全局PSF中的最亮点，所以可用中的最大值点作为强空域约束点，即：

${\mathrm{\Ω}}_s=\left\{(x_s,y_s)\right\}=\left\{(x,y)\right|{\hat{h}}_i(x,y)=max\left({\hat{h}}_i\right)\}---\left(20\right)$

式中，(x_s,y_s)表示强空域约束点的坐标。

弱空域约束霍夫变换，就是使霍夫变换提取出的直线尽量经过空间中的某些指定点(下称“弱空域约束点”)的方法。设弱空域约束点的集合为Ω_w，点数为n_w，也首先建立一幅与B_i(x,y)尺寸相同的弱空域约束二值图像B_w(x,y)：

$B_w(x,y)=\left\{\begin{array}{cc}1,& (x,y)\∈{\mathrm{\Ω}}_w\\ 0,& (x,y)\∈{\stackrel{\‾}{\mathrm{\Ω}}}_w\∩\mathrm{\Ω}\end{array}\right.---\left(21\right)$

然后，对B_w(x,y)进行霍夫变换，即得到弱空域约束霍夫变换并定义弱空域约束有效域∑_w为：

式中，ξ是一个可调节的松弛比例(本发明取0.9)。那么所有位于∑_w中的就对应尽量经过Ω_w的直线。其中Ω_w的确定方法为：将全局PSF其按灰度大到小排序，选择其中前n_w个点作为Ω_w，而n_w可取为：

$n_w=max(M_{{\hat{h}}_i},N_{{\hat{h}}_i})---\left(23\right)$

式中，和分别表示的行数和列数。

最后，只需在∑_s∩∑_w中检测的峰值，再反变换为直线，并取其中角度为中间值的直线，就得到了运动模糊PSF的粗估计。本发明的这种方法，基本可以屏蔽掉诸如二值化误差、车灯效应等大部分干扰，具有良好的抗噪性和稳定性。

步骤三、用双迭代法进行多轮粗估计，对运动模糊PSF的支持域和粗估计结果进行优化。

在通用迭代图像盲复原法中，一般都需要事先知道PSF的支持域。如果支持域的选择出现了较大的偏差，就会严重影响PSF的正确性和复原效果，因此有一定局限性。本发明则提供了一种双迭代法来打破这种局限性。本发明的双迭代法，就是先用符合模糊图像半宽的正方形(即边长为min(M,N)/2)作为步骤二中全局PSF的初始支持域，再用步骤二粗估计得到的直线尺寸为新的支持域，再次重复步骤二，如此循环，直到两轮粗估计的直线长度差小于指定阈值(可取1～5)时为止。由于该终止条件很容易达成，通常几轮以内就能收敛，故不存在收敛性问题。收敛后就会自动得到最优的支持域和粗估计结果。因为每轮粗估计中都要使用一次通用迭代图像盲复原法，所以称这种方法为“双迭代法”。该步骤可以进一步克服车灯效应、全局模糊方向掩盖掉运动模糊方向等问题，使粗估计结果更加稳定可靠。

本实施例仅进行了四轮粗估计就达成收敛，每轮的全局PSF以及提取的直线如图4所示，其中从左到右分别对应第一至第四轮的结果，全局PSF中的黑线即表示提取的直线。可见这种双迭代法很好地排除了全局模糊的干扰，辨识出了真正的运动模糊PSF。

步骤四、用对称性与功率联合度量搜索的方法，对运动方向进行精估计。

因为霍夫变换存在误差，所以以上步骤得到的直线还只是运动模糊PSF的粗估计，尚需利用全局PSF的全部信息进行更加精确的估计。准确的运动方向在全局PSF中表现为这样两个特征：一是该方向直线两侧附近的全局PSF非常对称，二是该方向直线的功率极大。于是，可以用这两个标准在步骤三最后得到直线的角度附近范围进行搜索。

设步骤三最后得到直线的角度为θ_i(这里θ_i表示直线与x轴正方向的夹角且θ_i∈[0,π))，霍夫变换中角度的量化分辨率为θ_r，待搜索的最大角度偏差为θ_d，角度自变量为θ，那么待搜索范围Θ则为：

Θ＝{θ|θ＝θ_i±nθ_r,|θ-θ_i|≤θ_d,n＝0,1,2,…}(24)

如上所述，搜索标准定义为：

J_a(θ)＝S(θ)+λ_a·P(θ)(25)

其中，J_a(θ)表示待搜索角度的评价函数，S(θ)表示经过强空域约束点且角度为θ的直线两侧附近的全局PSF对称性度量函数，P(θ)表示该直线的功率函数，λ_a是两种度量之间的加权调节参量。

对于S(θ)的计算，首先要确定直线“两侧附近”的范围。本发明用B_i(x,y)的1值区域作为直线“两侧附近”的基本限定，再从其中剔除对称点超出支持域、对称点在0值区域和对称点是其自身的无效点，也就是定义对称性的度量范围为：

${\mathrm{\Ω}}_{\mathrm{\θ}}^{+}=\left\{(x,y)\right|(x,y)\∈\mathrm{\Ω},(x_{\mathrm{\θ}}^{+},y_{\mathrm{\θ}}^{+})\∈\mathrm{\Ω},B_i(x,y)=1,B_i(x_{\mathrm{\θ}}^{+},y_{\mathrm{\θ}}^{+})=1,(x,y)\≠(x_{\mathrm{\θ}}^{+},y_{\mathrm{\θ}}^{+})\}---\left(26\right)$

式中，表示(x,y)关于经过(x_s,y_s)且角度为θ的直线的对称点，其计算公式为：

$x_{\mathrm{\θ}}^{+}=d_s\·cos(2\mathrm{\θ}-{\mathrm{\γ}}_s)+x_s---\left(27\right)$

$y_{\mathrm{\θ}}^{+}=d_s\·sin(2\mathrm{\θ}-{\mathrm{\γ}}_s)+y_s---\left(28\right)$

式中，d_s表示(x,y)到(x_s,y_s)的距离，γ_s表示(x,y)到(x_s,y_s)的连线与x轴正方向的夹角，且有：

$d_s=\sqrt{{(x-x_s)}^2+{(y-y_s)}^2}---\left(29\right)$

${\mathrm{\γ}}_s=arctan\left(\frac{y-y_s}{x-x_s}\right)+\mathrm{sgn}\left(\mathrm{sgn}\right(x-x_s)+1)\·\mathrm{sgn}(y-y_s-\mathrm{\ϵ})\·\mathrm{\π}---\left(30\right)$

式中，sgn(·)表示取符号，ε是一个极小的正数(0<ε＜＜θ_r)，γ_s∈(-π,π]。

于是，在中位于该直线某一侧的有效点的集合为：

${\mathrm{\&Omega;}}_{\mathrm{\&theta;}}^{++}=\left\{(x,y)\right|(x,y)\&Element;{\mathrm{\&Omega;}}_{\mathrm{\&theta;}}^{+},\mathrm{\&theta;}-\mathrm{\&pi;}<{\mathrm{\&gamma;}}_s<\mathrm{\&theta;}\}---\left(31\right)$

利用对称点之间的相关性，就可以进行对称性度量，其公式为：

$S\left(\mathrm{\&theta;}\right)=\frac{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_{\mathrm{\&theta;}}^{++}}{\&Sigma;}{\hat{h}}_i(x,y){\hat{h}}_i(x_{\mathrm{\&theta;}}^{+},y_{\mathrm{\&theta;}}^{+})}{\sqrt{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_{\mathrm{\&theta;}}^{++}}{\&Sigma;}{{\hat{h}}_{{}^i}}^2(x,y)}\sqrt{\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_{\mathrm{\&theta;}}^{++}}{\&Sigma;}{{\hat{h}}_i}^2(x_{\mathrm{\&theta;}}^{+},y_{\mathrm{\&theta;}}^{+})}}---\left(32\right)$

一般来说，使S(θ)最大的θ就是准确的运动方向。但因为量化误差的原因，有可能使S(θ)最大的θ不唯一，所以还需要加入另一个度量标准——直线的功率P(θ)。P(θ)就是中经过强空域约束点且角度为θ的直线的单位能量。如果令这些直线长度相等且点数相同，那么P(θ)就可以简化为灰度能量和。设表征这些直线的点集为Ω_θ，Ω_θ中的所有点到(x_s,y_s)的最大距离为ρ_M，则有：

Ω_θ＝{(x,y)|x＝ρcosθ+x_s,y＝ρsinθ+y_s,|ρ|＝0,1,2,…,ρ_M}(33)

其中ρ_M可取为：

${\mathrm{\&rho;}}_M=min(y_s-1,x_s-1,M_{{\hat{h}}_i}-y_s,N_{{\hat{h}}_i}-x_s)---\left(34\right)$

式中，min(·)表示取最小值。于是P(θ)的计算式为：

$P\left(\mathrm{\&theta;}\right)=\underset{(x,y)\&Element;{\mathrm{\&Omega;}}_{\mathrm{\&theta;}}}{\&Sigma;}\left|{\hat{h}}_i(x,y)\right|---\left(35\right)$

最后，将上述对称性S(θ)和功率P(θ)进行归一化并加权和，即可得到待搜索角度的评价函数J_a(θ)。为了便于选取合适的加权调节参量λ_a，还要先将S(θ)和P(θ)都分别归一到[0,1]中。从而，在Θ中使J_a(θ)最大的θ，就是运动方向的精估计，即：

$\widehat{\mathrm{\&theta;}}=arg\underset{\mathrm{\&theta;}\&Element;{\mathrm{\&Theta;}}_0}{max}\&lsqb;J_a\left(\mathrm{\&theta;}\right)\&rsqb;---\left(36\right)$

式中，argmax表示取指定区域中使表达式最大的变量值。

在本实施例中，各参数取值为θ_r＝0.5、θ_d＝10及λ_a＝1，待搜索角度的对称性、功率和评价函数曲线如图5所示，最后计算出的为-49.42°。

步骤五、用在去主线的扩展对数频谱中寻找全局转折处的方法，对运动距离进行精估计。

如前所述，运动模糊会使模糊图像的频谱存在一系列平行的暗条纹。但是当模糊图像的频谱被噪声强烈干扰时，暗条纹将部分或全部淹没在噪声中，很难进行分辨和计算。然而，就算暗条纹不可见，第一条暗条纹也会使形成一个明显的全局转折处。因此，利用该特性以及已经得到的运动方向，就可以对运动距离进行精估计。

首先，需要用运动方向的精估计值计算暗条纹的角度估计值，即将代入公式(5)，计算出然后，将去主线的扩展对数频谱顺时针旋转度，使可能存在的暗条纹呈竖直方向。令旋转后的频谱为再沿竖直的u′方向对求均值，得到的定向均值曲线C_d(v′)，即：

式中，E_u′[·]表示沿竖直的u′方向求均值。于是，C_d(v′)将会在第一条暗条纹处呈现一个全局转折点。设该全局转折点横坐标为v′₀，C_d(v′)的定义域为v′∈D_T＝[0,v′_T]，其中v′_T表示v′的最大取值，则v′₀的计算公式为：

$v_0^{\&prime;}=\arg \underset{v^{\&prime;}\&Element;D_T}{max}|\&lsqb;C_d\left(0\right)-C_d\left(v_T^{\&prime;}\right)\&rsqb;v^{\&prime;}+v_T^{\&prime;}{C}_d\left(v^{\&prime;}\right)-v_T^{\&prime;}{C}_d\left(0\right)|---\left(38\right)$

于是|v′₀|就是第一条暗条纹与频谱原点的距离。从而，利用它与运动距离的关系，即可得到运动距离的精估计值为：

$\hat{L}=\left|\frac{Mcos\widehat{\mathrm{\&phi;}}}{v_0^{\&prime;}cos\widehat{\mathrm{\&theta;}}}\right|---\left(39\right)$

本实施例的C_d(v′)曲线如图6所示，其中竖直的虚线表示v′＝v′₀所在位置。最后计算出的为24.9222。

步骤六、根据运动方向和运动距离的精估计结果，创建光学传输函数，并用常规图像复原方法对模糊图像进行复原。

在得到了运动方向和距离的精估计值之后，将其代入公式(3)即可得到运动模糊的OTF。从而，用常规图像复原方法(比如维纳滤波、约束最小平方滤波、连续边界慢演进滤波等)就可以对模糊图像进行复原。本实施例中选用基于中频的维纳滤波法进行复原，具体的公式可参考相关文献，此处不再赘述。需要说明的是，如果模糊图像是彩色图像，则要对每个色彩面分别复原，再合并成新的彩色图像。复原之后，还可以根据具体场合适当调节对比度，使复原图像更易观看。本实施例的复原结果如图7所示，可见其中的数字已经可以辨认。

本发明未详细阐述部分属于本领域技术人员的公知技术。

Claims (6)

1.一种运动模糊图像的双迭代混合盲复原方法，其特征在于实现步骤如下：

步骤一、计算待复原的运动模糊图像(以下简称“模糊图像”)的三种频谱，包括常规频谱、扩展对数频谱和去主线的扩展对数频谱；

步骤二、用通用迭代图像盲复原法对模糊图像进行初步复原，得到全局点扩散函数(PSF)，并用强空域约束霍夫变换和弱空域约束霍夫变换对运动模糊PSF进行粗估计；

步骤三、用双迭代法进行多轮粗估计，对运动模糊PSF的支持域和粗估计结果进行优化；

步骤四、用对称性与功率联合度量搜索的方法，对运动方向进行精估计；

步骤五、用在去主线的扩展对数频谱中寻找全局转折处的方法，对运动距离进行精估计；

步骤六、根据运动方向和运动距离的精估计结果，创建光学传输函数，并用常规图像复原方法对模糊图像进行复原。

2.根据权利要求1所述的一种运动模糊图像的双迭代混合盲复原方法，其特征在于：所述步骤一的实现过程为：

(1)用常规的快速傅里叶变换计算常规频谱；

(2)将模糊图像转换为灰度图像，再进行补0扩展，然后进行快速傅里叶变换，最后进行中心化和对数化得到扩展对数频谱；

(3)计算扩展对数频谱高频区域的竖直和水平投影曲线，再在该曲线中用适当的阈值找出频谱主线的对应点，最后将竖直和水平的主线合并成需要去除的区域，从而得到去主线的扩展对数频谱。

3.根据权利要求1所述的一种运动模糊图像的双迭代混合盲复原方法，其特征在于：所述步骤二的实现过程为：

(1)用最大似然估计算法对模糊图像进行初步复原，得到全局PSF；

(2)用最大类间方差法对全局PSF进行二值化；

(3)以全局PSF的最亮点为强空域约束点，计算全局PSF的强空域约束霍夫变换；

(4)以全局PSF中按灰度排序靠前的若干点为弱空域约束点，计算全局PSF的弱空域约束霍夫变换；

(5)用上述强空域约束霍夫变换和弱空域约束霍夫变换对全局PSF提取直线，作为运动模糊PSF的粗估计。

4.根据权利要求1所述的一种运动模糊图像的双迭代混合盲复原方法，其特征在于：所述步骤三的实现过程为：

(1)用符合模糊图像半宽的正方形作为步骤二中全局PSF的初始支持域；

(2)用步骤二粗估计得到的直线尺寸为新的支持域，再次重复步骤二，如此循环；

(3)当两轮粗估计的直线长度差小于指定阈值时终止循环。

5.根据权利要求1所述的一种运动模糊图像的双迭代混合盲复原方法，其特征在于：所述步骤四的实现过程为：

(1)以步骤三最后得到直线的角度附近范围作为待搜索范围；

(2)对于待搜索范围内的任一角度，计算其经过强空域约束点的直线两侧附近的全局PSF对称性；

(3)对于待搜索范围内的任一角度，计算其经过强空域约束点的直线的功率；

(4)将上述对称性和功率进行归一化并加权和，得到待搜索角度的评价函数，在待搜索范围中使该评价函数最大的角度就是运动方向的精估计。

6.根据权利要求1所述的一种运动模糊图像的双迭代混合盲复原方法，其特征在于：所述步骤五的实现过程为：

(1)用运动方向的精估计值计算暗条纹的角度估计值；

(2)根据暗条纹的角度估计值，将去主线的扩展对数频谱顺时针旋转，使可能存在的暗条纹呈竖直方向，并沿竖直方向求均值，得到其竖直均值曲线；

(3)计算上述竖直均值曲线的全局转折点，得到第一条暗条纹位置；

(4)利用第一条暗条纹位置与运动距离的关系，计算运动距离的精估计值。