CN103093431A - 基于pca字典和结构先验信息的压缩感知重构方法 - Google Patents

Abstract

本发明公开了一种基于PCA字典和结构先验信息的压缩感知重构方法，主要解决现有压缩感知重构方法纹理和边缘模糊的问题，其实现的步骤为：构造各方向黑白图像；将黑白图像分块得到训练样本；对训练样本PCA分解，获得PCA方向基；接收观测矩阵和测量向量；根据接收信息对稀疏表示模型求解，得到最优方向；通过最优方向得到最优稀疏表示图像，并对该稀疏表示图像进行优化；对优化后图像依次进行滤波和凸投影操作后输出。本发明能更加稀疏和自适应的表示任意方向的图像信号，提高了重构图像的纹理和边缘质量，可用于图像恢复时获得高质量的图像。

Description

基于PCA字典和结构先验信息的压缩感知重构方法

技术领域

本发明属于图像处理技术领域，更进一步涉及压缩感知图像重构方法，可用于在对原图像进行恢复时，获得高清晰质量的图像。

背景技术

在图像重构技术领域中，一种新的数据采集理论——压缩感知理论为数据采集技术带来了新的革命。该理论指出，信号可以进行低速采样和少量采样，并且可以精确重构，这样大大降低了设备存储限制和计算的复杂度。目前压缩感知已成为学术界研究的热点，并不断被应用在图像处理领域和无线传感领域中。压缩感知理论主要包括信号的稀疏表示、信号的观测和信号的重构等三个方面。在信号稀疏表示方面，常用字典有余弦字典、脊波字典等，在信号重构方面，通过求解l₀或l₁范数的优化问题来重构图像。

T.Blumensath等人在文献“T.Blumensath,M.Davies，Iterative HardThresholding for Compressed Sensing”中提出迭代硬阈值压缩感知重构方法。该方法对稀疏信号进行低采样的随机观测并使用正交基，进而通过选择重要的稀疏表示系数，舍弃不重要的稀疏表示系数来重构图像。该方法存在的不足是，由于在信号稀疏表示中使用正交基，因此不能保证对所有图像如纹理图像进行更加稀疏的表示，从而导致重构出的图像不够清晰。

西安电子科技大学的专利申请“压缩感知框架下基于非凸模型的图像压缩重构方法”（公开号：CN101877125A，申请号：201110001520.0，申请日：2011年1月6日）中公开了一种基于非凸模型的图像压缩重构方法。该方法对图像作小波变换，得到变换域的系数，对变化域的系数进行傅里叶变换并随机抽取获得压缩后的数据，然后对压缩后的数据采用梯度投影法，通过计算下降方向和下降步长来更新迭代及优化求解，重构变换域的系数，最终对重构后的变换域的系数做逆变换得到重构后的图像。该专利申请存在的不足是，由于采用非凸模型，计算复杂度比较高，并且小波变换具有方向性不足、非冗余性的缺点，在处理高频系数时不能很好的捕捉图像的线性奇异信息，最终导致重构后图像的边缘细节存在一定模糊。

综上所述，由于在信号稀疏表示中使用正交基，信号无法更加稀疏的表示，导致重构出的图像不够准确。因此，压缩感知重构问题的研究主要集中在如何构造更好的稀疏表示字典以及在字典下如何求解稀疏表示系数来精确地恢复原始信号。

发明内容

本发明的目的在于针对现有压缩感知重构技术中在观测数较少的情况下，正交基无法对图像信号进行有效的稀疏表示，导致图像纹理信息难以准确重构的缺点，提出一种基于PCA字典和结构先验信息的压缩感知重构方法，提高重构后图像的质量。

本发明借鉴了迭代滤波和凸集投影的方法，将滤波和凸集投影应用在框架中，可以去除块效应，重构出准确的图像的边缘和纹理信息。其实现步骤如下：

(1)构造黑白图像

分别过大小为21*21的全白图像中心点作直线，生成18个由不同斜率直线分割的图像，直线斜率依次取自角度集合{10*k|k=0,1,2…17}，在每幅分割图像中，将包含图像右下角顶点的一侧区域取值为1，另一侧区域取值为0，得到18个方向的黑白图像；

(2)分别对每个方向的黑白图像采用隔点法选取出所有8*8的块，得到每一个方向的训练样本；

(3)分别对每个方向的训练样本进行PCA分解，得到PCA正交基和特征值，每个方向k分别保留前16个最大特征值和其对应的基，得到特征值矩阵S_k和PCA方向基B_k；

(4)输入测试图像并分成8*8的不重叠块，利用随机高斯观测矩阵A分别对每一个块进行观测得到测量向量y，发送端发送观测矩阵A和每一块的测量向量y，接收端进行接收；

(5)根据观测矩阵A，k方向的PCA方向基B_k和观测向量y构建稀疏表示模型：

$(\widetilde{\mathrm{\α}},\tilde{k})=\underset{\mathrm{\α},k}{\arg \min }\left(\right||{\mathrm{AB}}_k\mathrm{\α}-y{\left|\right|}^2+{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|\mathrm{\α}\left(m\right)\right|}^2}{{\mathrm{\λ}}_k^m}+{\mathrm{\σ}}^2\log \underset{m=1}{\overset{N}{\mathrm{\Π}}}{\mathrm{\λ}}_k^m),$

其中，

为估计的最优块稀疏表示系数，

为估计的方向，α向量为待求解的块稀疏表示系数，σ为估计的噪声方差，α(m)为α系数向量的第m个系数值，

为k方向上第m个最大特征值；

(6)对接收到的每一块观测向量分别在18个PCA方向基上，利用EM算法对稀疏表示模型

求解，得到每一块观测向量的稀疏表示系数和最优方向；

(7)分别将每一块观测向量的最优方向对应的稀疏表示系数与此方向的PCA方向基相乘，得到每一观测向量对应的图像块，将这些图像块依次排列组成图像I⁽¹⁾；

(8)优化图像I⁽¹⁾

8a)以图像I⁽¹⁾中每一图像块为中心，统计其周围相邻8个图像块的方向，如果有7个以上的块方向一致，则保存此方向，查看中心块的方向与此记录的方向是否一致，如果不一致则需用此方向的PCA方向基根据稀疏表示模型来重新求解稀疏表示系数，将系数与其对应的PCA方向基相乘，得到修改的图像块；

8b)用修改的图像块替换步骤（7）中所得图像对应位置的图像块，得到优化的图像I⁽²⁾；

(9)对优化的图像I⁽²⁾做三维块匹配BM3D滤波操作，获得滤波后的图像I⁽³⁾；

(10)对滤波后的图像I⁽³⁾进行不重叠分块，并对这些分块进行凸投影求解，获得优化图像块，将优化的图像块依次排列组成图像I，并输出。

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

第一，本发明提出了在压缩感知领域中使用PCA学习任意方向的方向基，把所有方向学习获得的方向基集成起来就得到了PCA方向基字典。当该字典方向足够多时，它可以更加稀疏和自适应的表示任意方向的图像信号，克服了现有压缩感知重构技术中，正交基无法有效地稀疏表示图像信号的不足，提高了重构图像的质量。

第二,本发明在重构图像时除了利用图像在PCA方向基字典上具有的稀疏约束以外，还充分利用了图像本身存在的局部结构信息，从而提高了图像的重构效果。

第三,本发明将迭代滤波和凸集投影应用在框架中，克服了现有分块压缩感知中存在的块效应问题，改善了图像重构质量不高的缺点，使得本发明的重构图像质量得到了进一步提高。

附图说明

图1是本发明的流程图；

图2是采样率为45%时，本发明与现有技术的仿真对比图；

图3是本发明与现有技术重构出来的Barbara图的峰值信噪比PSNR随采样率变化的趋势图。

具体实施方式

下面结合附图对本发明做进一步描述。

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

步骤一，构造黑白图像。

分别过大小为21*21的全白图像中心点作直线，生成18个由不同斜率直线分割的图像，直线斜率依次取自角度集合{10*k|k=0,1,2…17}，在每幅分割图像中，将包含图像右下角顶点的一侧区域取值为1，另一侧区域取值为0，得到18个方向的黑白图像。

步骤二，得到训练样本。

分别对每个方向的黑白图像采用隔点法选取出所有8*8的块，得到每一个方向的训练样本。

步骤三，获得PCA方向基。

3.1)分别对每个方向的训练样本进行PCA分解，得到PCA正交基和特征值：

3.1a)设得到的k方向的训练样本集合为{f_i}，求出协方差矩阵∑_k为：

${\mathrm{\Σ}}_k=E\left[f_i{f}_i^T\right],$

其中，函数E的作用是求解自变量的数学期望，

为样本中的块f_i的转置；

3.1b)方差矩阵∑_k进行对角化，得到PCA正交基和特征值，即：

${\mathrm{\Σ}}_k=\widetilde{B_k}\widetilde{S_k}{\tilde{B}}_k^T,$

其中，

为k方向的PCA正交基，为k方向的特征值矩阵，

为k方向上第m个最大特征值；

3.2)对每个方向k分别保留前16个最大特征值和其对应的基，得到特征值矩阵S_k和PCA方向基B_k。

步骤四，接收信息。

输入测试图像，并将其分成8*8的不重叠块，利用随机高斯观测矩阵A分别对每一个块进行观测得到测量向量y，发送端发送观测矩阵A和每一块的测量向量y，接收端进行接收。

步骤五，求解最优方向。

5.1)构建稀疏表示模型：

$(\widetilde{\mathrm{\α}},\tilde{k})=\underset{\mathrm{\α},k}{\arg \min }({\left|\right|{\mathrm{AB}}_k\mathrm{\α}-y\left|\right|}^2+{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|\mathrm{\α}\left(m\right)\right|}^2}{{\mathrm{\λ}}_k^m}+{\mathrm{\σ}}^2\log \underset{m=1}{\overset{N}{\mathrm{\Π}}}{\mathrm{\λ}}_k^m),$

其中，

为估计的最优块稀疏表示系数，

为估计的方向，α向量为待求解的块稀疏表示系数，σ为估计的噪声方差，α(m)为α系数向量的第m个系数值，λ_k ^m为k方向上第m个最大特征值，B_k为k方向的PCA方向基，A为观测矩阵，y为观测向量；

5.2）对接收到的每一块观测向量分别在18个PCA方向基上，利用EM算法对稀疏表示模型

求解：

5.2a)固定PCA方向基，利用如下公式求解每一个方向基所在方向上的稀疏表示系数；

$\widetilde{{\mathrm{\α}}_k}=\underset{{\mathrm{\α}}_k}{\arg \min }({\left|\right|A{B}_k{\mathrm{\α}}_k-y\left|\right|}^2+{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|{\mathrm{\α}}_k\left(m\right)\right|}^2}{{\mathrm{\λ}}_k^m}),$

其中，A为观测矩阵，B_k为k方向的PCA方向基，y为块观测向量，σ为估计的噪声方差，向量为估计的k方向上的块稀疏表示系数，α_k向量为要求解的k方向上的块稀疏表示系数，α_k(m)为α_k系数向量的第m个系数值，λ_k ^m为k方向上第m个最大特征值；

5.2b)分别将每一个方向k上求解的稀疏表示系数向量代入如下公式，求解最优方向。

$\tilde{k}=\underset{k}{\arg \min }({\left|\right|{\mathrm{AB}}_k\widetilde{{\mathrm{\α}}_k}-y\left|\right|}^2{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{|\widetilde{{\mathrm{\α}}_k\left(m\right){|}^2}}{{\mathrm{\λ}}_k^m}+{\mathrm{\σ}}^2\log \underset{m=1}{\overset{N}{\mathrm{\Π}}}{\mathrm{\λ}}_k^m),$

其中，

为估计的PCA方向基的方向，

为

系数向量的第m个系数值。

步骤六，得到最优稀疏表示图像。

分别将每一块观测向量的最优方向对应的稀疏表示系数与此方向的PCA方向基相乘，得到每一观测向量对应的图像块，将这些图像块依次排列组成图像I⁽¹⁾。

步骤七，优化图像。

7.1）以图像I⁽¹⁾中每一图像块为中心，统计其周围相邻8个图像块的方向，如果有7个以上的块方向一致，则保存此方向，查看中心块的方向与此记录的方向是否一致，如果不一致则需用此方向的PCA方向基根据稀疏表示模型来重新求解稀疏表示系数，将系数与其对应的PCA方向基相乘，得到修改的图像块；

7.2）用修改的图像块替换步骤六中所得图像对应位置的图像块，得到优化的图像I⁽²⁾。

步骤八，滤波。

对优化的图像做三维块匹配BM3D滤波操作，获得滤波后的图像。对优化的图像I⁽²⁾做三维块匹配BM3D滤波操作，获得滤波后的图像I⁽³⁾。

步骤九，凸投影操作。

9.1)对滤波后的图像I⁽³⁾进行不重叠分块，并对这些分块进行凸投影求解：

9.1a)将滤波后的图像I⁽³⁾按序分成8*8的块，每个分块之间没有重叠的像素；

9.1b)利用如下公式对这些分块进行凸投影求解：

f=f+A^T(y-Af),

其中，f为图像块，A^T为观测矩阵A的转置，y为图像块f对应的观测向量。

9.2)获得优化图像块，将优化的图像块依次排列组成图像I。

步骤十，将步骤九中优化后的图像I输出。

本发明的效果可以通过以下仿真进一步说明。

1.仿真条件：

本发明的仿真在windows XP，SPI，CPU Pentium(R)4，基本频率2.4GHZ，软件平台为MatlabR2007上运行，仿真选用的是512×512的标准Barbara图像。

2.仿真内容与结果：

(1)仿真1：

本仿真中，分别用正交匹配追踪OMP、压缩采样匹配追踪CoSaMP、子空间匹配SP、迭代加权最小平方差IRLS这些现有压缩感知方法和本发明方法对大小为512×512的标准Barbara图像，在采样率均为45%的条件下进行图像重构，重构结果如图2所述。其中图2(a1)是Barbara原始图像，图2(b1)、图2(c1)、图2(d1)、图2(e1)分别是用OMP、CoSaMP、SP、IRLS得到的重构结果图，图2(f1)是本发明的重构结果图。

为了进一步观测图像的细节信息，图2(a2)、图2(b2)、图2(c2)、图2(d2)、图2(e2)、2(f2)分别展示了图2(a1)、图2(b1)、图2(c1)、图2(d1)、图2(e1)、图2(f1)中裤子纹理的局部放大图。上述现有压缩感知方法使用的稀疏表示字典是正交余弦字典，实验中这些方法和本发明方法都只使用了一次滤波和一次凸投影操作。

由图2(b1)、图2(c1)、图2(d1)、图2(e1)可以看出使用现有技术重构后的图像是模糊的，尤其是纹理部分重构效果有待提高。更进一步从裤子纹理局部放大图2(b2)、图2(c2)、图2(d2)、图2(e2)可以看出裤子纹理部分相当模糊，明显失去了原图像均匀纹理的效果。

由图2(f1)可以看出本发明的重构结果图接近原始图像的清晰度，进一步由其对应的裤子纹理局部放大图2(f2)可以清晰的看出本发明提出的重构方法能够得到清晰的纹理。

(2)仿真2：

本仿真中，用现有的OMP、CoSaMP、SP、IRLS方法和本发明方法分别在采样率为25%、30%、35%、40%、45%时对Barbara图像进行仿真，得到精确的峰值信噪比PSNR，如表1所示。

表1  各采样率下PSNR值

从表1数据可以看出，本发明的方法在采样率为25%、30%、35%、40%、45%下得到的结果图的峰值信噪比PSNR都要高于其他四种方法得到的PSNR，即本发明的方法的重构图像质量要比其他四种方法高。

根据表1数据得到OMP、CoSaMP、SP、IRLS方法和本发明方法重构出的Barbara图的PSNR随采样率变化的趋势图，其结果如图3所示，图3中的横坐标表示采样率，纵坐标表示峰值信噪比PSNR值。

由图3可以看出，本发明方法得到的重构结果图的PSNR值明显高于其他四种方法。

综上，本发明能够很好地重构图像的纹理和边缘部分，得到清晰的图像，与现有的其他重构方法相比，本发明提高了图像的重构质量。

Claims (4)

1.一种基于PCA字典和结构先验信息的压缩感知重构方法，包括如下步骤：

(1)构造黑白图像

分别过大小为21*21的全白图像中心点作直线，生成18个由不同斜率直线分割的图像，直线斜率依次取自角度集合{10*k|k=0,1,2…17}，在每幅分割图像中，将包含图像右下角顶点的一侧区域取值为1，另一侧区域取值为0，得到18个方向的黑白图像；

(2)分别对每个方向的黑白图像采用隔点法选取出所有8*8的块，得到每一个方向的训练样本；

(3)分别对每个方向的训练样本进行PCA分解，得到PCA正交基和特征值，每个方向k分别保留前16个最大特征值和其对应的基，得到特征值矩阵S_k和PCA方向基B_k；

(4)输入测试图像并分成8*8的不重叠块，利用随机高斯观测矩阵A分别对每一个块进行观测得到测量向量y，发送端发送观测矩阵A和每一块的测量向量y，接收端进行接收；

(5)根据观测矩阵A，k方向的PCA方向基B_k和观测向量y构建稀疏表示模型：

$(\widetilde{\mathrm{\α}},\tilde{k})=\underset{\mathrm{\α},k}{\arg \min }({\left|\right|{\mathrm{AB}}_k\mathrm{\α}-y\left|\right|}^2+{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|\mathrm{\α}\left(m\right)\right|}^2}{{\mathrm{\λ}}_k^m}{\mathrm{\σ}}^2\log \underset{m=1}{\overset{N}{\mathrm{\Π}}}{\mathrm{\λ}}_k^m),$

其中，

为估计的最优块稀疏表示系数，

为估计的方向，α向量为待求解的块稀疏表示系数，σ为估计的噪声方差，α(m)为α系数向量的第m个系数值，

为k方向上第m个最大特征值；

(6)对接收到的每一块观测向量分别在18个PCA方向基上，利用EM算法对稀疏表示模型

求解，得到每一块观测向量的稀疏表示系数和最优方向；

(7)分别将每一块观测向量的最优方向对应的稀疏表示系数与此方向的PCA方向基相乘，得到每一观测向量对应的图像块，将这些图像块依次排列组成图像I⁽¹⁾；

(8)优化图像I⁽¹⁾

8a)以图像I⁽¹⁾中每一图像块为中心，统计其周围相邻8个图像块的方向，如果有7个以上的块方向一致，则保存此方向，查看中心块的方向与此记录的方向是否一致，如果不一致则需用此方向的PCA方向基根据稀疏表示模型来重新求解稀疏表示系数，将系数与其对应的PCA方向基相乘，得到修改的图像块；

8b)用修改的图像块替换步骤（7）中所得图像对应位置的图像块，得到优化的图像I⁽²⁾；

(9)对优化的图像I⁽²⁾做三维块匹配BM3D滤波操作，获得滤波后的图像I⁽³⁾；

(10)对滤波后的图像I⁽³⁾进行不重叠分块，并对这些分块进行凸投影求解，获得优化图像块，将优化的图像块依次排列组成图像I，并输出。

2.根据权利要求1所述的基于PCA字典和结构先验信息的压缩感知重构方法，其中，步骤3)所述的对每个方向的训练样本进行PCA分解，按如下步骤进行：

3a)设得到的k方向的训练样本集合为{f_i}，求出协方差矩阵∑_k为：

${\mathrm{\Σ}}_k=E\left[f_i{f}_i^T\right],$

其中，函数E的作用是求解自变量的数学期望，

为样本中的块f_i的转置；

3b)对协方差矩阵∑_k进行对角化，得到PCA正交基和特征值，即：

${\mathrm{\Σ}}_k={\tilde{B}}_k{\tilde{S}}_k{\tilde{B}}_k^T,$

其中，为k方向的PCA正交基，

为k方向的特征值矩阵，

为k方向上第m个最大特征值。

3.根据权利要求1所述的基于PCA字典和结构先验信息的压缩感知重构方法，其中步骤（6）中所述的对接收到的每一块观测向量分别在18个PCA方向基上，利用EM算法对稀疏表示模型求解，按如下步骤进行：

6a)固定PCA方向基，利用如下公式求解每一个方向基所在方向上的稀疏表

示系数；

${\widetilde{\mathrm{\α}}}_k=\underset{{\mathrm{\α}}_k}{\arg \min }({\left|\right|{\mathrm{AB}}_k{\mathrm{\α}}_k-y\left|\right|}^2+{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|{\mathrm{\α}}_k\left(m\right)\right|}^2}{{\mathrm{\λ}}_k^m}),$

其中，A为观测矩阵，B_k为k方向的PCA方向基，y为块观测向量，σ为估计的噪声方差，

向量为估计的k方向上的块稀疏表示系数，α_k向量为要求解的k方向上的块稀疏表示系数，α_k(m)为α_k系数向量的第m个系数值，λ_k ^m为k方向上第m个最大特征值；

6b)分别将每一个方向k上求解的稀疏表示系数向量

代入如下公式，求解最优方向。

$\tilde{k}=\underset{k}{\arg \min }({\left|\right|{\mathrm{AB}}_k{\widetilde{\mathrm{\α}}}_k-y\left|\right|}^2+{\mathrm{\σ}}^2\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|{\widetilde{\mathrm{\α}}}_k\left(m\right)\right|}^2}{{\mathrm{\λ}}_k^m}+{\mathrm{\σ}}^2\log \underset{m=1}{\overset{N}{\mathrm{\Π}}}{\mathrm{\λ}}_k^m),$

其中，

为估计的PCA方向基的方向，

为

系数向量的第m个系数值。

4.根据权利要求1所述的基于PCA字典和结构先验信息的压缩感知重构方法，其中，步骤(10)所述的对滤波后的图像I⁽³⁾进行不重叠分块，并对这些分块进行凸投影求解，按如下步骤进行：

10a)将滤波后的图像I⁽³⁾按序分成8*8的块，每个分块之间没有重叠的像素；

10b)利用如下公式对这些分块进行凸投影求解：

f=f+A^T(y-Af)

其中，f为图像块，A^T为观测矩阵A的转置，y为图像块f对应的观测向量。

Images