基于非局部字典学习和双正则项的图像超分辨率重建方法

摘要

本发明公开了一种基于非局部字典学习和双正则项的图像超分辨率重建方法，主要解决现有字典学习方法无法有效利用待超低分辨率图像的先验信息，导致重建的图像不自然的问题。其主要步骤为：（1）得到初始的高分辨率图像（2）训练初始的残差字典集合d0和初始的期望字典集合D0；（3）在初始的高分辨率图像上计算初始的非局部正则权值矩阵W0和局部核回归正则权值矩阵K0；（4）对输入的初始高分辨率图像进行正则优化处理，得到优化图像（5）应用初始的残差字典集合d0和初始的期望字典集合D0重建优化图像得到重建后的图像本发明能够对遥感图像进行重建，并且能够有效的保持图像的边缘及纹理信息，可用于卫星监控和遥感成像。

说明书

技术领域

本发明属于图像处理技术领域，涉及一种图像的超分辨率重建方法，可用于自然图像、医学图像和遥感图像的超分辨率重建。

背景技术

图像超分辨率重建是利用单幅或多幅低分辨率图像重建出一幅高分辨率图像的逆问题。一般来说，图像超分辨率重建的方法分为三类：基于插值的方法，基于模型重建的方法以及基于学习的方法。其中，基于插值的方法包括经典的双线性插值法，三次样条插值法等。这类方法简单快速，但容易导致边缘模糊，达不到好的重建效果；基于模型重建的方法包括迭代反投影法，最大后验概率法，凸集投影法等，这类方法多是针对多幅低分辨率图像的重建，虽然能够产生较好的重建效果，但是重建模型的参数很难估计，提高图像分辨率的能力有限。随着机器学习技术的发展，大量的基于学习的超分辨率重建方法被提出。该类方法学习已知的高分辨率图像和低分辨率图像之间的关系，并将这种关系作为一种先验知识来指导低分辨率图像的超分辨重建。近年来比较流行的是Yang等人提出的基于稀疏表示的字典学习超分辨率重建方法，他们假设多数自然图像具有相似的结构表示基元，例如边缘、纹理等。自此，许多基于字典学习的重建方法得到了发展。

多数基于字典学习的重建方法虽然可以有效利用外来高分辨样本图像的先验信息，但该类方法忽略了待超低分辨率图像自身的先验知识，导致重建后的图像存在伪像和失真。

发明内容

本发明的目的在于针对基于字典学习的图像超分辨重建方法的缺点，提出一种基于非局部字典学习和双正则项的图像超分辨率重建方法，以利用待超低分辨率图像冗余的自相似性信息，引入改进的非局部正则项和改进的局部核回归正则项，提高图像重建的质量。

为实现上述目的，本发明的技术方案是首先从待超低分辨率图像中构造自适应的双字典，然后根据图像的降质模型设计目标优化函数，通过多次迭代得到高分辨率图像。具体步骤包括如下：

(1)输入低分辨率图像I_l，对该低分辨率图像I_l作双线性插值处理，得到初始的高分辨率图像

(2)对初始的高分辨率图像进行自适应的聚类字典训练，得到R个聚类中心C_center＝{C_i,i＝1,2,...,R}，初始的期望字典集合D₀和初始的残差字典集合d₀；

(3)在初始的高分辨率图像上计算初始的非局部正则权值矩阵W₀和初始的局部核正则权值矩阵K₀；

(4)设迭代次数为n，n＝0,1,2,...,299，初始迭代n＝0，重建超分辨率图像：

4a)输入初始迭代图像根据非局部正则权值矩阵W₀和初始的局部核正则权值矩阵K₀，对该初始迭代图像进行正则优化，得到优化图像

4b)应用初始的期望字典集D₀，初始的残差字典集d₀，R个聚类中心C_center＝{C_i,i＝1,2,...,R}和输入的初始迭代图像对优化图像进行重建，得到重建后的图像

4c)设误差精度ε＝2e-9，N为重建图像中总的图像块个数，根据所述参数ε,N,判断迭代是否终止，若则终止迭代，输出重建图像否则将重建图像作为新的迭代输入图像，迭代次数加1,即n＝n+1，并设指示参数p1＝120；

4d)根据参数p1,n判断是否更新初始的期望字典集合D₀，初始的残差典集合d₀，初始的非局部正则权值矩阵W₀以及初始的局部核正则权值矩阵K₀，若n为p1的整数倍，则将初始的期望字典集合D₀更新为D_n，将初始的残差字典集合d₀更新为d_n，将初始的非局部正则权值矩阵W₀更新为W_n，将初始的局部核正则权值矩阵K₀更新为K_n；否则不更新，返回步骤4a)。

本发明与现有的技术相比具有以下优点：

1)本发明充分利用了待超低分辨率图像自身的先验信息，克服了传统的基于字典学习的重建方法无法有效补充低分辨率图像缺失信息的缺点，能够更加高效和定向的补充低分辨率图像缺失的信息；

2)本发明在重建过程中引入改进的非局部双边正则约束项，在进行相似块匹配时，综合考虑了图像块的特征距离、亮度距离以及几何位置距离，克服了传统的非局部正则技术仅利用亮度距离衡量图像块相似性的缺点，有效提高了相似块匹配的精度；

3)本发明在重建过程中引入改进的局部核回归正则约束项，利用高频特征图求取权重，充分利用了图像的局部特征先验信息，克服了传统的核回归正则技术仅利用亮度信息计算权重的缺点，能够求得更加精确的局部权重系数，有效提高了重建图像的质量。

实验表明，本发明适用于遥感图像的超分辨率重建，较好的保持了图像的边缘以及纹理细节信息。

附图说明

图1是本发明的总流程图；

图2是本发明用到的两幅测试图像的原始高分辨率图；

图3是用本发明和现有两种方法对测试图像IM1的重建效果图；

图4是用本发明和现有两种方法对测试图像IM2的重建效果图。

具体实施方式

参照附图1，本发明的具体步骤包括：

步骤1.对初始的高分辨率图像进行自适应的聚类字典训练，得到R个聚类中心C_center＝{C_i,i＝1,2,...,R}，初始的期望字典集合D₀和初始的残差字典集合d₀。

1a)提取初始高分辨率图像的高频特征，得到高频特征图G；

1b)分别在初始高分辨率图像和高频特征图像G上取7×7的块，将该初始高分辨率图像上获取的所有图像块以列向量的形式依次排列，形成图像块集合将该高频特征图像G上获取的所有特征块以列向量的形式依次排列，形成特征块集合

1c)利用K均值聚类的方法对特征块集合进行聚类，得到R个聚类中心C_center＝{C_i,i＝1,2,...,R}；

1d)根据R个聚类中心C_center＝{C_i,i＝1,2,...,R}和特征块集合将图像块集合分为R个图像块子集：

1e)设图像块u_i是初始高分辨率图像上的第i个图像块，设该图像块u_i属于第X类图像块子集X∈{1,2,...,R}，利用非局部均值滤波的方法，获得图像块u_i的25个相似图像块u′₁,u'₂,...,u'₂₅，同时获得该25个相似图像块u′₁,u'₂,...,u'₂₅和图像块u_i的权重系数w₁,w₂,...,w₂₅；

1f)根据25个权重系数w₁,w₂,...,w₂₅和25个相似图像块u′₁,u'₂,...,u'₂₅，计算图像块u_i的估计图像块$>\overline{u_i}:\overline{u_i}=\underset{m=1}{\overset{25}{\Sigma }}{w}_m{u}_m^{\text{'}};$>

1g)根据25个相似图像块u′₁,u'₂,...,u'₂₅和估计图像块计算该25个相似图像块与估计图像块的残差图像块Δ₁,Δ₂,...,Δ₂₅，其中,$>{\Delta }_2=u_2^{\text{'}}-\overline{u_i},...,{\Delta }_{25}=u_{25}^{\text{'}}-\overline{u_i};$>

1h)重复步骤1e)到步骤1g)，获得图像块子集的所有相似图像块和图像块子集的所有残差图像块；将所有相似图像块顺次排列，形成相似块集合U，将所有残差图像块顺次排列，形成残差块集合Δ；

1i)利用主成分分析的方法对相似块集合U进行训练，得到图像块子集的期望字典利用主成分分析的方法对残差块集合Δ进行训练，得到图像块子集的残差字典

1j)重复步骤1e)到步骤1i)，获得R个期望字典和R个残差字典$>d_0^1,d_0^2,...,d_0^R;$>

1k)将这R个期望字典顺序排列，构成期望字典集合将这R个残差字典顺序排列，构成残差字典集合$>d_0=\{d_0^1,d_0^2,...,d_0^R\}.$>

上述步骤1c)中所用的聚类方法包括K均值聚类方法，模糊C均值聚类方法等，本实例采用但不局限于为K均值聚类方法，K均值聚类方法的程序调用自“WeishengDong,Guangming Shi,Lei Zhang,and Xiaolin Wu,Image Deblurring andSuper-Resolution by Adaptive Sparse Domain Selection and Adaptive Regularization[J].IEEE Transactions on image processing,vol.20,no.7,pp.1838-1857,July2011.”。

上述步骤1i)中所用的训练字典的方法包括K-奇异值分解法，主成分分析法等，本实例采用但不局限于为主成分分析法，主成分分析法的程序调用自“Weisheng Dong,Guangming Shi,Lei Zhang,and Xiaolin Wu,Image Deblurring and Super-Resolution byAdaptive Sparse Domain Selection and Adaptive Regularization[J].IEEE Transactionson image processing,vol.20,no.7,pp.1838-1857,July2011.”。

步骤2.在初始高分辨率图像上计算初始的非局部正则权值矩阵W₀。

2a)对初始高分辨率图像提取高频特征，得到高频特征图G；

2b)设图像块u_i为初始高分辨率图像的第i个图像块，设像素点c_i为图像块u_i的中心像素点；设图像块g_i为高频特征图G的第i个图像块，设像素点d_i为图像块g_i的中心像素点；

2c)在初始高分辨率图像上，以像素点c_i为中心取12×12大小的相似窗Q₁，该相似窗Q₁内排除像素点c_i，有143个像素点c′₁,c'₂,...,c′₁₄₃；在高频特征图G上，以像素点d_i为中心取12×12大小的相似窗Q₂，该相似窗Q₂内排除像素点d_i，有143个像素点d′₁,d'₂,...,d′₁₄₃；

2d)设图像块u_j是以相似窗Q₁内第j个像素点c'_j为中心的图像块，设图像块g_j是以相似窗Q₂内第j个像素点d'_j为中心的图像块，j＝1,2,...,143，设像素点c_i在初始高分辨率图像上的二维坐标为Z_i(x,y)，设像素点c'_j在初始高分辨率图像上的二维坐标为Z_j(x,y)；

2e)根据坐标Z_i(x,y)和坐标Z_j(x,y)，计算所述像素点c'_j和所述像素点c_i的亮度值欧氏距离几何位置欧氏距离特征欧式距离和总距离

$>d_{{1c}_j^{\text{'}}}={\left|\right|u_i-u_j\left|\right|}_2^2/49,$>

$>d_{{2c}_j^{\text{'}}}={\left|\right|Z_i(x,y)-Z_j(x,y)\left|\right|}_2^2/49,$>

$>d_{{3c}_j^{\text{'}}}={\left|\right|g_i-g_j\left|\right|}_2^2/49,$>

$>d_{c_j^{\text{'}}}=d_{{1c}_j^{\text{'}}}+d_{{2c}_j^{\text{'}}}+d_{{3c}_j^{\text{'}}};$>

2f)按照步骤2e)中所述方法，获得与像素点c_i距离最近的25个像素点{c'_pm,m＝1,2,...,25}，该25个像素点与像素点c_i的亮度值欧氏距离为该25个像素点与像素点c_i的几何位置欧氏距离为该25个像素点与像素点c_i的特征欧氏距离为

2g)设亮度平滑控制参数h₁＝300，设几何平滑控制参数h₂＝100，设特征平滑控制参数h₃＝100，根据亮度平滑控制参数h₁，几何平滑控制参数h₂和特征平滑控制参数h₃，将上述亮度值欧氏距离几何位置欧氏距离特征欧氏距离代入高斯核函数exp(.)，计算像素点c'_pm与像素点c_i的权重系数w(c_i,c'_pm)：$>w(c_i,c_{\mathrm{pm}}^{\text{'}})=\exp (-\frac{d_{{1c}_{\mathrm{pm}}^{\text{'}}}}{h_1})\exp (-\frac{d_{{2c}_{\mathrm{pm}}^{\text{'}}}}{h_2})\exp (-\frac{d_{3c_{\mathrm{pm}}^{\text{'}}}}{h_3});$>

2h)重复步骤2a)到步骤2g)，获得初始高分辨率图像中的全部像素点的权重系数，将所有权重系数按照坐标顺序排列，形成初始的正则权值矩阵W₀：$>W_0=\left(\begin{array}{c}w(c_i,c_{\mathrm{pm}}^{\text{'}}),c_i\in I_h^0,c_{\mathrm{pm}}^{\text{'}}\in I_h^0,m=\mathrm{1,2},...,25\\ 0\end{array}\right).$>

步骤3.在初始高分辨率图像上计算初始的局部核正则权值矩阵K₀。

3a)对初始高分辨率图像提取高频特征，得到高频特征图G；

3b)设像素点a_i为高频特征图G中的第i个像素点，设该像素点a_i的坐标向量为t，设该像素点a_i的亮度值为F(t)，以该像素点a_i为中心取7×7的窗口；

3c)该窗口内共有49个像素点，将该49个像素点顺次排列，形成列向量b＝[b₁,b₂,...,b₄₉]^T，设第j个像素点b_j的坐标向量为t_j，设像素点b_j的亮度值为q_j；

3d)设矩阵S_j为像素点b_j的局部协方差矩阵，设全局光滑参数h＝0.005，设局部采样密度μ＝0.5，根据上述参数，计算像素点b_j的旋转矩阵

3e)根据上述参数，计算像素点b_j的控制核矩阵$>W_{H_{b_j}^{\mathrm{steer}}}(t_i-t)=\frac{\sqrt{\det \left(S_j\right)}}{2\pi h^2{\mu }^2}\exp \{-\frac{{(t_j-t)}^T{S}_j(t_j-t)}{2h^2{\mu }^2}\},$>j＝1,...49，其中，exp(.)为高斯核函数，det(.)为求矩阵行列式的函数。

3f)根据上述参数，获得49个像素点{b_j,j＝1,2,...,49}的总控制核矩阵W_t：$>W_t=\mathrm{diag}[W_{H_{b_1}^{\mathrm{steer}}}(t_1-t),W_{H_{b_2}^{\mathrm{steer}}}(t_2-t)...,W_{H_{b_{49}}^{\mathrm{steer}}}(t_{49}-t)],$>其中，diag(.)是取对角操作；

3g)设总坐标矩阵$>E_t=\left(\begin{array}{cc}1{(t_1-t)}^T{\mathrm{ltr}}^T\left\{(t-t_1){(t_1-t)}^T\right\}& ···\\ 1{(t_2-t)}^T{\mathrm{ltr}}^T\left\{(t_2-t){(t_2-t)}^T\right\}& ···\\ ·& ·\\ ·& ·\\ ·& ·\\ 1{(t_{49}-t)}^T{\mathrm{ltr}}^T\left\{(t_{49}-t){(t_{49}-t)}^T\right\}& ···\end{array}\right),$>根据总坐标矩阵E_t和总控制核矩阵W_t，利用最小二乘的方法，计算亮度值F(t)的估计值其中，ltr(.)表示将一个矩阵的下三角元素排成一个列向量的操作，e₁表示第一个元素为1，其它元素都为0的列向量；

3h)设可变参数N'是泰勒展开阶数，将改写为线性加权的形式：$>\hat{F}\left(t\right)=\underset{j=1}{\overset{49}{\Sigma }}{W}_{b_j}(t;H_{b_j}^{\mathrm{steer}},N^{\text{'}})b_j,$>其中矩阵表示窗口内第j个像素点b_j和像素点a_i的权重系数：

3i)设可变参数N'＝2，重复步骤3a)到步骤3h)，获得初始高分辨率图像中的全部像素点的权重系数，将所有权重系数按照坐标顺序排列，形成初始的局部核函数正则权值矩阵K₀：

$>K_0=\left(\begin{array}{c}{W}_{b_j}(t;H_{b_j}^{\mathrm{steer}},2),b_j\in b\\ 0,b_j\notin b\end{array}\right).$>

上述步骤3e)中所用的核包括双边核，控制核等，本实例采用但不局限于为控制核，控制核的程序调用自“H.Takeda,S.Farsiu,and P.Milanfar,Kernel regression forimage processing and reconstruction.IEEE Transactions on image processing,vol.16,no.2,pp.349–366,2007.”。

步骤4.输入初始迭代图像根据非局部正则权值矩阵W₀和初始的局部核正则权值矩阵K₀，对该初始迭代图像进行正则优化，得到优化图像

4a)设矩阵C为采样矩阵，设低分辨图像为I_l，设未知高分辨率图像为I_h，设正则参数β₁＝0.14,β₂＝0.46，设矩阵E为单位矩阵；

4b)根据非局部正则权值矩阵W₀，初始的局部核正则权值矩阵K₀和步骤4a)中的参数，得到正则优化公式为：$>I_h=\underset{I_h}{\arg \min }({\left|\right|I_l-{\mathrm{CI}}_h\left|\right|}_2^2+{\beta }_1{\left|\right|(E-W_0)I_h\left|\right|}_2^2+{\beta }_2{\left|\right|(E-K_0)I_h\left|\right|}_2^2),$>其中是l₂范数求值函数；

4c)设常量矩阵$>{\tilde{I}}_l=\left(\begin{array}{c}{I}_l\\ 0\\ 0\end{array}\right),$>设操作矩阵$>\tilde{K}=\left(\begin{array}{c}C\\ {\beta }_1(E-W_0)\\ {\beta }_2(E-K_0)\end{array}\right),$>根据常量矩阵操作矩阵将步骤4b)中的正则优化公式简化为：

4d)根据输入图像和公式利用梯度下降法计算优化图像

$>I_h^{1/2}=I_h^0+K^T({\tilde{I}}_l-{\mathrm{KI}}_h^0)=I_h^0+C^T{I}_l-C^T{\mathrm{CI}}_h^0$>

$>-{\beta }_1^2{(E-W_0)}^T(E-W_0)I_h^0-{\beta }_2^2{(E-K_0)}^T(E-K_0)I_h^0.$>

该步骤4d)中l₂范数公式所用的求解方法包括梯度下降法，软阈值收缩法等，本实例采用但不局限于为梯度下降法，梯度下降法程序调用自“Weisheng Dong,Guangming Shi,Lei Zhang,and Xiaolin Wu,Image Deblurring and Super-Resolution byAdaptive Sparse Domain Selection and Adaptive Regularization[J].IEEE Transactionson image processing,vol.20,no.7,pp.1838-1857,July2011.”。

步骤5.应用初始的期望字典集D₀，初始的残差字典集d₀，R个聚类中心C_center＝{C_i,i＝1,2,...,R}和输入的初始迭代图像对优化图像进行重建，得到重建后的图像

5a)提取优化图像的高频特征，得到高频特征图像G；

5b)分别对优化图像和高频特征图像G取7×7的块，该优化图像上获取的所有图像块以列向量的形式依次排列，形成图像块集合该高频特征图G上获取的所有特征块以列向量的形式依次排列，形成特征块集合

5c)根据R个聚类中心C_center＝{C_i,i＝1,2,...,R}和特征块集合将图像块集合分为R个图像块子集：

5d)设图像块u_i是初始高分辨率图像上的第i个图像块，设该图像块u_i属于第X类图像块子集X∈{1,2,...,R}，根据类别X，在期望字典集合D₀中选定字典在残差字典集合d₀中选定字典

5e)利用非局部均值技术，获得图像块u_i的估计块和残差块

5f)利用估计块残差块Δ_i，期望字典残差字典计算该图像块u_i在期望字典上的稀疏表示系数α_D和图像块u_i在残差字典上的稀疏表示系数α_d：

$>{\alpha }_D={\left(D_0^X\right)}^T\overline{u_i},$>

$>{\alpha }_d={\left(d_0^X\right)}^T{\Delta }_i;$>

5g)设初始高分辨率图像中的第i个图像块为u′_i，按照步骤5a)到5f)中所述方法，设该图像块u′_i的期望字典稀疏表示系数为α′_D，设该图像块u′_i的残差字典稀疏表示系数α′_d；

5h)设阈值参数ε2＝0.35，根据阈值参数ε2，期望字典稀疏表示系数α′_D，残差字典稀疏表示系数α′_d，利用迭代加权的稀疏表示法，计算稀疏表示正则参数λ_D和λ_d：

$>{\lambda }_D=\frac{1}{{\alpha }_D^{\text{'}}+ϵ2},$>

$>{\lambda }_d=\frac{1}{{\alpha }_d^{\text{'}}+ϵ2};$>

5i)根据上述参数α_D，α_d，λ_D、λ_d，利用软阈值收缩法，将图像块u_i的期望字典稀疏表示系数α_D更新为将残差字典稀疏表示系数α_d更新为

$>{\alpha }_D^1=\mathrm{sign}\left({\alpha }_D\right)\times \max (\mathrm{abs}\left({\alpha }_D\right)-{\lambda }_D,0),$>

$>{\alpha }_d^1=\mathrm{sign}\left({\alpha }_d\right)\times \max (\mathrm{abs}\left({\alpha }_d\right)-{\lambda }_d,0),$>

其中，sign(.)表示符号函数，abs(.)表示取绝对值函数；

5j)根据上述参数期望字典残差字典计算图像块u_i重建后的图像块

$>\overline{u_i}=D_0^X{\alpha }_D^1+d_0^X{\alpha }_d^1;$>

5k)重复步骤5d）到步骤5j)，重建优化图像中所有的图像块，得到重建后的图像

上述步骤5h)中外字典正则参数λ_D和内字典正则参数λ_d所用的求解方法包括迭代加权的稀疏表示法，全局加权的稀疏表示法等。本实例采用但不局限于为迭代加权的稀疏表示法，迭代加权的稀疏表示法的程序调用自“E.J.Candes,M.B.Wakin,and S.Boyd,Enhancing Sparsity by Reweighted l₁Minimization[J],Journal of Fourier Analysisand Applications,vol.14,pp.877–905,2008.”。

上述步骤5i)中稀疏表示系数更新所用的求解方法包括梯度下降法，软阈值收缩法等。本实例采用但不局限于为软阈值收缩法，软阈值收缩法的程序调用自“I.Daubechies,M.Defriese,and C.DeMol,An iterative thresholding algorithm forlinearinverse problems with a sparsity constraint[J],Commun.Pure Appl.Math.,vol.57,pp.1413～1457,2004.”。

本发明的效果可以通过以下实验进一步说明：

1)实验条件

本实验采用的测试图像为IM1和IM2，如图2所示，采用MATLAB7.10.0软件作为仿真工具，计算机配置为Intel(R)Core(TM)2/2.33G/1.99G，本实验中，图像降质过程中的下采样因子设为3。

2)实验内容

实验1，利用Bicubic插值算法、SCDL算法和本发明方法，对IM1图像进行重建，结果如图3所示，其中图3(a)为IM1图像的高分辨率图像，图3(b)为输入的低分辨率图像，图3(c)为Bicubic插值算法的重建结果，图3(d)为SCDL算法的重建结果，图3(e)为本发明的重建结果；

实验2，利用Bicubic插值算法、SCDL算法和本发明方法，对IM2图像进行重建，结果如图3所示，其中图4(a)为IM2图像的高分辨率图像，图4(b)为输入的低分辨率图像，图4(c)为Bicubic插值算法的重建结果，图4(d)为SCDL算法的重建结果，图4(e)为本发明的重建结果。

实验结果分析

从图3—图4可以看出，本发明的重建图像的视觉效果要优于其它两种重建方法，对于输入的两幅测试图像IM1和IM2，图像的边缘和纹理细节信息都保持的比较好。