CN102609920B - 一种基于压缩感知的彩色数字图像修复方法 - Google Patents

Abstract

一种基于压缩感知的彩色数字图像修复方法，包括步骤一，对原始受损图像进行预处理，将图像破损处的值赋为0；步骤二，对图像进行分块，并分析生成各图像块破损处的退化矩阵；步骤三，将步骤二得到的图像块与退化矩阵相乘，得到待采集信号；然后将待采集信号利用CS稀疏基表示得到相应的稀疏信号，最后将稀疏信号与随机生成的高斯测量矩阵相乘，得到观测信号；步骤四，采用基于全变差的对数障碍函数法对步骤三中得到的CS观测信号进行恢复，得到修复图像。本发明不依附于图像的具体结构，针对不同的图像的恢复效果相似，可以在不同的环境中使用。

Description

一种基于压缩感知的彩色数字图像修复方法

技术领域

本发明涉及图像处理技术领域，尤其是涉及一种基于压缩感知的彩色数字图像修复方法。

背景技术

图像修复是图像复原研究中的一个重要内容，它的主要目的是改善图像质量。对于一幅退化的或者受到噪声污染的图像，其修复过程就是一个利用退化现象的某种先验知识来重建或恢复原有图像的处理过程。作为计算机图形学和计算机视觉的研究热点之一，图像修复技术在文物保护、影视特技制作、虚拟现实、多余物体剔除等方面具有重大的应用价值。

经文献检索发现，传统的图像修复算法都是基于某一类数学模型而提出的，常见的有逆滤波法、偏微分方程法及纹理合成等方法。Bertalmio、Sapiro和Caselles等人在论文“ImageInpainting(图像修复)”(Proceedings of the 27th annual conference on Computer graphics andinteractive techniques，p.417-424，July 2000)(第27届计算机图形与交互技术年度会议)中提出了一种基于数学模型的图像修复算法。该方法利用待修补区域的边缘信息，采用一种由粗到精的处理过程来估计等照度线的方向，利用传播机制将信息传播到待修补的区域内，得到了较好的修补效果。对于图像中的大块丢失信息，目前主要使用基于纹理结构的图像修复技术。Bertalmio等在“Simultaneous Structure and Texture Image Inpainting(基于结构和纹理的同步图像修复)”(Image Processing vol.12，no.8，pp.882-889，Aug.2003)(《图像处理》，2003年第12卷第8期)中首先对图像的结构部分和纹理部分进行分解；然后再采用不同的方法进行修补，结构部分用修补方法进行修补，纹理部分用合成方法填充；最后把这两部分修补的结果叠加起来，就得到了最终的修补图像。由于传统的图像修复算法都是基于某一类数学模型的，因此这类算法都易于理解。但是传统的图像修复算法需要根据原始图像的结构进行分析恢复，不同结构的图像恢复效果也不同，某些结构的图片通过传统的方法恢复出来的效果并不理想。

另一方面，作为一个新颖的信号采集与处理的理论，压缩感知(CS：Compressive Sensing)一经提出就得到了广泛关注。它由稀疏基变换、测量矩阵及信号重构三部分构成，其理论基础是如果信号是可压缩的或在某个变换域上是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得的高维信号投影到一个低维空间上，然后通过求解一个优化问题从这些低维空间的投影中以高概率重构出原始信号。压缩感知理论将信号的采样和编码放在同一步骤，利用信号的稀疏性，以远低于奈奎斯特率的采样频率对信号进行采样，并利用这些采样信号恢复出原始信号。由于图像破损区域相对于整个图像来说是稀疏的，因此将压缩感知理论应用于破损图像的修复是可行的。同时，由于使用压缩感知的图像修复技术不依靠于图像本身的结构进行恢复，与传统的图像修复方法相比，基于压缩感知的图像修复技术可以对各类图像具有相似的恢复效果，其效果不受具体图像结构的影响。

发明内容

本发明的目的是克服上述现有技术的不足，提供一种基于压缩感知的彩色数字图像修复方法，将压缩感知理论应用于图像修复领域。

本发明的技术解决方案如下：

1)对原始受损图像进行预处理，使图像破损处的值赋为0，从视觉效果上看，就是将受损部分图像均处理为黑色；

2)对由1)得到的图像进行分块，对每一图像块分别进行分析，生成该图像块破损处的退化矩阵；

3)将由2)得到的图像块B_ij与退化矩阵H_ij相乘得到待采集信号，然后将待采集信号利用CS稀疏基表示得到相应的稀疏信号，最后将稀疏信号与随机生成的高斯测量矩阵相乘，得到观测信号。

4)采用基于全变差(TV：Total Variation)的对数障碍函数法(Log-barrier algorithm)对3)中得到的观测信号进行恢复，得到修复图像。

所述的步骤1)具体如下：

11)观察原始受损图像，可直观地辨别、定位图像中的受损部分。

12)对于彩色图像，使用图像处理工具将图像受损部分的红色通道(R)、绿色通道(G)、蓝色通道(B)均赋为0，即使用工具将受损部分在视觉上涂成黑色区域；对于灰度图像，将受损区域的灰度置为0，即将图像的受损区域在视觉上涂成黑色区域。

所述的步骤2)具体如下：

21)对图像进行分块处理。经过1)处理后得到的图像尺寸为N＝n×n，将该图像分成互不重叠的尺寸固定为T＝m×m的图像块B_ij，其中n能够被m整除， $1\≤i,j\≤\frac{n}{m}.$

对所有图像块B_ij分别进行下述处理：

22)本发明适用于彩色图像和灰度图像。彩色图像的修复在R、G、B(红、绿、蓝)三通道上分别进行，而灰度图像的处理则仅在灰度上进行处理。对于彩色图像，将图像块B_ij按R、G、B(红、绿、蓝)三通道分为m×m的三个向量I1_ij、I2_ij、I3_ij。为了方便对其进行采样，将I1_ij、I2_ij、I3_ij分别重排为三个一维向量x1_ij、x2_ij、x3_ij，每个向量的维数为T×1，其中T＝m×m，是该图像块的像素值。对于灰度图像，将图像块图像块B_ij按灰色通道分为m×m的一个向量I_ij。

23)对x1_ij、x2_ij、x3_ij进行分析，生成图像块B_ij的退化矩阵H_ij，其步骤具体如下：

首先，将矩阵H_ij初始化为规模是T×1的一维零向量，即该向量内的向量元素均为0；同时，设置一个受损判定阈值 为可接受的误差值。

然后遍历x1_ij、x2_ij、x3_ij中所有的数据值，当x1_ij、x2_ij、x3_ij中对应位置的向量元素值均小于阈值时，算法即认为该点为受损点，将矩阵H_ij中该点值置为0，否则置为1。

重复上述步骤，最终获得的矩阵H_ij就是图像块B_ij破损处的退化矩阵。

所述的步骤3)是将第2)步中得到的的图像块B_ij基于压缩感知理论进行采样测量，具体的步骤如下：

首先，分别将22)中得到的三组向量x1_ij、x2_ij、x3_ij与23)中得到的退化矩阵H_ij相乘，以消除破损区域对压缩感知信号采集的影响，得到图像B_ij的三通道的向量组x1_ij、x2_ij、x3_ij在消除破损区域影响后的CS待采集信号

$\left\{\begin{array}{c}{x1}_{\mathrm{ij}}^{\′}=H_{\mathrm{ij}}\×{x1}_{\mathrm{ij}}\\ {x2}_{\mathrm{ij}}^{\′}=H_{\mathrm{ij}}\×{x2}_{\mathrm{ij}}\\ {x3}_{\mathrm{ij}}^{\′}=H_{\mathrm{ij}}\×{x3}_{\mathrm{ij}}\end{array}\right.$

然后，对上一步中得到的和分别使用快速傅里叶变换(FFT：FastFourier Transform)，得到分块图像B_ij三通道的稀疏信号a1_ij、a2_ij和a3_ij：

$\left\{\begin{array}{c}{a1}_{\mathrm{ij}}=\mathrm{\Ψ}\×{x1}_{\mathrm{ij}}^{\′}\\ {a2}_{\mathrm{ij}}=\mathrm{\Ψ}\×{x2}_{\mathrm{ij}}^{\′}\\ {a3}_{\mathrm{ij}}=\mathrm{\Ψ}\×{x3}_{\mathrm{ij}}^{\′}\end{array}\right.$

其中，Ψ是规模为T×1傅里叶变换基，即对向量分别进行傅里叶变换。

最后，将上一步中得到的稀疏信号a1_ij、a2_ij、a3_ij分别与随机生成的高斯测量矩阵相乘，得到最终的观测信号y1_ij、y2_ij和y3_ij：

$\left\{\begin{array}{c}{y1}_{\mathrm{ij}}=\mathrm{\Φ}\×{a1}_{\mathrm{ij}}=\mathrm{\Φ}\×\mathrm{\Ψ}\×{x1}_{\mathrm{ij}}\\ {y2}_{\mathrm{ij}}=\mathrm{\Φ}\×{a2}_{\mathrm{ij}}=\mathrm{\Φ}\×\mathrm{\Ψ}\×{x2}_{\mathrm{ij}}\\ {y3}_{\mathrm{ij}}=\mathrm{\Φ}\×{a2}_{\mathrm{ij}}=\mathrm{\Φ}\×\mathrm{\Ψ}\×{x3}_{\mathrm{ij}}\end{array}\right.$

其中，y1_ij、y2_ij和y3_ij的维数为M x 1，M为最终测量点的个数，Φ为高斯测量矩阵，大小为M×T，M＜＜T。

所述的高斯测量矩阵，是由随机函数随机生成的矩阵，其中所有的矩阵元素满足高斯分布。

对由整个图像分割得到的各个图像块进行上述操作，就可以获得整个图像的压缩感知测量值。由这些压缩感知测量值，按照下述步骤4)中的方法进行恢复，得到恢复图像。

所述的步骤4)具体如下：

对所有图像块B_ij分别执行下述操作步骤：

41)计算初始猜测值其步骤为：

首先，高斯测量矩阵的转置Φ^T与测量信号y1_ij、y2_ij、y3_ij相乘得到

$\left\{\begin{array}{c}{x1}_{{\mathrm{ij}}_0}^{\′}={\mathrm{\Φ}}^T\×{y1}_{\mathrm{ij}}\\ {x2}_{{\mathrm{ij}}_0}^{\′}={\mathrm{\Φ}}^T\×{y2}_{\mathrm{ij}}\\ {x3}_{{\mathrm{ij}}_0}^{\′}={\mathrm{\Φ}}^T\×{y3}_{\mathrm{ij}}\end{array}\right.$

然后，将快速傅里叶变换基的逆运算Ψ^-1与相乘，得到初始猜测值

$\left\{\begin{array}{c}{x1}_{{\mathrm{ij}}_0}={\mathrm{\Ψ}}^{-1}\×{x1}_{{\mathrm{ij}}_0}^{\′}={\mathrm{\Ψ}}^{-1}\×{\mathrm{\Φ}}^T\×{y1}_{\mathrm{ij}}\\ {x2}_{{\mathrm{ij}}_0}={\mathrm{\Ψ}}^{-1}\×{x2}_{{\mathrm{ij}}_0}^{\′}={\mathrm{\Ψ}}^{-1}\×{\mathrm{\Φ}}^T\×{y2}_{\mathrm{ij}}\\ {x3}_{{\mathrm{ij}}_0}={\mathrm{\Ψ}}^{-1}\×{x3}_{{\mathrm{ij}}_0}^{\′}={\mathrm{\Ψ}}^{-1}\×{\mathrm{\Φ}}^T\×{y3}_{\mathrm{ij}}\end{array}\right.$

其中，的维数均为T x 1。

在以后步骤中，就可以利用对数障碍函数法进行压缩感知重构。

42)设置压缩感知信号恢复迭代参数，包括：

a.恢复过程能够容忍的误差大小η，通常情况下，η≤1e-3

b.对数障碍法中每次迭代的障碍常数(barrier constant)的影响因子μ，通常情况下，μ＝10

c.初始的第一轮全变差参数τ¹＝0

d.初始化迭代计数s＝0

43)计算第一轮全变差参数τ¹：

${\mathrm{\τ}}^1=\frac{M}{{\mathrm{\Σ}}_{k=1}^3\sqrt[2]{({{\mathrm{Dhxk}}_{{\mathrm{ij}}_0}}^2+{{\mathrm{Dvxk}}_{{\mathrm{ij}}_0}}^2)}}$

其中：M为的维数，分别为初始猜测值(k＝1，2，3)的TV值。

所述的TV值的计算方法如下：对一幅大小为P×P的图像上任一点p_ij，其TV值为

$D_{h;\mathrm{ij}}p=\left\{\begin{array}{cc}{p}_{i+1,j}-p_{\mathrm{ij}}& i<P\\ 0& i=P\end{array}\right.,$ $D_{v;\mathrm{ij}}p=\left\{\begin{array}{cc}{p}_{i,j+1}-p_{\mathrm{ij}}& j<P\\ 0& j=P\end{array}\right.$

其中p_i+1，j和p_i，j+1分别是点p_ij在横向和纵向上的相邻点。

44)计算障碍对数法的迭代次数lbiter：

μ是对数障碍法中每次迭代的障碍常数(barrier constant)的影响因子。

45)使用牛顿迭代法(Newton’s Method)分别迭代求解在第s轮的猜测值为初始猜测值为的基于Φx＝y的最优化问题，得到下一轮(即第s+1轮)的猜测点

所述的牛顿迭代法求解最优化问题，具体的描述如下：

对于方程f(x)＝0根的求解，牛顿迭代法使用函数f(x)的泰勒级数的前面几项来寻找最优解。具体的步骤如下：

首先，选择一个接近函数f(x)零点的x₀，计算相应的f(x₀)和切线斜率f′(x₀)(这里f′表示函数f的导数)。然后我们计算穿过点(x₀，f(x₀))并且斜率为f′(x₀)的直线和x轴的交点的x坐标，也就是求如下方程的解：

x f′(x₀)+f(x₀)-x₀f′(x₀)＝0

我们将新求得的点的x坐标命名为x₁，通常x₁会比x₀更接近方程f(x)＝0的解。

因此我们现在可以利用x₁开始下一轮迭代。迭代公式可化简为如下所示：

$x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\&prime;}\left(x_n\right)}$

已经证明，如果f′是连续的，并且待求的零点x是孤立的，那么在零点x周围存在一个区域，只要初始值x₀位于这个邻近区域内，那么牛顿法必定收敛。并且，如果f′(x)不为0，那么牛顿法将具有平方收敛的性能。粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

46)如果说明x_s+1即为理想中的值，此处x_s+1代表即为恢复出来的结果，结束运算，否则进行47)

47)令τ^s+1＝μτ^s，当s＝lbiter时，说明x_s+1代表即为受损图像恢复出来的结果；否则，s＝s+1，并返回步骤45)。

与传统的图像修复方法相比，本发明利用了压缩感知算法的特点：不依附于图像的具体结构，针对不同的图像的恢复效果相似，可以在不同的环境中使用。此外，本发明方法中使用到了针对各个图像块的退化矩阵H_ij可以有效地去除破损区域数据对采样结果的影响，确保了最终的恢复效果。

附图说明

图1为本发明实施例中tree原图。

图2为图1受损后的示例图。

图3为图2经过修复之后的示例图。

具体实施方式

下面结合附图和实施例对本发明进行详细说明，但不应以此限制本发明的保护范围。

本实施例采用大小为256x256的Tree的彩色图像，参见图1，记为X∈R^n×n，其中n＝256。

一种基于压缩感知的彩色数字图像修复方法，包括如下步骤：

第一步，对原始受损图像进行预处理，本实施例利用图像处理工具Adobe Photoshop将图像破损处的R、G、B值赋为0，即受损部分图像均处理为黑色，以便于后续步骤进行处理，经过本步骤处理后得到彩色图像，如图2所示。

第二步，对图像进行分块，得到多个图像块，再对每个图像块进行分析，生成多个图像块破损处的退化矩阵。

首先，对第一步得到的图像进行分块处理。该图像尺寸为N＝256×256，将该图像分成不重叠的尺寸固定为T＝m x m，m＝64的图像块B_ij，其中1≤i，j≤4。接下来所有的步骤均针对该图像块B_ij进行处理。

然后，由于本实施例处理的Tree图像为彩色图像，因此将图像块B_ij分为R、G、B(红、绿、蓝)三通道处理，得到三个64×64的向量I1_ij、I2_ij、I3_ij。为了方便对其进行采样，将I1_ij、I2_ij、I3_ij分别重排为三个一维向量x1_ij、x2_ij、x3_ij，每个向量的维数为T×1，其中T＝4096。

最后，对x1_ij、x2_ij、x3_ij进行分析，生成图像块B_ij的退化矩阵H_ij，其步骤具体如下：

a.将H_ij初始化为规模为4096×1的一维零向量，即该向量内的向量元素均为0；同时，设置一个受损判定阈值

b.然后遍历x1_ij、x2_ij、x3_ij中所有的数据值，当x1_ij、x2_ij、x3_ij中对应位置的向量元素值均小于阈值时，算法即认为该点为受损点，将退化矩阵H_ij中该点值置为0，否则置为1。重复上述步骤，最终获得的矩阵H_ij就是图像块B_ij的退化矩阵。

第三步，将图像块B_ij与退化矩阵H_ij相乘消除破损处信号对采集信号的影响，得到待采集信号，然后将待采集信号利用CS稀疏基表示得到相应的稀疏信号，最后将稀疏信号与随机生成的高斯测量矩阵相乘，得到观测信号。

具体的步骤如下：

首先，将第二步得到的3组向量x1_ij、x2_ij、x3_ij与退化矩阵H_ij相乘，以消除破损区域对压缩感知信号采集的影响，得到图像块B_ij的三通道向量组x1_ij、x2_ij、x3_ij在消除破损区域影响后的CS待采集信号

$\left\{\begin{array}{c}{x1}_{\mathrm{ij}}^{\&prime;}=H_{\mathrm{ij}}\&times;{x1}_{\mathrm{ij}}\\ {x2}_{\mathrm{ij}}^{\&prime;}=H_{\mathrm{ij}}\&times;{x2}_{\mathrm{ij}}\\ {x3}_{\mathrm{ij}}^{\&prime;}=H_{\mathrm{ij}}\&times;{x3}_{\mathrm{ij}}\end{array}\right.$

其次，对上一步中得到的和分别使用快速傅里叶变换(FFT：Fast FourierTransform)，得到各通道的稀疏信号a1_ij、a2_ij和a3_ij：

$\left\{\begin{array}{c}{a1}_{\mathrm{ij}}=\mathrm{\&Psi;}\&times;{x1}_{\mathrm{ij}}^{\&prime;}\\ {a2}_{\mathrm{ij}}=\mathrm{\&Psi;}\&times;{x2}_{\mathrm{ij}}^{\&prime;}\\ {a3}_{\mathrm{ij}}=\mathrm{\&Psi;}\&times;{x3}_{\mathrm{ij}}^{\&prime;}\end{array}\right.$

其中，Ψ是规模为4096×1傅里叶变换基，即对向量分别进行傅里叶变换。

最后，将上一步中得到的稀疏信号a1_ij、a2_ij、a3_ij分别与随机生成的高斯测量矩阵相乘，得到最终的观测信号y1_ij、y2_ij和y3_ij：

$\left\{\begin{array}{c}{y1}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;{a1}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;\mathrm{\&Psi;}\&times;{x1}_{\mathrm{ij}}\\ {y2}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;{a2}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;\mathrm{\&Psi;}\&times;{x2}_{\mathrm{ij}}\\ {y3}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;{a2}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;\mathrm{\&Psi;}\&times;{x3}_{\mathrm{ij}}\end{array}\right.$

其中，y1_ij、y2_ij和y3_ij的维数为500 x 1，M＝500为最终测量点的个数，Φ为高斯测量矩阵，大小为M×T，M＜＜T，其中T＝4096。

第四步，根据观测信号，使用基于全变差(Total Variation)的对数障碍函数法(Log-barrieralgorithm)对观测信号结果进行恢复，得到原始图像。

41)计算初始猜测值其步骤为：

首先，将所述的高斯测量矩阵的转置Φ^T与观测信号y1_ij、y2_ij、y3_ij相乘得到

$\left\{\begin{array}{c}{x1}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Phi;}}^T\&times;{y1}_{\mathrm{ij}}\\ {x2}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Phi;}}^T\&times;{y2}_{\mathrm{ij}}\\ {x3}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Phi;}}^T\&times;{y3}_{\mathrm{ij}}\end{array}\right.$

其次，将傅里叶变换基的逆运算Ψ^-1与相乘，得到最终的初始猜测值

$\left\{\begin{array}{c}{x1}_{{\mathrm{ij}}_0}={\mathrm{\&Psi;}}^{-1}\&times;{x1}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Psi;}}^{-1}\&times;{\mathrm{\&Phi;}}^T\&times;{y1}_{\mathrm{ij}}\\ {x2}_{{\mathrm{ij}}_0}={\mathrm{\&Psi;}}^{-1}\&times;{x2}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Psi;}}^{-1}\&times;{\mathrm{\&Phi;}}^T\&times;{y2}_{\mathrm{ij}}\\ {x3}_{{\mathrm{ij}}_0}={\mathrm{\&Psi;}}^{-1}\&times;{x3}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Psi;}}^{-1}\&times;{\mathrm{\&Phi;}}^T\&times;{y3}_{\mathrm{ij}}\end{array}\right.$

其中，的维数均为4096 x 1。

在以后步骤中，就可以利用对数障碍函数法进行压缩感知重构。

42)设置的进行压缩感知信号恢复迭代的参数：

a.恢复过程中对数障碍法能够容忍的误差大小η＝1e-3

b.对数障碍法中每次迭代的障碍常数(barrier constant)的影响因子：μ＝10

c.初始的第一轮全变差参数τ¹＝0

d.初始迭代计数s＝0

43)计算第一轮全变差参数τ¹：

${\mathrm{\&tau;}}^1=\frac{M}{{\mathrm{\&Sigma;}}_{k=1}^3\sqrt[2]{({{\mathrm{Dhxk}}_{{\mathrm{ij}}_0}}^2+{{\mathrm{Dvxk}}_{{\mathrm{ij}}_0}}^2)}}$

其中：M＝4096，为的维数；分别为初始猜测值(k＝1，2，3)的TV值。

所述的TV值的计算方法如下：对一幅大小为P×P的图像上任一点p_ij，其TV值为

$D_{h;\mathrm{ij}}p=\left\{\begin{array}{cc}{p}_{i+1,j}-p_{\mathrm{ij}}& i<P\\ 0& i=P\end{array}\right.,$ $D_{v;\mathrm{ij}}p=\left\{\begin{array}{cc}{p}_{i,j+1}-p_{\mathrm{ij}}& j<P\\ 0& j=P\end{array}\right.$

其中p_i+1，j和p_i，j+1分别是点p_ij在横向和纵向上的相邻点。

44)计算障碍对数法的迭代次数lbiter：

45)使用牛顿迭代法(Newton’s Method)分别迭代求解在第s轮的猜测值为初始猜测值为的基于Φx＝y的最优化问题，得到第s+1轮的猜测值

所述的牛顿迭代法求解最优化问题，具体的描述如下：

对于方程f(x)＝0根的求解，牛顿迭代法使用函数f(x)的泰勒级数的前面几项来寻找最优解。具体的步骤如下：

首先，选择一个接近函数f(x)零点的x₀，计算相应的f(x₀)和切线斜率f′(x₀)(这里f′表示函数f的导数)。然后我们计算穿过点(x₀，f(x₀))并且斜率为f′(x₀)的直线和x轴的交点的x坐标，也就是求如下方程的解：

x f′(x₀)+f(x₀)-x₀f′(x₀)＝0

我们将新求得的点的x坐标命名为x₁，通常x₁会比x₀更接近方程f(x)＝0的解。因此我们现在可以利用x₁开始下一轮迭代。迭代公式可化简为如下所示：

$x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\&prime;}\left(x_n\right)}$

已经证明，如果f′是连续的，并且待求的零点x是孤立的，那么在零点x周围存在一个区域，只要初始值x₀位于这个邻近区域内，那么牛顿法必定收敛。并且，如果f′(x)不为0，那么牛顿法将具有平方收敛的性能。粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。

46)如果说明x_s+1即为理想中的值，此处x_s+1代表结束运算，否则进行47)

47)令τ^s+1＝μτ^s，当s＝lbiter时，说明x_s+1代表即为受损图像恢复出来的结果；否则，s＝s+1，并返回步骤45)。

最终，恢复出来的图像如图3所示，与原图像图1对比可以得到，本方法充分利用了压缩感知算法的特点：不依附于图像的具体结构，针对不同的图像的恢复效果相似，可以在不同的环境中使用。此外，本方法中使用到了针对各个分块的退化矩阵H_ij可以有效地去除破损区域数据对采样结果的影响，确保了最终的恢复效果。

本发明同样适用于灰度受损图像，其步骤基本相同，不同之处在于，采用图像处理工具将图像受损部分的灰度值赋为0；并将得到的图像进行分块，对每一图像块按灰度通道分为一个向量I_ij。

Claims (5)

1.一种基于压缩感知的彩色数字图像修复方法，其特征在于，该方法包括以下步骤：

步骤一，对原始受损图像进行预处理，将图像破损处的值赋为0；

步骤二，对由步骤一得到的图像进行分块，对每一图像块分别进行分析，生成各图像块破损处的退化矩阵；

步骤三，将步骤二得到的图像块与退化矩阵相乘，得到待采集信号；然后将待采集信号利用CS稀疏基表示得到相应的稀疏信号，最后将稀疏信号与随机生成的高斯测量矩阵相乘，得到观测信号，具体如下：

21)对图像进行分块处理：经过步骤一处理后得到的图像尺寸为N＝n×n，将该图像分成互不重叠的尺寸固定为T＝m×m的图像块B_ij，其中n能够被m整除，

对所有图像块B_ij分别进行下述处理：

22)当受损图像为彩色图像时，将图像块B_ij分为红、绿、蓝三通道处理，得到三个m×m的向量I1_ij、I2_ij、I3_ij；并将I1_ij、I2_ij、I3_ij分别重排为三个一维向量x1_ij、x2_ij、x3_ij，每个向量的维数为T×1；

23)对x1_ij、x2_ij、x3_ij进行分析，生成图像块B_ij的退化矩阵H_ij，其步骤具体如下：

首先，将H_ij初始化为规模为T×1的一维零向量；设置受损判定阈值；

然后，遍历x1_ij、x2_ij、x3_ij中所有的数据值，当x1_ij、x2_ij、x3_ij中对应位置的向量元素值均小于阈值时，判定该点为受损点，将矩阵H_ij中该点值置为0，否则置为1；

重复上述步骤，最终获得的矩阵H_ij就是图像块B_ij的退化矩阵;

步骤四，采用基于全变差的对数障碍函数法对步骤三中得到的CS观测信号进行恢复，得到修复图像。

2.根据权利要求1所述的基于压缩感知的彩色数字图像修复方法，其特征在于，所述的步骤一包括如下步骤：

11)获取受损图像，对原始受损图像中受损部分进行辨别；

12)对于彩色图像，使用图像处理工具将图像受损部分的红色通道(R)、绿色通道(G)、蓝色通道(B)均赋为0；对于灰度图像，将受损区域的灰度置为0。

3.根据权利要求1所述的基于压缩感知的彩色数字图像修复方法，其特征在于，所述的步骤三具体如下：

首先，分别将所述的一维向量x1_ij、x2_ij、x3_ij与所述的退化矩阵H_ij相乘，得到待采集信号x1′_ij、x2′_ij、x3′_ij：

$\left\{\begin{array}{c}{x1}_{\mathrm{ij}}^{\&prime;}=H_{\mathrm{ij}}\&times;{x1}_{\mathrm{ij}}\\ {x2}_{\mathrm{ij}}^{\&prime;}=H_{\mathrm{ij}}\&times;{x2}_{\mathrm{ij}}\\ {x3}_{\mathrm{ij}}^{\&prime;}=H_{\mathrm{ij}}\&times;{x3}_{\mathrm{ij}}\end{array}\right.$

其次，对上一步中得到的x1′_ij、x2′_ij和x3′_ij分别使用快速傅里叶变换，得到各通道的稀疏信号a1_ij、a2_ij和a3_ij：

$\left\{\begin{array}{c}{a1}_{\mathrm{ij}}=\mathrm{\&psi;}\&times;{x1}_{\mathrm{ij}}^{\&prime;}\\ {a2}_{\mathrm{ij}}=\mathrm{\&psi;}\&times;{x2}_{\mathrm{ij}}^{\&prime;}\\ {a3}_{\mathrm{ij}}=\mathrm{\&psi;}\&times;{x3}_{\mathrm{ij}}^{\&prime;}\end{array}\right.$

其中，Ψ是规模为T×1傅里叶变换基，即对向量x1′_ij、x2′_ij，x3′_ij分别进行傅里叶变换；

最后，将上一步中得到的稀疏信号a1_ij、a2_ij、a3_ij分别与随机生成的高斯测量矩阵相乘，得到最终的观测信号y1_ij、y2_ij和y3_ij：

$\left\{\begin{array}{c}{y1}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;{a1}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;\mathrm{\&psi;}\&times;{x1}_{\mathrm{ij}}\\ {y2}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;{a2}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;\mathrm{\&psi;}\&times;{x2}_{\mathrm{ij}}\\ {y3}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;{a2}_{\mathrm{ij}}=\mathrm{\&Phi;}\&times;\mathrm{\&psi;}\&times;{a3}_{\mathrm{ij}}\end{array}\right.$

其中，y1_ij、y2_ij和y3_ij的维数为M×1，M为最终测量点的个数，Φ为高斯测量矩阵，大小为M×T，M＜＜T。

4.根据权利要求3所述的基于压缩感知的彩色数字图像修复方法，其特征在于，所述的高斯测量矩阵，是由随机函数随机生成的矩阵，其中所有的矩阵元素满足高斯分布。

5.根据权利要求1所述的基于压缩感知的彩色数字图像修复方法，其特征在于，所述的步骤四具体如下：

41)计算初始猜测点其步骤为：

首先，将所述的高斯测量矩阵的转置Φ^T与测量信号y1_ij，y2_ij，y3_ij相乘得到

$\left\{\begin{array}{c}{x1}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Phi;}}^T\&times;{y1}_{\mathrm{ij}}\\ {x2}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Phi;}}^T\&times;{y2}_{\mathrm{ij}}\\ {x3}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&Phi;}}^T\&times;{y3}_{\mathrm{ij}}\end{array}\right.$

然后，再将傅里叶变换基的逆运算Ψ^-1与相乘，得到最终的初始猜测点

$\left\{\begin{array}{c}{x1}_{{\mathrm{ij}}_0}={\mathrm{\&psi;}}^{-1}\&times;{x1}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&psi;}}^{-1}\&times;{\mathrm{\&Phi;}}^T\&times;{y1}_{\mathrm{ij}}\\ {x2}_{{\mathrm{ij}}_0}={\mathrm{\&Psi;}}^{-1}\&times;{x2}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&psi;}}^{-1}\&times;{\mathrm{\&Phi;}}^T\&times;{y2}_{\mathrm{ij}}\\ {x3}_{{\mathrm{ij}}_0}={\mathrm{\&psi;}}^{-1}\&times;{x3}_{{\mathrm{ij}}_0}^{\&prime;}={\mathrm{\&psi;}}^{-1}\&times;{\mathrm{\&Phi;}}^T\&times;{y3}_{\mathrm{ij}}\end{array}\right.$

其中，的维数均为T×1;

在以后步骤中，就可以利用对数障碍函数法进行压缩感知重构；

42)设置压缩感知信号恢复迭代参数，包括：

a.恢复过程能够容忍的误差大小η，通常情况下，η≤1e-3

b.对数障碍法中每次迭代的障碍常数(barrier constant)的影响因子μ，通常情况下，μ＝10

c.初始的第一轮全变差参数τ¹＝0

d.初始化迭代计数s=0

43)计算第一轮全变差参数τ¹：

${\mathrm{\&tau;}}^1=\frac{M}{{\mathrm{\&Sigma;}}_{k=1}^3\sqrt[2]{\left({{\mathrm{Dhxk}}_{{\mathrm{ij}}_0}}^2+{{\mathrm{Dvxk}}_{{\mathrm{ij}}_0}}^2\right)}}$

其中：M为的维数，分别为初始猜测值k=1，2，3的TV值；

所述的TV值的计算方法如下：对一幅大小为P×P的图像上任一点p_ij，其TV值为

$D_{h;\mathrm{ij}}p=\left\{\begin{array}{cc}{p}_{i+1,j}-p_{\mathrm{ij}}& i<p\\ 0& i=p\end{array}\right.,D_{v;\mathrm{ij}}p=\left\{\begin{array}{cc}{p}_{i,j+1}-p_{\mathrm{ij}}& j<p\\ 0& j=p\end{array}\right.$

其中p_i+1，j和p_i，j+1分别是点p_ij在横向和纵向上的相邻点;

44)计算障碍对数法的迭代次数lbiter:

45)使用牛顿迭代法分别迭代求解在第s轮的猜测值为初始猜测值为的基于Φx＝y的最优化问题，得到第s+1轮的猜测值

46)如果说明x_s+1代表即为受损图像恢复出来的结果，结束运算，否则进行47)；

47)令τ^s+1＝μτ^s，当s=lbiter时，说明x_s+1代表即为受损图像恢复出来的结果；否则，s＝s+1，并返回步骤45)。