CN104952051A - 基于高斯混合模型的低秩图像修复方法 - Google Patents

Abstract

本发明公开了一种基于高斯混合模型的低秩图像修复方法,主要解决现有方法在图像修复中没有兼顾图像整体的低秩结构和局部的相似性，造成修复图像的模糊与失真问题。本发明将图像建模成低秩项和非局部相似性算子之和，其中：低秩项由核范数构成，用于捕捉修复图像整体的低秩结构；非局部相似性算子是基于高斯混合模型来构造，它能够根据相似图像块的统计特性来完善修复图像的局部细节。实验结果表明，与其它传统方法相比，本发明的修复结果无论在数值评价指标上还是在视觉效果上，都明显优于其它方法。可用于对具有低秩结构的受损自然图像修复。

Description

基于高斯混合模型的低秩图像修复方法

技术领域

本发明属于图像处理技术领域，特别是一种图像修复方法，可用于对具有低秩结构的受损自然图像修复。

背景技术

图像修复是图像处理领域中一个重要的分支，它试图利用受损图像仅有的信息来恢复出原始图像，这对于很多领域都有着重要的应用价值。近年来，随着低秩理论的进步，低秩图像修复问题也受到了越来越多的关注。

低秩图像修复往往被建模成秩最小化问题。假设原始图像用二维矩阵来表示，矩阵表示受损后的图像，则低秩图像修复可以建模成如下公式：

其中rank(·)是对矩阵X求秩算子，Ω是未受损像素点位置集合，是正交投影算子，公式如下：

由于求秩运算是NP-hard问题，所以求解秩最小问题是非凸优化问题。为了提高秩最小问题的可解性，Candès等人提出了利用核范数来代替求秩算子，并且证明了，核范数是求秩算子最紧致的凸逼近，则上述模型可转化为如下公式：

其中，||X||_*定义为X的核范数。解决核范数最小化的方法有很多，比如奇异值阈值法SVT、核范数最小二乘法NNLS、鲁棒主成份分解法RPCA等等。为了更好的逼近真实图像的秩，Hu等人提出了截断的核范数TNNR，并且给出了三种求解方法：TNNR-ADMM，TNNR-APGL和TNNR-ADMMAP。实验证明，TNNR比基于核范数的方法得到更精确的修复结果。

以上提到的方法都是从低秩矩阵填充的角度对图像修复问题进行建模的，但是自然图像很少是严格低秩的，运用秩最小化必然会损失掉图像内在的结构和丰富的纹理信息。所以一个鲁棒的低秩图像修复方法不仅要能够抓住图像整体的低秩结构，还要使得图像的细节信息得到保留。

近年来，非局部相似性的运用使得很多图像处理的性能有了质的提升。研究表明，图像中有很多结构和纹理相似的小图像块，利用这些图像块的相似性，对重构图像的特性加以约束，可以明显地提升图像重构的效果。Yu等人提出的分段线性估计器PLE就是典型的利用非局部相似性解决图像逆问题的方法。它采用高斯混合模型来对不同类型的图像块进行建模，通过一个线性估计器，来对相似的图像块进行聚类，之后，用属于同一类的图像块对该类的高斯分布进行参数的更新，交替以上两个步骤，直到得到满意结果。随着低秩技术的发展，很多非局部模型也嵌入了低秩的约束。Dong将相似图像块的系数矩阵替换成一个对角矩阵和一个右乘矩阵，通过稀疏性的假设，将系数求解问题转换成了低秩逼近的问题。Wang在解决图像超分辨问题时，利用相似图像块构造一个增广矩阵，在假设低分辨像素和高分辨像素具有相同线性表示系数的条件下，证明该增广矩阵具有低秩性，从而将图像超分辨问题转化成低秩矩阵求解问题。此外，Dong基于非局部相似性和低秩假设提出了一种压缩感知重构方法，叫做NLR-CS，NLR-CS，其效果优于目前许多优秀的压缩感知重构方法。虽然，以上提出的模型都将低秩与非局部相似性结合起来，但是，所有的方法都是利用低秩正则来约束相似的图像块。也就是说，这些模型利用低秩正则来增强图像块之间的非局部相似性，而没有考虑到整个图像潜在的低秩结构，从而导致修复结果的整体结构性不强，不能较好地保持图像中的低秩信息。

发明内容

本发明提出一种基于高斯混合模型的低秩图像修复方法，其目的在于构建一个鲁棒有效的图像修复模型，以充分发掘图像潜在的低秩结构，并保留图像的细节信息，获得更精确的修复结果。

实现本发明目的技术方案包括如下步骤：

(1)输入受损图像Y，构造如下修复图像模型：

其中，α是常数，X是待求解的修复图像，||X||_*是X的核范数，y_i是从受损图像Y中提取的第i个图像块，x_i是从待求解的修复图像X中提取的第i个图像块，其中，1≤i≤N，N是图像块数量，U_i是每个图像块对应的退化算子，是属于第j个高斯分布的所有图像块标号的集合，是图像块x_i服从的高斯分布，Ω是受损像素点位置的集合，是正交投影算子，其中，下标m,n表示矩阵的第m行、第n列的元素；

(2)设迭代计数器k＝0，利用迭代收缩阈值方法对待求解的修复图像X初始化，得到初始化修复结果X^k；

(3)初始化修复图像块的J个高斯分布其中为均值向量，为协方差矩阵；

(4)将修复结果X^k和高斯分布代入修复图像模型进行修复图像的3次迭代求解:

(4.1)利用分段线性估计器更新高斯分布为

(4.2)利用特征值收缩算子更新修复结果为X^k+1；

(5)判断迭代次数是否等于3，若是，则将更新后的修复结果X^k+1赋给待修复图像X，完成对受损图像的修复；否则，迭代计数器k＝k+1，返回步骤(4)。

本发明由于集成了低秩逼近和非局部相似性技术手段，因而具有如下优点：

(A)能够捕捉到图像潜在的整体低秩结构。

(B)充分抓住了图像局部的细节信息。

(C)修复结果优于传统的基于低秩逼近的算法和非局部模型。

实验证明，与传统的低秩逼近算法和非局部模型相比，本发明修复的图像具有更好的整体结构性和更精确的细节信息。

附图说明

图1是本发明的实现流程图；

图2是用本发明与现有方法对于图像Basketwork和Field修复结果的评价指标对比；

图3是用本发明与现有方法在采样率为30％时，对图像Basketwork的修复结果图；

图4是用本发明与现有方法在采样率为30％时，对图像Field的修复结果图；

图5是用本发明与现有方法对图像Sunset去除文字的修复结果图；

图6是用本发明与现有方法对图像Lena去除文字的修复结果图；

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1，输入受损图像Y，构造如下修复图像模型。

假设一幅自然图像Y中存在着许多被噪声或者文字遮挡的像素点，定义图像Y为受损图像，而受损图像Y中被遮挡住的像素点统称为受损像素点。

1.1)输入受损图像Y，设未受损像素点的位置集合为Ω；

1.2)对于待修复图像X，用核范数||X||_*来描述它潜在的低秩结构：

1.3)对于待修复图像X中的局部细节信息，用高斯混合模型来进行表示

1.4)联合步骤(1.2)和(1.3)中的公式，得到修复图像模型：

其中，α是常数，取值为30，X是待求解的修复图像，||X||_*是X的核范数，y_i是从受损图像Y中提取的第i个图像块，x_i是从待求解的修复图像X中提取的第i个图像块，其中，1≤i≤N，N是图像块数量，U_i是每个图像块对应的退化算子，是属于第j个高斯分布的所有图像块标号的集合，是图像块x_i服从的高斯分布，是正交投影算子，其中，下标m,n表示矩阵的第m行、第n列的元素。

步骤2,设迭代计数器k＝0，利用迭代收缩阈值方法对待求解的修复图像X初始化，得到初始化修复结果X^k。

2.1)生成一个尺寸与待修复图像X相等的图像X^k，使得X^k中所有像素值都为0，对图像X^k进行步骤(2.2)到(2.5)的迭代计算；

2.2)利用正交投影算子使图像X^k中的像素点尽量逼近受损图像Y中未受损的像素点：

2.3)对步骤(2.2)中得到的图像X^k做奇异值分解，将其分解为三个不同的矩阵，即：

X^k＝UΛV^T，

其中，U和V分别是左奇异矩阵和右奇异矩阵，Λ＝diag({λ₁,λ₂,...})是对角矩阵，其对角线元素{λ₁,λ₂,...}为图像X^k的特征值；

2.4)利用步骤(2.3)得到的左奇异矩阵U、右奇异矩阵V和对角矩阵Λ构造关于图像X^k的特征值收缩算子

其中，(λ_i-β)₊＝max(λ_i-β,0)，即取λ_i-β和0二者之中值较大者，实验中β取值为100；

2.5)利用步骤(2.4)构造的特征值收缩算子按如下公式对图像X^k进行计算：

2.6)若初始化数到达300次，则输出初始化修复结果Xk；否则，返回步骤(2.2)。

步骤3，初始化修复图像块的J个高斯分布为其中为均值向量，为协方差矩阵。

3.1)构造人工黑白图像：生成4幅大小为65×65像素的人工黑白图{I_j}_1≤j≤4，其中黑色像素点值为0，白色像素点值为1，每幅黑白图I_j的黑色区域与白色区域的分界线均经过图像的中心坐标(33,33)，这4幅人工黑白图分界线的角度分别为0°、45°、90°和135°；

3.2)从人工黑白图{I_j}_1≤j≤4中抽取图像块样本：在每幅人工黑白图I_j中，以尺寸为8×8像素的窗口取出那些接触到黑白分界线的图像块，从第j幅图像取出的图像块样本集合记为其中，表示从第j幅图像中取出的第i个图像块，表示集合中图像块的数量；

3.3)根据图像块样本集合按如下公式计算初始高斯分布的均值和协方差矩阵

步骤4，将修复结果X^k和高斯分布代入步骤1的修复图像模型中进行修复图像的3次迭代求解。

4.1)利用分段线性估计器更新高斯分布为

4.1.1)计算所有受损图像块{y_i}_1≤i≤N在每个高斯分布下的重构结果：

${\tilde{x}}_i^j={\mathrm{\Σ}}_j^k{U}_i^T{(U_i{\mathrm{\Σ}}_j^k{U}_i^T+{\mathrm{\σ}}^{2}I)}^{-1}{y}_i,$

其中，是第i个图像块在第j个高斯分布下的重构结果，σ是系统对噪声的容忍度，取值为3，I是一个单位矩阵；

4.1.2)利用每个高斯分布下的重构结果通过下式选取与第i个图像块最匹配的高斯分布：

${\tilde{j}}_i=\arg \underset{j}{\min }\left|\right|U_i{\tilde{x}}_i^j-y_i|{|}_2^2+{\mathrm{\σ}}^2{\left({\tilde{x}}_i^j\right)}^T{\left({\mathrm{\Σ}}_j^k\right)}^{-1}{\tilde{x}}_i^j+{\mathrm{\σ}}^{2}log|{\mathrm{\Σ}}_i^k|,$

其中，||·||是行列式算子；

4.1.3)根据步骤(4.1.2)求得最佳匹配的高斯分布标号将图像块分成J类，得到J个类别集合其中第j类集合按如下公式计算：

4.1.4)利用类别集合更新第j个高斯分布的均值和协方差矩阵

其中，是类别集合中元素的个数，是从X^k中抽取的第i个图像块；

4.2)利用特征值收缩算子更新修复结果为X^k+1；

4.2.1)对于所有的受损图像块{y_i}_1≤i≤N，利用更新后的高斯协方差矩阵按如下公式计算其重构结果：

${\tilde{x}}_i={\mathrm{\Σ}}_{{\tilde{j}}_i}^{k+1}{U}_i^T{(U_i{\mathrm{\Σ}}_{{\tilde{j}}_i}^{k+1}{U}_i^T+{\mathrm{\σ}}^{2}I)}^{-1}{y}_i;$

4.2.2)将所有重构的图像块拼接起来，构成完整的重构图像

4.2.3)按如下公式更新修复结果的每个像素

其中，是针对重构图像构造的特征值收缩算子；

根据对所有像素的更新，得到更新后的修复结果X^k+1。

步骤5，判断迭代次数是否等于3，若是，则将更新后的修复结果X^k+1赋给待修复图像X，完成对受损图像的修复；否则，迭代计数器k＝k+1，返回步骤(4)。

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

1.实验条件

实验所用计算机的处理器为Intel Core(TM)2Duo CPU 2.67GHz，内存4GB，编程平台是Matlab 2013。实验中采用的图像数据为Basketwork、Field、Sunset和Lena图像，大小为256×256像素。

仿真使用方法有：本发明方法NNR-PLE、奇异值阈值法SVT、固定点Bregman迭代法FPCA、非凸优化方法NonCVX、截断的核范数法TNNR-ADMM和分段线性估计法PLE。

2.实验内容

实验1：在不同采样率下用本发明方法和现有5种方法对图像Basketwork和Field进行修复，对比不同算法的峰值信噪比PSNR值，结果如图2。其中：

图2(a)是用本发明方法与现有5种方法对于图像Basketwork的修复结果的PSNR曲线；

图2(b)是用本发明方法与现有5种方法对于图像Field的修复结果的PSNR曲线。

从图2可以看出，本发明修复结果的PSNR值要明显高于其它方法。

实验2：在采样率为30％下，用本发明方法和现有5种方法对随机遮挡的图像Basketwork进行修复，结果如图3所示，其中：

图3(a)是Basketwork原图，

图3(b)是采样率为30％的Basketwork图像，

图3(c)是用本发明对图3(b)的修复结果，其峰值信噪比PSNR为26.32dB，

图3(d)是用SVT对图3(b)的修复结果，其峰值信噪比PSNR为21.82dB，

图3(e)是用FPCA对图3(b)的修复结果，其峰值信噪比PSNR为22.60dB，

图3(f)是用NonCVX对图3(b)的修复结果，其峰值信噪比PSNR为23.26dB，

图3(g)是用TNNR-ADMM对图3(b)的修复结果，其峰值信噪比PSNR为23.49dB，

图3(h)是用PLE对图3(b)的修复结果，其峰值信噪比PSNR为24.76dB。

实验3：在采样率为30％下，用本发明方法和现有5种方法对随机遮挡的图像Field进行修复，结果如图4所示，其中：

图4(a)是Field原图，

图4(b)是采样率为30％的Field图像，

图4(c)是用本发明对图4(b)的修复结果，其峰值信噪比PSNR为29.26dB，

图4(d)是用SVT对图4(b)的修复结果，其峰值信噪比PSNR为23.77dB，

图4(e)是用FPCA对图4(b)的修复结果，其峰值信噪比PSNR为23.94dB，

图4(f)是用NonCVX对图4(b)的修复结果，其峰值信噪比PSNR为25.21dB，

图4(g)是用TNNR-ADMM对图4(b)的修复结果，其峰值信噪比PSNR为25.09dB，

图4(h)是用PLE对图4(b)的修复结果，其峰值信噪比PSNR为27.81dB。

从图3和图4中可以看出，本发明的修复结果不仅反映出了图像中潜在的低秩结构，而且对于图像的细节信息也能够很好的还原出来，从峰值信噪比指标来看，本发明都优于其它方法。

实验4：用本发明方法和现有5种方法对图像Sunset进行去除文字修复，结果如图5所示，其中：

图5(a)是Sunset原图，

图5(b)是被文字污染的Sunset图像，

图5(c)是用本发明对图5(b)的修复结果，其峰值信噪比PSNR为42.80dB，

图5(d)是用SVT对图5(b)的修复结果，其峰值信噪比PSNR为30.53dB，

图5(e)是用FPCA对图5(b)的修复结果，其峰值信噪比PSNR为31.12dB，

图5(f)是用NonCVX对图5(b)的修复结果，其峰值信噪比PSNR为35.49dB，

图5(g)是用TNNR-ADMM对图5(b)的修复结果，其峰值信噪比PSNR为25.82dB，

图5(h)是用PLE对图5(b)的修复结果，其峰值信噪比PSNR为41.20dB，

实验5：用本发明方法和现有5种方法对图像Lena进行去除文字修复，结果如图6所示，其中：

图6(a)是Lena原图，

图6(b)是被文字污染的Lena图像，

图6(c)是用本发明对图6(b)的修复结果，其峰值信噪比PSNR为35.91dB，

图6(d)是用SVT对图6(b)的修复结果，其峰值信噪比PSNR为26.14dB，

图6(e)是用FPCA对图6(b)的修复结果，其峰值信噪比PSNR为26.79dB，

图6(f)是NonCVX对图6(b)的修复结果，其峰值信噪比PSNR为29.81dB，

图6(g)是用TNNR-ADMM对图6(b)的修复结果，其峰值信噪比PSNR为30.53dB，

图6(h)是用PLE对图6(b)的修复结果，其峰值信噪比PSNR为35.63dB。

从图5和图6中可以看出，本发明方法和PLE方法的修复结果更加真实，清晰，文字去除的非常干净，而其它算法有不同程度的失真和模糊。

综合以上实验结果，在多种图像修复应用中，本发明无论在视觉效果上还是峰值信噪比评价指标上，都优于其它现有5种方法。

Claims (4)

1.一种基于高斯混合模型的低秩图像修复方法，包括如下步骤：

(1)输入受损图像Y，构造如下修复图像模型：

其中，α是常数，X是待求解的修复图像，||X||_*是X的核范数，y_i是从受损图像Y中提取的第i个图像块，x_i是从待求解的修复图像X中提取的第i个图像块，其中，1≤i≤N，N是图像块数量，U_i是每个图像块对应的退化算子，是属于第j个高斯分布的所有图像块标号的集合，是图像块x_i服从的高斯分布，Ω是未受损像素点位置的集合，是正交投影算子，其中，下标m,n表示矩阵的第m行、第n列的元素；

(2)设迭代计数器k＝0，利用迭代收缩阈值方法对待求解的修复图像X初始化，得到初始化修复结果X^k；

(3)初始化修复图像块的J个高斯分布其中为均值向量，为协方差矩阵；

(4)将修复结果X^k和高斯分布代入修复图像模型进行修复图像的3次迭代求解:

(4.1)利用分段线性估计器更新高斯分布为

(4.2)利用特征值收缩算子更新修复结果为X^k+1；

(5)判断迭代次数是否等于3，若是，则将更新后的修复结果X^k+1赋给待修复图像X，完成对受损图像的修复；否则，迭代计数器k＝k+1，返回步骤(4)。

2.根据权利要求1所述的基于高斯混合模型的低秩图像修复方法，其中所述步骤(1)的修复图像模型，按如下步骤构造：

(1.1)对于待修复图像X，用核范数||X||_*来描述它潜在的低秩结构：

(1.2)对于待修复图像X中的局部细节信息，用高斯混合模型来进行表示：

(1.3)联合步骤(1.1)和(1.2)中的公式，得到修复图像模型：

3.根据权利要求1所述的基于高斯混合模型的低秩图像修复方法，其中所述步骤(4.1)中利用分段线性估计器更新高斯分布为按如下步骤进行：

(4.1.1)计算所有受损图像块{y_i}_1≤i≤N在每个高斯分布下的重构结果：

${\tilde{x}}_i^j={\mathrm{\Σ}}_j^k{U}_i^T{(U_i{\mathrm{\Σ}}_j^k{U}_i^T+{\mathrm{\σ}}^{2}I)}^{-1}{y}_i,$

其中，是第i个图像块在第j个高斯分布下的重构结果，σ是模型对噪声的容忍度，取值为3，I是一个单位矩阵；

(4.1.2)利用每个高斯分布下的重构结果通过下式选取与第i个图像块最匹配的高斯分布：

${\tilde{j}}_i=\arg \underset{j}{\min }\left|\right|U_i{\tilde{x}}_i^j-y_i|{|}_2^2+{\mathrm{\σ}}^2{\left({\tilde{x}}_i^j\right)}^T{\left({\mathrm{\Σ}}_j^k\right)}^{-1}{\tilde{x}}_i^j+{\mathrm{\σ}}^{2}log\left|{\mathrm{\Σ}}_j^k\right|,$

其中，|·|是行列式算子；

(4.1.3)利用步骤(4.1.2)求得的最匹配的高斯分布将图像块分成J类，得到J个类别标号集合其中第j类的标号按如下公式计算：

(4.1.4)按如下公式更新第j个高斯分布的均值和协方差矩阵

其中，是集合中元素的个数，是从X^k中抽取的第i个图像块。

4.根据权利要求1所述的基于高斯混合模型的低秩图像修复方法，其中所述步骤(4.2)中利用特征值收缩算子更新待修复图像为X^k+1，按如下步骤进行：

(4.2.1)对于所有的受损图像块{y_i}_1≤i≤N，利用更新后的高斯协方差矩阵计算其新的重构结果，其中第i个图像块的计算公式如下：

${\tilde{x}}_i={\mathrm{\Σ}}_{{\tilde{j}}_i}^{k+1}{U}_i^T{(U_i{\mathrm{\Σ}}_{{\tilde{j}}_i}^{k+1}{U}_i^T+{\mathrm{\σ}}^{2}I)}^{-1}{y}_i;$

(4.2.2)将所有的重构图像块拼接起来，构成完整的重构图像

(4.2.3)对重构图像做奇异值分解，将其分解为三个不同的矩阵，即：

$\tilde{X}={\mathrm{U\ΛV}}^T,$

其中，U和V分别是左奇异矩阵和右奇异矩阵，Λ＝diag({λ₁,λ₂,...})是对角矩阵，其对角线元素{λ₁,λ₂,...}为重构图像的奇异值，diag(·)为对角矩阵构造算子；

(4.2.4)利用左奇异矩阵U、右奇异矩阵V和对角矩阵Λ构造关于重构图像的奇异值收缩算子

其中，(λ_i-α)₊＝max(λ_i-α,0)。

(4.2.4)按如下公式更新修复结果的每个像素

其中，是针对重构图像构造的特征值收缩算子，

根据对所有像素的更新，得到更新后的修复结果X^k+1。