一种基于中频的离焦模糊图像盲复原方法

摘要

本发明提供一种基于中频的离焦模糊图像盲复原方法，该盲复原方法步骤为：计算离焦模糊图像的频谱；计算离焦模糊图像频谱的椭圆轨迹均值函数；对椭圆轨迹均值函数进行“变跨度平滑滤波”得到一个基准函数；用椭圆轨迹均值函数计算离焦模糊图像频谱的中频域；在中频域中用椭圆轨迹均值函数及其基准函数估计第一暗环的位置；用估计的第一暗环位置计算离焦半径并生成离焦光学传输函数；用基于中频的维纳滤波器复原图像。本发明与现有技术相比的有益效果在于：能够从离焦模糊图像中自动辨识出离焦半径等参数，实现离焦模糊图像的有效复原，并且具有抑噪能力强、交互参数少、复原速度快的优点。

说明书

技术领域

本发明属于图像盲复原技术领域，具体涉及一种基于中频的离焦模糊图像盲复原方法。

背景技术

图像复原是近年来的一个热点技术，已经广泛地用于天文、军事、医学、遥感、电视等 领域。它的目的是从降质的观测图像中恢复出原来的清晰场景。数学上降质模型可以描述为 原始图像与点扩展函数(point spread function，PSF)在有加性噪声的情况下的卷积，如下式：

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

其中x,y是平面的二维坐标，g(x,y)是降质图像的灰度，f(x,y)原始清晰图像的灰度，h(x,y)是 降质PSF的灰度，z(x,y)是加性噪声，代表卷积运算。通过傅立叶变换，该模型可在频域 表示为：

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

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

在各种成像系统中，离焦模糊是一种常见的降质因素。只要成像系统的镜头尚未对所拍 对象聚焦，就会产生离焦模糊。一般来说，在镜头聚焦不良的情况下，一个点源的影像可以 近似为一个圆盘，系统的点扩展函数可以用数学模型表示为：

其中，r为离焦半径。对该点扩展函数进行离散傅里叶变换，则得到：

$H(u,v)=K\frac{J_1\left(2\mathrm{\pi r\rho }\right)}{2\mathrm{\pi r\rho }}---\left(4\right)$

其中，J₁(·)为一阶第一类Bessel函数，K为归一化常数，表示距离量， M为图像的行数，N为图像的列数，且u＝0,1,…,M-1和v＝0,1,…,N-1。根据一阶第一类 Bessel函数的性质，J₁(·)的零点在以原点为圆心的同心圆周上，因而理论上离焦模糊图像的 频谱存在许多同心暗环，其中第一暗环的轨迹方程为：

2πrρ₁＝3.85       (5)

式中，ρ₁表示第一暗环的距离量。在噪声很小时，离焦模糊图像的频谱中往往可以直接看到 这些同心暗环。但是当噪声较大时，暗环将会被噪声淹没，退化成一个圆斑。

目前离焦模糊图像盲复原方法大致有三类，一是通过在变换域寻找零点位置确定参数， 二是通过直线或直线边缘区域的模糊图像在空间域直接计算参数，三是根据复原结果的评价 (如峰度)来对参数进行搜索。例如，刘克等人首先利用梯度图像估计离焦半径的初值，然 后利用最小二乘图像复原技术和斐波拉契最优搜索方法确定精确的离焦半径；L.Dalong等 人提出一种基于最小峰度的模糊辨识技术，根据最小峰度的复原图像来确定模糊参数；P. Marziliano等人基于模糊边缘的平滑效应，提出一种在空间域中衡量边缘扩散从而辨识模糊 参数的方法；E.Ong等人改进为计算与边缘垂直的梯度方向及反方向上的边缘扩散；C. Yunchung通过获取边缘梯度幅度信息，以其标准差来表示模糊边缘宽度，结合加权的边缘 梯度幅值，使其结果更加稳定，等等。另外还有一类并非专门针对离焦模糊的一般图像盲复 原方法，比如迭代盲反卷积、非负性支持域约束递归逆滤波法、RL算法、极大似然估计算 法、EM算法、模拟退火法、ARMA模型参数估计法、各向异性扩散方程法等等，也可适用 于离焦模糊图像的复原，但效果和速度上不如专门针对离焦模糊图像的复原方法优越。

这些现有方法虽然能够在某些场合下取得一定的成效，但存在一些制约其广泛应用的通 病。首先，这些方法大多属于模糊检测，所得结果不够准确，尤其在强噪情况下可靠性较差， 复原效果不佳。其次，这些方法常涉及到迭代运算，计算量很大，导致复原效率很低，复原 一幅图像往往要耗时几秒至几十分钟不等。另外，这些方法几乎都需要交互式地调节繁多的 参数，使得人们难以操作和使用。以上种种因素，都限制了现有离焦模糊图像盲复原方法的 推广和应用。

发明内容

本发明要解决的技术问题为：克服现有方法对离焦模糊图像的模糊参数辨识精度低、复 原效果不佳、复原速度慢的问题，利用图像中频域的原理和方法，提供一种离焦模糊图像的 快速盲复原方法。

本发明采用的技术方案为：一种基于中频的离焦模糊图像的中频盲复原法，实现步骤如 下：

步骤(1)、计算离焦模糊图像的频谱，方法为：计算离焦模糊图像的离散傅里叶变换并 进行中心化和对数化。

步骤(2)、计算离焦模糊图像频谱的椭圆轨迹均值函数，方法为：从离焦模糊图像频谱 中心开始，每隔一定距离就在频谱中生成一个椭圆轨迹，该椭圆的双轴之比等于图像的高宽 比，然后计算该椭圆轨迹下所有频谱值的均值，从而构造一个以上述距离为自变量的椭圆轨 迹均值函数。

步骤(3)、对椭圆轨迹均值函数进行“变跨度平滑滤波”得到一个基准函数，方法为：用 指定的跨度对椭圆轨迹均值函数进行移动窗平均平滑滤波，但在边界附近如果窗口超出边界 则相应地减小跨度，使窗口刚好包含边界，平滑后即得到上述基准函数。

步骤(4)、用椭圆轨迹均值函数计算离焦模糊图像频谱的中频域，方法为：先对椭圆轨 迹均值函数进行“移动右平均”平滑，然后用直线连接该平滑曲线的两个端点并计算到该直线 最大距离的点，即得到中频域和高频域之间的分界点；再在中低频域上对椭圆轨迹均值函数 的自然指数函数进行“移动左平均”平滑，然后用直线连接该平滑曲线的两个端点并计算到该 直线最大距离的点，即得到低频域和中频域之间的分界点。

步骤(5)、在中频域中用椭圆轨迹均值函数及其基准函数估计第一暗环的位置，方法为： 在中频域中计算基准函数与椭圆轨迹均值函数的差值函数，并以该差值函数第一个凸包的全 局最尖锐处作为第一暗环的估计位置。

步骤(6)、用估计的第一暗环位置计算离焦半径并生成离焦光学传输函数，方法为：用 第一暗环的轨迹方程求解离焦半径，再代入离焦光学传输函数。

步骤(7)、用基于中频的维纳滤波器复原图像，方法为：将上述中频域和高频域之间的 分界点代入基于中频的维纳滤波器，得到复原图像的离散傅里叶变换，再进行离散傅里叶逆 变换，即可得到复原图像。

本发明与现有技术相比的有益效果在于：能够从离焦模糊图像中自动辨识出离焦半径等 参数，实现离焦模糊图像的有效复原，并且具有抑噪能力强、交互参数少、复原速度快的优 点。图1是本发明的一个典型的实施例，(a)是一幅噪声很大的真实离焦模糊图像，在光线 较暗且尚未完成对焦时拍下，原图中完全无法辨认任何信息。经本发明复原后，里面的文字 已经清晰可辨，复原时间仅为毫秒级，复原结果如图(b)所示。

附图说明

图1是本发明的一个典型的实施例，其中(a)是一幅噪声很大的真实离焦模糊图像，(b) 是本发明的复原结果。

图2是该实施例中离焦模糊图像的频谱示意图。

图3是该实施例中离焦模糊图像频谱的椭圆轨迹均值函数曲线、基准函数曲线、中频域 和高频域之间的分界线和低频域和中频域之间的分界线示意图。

图4是该实施例中基准函数与椭圆轨迹均值函数的差值函数曲线和第一暗环估计位置示 意图。

具体实施方式

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

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

(1)计算离焦模糊图像的频谱，方法为：计算离焦模糊图像的离散傅里叶变换并进行 中心化和对数化。其中，离散傅里叶变换及其中心化用常规的计算公式即可实现，此处不再 赘述，对数化可由以下公式实现：

G_ln(u,v)＝ln[|G(u,v)|+ε]    (6)

式中，G_ln(u,v)是对数化后的频谱，G(u,v)是中心化后的离散傅里叶变换，ε是一个很小的 正数，用来防止出现负无穷值。下文中如无特别说明，“频谱”都是指对数化后的频谱G_ln(u,v)。 本实施例的离焦模糊图像如图1(a)所示，其频谱如图2所示。

(2)计算离焦模糊图像频谱的椭圆轨迹均值函数，方法为：从离焦模糊图像频谱中心 开始，每隔一定距离就在频谱中生成一个椭圆轨迹，该椭圆的双轴之比等于图像的高宽比， 然后计算该椭圆轨迹下所有频谱值的均值，从而构造一个以上述距离为自变量的椭圆轨迹均 值函数，其计算公式为：

C(ρ)＝E[G_ln(u,v)],(u,v)∈{(u,v)|(u/M)²+(v/N)²＝ρ²}     (7)

式中，ρ是距离量，C(ρ)是椭圆轨迹均值函数，E[·]表示求均值。即指定一个ρ时，用轨 迹(u/M)²+(v/N)²＝ρ²上所有G_ln(u,v)的均值，作为C(ρ)的函数值。因频谱为中心对称， 故ρ的取值范围为0≤ρ≤0.5。本实施例的椭圆轨迹均值函数曲线如图3所示。

(3)对椭圆轨迹均值函数进行“变跨度平滑滤波”得到一个基准函数，方法为：用指定 的跨度对椭圆轨迹均值函数进行移动窗平均平滑滤波，但在边界附近如果窗口超出边界则相 应地减小跨度，使窗口刚好包含边界，平滑后即得到上述基准函数。设ρ的量化间隔为Δρ， 则基准函数的计算公式为：

其中，表示基准函数，s_ρ是移动窗的跨度(奇数)，i是量化间隔的计数变量，表示 向下取整。移动窗平均平滑在边界附近通常是采用左右扩展的方法来进行处理，但这样会一 定程度上偏离原频谱。因此，本发明采用“变跨度”的方法来使边界附近的和C(ρ)保持 基本一致，也就是让跨度s_ρ随ρ变化：在中间部分用指定的跨度进行处理，而在边界附近用 不大于ρ到边界距离的跨度进行处理。又因为ρ的取值范围为0≤ρ≤0.5，故取

其中，ξ是跨度与总宽度的比例，通常用ξ＝0.1～0.2，min(·)表示取最小值。本实施例的基 准函数函数曲线如图3所示。

(4)用椭圆轨迹均值函数计算离焦模糊图像频谱的中频域，方法为：先对椭圆轨迹均 值函数进行“移动右平均”平滑，然后用直线连接该平滑曲线的两个端点并计算到该直线最大 距离的点，即得到中频域和高频域之间的分界点；再在中低频域上对椭圆轨迹均值函数的自 然指数函数进行“移动左平均”平滑，然后用直线连接该平滑曲线的两个端点并计算到该直线 最大距离的点，即得到低频域和中频域之间的分界点。

根据中频域的特性，模糊图像的频谱在中频域中存有点扩散函数的信息，而低频域和高 频域都会受到不同程度的污染。因此如果能定位中频域，则更能准确的找出第一暗环的位置。 假设中频域和高频域之间的分界点为ρ_MH，低频域和中频域之间的分界点为ρ_LM，则 0≤ρ＜ρ_LM为低频域，ρ_LM≤ρ≤ρ_MH为中频域，ρ_MH＜ρ≤0.5为高频域。因为位于频域以 右的高频域是用转折特性来界定，所以可以直接对C(ρ)进行“移动右平均”平滑，即：

${\tilde{C}}_r\left(\rho \right)=\frac{1}{s_r}\underset{i=0}{\overset{s_r-1}{\Sigma }}C(\rho +i·\mathrm{\Delta \rho }),0\le \rho \le 0.5---\left(10\right)$

式中，是C(ρ)的移动右平均平滑结果，s_r是移动右平均的跨度。从而，用直线连接的端点并计算到该直线最大距离的点，即得到中频域和高频域之间的分界点，计算公式为：

${\rho }_{\mathrm{MH}}=\arg \underset{0<\rho <0.5}{\max }\left|\right[{\tilde{C}}_r\left(0\right)-{\tilde{C}}_r\left(0.5\right)]\rho +0.5{\tilde{C}}_r\left(\rho \right)-0.5{\tilde{C}}_r\left(0\right)|---\left(11\right)$

式中，argmax表示取指定区域中使表达式最大的变量值。另一方面，由于0＜ρ_LM＜ρ_MH， 故计算ρ_LM时只需考虑中低频域0≤ρ≤ρ_MH。因为位于频域以左的低频域是用图像能量来界 定，所以需要对exp[C(ρ)]进行“移动左平均”平滑，即：

${\tilde{C}}_l\left(\rho \right)=\frac{1}{s_l}\underset{i=0}{\overset{s_l-1}{\Sigma }}\exp \left[C(\rho -i·\mathrm{\Delta \rho })\right],0\le \rho \le {\rho }_{\mathrm{MH}}---\left(12\right)$

式中，是exp[C(ρ)]的移动左平均平滑结果，s_l是移动左平均的跨度。类似的，用直线 连接exp[C(ρ)]的端点并计算到该直线最大距离的点，即得到低频域和中频域之间的分界点， 计算公式为：

${\rho }_{\mathrm{LM}}=\arg \underset{0<\rho <{\rho }_{\mathrm{MH}}}{\max }\left|\right[{\tilde{C}}_l\left(0\right)-{\tilde{C}}_l\left({\rho }_{\mathrm{MH}}\right)]\rho +{\rho }_{\mathrm{MH}}{\tilde{C}}_l\left(\rho \right)-{\rho }_{\mathrm{MH}}{\tilde{C}}_l\left(0\right)|---\left(13\right)$

本实施例的中频域和高频域之间的分界线以及低频域和中频域之间的分界线都示于如 图3之中。

(5)在中频域中用椭圆轨迹均值函数及其基准函数估计第一暗环的位置，方法为：在 中频域中计算基准函数与椭圆轨迹均值函数的差值函数，并以该差值函数第一个凸包的全局 最尖锐处作为第一暗环的估计位置。

理论上，第一暗环的零点会使基准函数与椭圆轨迹均值函数的差值函数产生一个凸包， 但受图像能量和噪声影响，实际中该凸包不一定是最大，而且第一暗环也不一定位于凸包的 峰值处。但是第一暗环几乎都位于第一个凸包中，因此需要用最大值的一定容差来找出第一 个凸包，并且以变化率最大的原则来确定第一暗环，即凸包中的全局最尖锐处。但是噪声的 存在可能会导致个别局部尖锐点，因此还需要首先对差值函数进行平滑处理，从而差值函数 的计算公式为：

$D\left(\rho \right)=\mathrm{smooth}[{\tilde{C}}_r\left(\rho \right)-C_r\left(\rho \right)],{\rho }_{\mathrm{LM}}\le \rho \le {\rho }_{\mathrm{MH}}---\left(14\right)$

式中，D(ρ)表示上述差值函数，smooth[·]是平滑算子。该差值函数定义在中频域中，这是 因为只有中频域才较好地保存了点扩散函数的信息，并且本发明的中频域计算方法会自动地 将第一个凸包包含在内。本实施例的差值函数曲线如图4所示。

然后，以该差值函数最大值的一定比例为阈值，找出差值函数的第一个凸包的区域，即 找出D(ρ)在ρ_LM≤ρ≤ρ_MH中第一个满足D(ρ)≥λ·max[D(ρ)]的连续区域(用Ω₁表示)，其 中λ为上述用来确定阈值的比例，max[·]表示取最大值。再用Ω₁中各点的夹角来量化尖锐的 程度，其计算公式为：

$A\left(\rho \right)=\pi -\arctan \left\{\frac{\mathrm{\Delta \rho }·[2D\left(\rho \right)-D(\rho -\mathrm{\Delta \rho })-D(\rho +\mathrm{\Delta \rho })]}{{\mathrm{\Delta \rho }}^2+[D\left(\rho \right)-D(\rho -\mathrm{\Delta \rho })][D(\rho +\mathrm{\Delta \rho })-D\left(\rho \right)]}\right\},\rho \in {\Omega }_1---\left(15\right)$

式中，A(ρ)表示D(ρ)曲线在ρ处的夹角。于是，A(ρ)最小的位置就是Ω₁中的最尖锐处及 第一暗环的估计位置即：

${\hat{\rho }}_1=\arg \underset{\rho \in {\Omega }_1}{\min }\left[A\left(\rho \right)\right]---\left(16\right)$

式中，argmin表示取指定区域中使表达式最小的变量值。本实施例的第一暗环估计位置如 图4所示。

(6)用估计的第一暗环位置计算离焦半径并生成离焦光学传输函数，方法为：用第一 暗环的轨迹方程求解离焦半径，再代入离焦光学传输函数。用第一暗环的轨迹方程可推出离 焦半径的计算公式为：

$\hat{r}=\frac{3.85}{2\pi {\hat{\rho }}_1}---\left(17\right)$

式中，表示估计的离焦半径。本实施例中的计算结果为再将代入离焦光学 传输函数，即可生成估计的离焦光学传输函数

$\hat{H}(u,v)=K\frac{J_1\left(2\pi \hat{r}\sqrt{{(u/M)}^2+{(v/N)}^2}\right)}{2\pi \hat{r}\sqrt{{(u/M)}^2+{(v/N)}^2}}---\left(18\right)$

(7)用基于中频的维纳滤波器复原图像，方法为：将上述中频域和高频域之间的分界 点代入基于中频的维纳滤波器，得到复原图像的离散傅里叶变换，再进行离散傅里叶逆变换， 即可得到复原图像。本发明针对离焦模糊图像的复原，推导出基于中频的维纳滤波器的计算 公式为：

$\hat{F}(u,v)=\frac{{\hat{H}}^{*}(u,v)G(u,v)}{{\left|\hat{H}(u,v)\right|}^2+{\left|\tilde{H}\left({\rho }_{\mathrm{MH}}\right)\right|}^{2\eta }·{\left[\max \right|\tilde{H}\left(\rho \right)\left|\right]}^{2(1-\eta )}-{\left[\min \right|\tilde{H}\left(\rho \right)\left|\right]1}^2}---\left(19\right)$

式中，是复原图像的离散傅里叶变换，η是噪声抑制参数，是在v＝0截 面上的曲线的一个平滑版本，即：

$\tilde{H}\left(\rho \right)=\mathrm{smooth}\left[\hat{H}(\mathrm{\rho M},0)\right]---\left(20\right)$

其中平滑算子smooth[·]的目的是使基本单调，消除其零点以便取值。另外，噪声 抑制参数η用于调和噪声与图像细节的矛盾：η接近1时噪声和细节都较为清晰，η减小时 噪声和细节同时有所减弱，从中选取一个平衡值以实现图像的最佳复原。最后，计算的离散傅里叶逆变换并标准化(必要时可调整对比度)，即得复原图像。本实施例的复原图 像如图1(b)所示。

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