CN104867119A - 基于低秩矩阵重建的结构性缺失图像填充方法 - Google Patents

Abstract

本发明属于计算机视觉领域，为实现对像素结构性缺失的图片进行填充。本发明采取的技术方案是，基于低秩矩阵重建的结构性缺失图像填充方法，包括下列步骤：1)将图像看作矩阵，则原图像用矩阵A表示，解决像素结构性缺失的图像填充问题为求解优化方程；2)训练字典Φ；3)使用交替方向法ADM将序列转换成序列进行求解，然后按照步骤进行迭代求解得到最终的结果。本发明主要应用于计算机图像处理。

Description

基于低秩矩阵重建的结构性缺失图像填充方法

技术领域

本发明属于计算机视觉领域。特别涉及基于低秩矩阵重建的结构性缺失图像填充方法。

背景技术

矩阵重建问题包括矩阵填充和矩阵恢复，自提出以来就一直倍受关注，且具有很强的生命力。尤其在大数据的时代背景下更加成为了数学和计算机领域的研究热点。近年来，关于解决矩阵重建问题的算法已经有很多研究成果。这些算法主要是利用迭代的奇异值分解的方法来逼近原来模型的解。如SVT(奇异值阈值)算法，APG(加速近邻梯度)算法，ALM(增广拉格朗日乘子)算法等。在现有的算法中，解决矩阵填充问题时，SVT算法在编程实现过程中需要的内存很小，所以比较适合大规模矩阵的情况使用，但是SVT只适用于矩阵秩非常低的情况。APG是从FISTA(快速迭代收缩阈值)算法变化得到的，其收敛速度比SVT等算法快很多。在解决矩阵恢复问题时，SVT和APG算法依然可以有着不错的表现，但是从收敛速度来看，这些算法都是次线性的。相比之下ALM算法则有着较快的收敛速度。

图像填充是计算机视觉领域中一类很重要的问题，目前已有的矩阵重建算法如SVT算法、APG算法和ALM算法等都是利用图像低秩特性将缺失的像素填充好。但是当图像中像素为结构性缺失如线状缺失，甚至是整行缺失时，已有算法则不能解决这种图像填充问题。因为它们所作用的矩阵都有个共同的特点，即观测矩阵的元素空缺位置是稀疏且随机的，这对于诸如图像去噪等应用具有一定效果，但对于空缺整行元素的情况则完全不起作用。因为大量整行元素空缺的矩阵恢复仅仅应用低秩条件进行限制是无法求解的。随着大数据时代的到来，信息量急剧增长，图像在传输过程中很可能会遭遇结构性缺失的污染。因此设计出一种能够解决像素结构性缺失的图像填充的算法是很有必要的。

发明内容

本发明意在弥补现有技术的不足，即实现对像素结构性缺失的图片进行填充。本发明采取的技术方案是，基于低秩矩阵重建的结构性缺失图像填充方法，将矩阵重建理论与稀疏表示理论相结合，在传统的矩阵重建模型基础上引入字典学习模型，从而解决已有技术无法处理的问题。本发明包括下列步骤：

1)将图像看作矩阵，则原图像用矩阵A表示，解决像素结构性缺失的图像填充问题为求解如下优化方程：

min||A||_*+λ||B||₁

              (1)

约束条件A＝ΦB,A+E＝D,P_Ω(E)＝0

其中||A||_*表示矩阵A的核范数，||·||₁表示矩阵的一范数，Ω是观测空间，p_Ω(·)是投影算子，表示变量投射到空间域Ω内的值，λ为权重系数，约束条件中的Φ为训练好的字典，B是字典对应的系数，D为已知的降质观测矩阵即缺失的受损图像，E代表受损图像中缺失的像素；

求解此方程时，本发明采用增广拉格朗日乘子法，方程如下：

$\begin{array}{c}L(A,B,E,Y_1,Y_2)={\left|\right|A\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1,A-\mathrm{\&Phi;B}>\\ +<Y_2,D-A-E>+\frac{{\mathrm{\&mu;}}_1}{2}{\left|\right|A-\mathrm{\&Phi;B}\left|\right|}_F^2+\frac{{\mathrm{\&mu;}}_2}{2}{\left|\right|D-A-E\left|\right|}_F^2\end{array}---\left(2\right)$

L(A,B,E,Y₁,Y₂)即为增广拉格朗日函数，其中μ₁和μ₂是惩罚因子；Y₁、Y₂均为拉格朗日乘子矩阵；<·,·>表示两个矩阵的内积；||·||_F表示矩阵的弗罗贝尼乌斯Frobenius范数；

(2)式的迭代求解法方程如下：

$\left\{\begin{array}{c}(A^{k+1},B^{k+1},E^{k+1})=\arg \underset{A,B,E}{\min }\left\{L(A,B,E,Y_1^k,Y_2^k)\right\}\\ Y_1^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ Y_2^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\\ \end{array}\right.---\left(3\right)$

上式中的表示使目标函数取最小值时的变量A,B,E的值，ρ₁和ρ₂为倍数因子，k是迭代次数；

2)训练字典Φ：在高质量的图像数据集上使用在线学习算法训练出字典Φ；

3)使用交替方向法ADM将序列(3)转换成如下序列进行求解：

$\left\{\begin{array}{c}{B}^{k+1}=\arg \underset{B}{\min }\left\{L(A^k,B,E^k,Y_1^k,Y_2^k)\right\}\\ A^{k+1}=\arg \underset{A}{\min }\left\{L(A,B^{k+1},E^k,Y_1^k,Y_2^k)\right\}\\ E^{k+1}=\arg \underset{E}{\min }\left\{L(A^{k+1},B^{k+1},E,Y_1^k,Y_2^k)\right\}\\ Y_1^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ Y_2^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\\ \end{array}\right.---\left(4\right)$

然后按照步骤4)、5)、6)的方法进行迭代求解得到最终的结果；

4)求解B^k+1：使用加速近邻梯度算法求得B^k+1；

去掉式子(4)中求解B的目标函数里与B无关的项，得到如下方程:

$B^{k+1}=\arg \underset{B}{\min }\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2---\left(5\right)$

使用泰勒展开的方法，构造出一个二阶函数来逼近上式，然后针对这个二阶的函数来求解原方程，令 $f\left(B\right)=<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2,$ 最终可以解得：

$B^{k+1}=B_{j+1}^k=\mathrm{soft}(U_{j+1},\mathrm{\&lambda;}/L_f)---\left(6\right)$

其中的soft(x,α)为收缩算子， 代表函数f对Z_j+1的Frechet梯度，这里f即为f(B)且Z_j+1代替B，L_f是一个常数，取值为变量Z_j的更新规则如下：

$\left\{\begin{array}{c}{Z}_{j+1}=B_{j+1}^k+\frac{t_{j+1}-1}{t_{j+1}}(B_{j+1}^k-B_j^k)\\ t_{j+1}=(1+\sqrt{4t_j^2+1})/2\end{array}\right.---\left(7\right)$

t_j是一组序列，j是变量迭代次数；

5)求解A^k+1：使用SVT算法求解A^k+1；

去掉式子(4)中求解A的目标函数里与A无关的项，并且通过配方得到：

$A^{k+1}=\arg \underset{A}{\min }{\left|\right|A\left|\right|}_{*}+\frac{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}{2}{\left|\right|A-W^k\left|\right|}_F^2---\left(8\right)$

其中：

$W^k=(Y_2^k-Y_1^k+{\mathrm{\&mu;}}_1^k\mathrm{\&Phi;}{B}^{k+1}-{\mathrm{\&mu;}}_2^k{E}^k+{\mathrm{\&mu;}}_2^{k}D)/({\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k)---\left(9\right)$

对于式子(9)使用奇异值阈值法解得：

$A^{k+1}=U^k\mathrm{soft}({\mathrm{\&Sigma;}}^k,\frac{1}{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}){\left(V^k\right)}^T---\left(10\right)$

其中U^k，V^k分别是W^k的左奇异矩阵和右奇异矩阵；

6)求解E^k+1：将E^k+1分为两部分求解，对于空间域Ω内部，E的值为0，对于在空间域Ω以外的部分即中，使用一阶求导来求解，将两部分合起来即为E的最终解：

$E^{k+1}=P_{\mathrm{\&Omega;}}\left(0\right)+P_{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}(\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}+D-A^{k+1})---\left(11\right)$

7)重复上述步骤4)、5)、6)直到算法收敛，这时迭代的结果A^k+1、B^k+1和E^k+1就是原问题的最终解A、B和E。

本发明的技术特点及效果：

本发明方法针对像素结构性缺失的图像填充问题,通过引入字典将图像阵稀疏表示,实现了对结构性缺失的图像填充问题的求解。本发明具有以下特点：

1、运用了ALM、APG、SVT等算法求解子问题，整合了已有算法的优点。

2、使用了列字典，与传统的块字典相比更加高效。

3、将矩阵重建理论和稀疏表示理论相结合，在传统的矩阵重建模型中引入字典学习，这样使得整行元素缺失的矩阵也能够被重建。

附图说明

图1是算法流程图；

图2是原始的Lena图；

图3是缺失后的受损图；黑色为缺失像素(左上：10％缺失；右上：20％缺失；左下：30％缺失；右下：50％缺失)。

图4是填充后的结果图；(左上：10％缺失填充结果，PSNR＝42.09；右上：20％缺失填充结果PSNR＝37.98；左下：30％缺失填充结果PSNR＝35.05；右下：50％缺失填充结果PSNR＝31.20)。

具体实施方式

是在传统矩阵重建模型的基础上引入字典学习模型，使得能够重建出结构性缺失的低秩矩阵从而得到填充后的图像，即基于低秩矩阵重建的结构性缺失图像填充方法，从而解决已有技术无法处理的问题。下面结合附图和实施例对本发明作详细说明。

1)将图像看作矩阵，则原图像可以用矩阵A表示，解决像素结构性缺失的图像填充问题为求解如下优化方程：

min||A||_*+λ||B||₁

           (1)

约束条件A＝ΦB,A+E＝D,P_Ω(E)＝0

||A||_*表示矩阵A的核范数。||·||₁表示矩阵的一范数。Ω是观测空间，p_Ω(·)是投影算子，表示变量投射到空间域Ω内的值。λ为权重系数。约束条件中的Φ为训练好的字典，B是字典对应的系数矩阵。D为已知的降质观测矩阵即缺失的受损图像。E代表受损图像中缺失的像素。实验中使用Lena图像作为测试数据。

11)求解此方程时，本发明采用增广拉格朗日乘子法，方程如下：

$\begin{array}{c}L(A,B,E,Y_1,Y_2)={\left|\right|A\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1,A-\mathrm{\&Phi;B}>\\ +<Y_2,D-A-E>+\frac{{\mathrm{\&mu;}}_1}{2}{\left|\right|A-\mathrm{\&Phi;B}\left|\right|}_F^2+\frac{{\mathrm{\&mu;}}_2}{2}{\left|\right|D-A-E\left|\right|}_F^2\end{array}---\left(2\right)$

L(A,B,E,Y₁,Y₂)即为增广拉格朗日函数。其中的μ₁和μ₂是惩罚因子；Y₁、Y₂均为拉格朗日乘子矩阵；<·,·>表示两个矩阵的内积，如<M,N>＝trace(M^TN)，M^T表示矩阵M的转置，trace(·)表示求矩阵的迹；||·||_F表示矩阵的Frobenius(弗罗贝尼乌斯)范数。

12)用迭代求解法求解(2)式，迭代的方程如下：

$\left\{\begin{array}{c}(A^{k+1},B^{k+1},E^{k+1})=\arg \underset{A,B,E}{\min }\left\{L(A,B,E,Y_1^k,Y_2^k)\right\}\\ Y_1^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ Y_2^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\\ \end{array}\right.---\left(3\right)$

上式中的表示使目标函数取最小值时的变量A,B,E的值。ρ₁和ρ₂为倍数因子。k是迭代次数。

2)使用Online Learning字典学习方法训练出字典Φ。

21)构造字典Φ使得矩阵A能够由字典稀疏表示，即满足A＝ΦB，其中B是系数矩阵且是稀疏的。本发明使用Online Learning算法在Kodak图像集上训练出字典Φ。

22)在训练字典时，我们设定字典相关参数如下：待重建矩阵A的行数与字典Φ中元素的维数m相等，即A的行数与Φ的行数均为m。训练出的字典Φ是一个过完备的字典，即字典的列数必须大于其行数。

3)本发明使用交替方向法来求解序列(3)，即将序列(3)转换成下面的序列进行求解：

$\left\{\begin{array}{c}{B}^{k+1}=\arg \underset{B}{\min }\left\{L(A^k,B,E^k,Y_1^k,Y_2^k)\right\}\\ A^{k+1}=\arg \underset{A}{\min }\left\{L(A,B^{k+1},E^k,Y_1^k,Y_2^k)\right\}\\ E^{k+1}=\arg \underset{E}{\min }\left\{L(A^{k+1},B^{k+1},E,Y_1^k,Y_2^k)\right\}\\ Y_1^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ Y_2^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\\ \end{array}\right.---\left(4\right)$

设定好各参数初值，然后按照步骤4)、5)、6)的方法进行迭代求解最终结果。

4)使用加速近邻梯度算法求解第一个未知元B^k+1。

41)去掉式子(4)中求解B的目标函数里与B无关的项后得到如下方程：

$B^{k+1}=\arg \underset{B}{\min }\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2---\left(5\right)$

使用泰勒展开的方法，构造出一个二阶函数来逼近上式，然后针对这个二阶的函数来求解原方程。

令 $f\left(B\right)=<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2,F\left(B\right)=\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2$ 再引入变量Z，定义如下函数：

$Q(B,Z)=\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+f\left(Z\right)+<\&dtri;f\left(Z\right),B-Z>+\frac{L_f}{2}{\left|\right|B-Z\left|\right|}_F^2---\left(6\right)$

其中，代表函数f对Z的Frechet梯度，这里f即为f(B)且Z_j+1代替B，L_f是一个常数，取值为用来保证对于所有的Z，都有F(Z)≤Q(B,Z)。

42)经过上步转化，式子(5)转化成了求解Q(B,Z_j)的最小值来求解，再通过配方得到如下形式：

$\begin{array}{c}{B}^{k+1}=B_{j+1}^k=\arg \underset{B}{\min }Q(B,Z)\\ =\arg \underset{B}{\min }\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+\frac{L_f}{2}{\left|\right|B-U_{j+1}\left|\right|}_F^2\end{array}---\left(7\right)$

其中变量Z_j的更新规则如下：

$\left\{\begin{array}{c}{Z}_{j+1}=B_{j+1}^k+\frac{t_{j+1}-1}{t_{j+1}}(B_{j+1}^k-B_j^k)\\ t_{j+1}=(1+\sqrt{4t_j^2+1})/2\end{array}\right.---\left(8\right)$

t_j是一组序列，j是变量迭代次数。使用收缩算子解得：

$B^{k+1}=B_{j+1}^k=\mathrm{soft}(U_{j+1},\mathrm{\&lambda;}/L_f)---\left(9\right)$

其中soft(x,α)为收缩算子。

5)求解A^k+1，去掉式子(4)中求解A的目标函数里与A无关的项得到：

$\begin{array}{c}{A}^{k+1}=\arg \underset{A}{\min }{\left|\right|A\left|\right|}_{*}+<Y_1^k,A-\mathrm{\&Phi;}{B}^{k+1}>+<Y_2^k,D-A-E^k>\\ +\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A-\mathrm{\&Phi;}{B}^{k+1}\left|\right|}_F^2+\frac{{\mathrm{\&mu;}}_2^k}{2}{\left|\right|D-A-E^k\left|\right|}_F^2\end{array}---\left(10\right)$

使用配方法将上式改写成：

$A^{k+1}=\arg \underset{A}{\min }{\left|\right|A\left|\right|}_{*}+\frac{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}{2}{\left|\right|A-W^k\left|\right|}_F^2---\left(11\right)$

其中：

$W^k=(Y_2^k-Y_1^k+{\mathrm{\&mu;}}_1^k\mathrm{\&Phi;}{B}^{k+1}-{\mathrm{\&mu;}}_2^k{E}^k+{\mathrm{\&mu;}}_2^{k}D)/({\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k)---\left(12\right)$

对于式子(11)使用奇异值阈值法解得：

$A^{k+1}=U^k\mathrm{soft}({\mathrm{\&Sigma;}}^k,\frac{1}{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}){\left(V^k\right)}^T---\left(13\right)$

其中U^k，V^k分别是W^k的左奇异矩阵和右奇异矩阵。

6)求解E^k+1，关于E^k+1的方程如下式所示：

$E^{k+1}=\arg \underset{E}{\min }\frac{{\mathrm{\&mu;}}_2^k}{2}{\left|\right|D-A^{k+1}-E+\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}\left|\right|}_F^2---\left(14\right)$

将E^k+1分为两部分求解，对于在空间域Ω中的部分，已知E值为0。对于空间Ω以外的地方即中，使用一阶求导求解，将两部分合起来即为E的最终解。

61)在空间域Ω外部的地方。可以一阶求导求得式子(14)的解为：

$E^{k+1}=\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}+D-A^{k+1}---\left(15\right)$

62)再将空间域Ω内部和外部的解联合起来就是E的最终解：

$E^{k+1}=P_{\mathrm{\&Omega;}}\left(0\right)+P_{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}(\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}+D-A^{k+1})---\left(16\right)$

7)重复上述步骤4)、5)、6)直到算法收敛。这时迭代的结果A^k+1、B^k+1和E^k+1就是原问题的最终解A、B和E。

本发明方法将矩阵重建理论与稀疏表示理论相结合，在传统的矩阵重建模型基础上引入字典学习模型，从而解决已有技术无法处理的问题，即实现对像素结构性缺失的图片进行填充(实验流程图如图1所示)。结合附图和实施例的详细说明如下：

1)实验中使用512×512像素的Lena图片(如图2所示)作为原始图，在其上构造了4种不同程度缺失率(10％，20％，30％，50％)的受损图像进行测试(如图3所示)。本发明采用固定大小的字典，所以先将待填充图按照从上向下滑窗的方式分成若干个m×512的图像块。m表示字典中码字的维数，滑窗的步长为5像素。将这若干个m×512的图像块依次填充好，最终再组合起来即可得到原始尺寸512×512的恢复图。填充第一个图像块时，则将其用矩阵A表示，填充当前处理图像块问题为求解如下优化方程：

min||A||_*+λ||B||₁

                 (1)

约束条件A＝ΦB,A+E＝D,P_Ω(E)＝0

||A||_*表示矩阵A的核范数。||·||₁表示矩阵的一范数。Ω是观测空间，p_Ω(·)是投影算子，表示变量投射到空间域Ω内的值。λ为权重系数，在实验中将λ取值为0.01。约束条件中的Φ为训练好的字典，B是字典对应的系数矩阵。D为已知的降质观测矩阵即缺失的受损图像。E代表受损图像中缺失的像素。

11)求解此方程时，采用增广拉格朗日乘子法，方程如下：

$\begin{array}{c}L(A,B,E,Y_1,Y_2)={\left|\right|A\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1,A-\mathrm{\&Phi;B}>\\ +<Y_2,D-A-E>+\frac{{\mathrm{\&mu;}}_1}{2}{\left|\right|A-\mathrm{\&Phi;B}\left|\right|}_F^2+\frac{{\mathrm{\&mu;}}_2}{2}{\left|\right|D-A-E\left|\right|}_F^2\end{array}---\left(2\right)$

L(A,B,E,Y₁,Y₂)即为增广拉格朗日函数。其中的μ₁和μ₂是惩罚因子；Y₁、Y₂均为拉格朗日乘子矩阵；<·,·>表示两个矩阵的内积，<M,N>＝trace(M^TN)，M^T表示矩阵M的转置，trace(·)表示求矩阵的迹；||·||_F表示矩阵的Frobenius(弗罗贝尼乌斯)范数。

12)用迭代求解法求解(2)式，迭代的方程如下：

$\left\{\begin{array}{c}(A^{k+1},B^{k+1},E^{k+1})=\arg \underset{A,B,E}{\min }\left\{L(A,B,E,Y_1^k,Y_2^k)\right\}\\ Y_1^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ Y_2^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\\ \end{array}\right.---\left(3\right)$

上式中的表示使目标函数取最小值时的变量A,B,E的值。ρ₁和ρ₂为倍数因子。k是迭代次数。

2)使用Online Learning字典学习方法训练出字典Φ。

21)构造字典Φ使得矩阵A能够由字典稀疏表示，即满足A＝ΦB，其中B是系数矩阵且是稀疏的。实验中使用Online Learning算法在Kodak图像集上训练出字典Φ。在Kodak图像集中所有图像上共随机选取50000个大小为30×1的像素列作为训练数据。

22)在训练字典时，我们设定字典相关参数如下：待重建矩阵A的行数与字典Φ中元素的维数m相等，即A的行数与Φ的行数均为m，实验中取m＝30。训练出的字典Φ是一个过完备的字典，即字典的列数必须大于其行数。实验中取字典列数为300，则字典Φ的规格为30×300。

3)在求解序列(3)时，本发明使用交替方向法将序列(3)转换成如下序列进行求解：

$\left\{\begin{array}{c}{B}^{k+1}=\arg \underset{B}{\min }\left\{L(A^k,B,E^k,Y_1^k,Y_2^k)\right\}\\ A^{k+1}=\arg \underset{A}{\min }\left\{L(A,B^{k+1},E^k,Y_1^k,Y_2^k)\right\}\\ E^{k+1}=\arg \underset{E}{\min }\left\{L(A^{k+1},B^{k+1},E,Y_1^k,Y_2^k)\right\}\\ Y_1^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ Y_2^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\\ \end{array}\right.---\left(4\right)$

设定好各参数初值，然后按照步骤4)、5)、6)的方法进行迭代求解。实验中设定初值为：K＝1；ρ₁＝ρ₂＝1.1； ${\mathrm{\&mu;}}_1^1={\mathrm{\&mu;}}_2^1=0.05;Y_1^1=Y_2^1=0;$ A¹＝E¹＝B¹＝0。

4)使用APG方法求解第一个未知元B^k+1。

41)去掉式子(4)中求解B的目标函数里与B无关的项后得到如下方程：

$B^{k+1}=\arg \underset{B}{\min }\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2---\left(5\right)$

使用泰勒展开的方法，构造出一个二阶函数来逼近上式，然后针对这个二阶的函数来求解原方程。

令 $f\left(B\right)=<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2,F\left(B\right)=\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2$ 再引入变量Z，定义如下函数：

$Q(B,Z)=\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+f\left(Z\right)+<\&dtri;f\left(Z\right),B-Z>+\frac{L_f}{2}{\left|\right|B-Z\left|\right|}_F^2---\left(6\right)$

其中，代表函数f对Z的Frechet梯度；这里f即为f(B)且Z_j+1代替B，L_f是一个常数，取值为用来保证对于所有的Z，都有F(Z)≤Q(B,Z)。

42)经过上步转化，式子(5)转化成了求解Q(B,Z_j)的最小值来求解，再通过配方得到如下形式：

$\begin{array}{c}{B}^{k+1}=B_{j+1}^k=\arg \underset{B}{\min }Q(B,Z)\\ =\arg \underset{B}{\min }\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+\frac{L_f}{2}{\left|\right|B-U_{j+1}\left|\right|}_F^2\end{array}---\left(7\right)$

其中变量Z_j的更新规则如下：

$\left\{\begin{array}{c}{Z}_{j+1}=B_{j+1}^k+\frac{t_{j+1}-1}{t_{j+1}}(B_{j+1}^k-B_j^k)\\ t_{j+1}=(1+\sqrt{4t_j^2+1})/2\end{array}\right.---\left(8\right)$

t_j是一组序列，j是变量迭代次数。经过上述转化，设定各参数初始值如下：j＝1；t₁＝1；Z₁＝0。收敛时可以解得：

$B^{k+1}=B_{j+1}^k=\mathrm{soft}(U_{j+1},\mathrm{\&lambda;}/L_f)---\left(9\right)$

其中soft(x,α)为收缩算子。

5)求解A^k+1，去掉式子(4)中求解A的目标函数里与A无关的项得到：

$\begin{array}{c}{A}^{k+1}=\arg \underset{A}{\min }{\left|\right|A\left|\right|}_{*}+<Y_1^k,A-\mathrm{\&Phi;}{B}^{k+1}>+<Y_2^k,D-A-E^k>\\ +\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A-\mathrm{\&Phi;}{B}^{k+1}\left|\right|}_F^2+\frac{{\mathrm{\&mu;}}_2^k}{2}{\left|\right|D-A-E^k\left|\right|}_F^2\end{array}---\left(10\right)$

使用配方法将上式改写成：

$A^{k+1}=\arg \underset{A}{\min }{\left|\right|A\left|\right|}_{*}+\frac{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}{2}{\left|\right|A-W^k\left|\right|}_F^2---\left(11\right)$

其中：

$W^k=(Y_2^k-Y_1^k+{\mathrm{\&mu;}}_1^k\mathrm{\&Phi;}{B}^{k+1}-{\mathrm{\&mu;}}_2^k{E}^k+{\mathrm{\&mu;}}_2^{k}D)/({\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k)---\left(12\right)$

对于式子(11)使用奇异值阈值法解得：

$A^{k+1}=U^k\mathrm{soft}({\mathrm{\&Sigma;}}^k,\frac{1}{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}){\left(V^k\right)}^T---\left(13\right)$

其中U^k，V^k分别是W^k的左奇异矩阵和右奇异矩阵。

6)求解E^k+1，关于E^k+1的方程如下式所示：

$E^{k+1}=\arg \underset{E}{\min }\frac{{\mathrm{\&mu;}}_2^k}{2}{\left|\right|D-A^{k+1}-E+\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}\left|\right|}_F^2---\left(14\right)$

将E^k+1分为两部分求解，对于空间Ω以外的地方即中，使用一阶求导即可求解。对于在空间域Ω中的部分，使用一阶求导来求解，将两部分合起来即为E的最终解。

61)在空间域Ω外部的地方。可以使用一阶求导求得式子(14)的解为：

$E^{k+1}=\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}+D-A^{k+1}---\left(15\right)$

62)再将空间域Ω内部和外部的解联合起来就是E的最终解：

$E^{k+1}=P_{\mathrm{\&Omega;}}\left(0\right)+P_{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}(\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}+D-A^{k+1})---\left(16\right)$

7)重复上述步骤4)、5)、6)直至收敛。这时迭代的结果A^k+1、B^k+1和E^k+1就是原问题的最终解A、B和E。

8)依次处理步骤1)中得到的其余若干个图像块直至全部填充好，再将这些图像块组合成最终的填充图(如图4所示)。组合时，被多次填充的像素点取多次填充的均值作为最终值。

实验结果：本发明采用PSNR(峰值信噪比)作为图像填充结果的度量测度：

$\mathrm{PSNR}=10\&times;{\log }_{10}\left(\frac{{(2^n-1)}^2}{\mathrm{MSE}}\right)---\left(17\right)$

$\mathrm{MSE}=\frac{1}{\mathrm{wh}}\underset{y=0}{\overset{w-1}{\mathrm{\&Sigma;}}}\underset{x=0}{\overset{h-1}{\mathrm{\&Sigma;}}}{|I(x,y)-I_0(x,y)|}^2$

其中I₀代表没有受损的真实图像，I为填充后的图像，h为图像的高度，w为图像的宽度，(x,y)为第x行第y列的像素值，Σ表示求和运算，|·|为绝对值。本实验取n＝8，实验中用于测试的4张不同程度缺失的图片恢复的结果见图4标注。

Claims (1)

1.一种基于低秩矩阵重建的结构性缺失图像填充方法，其特征是，包括如下步骤：

1)将图像看作矩阵，则原图像用矩阵A表示，解决像素结构性缺失的图像填充问题为求解如下优化方程：

min||A||_*+λ||B||₁

                                (1)

约束条件 A＝ΦB,A+E＝D,P_Ω(E)＝0

其中||A||_*表示矩阵A的核范数，||·||₁表示矩阵的一范数，Ω是观测空间，p_Ω(·)是投影算子，表示变量投射到空间域Ω内的值，λ为权重系数，约束条件中的Φ为训练好的字典，B是字典对应的系数，D为已知的降质观测矩阵即缺失的受损图像，E代表受损图像中缺失的像素；

求解此方程时，本发明采用增广拉格朗日乘子法，方程如下：

$\begin{array}{c}L(A,B,C,E,Y_1,Y_2)={\left|\right|A\left|\right|}_{*}+\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<Y_1,A-\mathrm{\&Phi;B}>\\ +<Y_2,D-A-E>+\frac{{\mathrm{\&mu;}}_1}{2}{\left|\right|A-\mathrm{\&Phi;B}\left|\right|}_F^2+\frac{{\mathrm{\&mu;}}_2}{2}{\left|\right|D-A-E\left|\right|}_F^2\end{array}---\left(2\right)$

L(A,B,E,Y₁,Y₂)即为增广拉格朗日函数，其中μ₁和μ₂是惩罚因子；Y₁、Y₂均为拉格朗日乘子矩阵；表示两个矩阵的内积；||·||_F表示矩阵的弗罗贝尼乌斯Frobenius范数；

(2)式的迭代求解法方程如下：

$\left\{\begin{array}{c}(A^{k+1},B^{k+1},E^{k+1})=\arg \underset{A,B,E}{\min }\left\{L(A,B,E,Y_1^k,Y_2^k)\right\}\\ {Y_1}^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ {Y_2}^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\end{array}\right.---\left(3\right)$

上式中的表示使目标函数取最小值时的变量A,B,E的值，ρ₁和ρ₂为倍数因子，k是迭代次数；

2)训练字典Φ：在高质量的图像数据集上使用在线学习算法训练出字典Φ；

3)使用交替方向法ADM将序列(3)转换成如下序列进行求解：

$\left\{\begin{array}{c}{B}^{k+1}=\arg \underset{B}{\min }\left\{L(A^k,B,E^k,Y_1^k,Y_2^k)\right\}\\ A^{k+1}=\arg \underset{A}{\min }\left\{L(A,B^{k+1},E^k,Y_1^k,Y_2^k)\right\}\\ E^{k+1}=\arg \underset{E}{\min }\left\{L(A^{k+1},B^{k+1},E,Y_1^k,Y_2^k)\right\}\\ {Y_1}^{k+1}=Y_1^k+{\mathrm{\&mu;}}_1^k(A^{k+1}-\mathrm{\&Phi;}{B}^{k+1})\\ {Y_2}^{k+1}=Y_2^k+{\mathrm{\&mu;}}_2^k(D-A^{k+1}-E^{k+1})\\ {\mathrm{\&mu;}}_1^{k+1}={\mathrm{\&rho;}}_1{\mathrm{\&mu;}}_1^k\\ {\mathrm{\&mu;}}_2^{k+1}={\mathrm{\&rho;}}_2{\mathrm{\&mu;}}_2^k\end{array}\right.---\left(4\right)$

然后按照步骤4)、5)、6)的方法进行迭代求解得到最终的结果；

4)求解B^k+1：使用加速近邻梯度算法求得B^k+1；

去掉式子(4)中求解B的目标函数里与B无关的项，得到如下方程:

$B^{k+1}=\arg \underset{B}{\min }\mathrm{\&lambda;}{\left|\right|B\left|\right|}_1+<{Y_1}^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2---\left(5\right)$

使用泰勒展开的方法，构造出一个二阶函数来逼近上式，然后针对这个二阶的函数来求解原方程，令 $f\left(B\right)=<{Y_1}^k,A^k-\mathrm{\&Phi;B}>+\frac{{\mathrm{\&mu;}}_1^k}{2}{\left|\right|A^k-\mathrm{\&Phi;B}\left|\right|}_F^2,$ 最终可以解得：

$B^{k+1}=B_{j+1}^k=\mathrm{soft}(U_{j+1},\mathrm{\&lambda;}/L_f)---\left(6\right)$

其中的soft(x,α)为收缩算子， $U_{j+1}=Z_{j+1}-\&dtri;f\left(Z_{j+1}\right)/L_f,$ 代表函数f对Z_j+1的Frechet梯度，这里f即为f(B)且Z_j+1代替B，L_f是一个常数，取值为变量Z_j的更新规则如下：

$\left\{\begin{array}{c}{Z}_{j+1}=B_{j+1}^k+\frac{t_{j+1}-1}{t_{j+1}}(B_{j+1}^k-B_j^k)\\ t_{j+1}=(1+\sqrt{4t_j^2+1})/2\end{array}\right.---\left(7\right)$

t_j是一组序列，j是变量迭代次数；

5)求解A^k+1：使用SVT算法求解A^k+1；

去掉式子(4)中求解A的目标函数里与A无关的项，并且通过配方得到：

$A^{k+1}=\arg \underset{A}{\min }{\left|\right|A\left|\right|}_{*}+\frac{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}{2}{\left|\right|A-W^k\left|\right|}_F^2---\left(8\right)$

其中：

$W^k=(Y_2^k-{Y_1}^k+{\mathrm{\&mu;}}_1^k\mathrm{\&Phi;}{B}^{k+1}-{\mathrm{\&mu;}}_2^k{E}^k+{\mathrm{\&mu;}}_2^{k}D)/({\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k)---\left(9\right)$

对于式子(9)使用奇异值阈值法解得：

$A^{k+1}=U^k\mathrm{soft}({\mathrm{\&Sigma;}}^k,\frac{1}{{\mathrm{\&mu;}}_1^k+{\mathrm{\&mu;}}_2^k}){\left(V^k\right)}^T---\left(10\right)$

其中U^k，V^k分别是W^k的左奇异矩阵和右奇异矩阵；

6)求解E^k+1：将E^k+1分为两部分求解，对于空间域Ω内部，E的值为0，对于在空间域Ω以外的部分即中，使用一阶求导来求解，将两部分合起来即为E的最终解：

$E^{k+1}=P_{\mathrm{\&Omega;}}\left(0\right)+P_{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}(\frac{Y_2^k}{{\mathrm{\&mu;}}_2^k}+D-A^{k+1})---\left(11\right)$

7)重复上述步骤4)、5)、6)直到算法收敛，这时迭代的结果A^k+1、B^k+1和E^k+1就是原问题的最终解A、B和E。