CN111815537B - 一种新型图像盲解去模糊方法 - Google Patents

Abstract

本发明公开了一种新型图像盲解去模糊方法，首先根据模糊图像的信息估计卷积核的近似值，然后通过构造目标函数求解理想清晰图像的近似值，直接采用无参考图像质量评价函数作为优化函数，最后利用模拟退火算法对优化函数进行优化，对初步估计卷积核和初步理想清晰图像做进一步迭代求解，得到最终的理想清晰图像和卷积核。本发明直接将图像质量评价函数作为优化函数，只要选取的图像质量评价函数可靠，可以保证最终得到的去模糊图像质量。

Description

一种新型图像盲解去模糊方法

技术领域

本发明属于图像处理领域，具体涉及一种图像去模糊方法。

背景技术

模糊图像可视为清晰图像与一卷积核卷积后叠加噪声后形成的。传统的图像去卷积方法包括维纳滤波算法、RL滤波算法、总变分算法，这些算法未能充分利用图像的先验信息故而效果有待提升。采用正则化的方法即构造目标函数对待求解的理想清晰图像和卷积核进行约束对模糊图像去模糊取得了不错的效果；利用图像的先验信息采用最大后验概率的方法对图像进行去模糊也有很大的进展，但是这些方法缺乏理论基础，在有些情况下甚至无法得到最优解。

发明内容

为了克服现有技术的不足，本发明提供一种新型图像盲解去模糊方法，直接采用无参考图像质量评价函数作为优化函数，结合模拟退火算法，对所给模糊退化的图像进行去模糊，可同时得到卷积核和理想清晰图像。通过选取可靠的图像质量评价函数，最终一定可以得到理想的清晰图像。

为达到上述目的，本发明提供的一种新型图像盲解去模糊方法，包括以下步骤：

步骤1：建立图像盲解去模糊目标函数；

定义模糊图像生成的数学模型为：

Y＝X*h+n (1)

式中，X表示理想清晰图像，Y表示生成的模糊图像，h表示卷积核，n表示加性噪声；

构造目标函数为：

式中，s表示拉普拉斯算子，λ₁，λ₂，λ₃，λ₄表示各项的惩罚因子；tv(X)表示图像X的全变分约束：

g_h，g_v表示水平和垂直微分算子；

固定h，将目标函数转化到频域为：

固定H_w，进行求解得到：

X_w，H_w，Y_w，S，G_h，G_v分别为X，h，Y，s，g_h，g_v的傅里叶变换，F^-1(*)为傅里叶反变换，I为单位矩阵；

计算结果

记为X⁰，为初始理想清晰图像；

步骤2：构造优化函数f(X)；

将图像X转换为灰度图像，再分为若干块大小相同的矩形小块，其中任一小块记为x，大小为n₁×n₂，构造三个中间函数S₁(x)，S₂(x)，S₃(x)；

步骤2-1：基于频域构造S₁(x)；

当max(con)-min(con)≤th1或者μ(x)≤th2时，此时称x的对比度为0，当x的对比度为0时记为S₁(x)＝0；式中，

表示x的照度值，φ，

γ均为照度参数，th1、th2均为阈值；μ(x)表示图像小块x的像素灰度算术均值；

当max(con)-min(con)>th1且μ(x)>th2时，此时称x的对比度大于0；

记y_x(ω,θ)为x的傅里叶变换，ω为圆频率，θ为方位角，计算

z_x(ω)为x的一维幅度谱；

用-αlogω+logβ拟合z_x(ω)，α、β均为为拟合参数，求解x的一维幅度谱的斜率α_x：

则

其中τ₁、τ₂为均为中间参数；

步骤2-2：基于空间域构造S₂(x)：

定义x的全变分为：

构造

式中，ξ为图像小块x中2×2的小块，x_i，x_j为2×2小块中任意两个不同位置的像素；

步骤2-3：构造S₃(x)：

S₃(x)＝S₁(x)^η×S₂(x)^1-η (7)

式中，η为平衡参数；

将图像X中划分的所有图像小块代入S₃(x)计算，得到向量S₃(X)，对S₃(X)中的元素S₃(x)从大到小排序，取前N个元素求平均值：

式中，

为S₃(X)中的第g个元素；

令

为模拟退火算法中的优化函数；

步骤3：确定迭代初始卷积核：

步骤3-1：定义矩阵向量化表达式：

式中，v(*)表示对矩阵进行向量化，

表示将参与卷积运算的矩阵展开形成的卷积矩阵，m₁，m₂为h的长和宽；

定义卷积特征值：

式中，s₁、s₂为采样尺寸，l为一个卷积核，当‖E(X)*l‖_F取得最大值时，此时的l是E(X)的第i个卷积特征向量，记为

此时‖E(X)*l‖_F的值称为E(X)的第i个卷积特征值，记为

E(X)＝X*s，表示提取图像的边缘特征；<*>表示内积运算；约束条件||l||_F＝1,

表示要求卷积特征向量之间互相正交且模为1；

对式(1)，忽略加性噪声n，则

Y＝X*h (11)

对等式两边卷积拉普拉斯算子，得到：

E(Y)＝E(X)*h (12)

根据卷积特征值与卷积特征向量的定义，有如下性质：

其中

表示图像的最小卷积特征值；

步骤3-2：推导卷积核的初始迭代值h⁰；

根据模拟退火算法定义有：

进一步有：

因为

且

则可推导出：

等价为如下函数：

令

可进一步化简为：

w^E(Y)(h)＝(v(h))^TWv(h)

待求解的

记为迭代初始卷积核h⁰；对W进行特征值分解，最小特征值对应的特征向量展成m₁×m₂大小的矩阵即为h⁰；

步骤4：对模糊图像去模糊；

步骤4-1：求迭代初始值：

根据步骤1中的式(4)求解初始理想图像X⁰，作为迭代的初始值；

根据步骤3求解迭代初始值h⁰；

步骤4-2：用模拟退火算法求解最终的理想解X^*，h^*，迭代过程为：

步骤4-2-1：令i等于0，k等于0，预设初始温度T₀；

定义δ＝a_i·randn(m₁,m₂)，a_i为控制δ中元素范围的因子，randn(*)为随机数产生矩阵；

步骤4-2-2：根据h^k+1＝h⁰+δ得到h^k+1，按下式计算接受概率：

式中，T_i为模拟退火算法中第i次迭代时的温度；

若p(h^k→h^k+1)>u，则h^k+1＝h⁰+δ，利用h^k+1和步骤1中的式(4)，重新求解更新X^k，并重新生成δ；若p(h^k→h^k+1)≤u，则h^k+1＝h^k；u为在预设区间范围内生成的随机数；

步骤4-2-3：如果k小于迭代阈值，则重复步骤4-2-2；如果k大于或等于迭代阈值，则进行下一步；令k加1；

步骤4-2-4：如果T_i≤thresh，则结束迭代，否则继续迭代过程，thresh为温度阈值；

步骤4-2-5：计算T_i+1＝ρT_i，a_i+1＝ρa_i，ρ为温度衰减系数且ρ＜1，令i加1，转到步骤4-2-2；

步骤4-3：当迭代结束时，得到X^*＝X^k为最终理想清晰图像，h^*＝h^k为最终的卷积核。

由于采用了本发明提出的一种新型图像盲解去模糊方法，得到了如下有益效果：

1、其他图像去模糊算法的基本思路是对待求解的图像进行约束后构造目标函数然后求解，对得到的解进行图像质量评价从而说明算法的有效性；本发明直接将图像质量评价函数作为优化函数，只要选取的图像质量评价函数可靠，最终得到的去模糊图像是有保证的。

2、求解去模糊图像时，用卷积核和原图像表示待求解的去模糊图像，求解去模糊图像是通过对卷积核采样更新实现的，所以最终可同时得到去模糊图像和卷积核。

附图说明

图1是一幅原始图像。

图2是图1经过卷积核处理之后的模糊图像。

图3是采用本发明方法对图2去模糊后得到的理想清晰图像。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

本发明提供一种新型图像盲解去模糊方法，包括以下步骤：

步骤1：建立图像盲解去模糊目标函数；

定义模糊图像生成的数学模型为：

Y＝X*h+n (1)

式中，X表示理想清晰图像，Y表示生成的模糊图像，h表示卷积核，n表示加性噪声；

构造目标函数为：

式中，s表示拉普拉斯算子，λ₁，λ₂，λ₃，λ₄表示各项的惩罚因子；tv(X)表示图像X的全变分约束：

g_h，g_v表示水平和垂直微分算子；

固定h，将目标函数转化到频域为：

固定H_w，进行求解得到：

X_w，H_w，Y_w，S，G_h，G_v分别为X，h，Y，s，g_h，g_v的傅里叶变换，F^-1(*)为傅里叶反变换，I为单位矩阵；

计算结果

记为X⁰，为初始理想清晰图像；

步骤2：构造优化函数f(X)；

将图像X转换为灰度图像，再分为若干块大小相同的矩形小块，其中任一小块记为x，大小为n₁×n₂，构造三个中间函数S₁(x)，S₂(x)，S₃(x)；

步骤2-1：基于频域构造S₁(x)；

当max(con)-min(con)≤th1或者μ(x)≤th2时，此时称x的对比度为0，当x的对比度为0时记为S₁(x)＝0；式中，

表示x的照度值，φ，

γ均为照度参数，th1、th2均为阈值；μ(x)表示图像小块x的像素灰度算术均值；

当max(con)-min(con)>th1且μ(x)>th2时，此时称x的对比度大于0；

记y_x(ω,θ)为x的傅里叶变换，ω为圆频率，θ为方位角，计算

z_x(ω)为x的一维幅度谱；

用-αlogω+logβ拟合z_x(ω)，α、β均为为拟合参数，求解x的一维幅度谱的斜率α_x：

则

其中τ₁、τ₂为均为中间参数；

步骤2-2：基于空间域构造S₂(x)：

定义x的全变分为：

构造

式中，ξ为图像小块x中2×2的小块，x_i，x_j为2×2小块中任意两个不同位置的像素；

步骤2-3：构造S₃(x)：

S₃(x)＝S₁(x)^η×S₂(x)^1-η (7)

式中，η为平衡参数；

将图像X中划分的所有图像小块代入S₃(x)计算，得到向量S₃(X)，对S₃(X)中的元素S₃(x)从大到小排序，取前N个元素求平均值：

式中，

为S₃(X)中的第g个元素；

令

为模拟退火算法中的优化函数；

步骤3：确定迭代初始卷积核：

步骤3-1：定义矩阵向量化表达式：

式中，v(*)表示对矩阵进行向量化，

表示将参与卷积运算的矩阵展开形成的卷积矩阵，m₁，m₂为h的长和宽；

定义卷积特征值：

式中，s₁、s₂为采样尺寸，l为一个卷积核，当‖E(X)*l‖_F取得最大值时，此时的l是E(X)的第i个卷积特征向量，记为

此时‖E(X)*l‖_F的值称为E(X)的第i个卷积特征值，记为

E(X)＝X*s，表示提取图像的边缘特征；<*>表示内积运算；约束条件||l||_F＝1,

表示要求卷积特征向量之间互相正交且模为1；

对式(1)，忽略加性噪声n，则

Y＝X*h (11)

对等式两边卷积拉普拉斯算子，得到：

E(Y)＝E(X)*h (12)

根据卷积特征值与卷积特征向量的定义，有如下性质：

其中

表示图像的最小卷积特征值；

步骤3-2：推导卷积核的初始迭代值h⁰；

根据模拟退火算法定义有：

进一步有：

因为

且

则可推导出：

等价为如下函数：

令

可进一步化简为：

w^E(Y)(h)＝(v(h))^TWv(h)

待求解的

记为迭代初始卷积核h⁰；对W进行特征值分解，最小特征值对应的特征向量展成m₁×m₂大小的矩阵即为h⁰；

步骤4：对模糊图像去模糊；

步骤4-1：求迭代初始值：

根据步骤1中的式(4)求解初始理想清晰图像X⁰，作为迭代的初始值；

根据步骤3求解迭代初始值h⁰；

步骤4-2：用模拟退火算法求解最终的理想解X^*，h^*，迭代过程为：

步骤4-2-1：令i等于0，k等于0，预设初始温度T₀；

定义δ＝a_i·randn(m₁,m₂)，a_i为控制δ中元素范围的因子，randn(*)为随机数产生矩阵；

步骤4-2-2：根据h^k+1＝h⁰+δ得到h^k+1，按下式计算接受概率：

式中，T_i为模拟退火算法中第i次迭代时的温度；

若p(h^k→h^k+1)>u，则h^k+1＝h⁰+δ，利用h^k+1和步骤1中的式(4)，重新求解更新X^k，并重新生成δ；若p(h^k→h^k+1)≤u，则h^k+1＝h^k；u为在预设区间范围内生成的随机数；

步骤4-2-3：如果k小于迭代阈值，则重复步骤4-2-2；如果k大于或等于迭代阈值，则进行下一步；令k加1；

步骤4-2-4：如果T_i≤thresh，则结束迭代，否则继续迭代过程，thresh为温度阈值；

步骤4-2-5：计算T_i+1＝ρT_i，a_i+1＝ρa_i，ρ为温度衰减系数且ρ＜1，令i加1，转到步骤4-2-2；

步骤4-3：当迭代结束时，得到X^*＝X^k为最终理想清晰图像，h^*＝h^k为最终的卷积核。

实施例：

对于一幅模糊图像，采用本发明方法按如下步骤去除模糊。

1、建立目标函数求解初始理想清晰图像；

构造目标函数为：

固定h，目标函数变为：

由于图像较大时运算效率低为求解方便，再将其转换到频域进行求解，最终可以得到初始理想清晰图像X⁰。

2、构造优化函数f(X)；

将图像X转换为灰度图像，再分为若干块大小相同的矩形小块，其中任一小块记为x，大小为15×15，构造三个中间函数S₁(x)，S₂(x)，S₃(x)；

(1)设定th1＝5，th2＝2，φ＝0.7656，

γ＝2.2。

当max(con)-min(con)≤5或者μ(x)≤2时，此时称x的对比度为0，记S₁(x)＝0；con＝(0.7656+0.0364)^2.2，代表x的照度值，图像的数据类型为八位无符号整型。

当max(con)-min(con)>5且μ(x)>2时，此时称x的对比度大于0；

取τ₁＝-3，τ₂＝2，则得到

(2)

(3)取η＝0.5，N＝100；

S₃(x)＝S₁(x)^0.5×S₂(x)^0.5

将图像X中划分的所有图像小块代入S₃(x)计算，得到向量S₃(X)，对S₃(X)中的元素S₃(x)从大到小排序，取前100个元素求平均值：

最终得到优化函数：

3、根据卷积特征向量的定义易知，求解E(X)的卷积特征值和卷积特征向量可通过对矩阵

进行奇异值分解，得到S、V、D，矩阵V的对角线元素即为卷积特征值，矩阵D中列向量为相应的卷积特征向量。

卷积核的尺寸设定m₁＝m₂＝15，采样尺寸s₁＝s₂＝17；

根据步骤3的方法得到迭代初始卷积核h⁰。

4、对模糊图像去模糊；

(1)根据步骤3求解迭代初始值h⁰。

根据步骤1中的目标函数和输入的模糊图像Y，求解

记为X⁰，作为迭代的初始值，Y如图2所示；

(2)首先选取T_i+1＝ρT_i，ρ＝0.9；其次选取h^k+1＝h⁰+δ，引入与h^k维度相同的随机向量δ＝a_i·randn(15,15)，randn(*)的输入参数为待生成随机数矩阵的大小，δ在每次迭代时均需重新生成；卷积核中各元素互相独立且满足高斯分布，其均值为h⁰中对应的元素，标准差为δ中对应元素；a₀＝0.01，a_i与T_i同步且采用相同的机制更新，即a_i+1＝ρa_i；在每个固定的温度T_i下迭代1500次认为达到平衡状态，当T_i≤0.1则停止迭代，具体迭代过程为：

S1：令i＝0，k＝0，给出初始温度T₀＝1，δ＝0.01randn(15,15)，估计h⁰，结合Y求解X⁰；

S2：根据h^k+1＝h⁰+δ得到h^k+1，按下式计算接受概率：

若p(h^k→h^k+1)>u，则h^k+1＝h⁰+δ，重新求解更新X^k，重新生成δ；

若p(h^k→h^k+1)≤u，则h^k+1＝h^k；u为在区间[0.7,1]范围内生成的随机数；

S3：如果k小于1500，则重复步骤S2；如果k大于或等于1500，则进行下一步；令k加1；

S4：如果T_i≤0.1，则结束迭代，否则继续迭代过程；

S5：计算T_i+1＝0.9T_i和a_i+1＝0.9a_i，更新T_i和a_i，令i加1，转到步骤S2；

(3)当迭代结束时，得到X^*＝X^k为最终理想清晰图像，如图3所示，h^*＝h^k为最终的卷积核。

如图1是原始图像，对比图2和图3可以看出，对于图2模糊甚至白化的区域，在图3中有较好的恢复。

Claims (1)

1.一种新型图像盲解去模糊方法，其特征在于，包括以下步骤：

步骤1：建立图像盲解去模糊目标函数；

定义模糊图像生成的数学模型为：

Y＝X*h+n (1)

式中，X表示理想清晰图像，Y表示生成的模糊图像，h表示卷积核，n表示加性噪声；

构造目标函数为：

式中，s表示拉普拉斯算子，λ₁，λ₂，λ₃，λ₄表示各项的惩罚因子；tv(X)表示图像X的全变分约束：

g_h，g_v表示水平和垂直微分算子；

固定h，将目标函数转化到频域为：

固定H_w，进行求解得到：

X_w，H_w，Y_w，S，G_h，G_v分别为X，h，Y，s，g_h，g_v的傅里叶变换，F^-1(*)为傅里叶反变换，I为单位矩阵；

计算结果

记为X⁰，为初始理想清晰图像；

步骤2：构造优化函数f(X)；

将图像X转换为灰度图像，再分为若干块大小相同的矩形小块，其中任一小块记为x，大小为n₁×n₂，构造三个中间函数S₁(x)，S₂(x)，S₃(x)；

步骤2-1：基于频域构造S₁(x)；

当max(con)-min(con)≤th1或者μ(x)≤th2时，此时称x的对比度为0，当x的对比度为0时记为S₁(x)＝0；式中，

表示x的照度值，φ，

γ均为照度参数，th1、th2均为阈值；μ(x)表示图像小块x的像素灰度算术均值；

当max(con)-min(con)>th1且μ(x)>th2时，此时称x的对比度大于0；

记y_x(ω,θ)为x的傅里叶变换，ω为圆频率，θ为方位角，计算

z_x(ω)为x的一维幅度谱；

用-αlogω+logβ拟合z_x(ω)，α、β均为为拟合参数，求解x的一维幅度谱的斜率α_x：

则

其中τ₁、τ₂为均为中间参数；

步骤2-2：基于空间域构造S₂(x)：

定义x的全变分为：

构造

式中，ξ为图像小块x中2×2的小块，x_i，x_j为2×2小块中任意两个不同位置的像素；

步骤2-3：构造S₃(x)：

S₃(x)＝S₁(x)^η×S₂(x)^1-η (7)

式中，η为平衡参数；

将图像X中划分的所有图像小块代入S₃(x)计算，得到向量S₃(X)，对S₃(X)中的元素S₃(x)从大到小排序，取前N个元素求平均值：

式中，

为S₃(X)中的第g个元素；

令

为模拟退火算法中的优化函数；

步骤3：确定迭代初始卷积核：

步骤3-1：定义矩阵向量化表达式：

式中，v(*)表示对矩阵进行向量化，

表示将参与卷积运算的矩阵展开形成的卷积矩阵，m₁，m₂为h的长和宽；

定义卷积特征值：

式中，s₁、s₂为采样尺寸，l为一个卷积核，当‖E(X)*l‖_F取得最大值时，此时的l是E(X)的第i个卷积特征向量，记为

此时‖E(X)*l‖_F的值称为E(X)的第i个卷积特征值，记为

E(X)＝X*s，表示提取图像的边缘特征；<*>表示内积运算；约束条件||l||_F＝1,

表示要求卷积特征向量之间互相正交且模为1；

对式(1)，忽略加性噪声n，则

Y＝X*h (11)

对等式两边卷积拉普拉斯算子，得到：

E(Y)＝E(X)*h (12)

根据卷积特征值与卷积特征向量的定义，有如下性质：

其中

表示图像的最小卷积特征值；

步骤3-2：推导卷积核的初始迭代值h⁰；

根据模拟退火算法定义有：

进一步有：

因为

且

则可推导出：

等价为如下函数：

令

可进一步化简为：

w^E(Y)(h)＝(v(h))^TWv(h)

待求解的

记为迭代初始卷积核h⁰；对W进行特征值分解，最小特征值对应的特征向量展成m₁×m₂大小的矩阵即为h⁰；

步骤4：对模糊图像去模糊；

步骤4-1：求迭代初始值：

根据步骤1中的式(4)求解初始理想图像X⁰，作为迭代的初始值；

根据步骤3求解迭代初始值h⁰；

步骤4-2：用模拟退火算法求解最终的理想解X^*，h^*，迭代过程为：

步骤4-2-1：令i等于0，k等于0，预设初始温度T₀；

定义δ＝a_i·randn(m₁,m₂)，a_i为控制δ中元素范围的因子，randn(*)为随机数产生矩阵；

步骤4-2-2：根据h^k+1＝h⁰+δ得到h^k+1，按下式计算接受概率：

式中，T_i为模拟退火算法中第i次迭代时的温度；

若p(h^k→h^k+1)>u，则h^k+1＝h⁰+δ，利用h^k+1和步骤1中的式(4)，重新求解更新X^k，并重新生成δ；若p(h^k→h^k+1)≤u，则h^k+1＝h^k；u为在预设区间范围内生成的随机数；

步骤4-2-3：如果k小于迭代阈值，则重复步骤4-2-2；如果k大于或等于迭代阈值，则进行下一步；令k加1；

步骤4-2-4：如果T_i≤thresh，则结束迭代，否则继续迭代过程，thresh为温度阈值；

步骤4-2-5：计算T_i+1＝ρT_i，a_i+1＝ρa_i，ρ为温度衰减系数且ρ＜1，令i加1，转到步骤4-2-2；

步骤4-3：当迭代结束时，得到X^*＝X^k为最终理想清晰图像，h^*＝h^k为最终的卷积核。

Images