CN105046659A - 一种基于稀疏表示的单透镜计算成像psf估算方法 - Google Patents

Abstract

本发明公开了一种基于稀疏表示的单透镜计算成像PSF估算方法。针对单透镜计算成像中对PSF估计速度和精度的要求，尝试用稀疏表示的方法来估计PSF。首先将目标函数中的清晰图像表示为过完备字典和稀疏系数的乘积，并对稀疏系数进行约束。然后通过迭代优化算法依次交替估计出模糊核、过完备字典和稀疏系数。即可得到单透镜计算成像所需的模糊核。在训练过完备字典时，只需知道由单透镜成像系统得到的模糊图像，无需更多的额外信息，而且在计算过程中避免了多次卷积操作，减小了对复原图像边缘信息的影响。本方法操作简单，在图像处理和相机设计领域都具有非常重要的意义。

Description

一种基于稀疏表示的单透镜计算成像PSF估算方法

技术领域

本发明主要涉及到数字图像处理领域，特指一种基于稀疏表示的单透镜计算成像PSF估算方法。

背景技术

目前，单反相机以其高清的成像质量、丰富的镜头选择、迅捷的响应速度、卓越的手控能力等优势在人们的日常生活中发挥着越来越重要的作用。然而，为弥补单反镜头中镜片的几何畸变和像差，进一步提高成像质量，单反镜头的设计日益复杂，甚至包含数十个独立的光学器件。复杂的镜头在提高成像质量的同时，无疑也会增加镜头的体积和重量，也导致镜头的成本大大提高。镜头体积和重量的增加给用户的日常使用带来了不便，成本的提高也不便单反向大面积用户推广使用。因此，在尽量消除镜片像差，增加成像质量的同时，如何降低镜头成本，使其更为轻便，也成为目前单反相机设计的重要需求之一。近年来，随着图像复原技术的快速发展，图像去模糊等方法越来越成熟，镜头中某些消除像差和修正几何畸变的镜片可由去模糊等计算摄影技术代替，因此，单透镜成像(如图1所示)与图像复原技术的结合也逐渐成为单反相机设计的一个新的研究方向。

单透镜成像的关键在于必须首先对单透镜的模糊核进行标定，然后基于所得到的模糊核，采用一定的图像复原算法得到清晰图像。

单透镜的模糊核标定通常是在尺度空间中进行的，如图2所示，在尺度空间中，一般取3×3的高斯函数或者delta函数作为模糊核的初始值，通过依次在不同层次尺度空间中的迭代逐步求得最终的理想模糊核。而且在每一层尺度空间，首先将上一层次尺度空间中所求得的模糊核作为初始值，结合模糊图像求得潜在的清晰图像，然后把潜在清晰图像和模糊核作为已知项，再求出清晰图像。这个过程也需要经过多次迭代才能得到本层次尺度空间中较为理想的模糊核，如图3所示。

目前主流的图像复原算法中，求模糊核都要通过卷积运算，而卷积运算在一定程度上会影响模糊核估计的精度，导致图像存在振铃效应。

过完备稀疏表示作为一种有效的表示模型，广泛地应用于各种信号和图像处理任务中。过完备稀疏表示采用过完备基来代替传统的正交基函数，在过完备基里寻求少量元素来表示图像的特征，由于所采用基函数的冗余性，因此具有更好的稀疏表示能力。众所周知，图像或信号可以用过完备字典中的某几个基来线性表示，如图4所示，y表示给定信号，D表示过完备字典，x表示信号对应的稀疏系数，如果给定过完备字典D，则信号y可近似由字典中的某些原子线性组合表示，即其中x_i表示向量的每个分量，d_i表示字典中的一个原子。

稀疏表示在图像去噪、图像复原等领域都取得了很好的效果，而单透镜计算成像本质也是图像复原问题，因此，可以尝试使用稀疏表示来计算单透镜光学成像系统的模糊核，如果从已知模糊图像中训练过完备字典，则可以通过求解稀疏系数得到单透镜成像系统对应的模糊核，这个过程无需进行卷积操作，在一定程度上可以提高模糊核的精确度，从而提高最终的图像复原效果。

发明内容

针对现有单透镜成像方法中模糊核初始值与最终理想值差别较大，迭代估计过程耗时太长，不便于实际操作等问题，本发明提出一种基于稀疏表示的单透镜计算成像PSF估算方法。充分利用单透镜模糊核的空间对称性特点，首先估计出部分区域内的图像块的模糊核，然后将这些模糊核经过适当的旋转调整，作为其他区域内图像块模糊核估计的初始值。相比于高斯函数和delta函数，利用对称性所得到的模糊核更加接近最终的理想模糊核，所需的迭代次数会大大减少，从而在整体上加快模糊核的估计速度。

为解决上述技术问题，本发明的技术方案是，

一种基于稀疏表示的单透镜计算成像PSF估算方法，包括以下步骤：

步骤一：利用单透镜相机获取模糊图像，模糊图像是在正常光圈大小下，由单透镜相机得到的模糊图像。

步骤二：将模糊图像去模糊问题转化为联合优化问题，将目标函数中的清晰图像采用过完备字典D和稀疏系数α的乘积表示，并对稀疏系数的稀疏性加以约束，最终目标函数表示为：

$\underset{k,\mathrm{\α}}{min}{||b-k*D\·A||}_2^2+\mathrm{\λ}{||\mathrm{\α}||}_0---\left(1\right)$

其中，b表示模糊图像，k表示模糊核，D表示过完备字典，A表示稀疏系数矩阵，表示稀疏系数矩阵中的每一列，λ为控制权重的参数。

步骤三：目标函数中有三个未知变量：模糊核k，过完备字典D和稀疏系数α，需要通过交替迭代优化算法分别估计出这三个标量；

首先固定k和D来估计α，给定初始的k₀和D₀，通过OMP(OrthogonalMatchingPursuit)优化算法来估计α。

所采用OMP算法是：以贪婪迭代的方法选择字典D的列，使得在每次迭代中所选择的列与当前列的冗余向量最大程度地相关，通过反复迭代得到满足稀疏程度的稀疏系数并停止迭代，步骤如下：

输入：过完备字典D，清晰图像的图像块对应的向量y，稀疏度l。

输出：稀疏系数α的l-稀疏的逼近α'；

初始化：残差r₀＝y，索引集t＝1；

依次循环步骤1)-5)：

1)找出残差r和过完备字典的列d_i积中最大值所对应的脚标λ，即

λ_i＝argmax_j＝1...N|<r_t-1，d_i>|；

2)更新索引集Λ_t＝Λ_t-1∪{λ_i}，记录找到的过完备字典中重建原子集合

3)由最小二乘得到α_t'＝argmin||y-D_tα_t'||₂；

4)更新残差r_t＝y-D_tα_t'，t＝t+1；

5)判断是否满足t＞l，若满足，则停止迭代；若不满足，则执行1)。

步骤四：固定k和α来更新过完备字典D。模糊图像b不能直接用来构造，首先通过非盲卷积算法对模糊图像进行处理得到对应的潜在清晰图像。

求潜在清晰图像的非盲卷积算法是Levin在2007年所提出的的非盲卷积算法。使用如下概率模型：

$P_k\left(x\right|b)\&Proportional;\exp (-(\frac{1}{{\mathrm{\&eta;}}^2}{\left|C_{k}x-b\right|}^2+\mathrm{\&beta;}{\left|C_{g_x}x\right|}^2+\mathrm{\&beta;}{\left|C_{g_b}x\right|}^2\left)\right)---\left(2\right)$

其中b表示模糊图像，x表示清晰图像，k表示模糊核，g_x＝[1-1]，g_b＝[1-1]^T，C表示卷积操作，对于函数f，C_fx表示f与x的卷积操作，即C_fx≡f*x，β和η是控制权重的参数；则式(2)的最大值即为所求得清晰图像，即x^*＝argmaxP_k(x|b)；

步骤五：固定α和D估计出模糊核k。

清晰图像块I_p可由清晰图像对应的过完备字典中一些原子与稀疏系数线α性组合表示，而模糊图像块可由模糊图像对应的过完备字典中一些原子与稀疏系数线α'性组合表示，因为模糊图像的字典可由清晰图像的字典与模糊核卷积得到，所以不妨假设稀疏系数满足α＝α'。则潜在清晰图像可由下式得到：

$I_p^{(n+1)}=D^{(n+1)}{\mathrm{\&alpha;}}^{(n+1)}---\left(3\right)$

将卷积操作k⁽ⁿ⁺¹⁾*I⁽ⁿ⁺¹⁾表示为矩阵相乘的形式K·X，其中K和X分别为k⁽ⁿ⁺¹⁾和I⁽ⁿ⁺¹⁾对应的矩阵形式，在估计模糊核的目标函数中加入L₁范数图像先验：

$\min {||b-Ax||}_2^2+{||x||}_1---\left(4\right)$

通过EM迭代优化算法来求解目标函数公式(4)即能得到所需的模糊核。

EM算法的具体步骤如下：

(1)E-step：令Q(x)＝P(x|b,k)，并计算出Q(x)的均值u和方差C，其中P(x|b,k)表示在已知模糊图像b和模糊核k的条件下，对应的清晰图像为x的概率；均值u就代表给定当前模糊核时潜在的清晰图像，C代表清晰图像周围的方差；

(2)M-step：找到使k*x-b达到最小值的K值，即式(5)达到最小值；

E_Q||k*x-b||²(5)

在M-step中使用的是求导函数的方法，即找到式(5)的极值点，因为式(5)包含一个二次项，而在E-step中所求出的均值和方差足以满足式(5)的求解条件，此时所求出的是局部最优解，然后多次交替迭代E-step和M-step，直到得到最终的全局最优解，即最佳的模糊核。

如上所述，本发明基于单透镜计算成像系统，将单透镜PSF的估计通过稀疏表示来计算。这只是预处理过程，如果已经求得模糊核，则对同种类型的单透镜，可以直接利用该模糊核，或者将其作为估计模糊核的初始值，以期望估计出更加精确的模糊核。这种方法的输入参数只是由单透镜直接得到的模糊图像，而且计算过程避免了多次的卷积操作，能在一定程度上提高结果图像的边缘细节信息。这种方法在图像处理和相机设计领域都具有非常重要的意义。

附图说明

图1为单透镜成像示意图；

图2为尺度空间中估计PSF的迭代过程示意图；

图3为盲卷积图像复原算法的迭代过程；

图4为稀疏表示示意图；

图5为基于稀疏表示的单透镜计算成像PSF估算方法流程图；

图6为本实验中的单透镜相机以及由相机得到的模糊图像；

图7为最终求得的PSF。

具体实施方式

下面，将结合附图和具体实施例对本发明作进一步说明。

如图5所示，本实施例提供的一种基于稀疏表示的单透镜计算成像PSF估算方法，包括如下步骤：

步骤一：利用单透镜相机获取模糊图像，本实验所做的单透镜相机以及得到的模糊图像如图6所示；

步骤二：将图像去模糊问题转化为联合优化问题，将目标函数中的清晰图像采用过完备字典D和稀疏系数α的乘积表示，并对稀疏系数的稀疏性加以约束。最终目标函数可以表示为：

$\underset{k,\mathrm{\&alpha;}}{min}{||b-k*D\&CenterDot;A||}_2^2+\mathrm{\&lambda;}{||\mathrm{\&alpha;}||}_0---\left(1\right)$

其中，b表示模糊图像，k表示模糊核，D表示过完备字典，A表示稀疏系数矩阵，表示稀疏系数矩阵中的每一列，λ为控制权重的参数。

在具体实现过程中，取权重参数λ＝0.65。

步骤三：目标函数中有三个未知变量：模糊核k，过完备字典D和稀疏系数α，需要通过交替迭代优化算法分别估计出这三个标量。

首先固定k和D来估计α，给定初始的k₀和D₀，通过OMP(OrthogonalMatchingPursuit)优化算法来估计α。

所采用OMP算法的本质思想是：以贪婪迭代的方法选择字典D的列，使得在每次迭代中所选择的列与当前列的冗余向量最大程度地相关，通过反复迭代得到满足稀疏程度的稀疏系数并停止迭代。主要步骤如下：

输入：过完备字典D，清晰图像的图像块对应的向量y，稀疏度l。

输出：稀疏系数α的l-稀疏的逼近α'；

初始化：残差r₀＝y，索引集t＝1；

依次循环步骤1)-5)：

1)找出残差r和过完备字典的列d_i积中最大值所对应的脚标λ，即λ_i＝argmax_j=1...N|<r_t-1，d_i>|

2)更新索引集Λ_t＝Λ_t-1∪{λ_i}，记录找到的过完备字典中重建原子集合

3)由最小二乘得到α_t'＝argmin||y-D_tα_t'||₂；

4)更新残差r_t＝y-D_tα_t'，t＝t+1；

5)判断是否满足t＞l，若满足，则停止迭代；若不满足，则执行1)。

在具体实行过程中，取稀疏度l＝4。

步骤四：固定k和α来更新过完备字典D。模糊图像b不能直接用来构造，首先通过非盲卷积算法对模糊图像进行处理得到对应的潜在清晰图像。

求潜在清晰图像的非盲卷积算法是Levin在2007年所提出的的非盲卷积算法。使用如下概率模型：

$P_k\left(x\right|b)\&Proportional;\exp (-(\frac{1}{{\mathrm{\&eta;}}^2}{\left|C_{k}x-b\right|}^2+\mathrm{\&beta;}{\left|C_{g_x}x\right|}^2+\mathrm{\&beta;}{\left|C_{g_b}x\right|}^2\left)\right)---\left(2\right)$

其中b表示模糊图像，x表示清晰图像，k表示模糊核，g_x＝[1-1]，g_b＝[1-1]^T，C表示卷积操作，对于函数f，C_fx表示f与x的卷积操作，即C_fx≡f*x，β和η是控制权重的参数；则式(2)的最大值即为所求得清晰图像，即x^*＝argmaxP_k(x|b)；

在具体实行过程中，取权重控制参数β＝0.65，η＝0.45。

步骤五：固定α和D估计出模糊核k。清晰图像块I_p可由清晰图像对应的过完备字典中一些原子与稀疏系数线α性组合表示，而模糊图像块可由模糊图像对应的过完备字典中一些原子与稀疏系数线α'性组合表示，因为模糊图像的字典可由清晰图像的字典与模糊核卷积得到，所以不妨假设稀疏系数满足α＝α'。则潜在清晰图像可由下式得到：

$I_p^{(n+1)}=D^{(n+1)}{\mathrm{\&alpha;}}^{(n+1)}---\left(3\right)$

将卷积操作k⁽ⁿ⁺¹⁾I⁽ⁿ⁺¹⁾表示为矩阵相乘的形式K·X，其中K和X分别为k⁽ⁿ⁺¹⁾和I⁽ⁿ⁺¹⁾对应的矩阵形式，在估计模糊核的目标函数中加入L₁范数图像先验：

$\min {||b-Ax||}_2^2+{||x||}_1---\left(4\right)$

通过EM迭代优化算法来求解目标函数公式(4)即能得到所需的模糊核。

EM算法的具体步骤如下：

(1)E-step：令Q(x)＝P(x|b,k)，并计算出Q(x)的均值u和方差C，其中P(x|b,k)表示在已知模糊图像b和模糊核k的条件下，对应的清晰图像为x的概率；均值u就代表给定当前模糊核时潜在的清晰图像，C代表清晰图像周围的方差；

(2)M-step：找到使k*x-b达到最小值的K值，即式(5)达到最小值；

E_Q||k*x-b||²(5)

在M-step中使用的是求导函数的方法，即找到式(5)的极值点，因为式(5)包含一个二次项，而在E-step中所求出的均值和方差足以满足式(5)的求解条件，此时所求出的是局部最优解，然后多次交替迭代E-step和M-step，直到得到最终的全局最优解，即最佳的模糊核。

最终所得到的模糊核如图7所示，所求得的模糊核是空间不变的模糊核。

如上所述，本发明基于单透镜计算成像系统，将单透镜PSF的估计通过稀疏表示来计算。这只是预处理过程，如果已经求得模糊核，则对同种类型的单透镜，可以直接利用该模糊核，或者将其作为估计模糊核的初始值，以期望估计出更加精确的模糊核。这种方法的输入参数只是由单透镜直接得到的模糊图像，而且计算过程避免了多次的卷积操作，能在一定程度上提高结果图像的边缘细节信息。

以上包含了本发明优选实施例的说明，这是为了详细说明本发明的技术特征，并不是想要将发明内容限制在实施例所描述的具体形式中，依据本发明内容主旨进行的其他修改和变型也受本专利保护。本发明内容的主旨是由权利要求书所界定，而非由实施例的具体描述所界定。

Claims (2)

1.一种基于稀疏表示的单透镜计算成像PSF估算方法，其特征在于，包括以下步骤：

步骤一：利用单透镜相机获取模糊图像，模糊图像是在正常光圈大小下，由单透镜相机得到的模糊图像；

步骤二：将模糊图像去模糊问题转化为联合优化问题，将目标函数中的清晰图像采用过完备字典D和稀疏系数α的乘积表示，并对稀疏系数的稀疏性加以约束，最终目标函数表示为：

其中，b表示模糊图像，k表示模糊核，D表示过完备字典，A表示稀疏系数矩阵，表示稀疏系数矩阵中的每一列，λ为控制权重的参数；

步骤三：目标函数中有三个未知变量：模糊核k，过完备字典D和稀疏系数α，需要通过交替迭代优化算法分别估计出这三个标量；

首先固定k和D来估计α，给定初始的k₀和D₀，通过OMP优化算法来估计α；

所采用OMP算法是以贪婪迭代的方法选择字典D的列，使得在每次迭代中所选择的列与当前列的冗余向量最大程度地相关，通过反复迭代得到满足稀疏程度的稀疏系数并停止迭代，具体步骤如下：

输入：过完备字典D，清晰图像的图像块对应的向量y，稀疏度l；

输出：稀疏系数α的l-稀疏的逼近α'；

初始化：残差r₀＝y，索引集t＝1；

依次循环以下步骤1)-5)：

1)找出残差r和过完备字典的列d_i积中最大值所对应的脚标λ，即

λ_i＝argmax_j＝1...N|〈r_t-1，d_i〉|；

2)更新索引集Λ_t＝Λ_t-1∪{λ_i}，记录找到的过完备字典中重建原子集合

3)由最小二乘得到α_t'＝argmin||y-D_tα_t'||₂；

4)更新残差r_t＝y-D_tα_t'，t＝t+1；

5)判断是否满足t＞l，若满足，则停止迭代；若不满足，则执行1)；

步骤四：固定k和α来更新过完备字典D；

模糊图像b不能直接用来构造，首先通过非盲卷积算法对模糊图像进行处理得到对应的潜在清晰图像；

求潜在清晰图像的非盲卷积算法，使用如下概率模型：

其中b表示模糊图像，x表示清晰图像，k表示模糊核，g_x＝[1-1]，g_b＝[1-1]^T，C表示卷积操作，对于函数f，C_fx表示f与x的卷积操作，即C_fx≡f*x，β和η是控制权重的参数；则式(2)的最大值即为所求得清晰图像，即x^*＝argmaxP_k(x|b)；

步骤五：固定α和D估计出模糊核k；

清晰图像块I_p可由清晰图像对应的过完备字典中一些原子与稀疏系数线α性组合表示，而模糊图像块可由模糊图像对应的过完备字典中一些原子与稀疏系数线α'性组合表示，因为模糊图像的字典可由清晰图像的字典与模糊核卷积得到，所以不妨假设稀疏系数满足α＝α'。则潜在清晰图像可由下式得到：

将卷积操作k⁽ⁿ⁺¹⁾*I⁽ⁿ⁺¹⁾表示为矩阵相乘的形式K·X，其中K和X分别为k⁽ⁿ⁺¹⁾和I⁽ⁿ⁺¹⁾对应的矩阵形式，在估计模糊核的目标函数中加入L₁范数图像先验：

通过EM迭代优化算法来求解目标函数公式(4)即能得到所需的模糊核。

2.根据权利要求1所述的基于稀疏表示的单透镜计算成像PSF估算方法，其特征在于，步骤五中，EM迭代优化算法的具体步骤如下：

(1)E-step：令Q(x)＝P(x|b,k)，并计算出Q(x)的均值u和方差C，其中P(x|b,k)表示在已知模糊图像b和模糊核k的条件下，对应的清晰图像为x的概率；均值u就代表给定当前模糊核时潜在的清晰图像，C代表清晰图像周围的方差；

(2)M-step：找到使k*x-b达到最小值的K值，即式(5)达到最小值；

E_Q||k*x-b||²(5)

在M-step中使用的是求导函数的方法，即找到式(5)的极值点，因为式(5)包含一个二次项，而在E-step中所求出的均值和方差足以满足式(5)的求解条件，此时所求出的是局部最优解，然后多次交替迭代E-step和M-step，直到得到最终的全局最优解，即最佳的模糊核。