CN103400349A - 基于盲压缩感知的图像重构方法 - Google Patents

Abstract

本发明公开了一种基于盲压缩感知模型的图像重构方法，主要解决传统压缩感知只能对稀疏信号进行观测，且重构图像质量较差的问题。其实现步骤是：(1)对输入图像进行冗余变换，获得冗余矩阵；(2)对冗余矩阵在观测矩阵下进行压缩观测；(3)根据压缩观测的结果，利用OMP算法进行自适应更新稀疏矩阵；（4）根据更新后的稀疏矩阵，利用奇异值分解方法更新稀疏基；（5）将更新后的稀疏矩阵和稀疏基相乘得到重构后的图像冗余矩阵；（6）将重构后的图像冗余矩阵通过冗余逆变换获得重构图像；通过图像的峰值信噪比对重构图像进行评估。本发明具有重构图像质量高，对噪声抑制效果好的优点，可应用于图像去噪和图像压缩。

Description

基于盲压缩感知的图像重构方法

技术领域

本发明属于图像处理技术领域，特别涉及压缩感知图像重构方法，可用于图像去噪和图像压缩。

背景技术

压缩感知理论能够从少量的观测数据中恢复信号，已经被应用到图像重构领域，而在实际应用中，由于处理的图像信号往往不具有稀疏性，传统的压缩感知理论在稀疏化时常常需要选择固定的稀疏基，现有的方法如基于小波变换的压缩感知技术、基于离散余弦变换的压缩感知技术，即先在固定的稀疏基下稀疏化，再对稀疏化信号进行观测，而且观测信号的信号越稀疏，恢复重构的信号效果越好，但这种固定稀疏基的适用范围有限，即无法满足任何对象压缩感知的稀疏表示，而且对噪声的抑制作用有限。因此，如何选择一种全自动的稀疏表达是压缩感知图像重建面临的重要问题，该项研究具有重要的实际应用价值。

盲压缩感知（Blind Compressed Sensing，BCS）是一种新的压缩感知模型，它可以直接对非稀疏信号进行观测，重构时能对稀疏基进行自适应更新，这一思想的先进性在于信号的非稀疏观测以及信号重构时稀疏基的自动更新，但该模型的初始稀疏基的选择形式要求严格，必须具有分块对角形式，无法直接对含有复杂信息的图像信号进行重构。

发明内容

本发明的目的在于提出一种基于盲压缩感知的图像重构方法，以解决压缩观测时信号必须是稀疏的问题和重构时初始基的严格限制问题，将盲压缩感知直接用于对含有复杂信息的图像重构。

实现本发明目的的技术方案是：利用盲压缩感知模型，将合成的离散余弦变换基作为初始基，利用正交匹配追踪OMP算法挑选图像的稀疏系数，利用奇异值分解的方法更新初始合成的离散余弦变换基，进而得到重构图像，具体步骤包括如下：

（1）设输入图像矩阵为I、冗余矩阵为X，对输入图像I进行冗余变换,获得冗余矩阵X；

（2）设观测矩阵为A，根据盲压缩感知模型得出对图像信号重构的模型为：

$\min {\left|\right|B-\mathrm{AX}\left|\right|}_F^2=\underset{P,S}{\min }{\left|\right|B-\mathrm{APS}\left|\right|}_F^2$

其中P、S分别为X的稀疏基、稀疏矩阵，||||_F表示求矩阵的F-范数，利用A对冗余矩阵X进行压缩观测，得到观测后的矩阵B＝A*X；

（3）设冗余矩阵X的初始化稀疏基为P、稀疏矩阵为S，利用OMP算法更新稀疏矩阵S，使其满足：

$\underset{S}{\min }{\left|\right|B-\mathrm{DS}\left|\right|}_F^2,s.t{\left|\right|s_i\left|\right|}_0\≤k$

其中s_i为S的第i列，且i取S的任意一列，k为s_i的稀疏度，D＝A*P，*为乘号，||||_F表示求矩阵的F-范数；

（4）利用奇异值分解方法更新稀疏基P，使其满足：

$\underset{P}{\min }{\left|\right|B-\mathrm{APS}\left|\right|}_F^2;$

（5）将更新后的稀疏基 $\left[\begin{array}{ccccc}{P}_1& & & & \\ & \·\·\·& & & \\ & & P_i^{\′}& & \\ & & & \·\·\·& \\ & & & & P_{2L}\end{array}\right]$ 和更新后的稀疏矩阵S相乘，得到重构后的冗余矩阵X_r，即X_r＝P*S，其中P_i'为更新后的P的子块，i=1…2L；

（6）对重构后的冗余矩阵X_r进行冗余逆变换，求出重构后的图像矩阵I_r；

（7）采用图像的峰值信噪比PSNR作为评价指标对重构图像I_r进行评估，若重构图像的PSNR大于等于峰值信噪比阈值v，则输出重构图像I_r，若不满足返回（3）重复执行，该峰值信噪比阈值v根据不同图像和加入噪声值的不同设定。

本发明与传统压缩感知（Compressed Sensing，CS）方法相比具有如下优点：

1.本发明采用稀疏矩阵、稀疏系数交替更新的思想，在优化迭代中自适应调节稀疏矩阵，使重构图像对噪声具有更强的鲁棒性，避免了传统CS图像重构中对稀疏化要求高的缺陷；

2.本发明将盲压缩感知模型应用到图像观测重构中，解决了传统压缩感知不能对非稀疏图像信号进行观测重构问题，能直接对非稀疏的图像信号进行观测重构。

以下结合附图和实施例对本发明作进一步详细描述。

附图说明

图1是BCS图像重构流程图；

图2是本发明方法与现有Wavelet_OMP方法对Barbara图像处理的结果对比图；

图3是本发明方法与现有Wavelet_OMP方法对Desert图像处理的结果对比图；

图4是本发明方法与现有Wavelet_OMP方法对人脑图像处理的结果对比图。

具体实施方式

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

步骤1，对输入图像I进行冗余变换,获得冗余矩阵X：

（1a）从输入图像I的第i行第j列的元素开始，按上到下从左到右顺序取出p*q的图像块，每次取块时i依次加1而j不变，等到i加到c-p+1时，j再加1，此时i变为1，如此反复执行，直到i、j分别加到c-p+1、d-q+1为止，其中c为输入图像I的行数，d为输入图像I的列数；

（1b）将取出的图像块的每一列按顺序首尾相接得到列向量，并将所有的列向量再进行顺序组合，得到冗余矩阵X。

步骤2，对冗余矩阵X进行压缩观测：

设观测矩阵为A，并将该观测矩阵A与冗余矩阵X相乘，得到压缩观测后的矩阵B为：B＝A*X，

其中A＝[A₁,...,A_j,A_j+1...,A_2L]且A∈R^n*nL，R^n*nL表示n*nL实矩阵，n表示A的行数，L为A的压缩比，即A的列数与行数之比，[A_j A_j+1]为n*n维服从高斯分布的正交单位化矩阵，其中A_j、A_j+1分别为观测矩阵A的第j、第j+1块，j取1到2L的任意奇数。

步骤3，设冗余矩阵X的初始化稀疏基为P、稀疏矩阵为S，利用OMP算法更新稀疏矩阵S，使其满足：

$\underset{S}{\min }{\left|\right|B-\mathrm{DS}\left|\right|}_F^2,s.t{\left|\right|s_i\left|\right|}_0\≤k$

其中s_i为S的第i列，且i取S的任意一列，k为s_i的稀疏度，D＝A*P,*为乘号，||||_F表示求矩阵的F-范数，||||₀表示求向量的零范数。

本步骤除了可以采用OMP算法更新稀疏矩阵S外，还可以采用MP算法、BP算法对稀疏矩阵S进行更新，这些算法都是现有技术中成熟的常用算法。

步骤4，利用奇异值分解方法更新稀疏基P，使其满足

（4a）将稀疏基P分为大小为

的子块，即 $P=\left[\begin{array}{ccc}{P}_1& & \\ & \·\·\·& \\ & & P_{2L}\end{array}\right],$ P∈R^nL*nL，且令P_i为稀疏基P的任意子块，i=1…2L，将稀疏矩阵S分为大小为

的子块，即 $S=\left(\begin{array}{c}{S}_1\\ \·\·\·\\ S_{2L}\end{array}\right),$ S∈R^nL*e，其中n为观测矩阵A的列数，L为观测矩阵A的压缩比，e为冗余矩阵X的列数；

（4b）对P中的子块P_i初始化：利用离散余弦基公式DCT产生m₁*m₁的第一变换阵DC1和m₂*m₂的第二变换阵DC2，其中

并对两个变换阵求K积，得到初始化后的稀疏基P的子块P_i，其中i=1…2L，L为A的压缩比，

所述的离散余弦基公式DCT为：

$\mathrm{DCT}=\left\{\cos \right[\left(\begin{array}{c}0\\ \·\·\·\\ m-1\end{array}\right)*\mathrm{\π}*\left(\begin{array}{ccc}0& ...& m-1\end{array}\right){/m\left]\right\}}^T*\mathrm{diag}{(1./\mathrm{diag}(\left\{\cos \right[\left(\begin{array}{c}0\\ \·\·\·\\ m-1\end{array}\right)*\mathrm{\π}*\begin{array}{}\end{array}\left(\begin{array}{ccc}0& ...& m-1\end{array}\right)/m\left]\right\}}^T*\left\{\cos \right[\left(\begin{array}{c}0\\ \·\·\·\\ m-1\end{array}\right)*\mathrm{\π}*\begin{array}{}\end{array}\left(\begin{array}{ccc}0& ...& m-1\end{array}\right)/m\left]\right\}\left)\right)$

其中，diag()表示取矩阵的对角戏元素，()^T表示矩阵或向量的转置，m表示最后产生的矩阵维数；

（4c）依次对初始化后的稀疏基P的子块P_i进行更新，更新后的P的子块为P_i'，其中i=1…2L：

（4c1）令

其中，S_j为稀疏矩阵S的第j块，A_j为观测矩阵A的第j块，P_j为稀疏基P的第j块，j为1到2L，j不等于i，B为压缩观测后的矩阵；

（4c2）对S_i(B_i)^TA_i进行奇异值分解，得到分解后的结果为：S_i(B_i)^TA_i＝U*M*V^T，其中S_i为稀疏矩阵S的第i块，A_i为A的第i块，()^T为求矩阵的转置，U、M、V为分解得到的三个矩阵；

（4c3）将稀疏基P的第i个子块更新为：P_i'＝V*U^T。

步骤5，将更新后的稀疏基 $\left[\begin{array}{ccccc}{P}_1& & & & \\ & \·\·\·& & & \\ & & P_i^{\′}& & \\ & & & \·\·\·& \\ & & & & P_{2L}\end{array}\right]$ 和更新后的稀疏矩阵S相乘，得到重构后的冗余矩阵X_r，即X_r＝P*S，其中P_i'为更新后的P的子块，i=1…2L。

步骤6，对重构后的冗余矩阵X_r进行冗余逆变换，求出重构后的图像矩阵I_r：

（6a）将重构后的冗余矩阵X_r从第一列开始顺序取出X_r的列，将取出的列每隔p个元素截断，这样就将X_r一列变为q个列向量，每个列向量含有p个元素；再将q个列向量按截断出的先后顺序组合成p*q的图像块矩阵Q_h，h等于1到e，e为冗余矩阵X_r的列数；

（6b）建立一个c*d的空矩阵I_r0，将变换得到的图像块矩阵Q_h按从上到下、从左到右的顺序组合后，放入矩阵I_r0中，其中每次放入的行号加1，且保持列号不变，当行号加到d-q+1时，再从第一行开始，此时列号加1，如此反复执行，直到行号加到c-p+1和列号加到d-q+1停止；执行过程中每次放入p*q图像块矩阵Q_h根据p、q取值不同可能会与已经放入的元素在位置上重合，对于位置上重合的元素按加和取平均计算得到的c*d矩阵，即为重构后的图像矩阵I_r。

步骤7，根据不同图像和加入噪声值设定不同设定峰值信噪比阈值v，将图像的峰值信噪比PSNR作为评价指标对重构图像I_r进行评估，若重构图像的PSNR大于等于峰值信噪比阈值v，则输出重构图像I_r，反之返回步骤3重复执行。

本发明的优点由以下仿真数据和图像进一步说明：

1．仿真条件

仿真中选用256*256的图像为实验对象，块变换选取块大小为8*8，这样，观测阵A选为32*64，稀疏基矩阵P为64*64，其中用到的DCT基为2*2与8*8的矩阵做K式积，L=2为压缩率,A、P、S相应分为4块，循环迭代评价指标为图像的峰值信噪比PSNR，且阈值为20dB。

2．仿真内容

本实验主要目的是对比本发明与现有Wavelet_OMP方法的图像重构效果，对输入图像分别加入不同程度的噪声，即加入不同的噪声方差，选用Barbara图像、Desert图像、人脑图像做了仿真实验：

仿真1，分别用本发明方法和现有的Wavelet_OMP方法对Barbara图像进行重构，如图2，其中：

图2a为噪声方差sigma为10时的重构前的图像，

图2b为为噪声方差sigma为20时的重构前的图像，

图2c为噪声方差sigma为30时的重构前的图像，

图2d为噪声方差sigma为40时的重构前的图像，

图2e为噪声方差sigma为10时用本发明方法重构出的图像，

图2f为噪声方差sigma为20时用本发明方法重构出的图像，

图2g为噪声方差sigma为30时用本发明方法重构出的图像，

图2h为噪声方差sigma为40时用本发明方法重构出的图像，

图2m为噪声方差sigma为10时用现有Wavelet_OMP方法重构出的图像，

图2n为噪声方差sigma为20时用现有Wavelet_OMP方法重构出的图像，

图2o为噪声方差sigma为30时用现有Wavelet_OMP方法重构出的图像，

图2p为噪声方差sigma为40时用现有Wavelet_OMP方法重构出的图像。

表1为分别用本发明方法和现有Wavelet_OMP方法，在噪声方差sigma分别取10、20、30、40时，对Barbara图像的处理前后的图像PSNR值对比结果。

表1    Barbara实验

仿真2，分别用本发明方法和Wavelet_OMP方法对Desert图像进行图像重构，如图3，其中：

图3a为噪声方差sigma为10时的重构前的图像，

图3b为为噪声方差sigma为20时的重构前的图像，

图3c为噪声方差sigma为30时的重构前的图像，

图3d为噪声方差sigma为40时的重构前的图像，

图3e为噪声方差sigma为10时用本发明方法重构出的图像，

图3f为噪声方差sigma为20时用本发明方法重构出的图像，

图3g为噪声方差sigma为30时用本发明方法重构出的图像，

图3h为噪声方差sigma为40时用本发明方法重构出的图像，

图3m为噪声方差sigma为10时用现有Wavelet_OMP方法重构出的图像，

图3n为噪声方差sigma为20时用现有Wavelet_OMP方法重构出的图像，

图3o为噪声方差sigma为30时用现有Wavelet_OMP方法重构出的图像，

图3p为噪声方差sigma为40时用现有Wavelet_OMP方法重构出的图像

表2为分别用本发明方法和Wavelet_OMP方法，在噪声方差sigma分别取10、20、30、40下，对Desert图像的处理前后的图像PSNR值对比。

表2    Desert实验

仿真3，分别用本发明方法和Wavelet_OMP方法对人脑图像进行图像重构，结果如图4，其中：

图4a为噪声方差sigma为10时的重构前的图像，

图4b为为噪声方差sigma为20时的重构前的图像，

图4c为噪声方差sigma为30时的重构前的图像，

图4d为噪声方差sigma为40时的重构前的图像，

图4e为噪声方差sigma为10时用本发明方法重构出的图像，

图4f为噪声方差sigma为20时用本发明方法重构出的图像，

图4g为噪声方差sigma为30时用本发明方法重构出的图像，

图4h为噪声方差sigma为40时用本发明方法重构出的图像，

图4m为噪声方差sigma为10时用现有Wavelet_OMP方法重构出的图像，

图4n为噪声方差sigma为20时用现有Wavelet_OMP方法重构出的图像，

图4o为噪声方差sigma为30时用现有Wavelet_OMP方法重构出的图像，

图4p为噪声方差sigma为40时用现有Wavelet_OMP方法重构出的图像

表3为分别用本发明方法和Wavelet_OMP方法对磁共振成像MRI人脑图像的处理前后的图像PSNR值对比表，其噪声方差sigma分别取10、20、30、40。

表3    人脑图像实验

3.仿真结果分析

从表1和表2都可以看出，用本发明方法重构图像质量比现有Wavelet_OMP方法在PSNR上有至少2dB的提高。

从表3可以看出，用本发明方法重构图像质量比现有Wavelet_OMP方法在PSNR上有至少0.3dB的提高。

图2、图3、图4表明，Wavelet_OMP方法对抗噪声能力较弱，而本发明方法有一定的去噪功能。并且从视觉效果上看，在不同的噪声背景下，本发明方法的重构质量噪声很少，而Wavelet OMP方法噪声明显。

综上，本发明方法结合盲压缩感知模型稀疏矩阵、稀疏系数交替更新的思想，使得本发明方法比Wavelet_OMP方法重构图像质量更好，且对噪声具有更好的抵抗能力。

Claims (6)

1.一种基于盲压缩感知的图像重构方法，包括以下步骤：

（1）设输入图像矩阵为I、冗余矩阵为X，对输入图像I进行冗余变换,获得冗余矩阵X

（2）设观测矩阵为A，根据盲压缩感知模型得出对图像信号重构的模型为：

$\min {\left|\right|B-\mathrm{AX}\left|\right|}_F^2=\underset{P,S}{\min }{\left|\right|B-\mathrm{APS}\left|\right|}_F^2$

其中P、S分别为X的稀疏基、稀疏矩阵，||||_F表示求矩阵的F-范数，利用A对冗余矩阵X进行压缩观测，得到观测后的矩阵B＝A*X；

（3）设冗余矩阵X的初始化稀疏基为P、稀疏矩阵为S，利用OMP算法更新稀疏矩阵S，使其满足：

$\underset{S}{\min }{\left|\right|B-\mathrm{DS}\left|\right|}_F^2,s.t{\left|\right|s_i\left|\right|}_0\≤k$

其中s_i为S的第i列，且i取S的任意一列，k为s_i的稀疏度，D＝A*P,*为乘号，||||_F表示求矩阵的F-范数，||||₀表示求向量的零范数;

（4）利用奇异值分解方法更新稀疏基P，使其满足：

$\underset{P}{\min }{\left|\right|B-\mathrm{APS}\left|\right|}_F^2;$

（5）将更新后的稀疏基 $\left[\begin{array}{ccccc}{P}_1& & & & \\ & \·\·\·& & & \\ & & P_i^{\′}& & \\ & & & \·\·\·& \\ & & & & P_{2L}\end{array}\right]$ 和更新后的稀疏矩阵S相乘，得到重构后的冗余矩阵X_r，即X_r＝P*S，其中P_i'为更新后的P的子块，i=1…2L；

（6）对重构后的冗余矩阵X_r进行冗余逆变换，求出重构后的图像矩阵I_r；

（7）采用图像的峰值信噪比PSNR作为评价指标对重构图像I_r进行评估，若重构图像的PSNR大于等于峰值信噪比阈值v，则输出重构图像I_r，若不满足返回（3）重复执行，该峰值信噪比阈值v根据不同图像和加入噪声值的不同设定。

2.根据权利要求1所述的基于盲压缩感知的图像重构方法，其中步骤（1）所述的对输入图像I进行冗余变换得到冗余矩阵X，按如下步骤进行：

（1a）从输入图像I的(i,j)像素点开始，取p*q的图像块，即取出从i到i+p-1行、从j到j+q-1列的图像块，其中i从1到c-p+1，j从1到d-q+1，c为输入图像I的行数，d为输入图像I的列数；

（1b）将取出的图像块的每一列按顺序首尾相接得到列向量，并将所有的列向量再进行顺序组合，得到冗余矩阵X。

3.根据权利要求1所述的基于盲压缩感知的图像重构方法，其中所述步骤（2）中对观测矩阵A进行分块，表示为：

A＝[A₁,...,A_j,A_j+1...,A_2L]，其中A∈R^n*nL，R^n*nL表示A属于n*nL实矩阵，[A_j A_j+1]为n*n的服从高斯分布的正交单位化矩阵，其中A_j、A_j+1分别为观测矩阵A的第j、j+1块，j取1到2L的任意奇数，L为A的压缩比，即A的块数。

4.根据权利要求1所述的基于盲压缩感知的图像重构方法，其中步骤（4）所述的利用奇异值分解方法更新稀疏基P，按如下步骤进行：

（4a）将稀疏基P分为大小为

的子块，即 $P=\left[\begin{array}{ccc}{P}_1& & \\ & \·\·\·& \\ & & P_{2L}\end{array}\right],$ P∈R^nL*nL；将稀疏矩阵S分为大小为

的子块，即 $S=\left(\begin{array}{c}{S}_1\\ \·\·\·\\ S_{2L}\end{array}\right),$ S∈R^nL*d，其中n为观测矩阵A的列数，L为观测矩阵A的压缩比，e为冗余矩阵X的列数；

（4b）对P中的子块P_i初始化，即利用离散余弦基公式DCT产生m₁*m₁的第一变换阵DC1和m₂*m₂的第二变换阵DC2，其中

并对两个变换阵求K积，得到P_i，其中i=1…2L，L为A的压缩比：

（4c）依次对P_i进行更新：

首先，令

其中，S_j为稀疏矩阵S的第j块，A_j为观测矩阵A的第j块，P_j为稀疏基P的第j块，j为1到2L，j不等于i；

然后，对S_i(B_i)^TA_i进行奇异值分解，得到分解后的结果为：S_i(B_i)^TA_i＝U*M*V^T，其中S_i为稀疏矩阵S的第i块，A_i为A的第i块，()^T为求矩阵的转置，U、M、V为分解得到的三个矩阵；

最后，将P_i更新为：P_i＝V*U^T。

5.根据权利要求4所述的基于盲压缩感知的图像重构方法，其中所述步骤（4b）中的离散余弦基公式DCT，表示为：

$\mathrm{DCT}=\left\{\cos \right[\left(\begin{array}{c}0\\ \·\·\·\\ m-1\end{array}\right)*\mathrm{\π}*\left(\begin{array}{ccc}0& ...& m-1\end{array}\right){/m\left]\right\}}^T*\mathrm{diag}{(1./\mathrm{diag}(\left\{\cos \right[\left(\begin{array}{c}0\\ \·\·\·\\ m-1\end{array}\right)*\mathrm{\π}*\begin{array}{}\end{array}\left(\begin{array}{ccc}0& ...& m-1\end{array}\right)/m\left]\right\}}^T*\left\{\cos \right[\left(\begin{array}{c}0\\ \·\·\·\\ m-1\end{array}\right)*\mathrm{\π}*\begin{array}{}\end{array}\left(\begin{array}{ccc}0& ...& m-1\end{array}\right)/m\left]\right\}\left)\right)$

其中，diag()表示取矩阵的对角线元素，()^T表示矩阵或向量的转置，m表示最后产生的矩阵维数。

6.根据权利要求1所述的一种基于盲压缩感知的图像重构方法，其中步骤（6）所述的对重构后的冗余矩阵X_r进行冗余逆变换求出重构后的图像矩阵I_r，具体操作如下：

（6a）将重构后的冗余矩阵X_r从第一列开始顺序取出X_r的列，将取出的列每隔p个元素截断，这样就将X_r一列变为q个列向量，每个列向量含有p个元素，再将q个列向量按截断出的先后顺序组合成p*q的图像块；

（6b）建立一个c*d的空矩阵I_r0，将变换得到的图像块按从上到下、从左到右的顺序组合后，放入矩阵I_r0中，其中每次放入的行号加一，且保持列号不变，当行号加到d-q+1时，再从第一行开始，此时列号加一，如此反复执行，直到行号加到c-p+1和列号加到d-q+1停止，执行过程中行列重合部分的元素按加和取平均计算，最后得到的c*d矩阵，即为重构后的图像矩阵I_r。

Images