CN104299201A - 基于遗传稀疏优化和贝叶斯估计模型的图像重构方法 - Google Patents

Abstract

本发明公开了一种基于遗传稀疏优化和贝叶斯估计模型的图像重构方法,主要解决现有方法对压缩感知图像的分块重构中，存在的边界模糊和块效应明显的问题。本发明将图像块分为光滑和非光滑两类，分别对其建模重构；根据光滑块的统计特性，直接采用伪逆解对其直流分量和变化分量进行快速重构；对于非光滑块，利用遗传算法在PCA字典中选择一组原子来对其进行最优重构。实验结果表明，利用本发明重构的图像比传统的正交匹配追踪方法OMP和统计压缩感知方法SCS具有更好的边界和区域一致性，细节信息更加清晰，并且块效应明显减少，可用于低采样率下图像获取中的重构。

Description

基于遗传稀疏优化和贝叶斯估计模型的图像重构方法

技术领域

本发明属于图像处理技术领域，特别是一种图像重构方法，可用于解决自然图像的分块压缩感重构问题。

背景技术

在图像处理的发展历程中，压缩感知的提出使得图像获取技术有了长足的进步，它仅通过少量的观测值即可重构出原始的图像，使得图像获取的效率大大提高。这对于国民建设和军事发展有着重要的意义。

近年来，解决压缩感知图像重构的算法层出不穷，最常用的是基于贪婪追踪技术的迭代方法。其基本思想是，在每次迭代中通过选择与残差相关性最大的原子，使观测误差快速收敛。代表性的算法有匹配追踪、正交匹配追踪、子空间追踪和正则化正交匹配追踪等。但是由于其方法过于贪婪，并且回溯能力有限，往往容易陷入局部最优解，使重构结果失真。

图像重构的另一个重要的方面是字典的设计，根据结构和方向的特性，许多经典的字典被设计出来，例如小波基字典、Ridgelets、Curvelets、Bandelets和Contourlets等等。随着需求的日益复杂，利用传统字典对图像建模的局限性也越见明显，为了打破这种限制，字典学习的概念被提出。许多字典学习的策略应运而生，比如K-SVD，MOD等等。通过对训练样本的学习，使得字典能够自适应地获得图像的结构信息，并且使字典稀疏表示图像的能力更强。

近年来，Yu等人在文献“G.Yu and G.Sapiro,Statistical Compressed Sensing of GaussianMixture Models”中提出了基于高斯混合模型的统计压缩感知，从统计学的角度求解图像重构的逆问题。与传统压缩感知的稀疏性假设不同，该方法旨在重构一组服从某一特定分布的信号，所以信号分布的估计是该方法的关键问题。在基于混合高斯模型的统计压缩感知中，假设每个信号服从多个高斯分布中的一个，在每一次迭代中，该方法首先假定信号的高斯分布已知，对所有信号进行估计和分类；然后，假定信号已知，用属于同一类高斯分布的所有信号来估计对应高斯分布的参数，通过这样一个交替的迭代过程，使得信号可以高概率的重构出来。从PCA域分析，该方法其实也包含字典更新的过程。但是，该方法在实际应用中面临着一个重要的问题：当样本数量不足时，这种更新会破坏PCA字典的结构，使重构结果变差。

发明内容

本发明提出一种基于遗传稀疏优化和贝叶斯估计模型的图像重构方法，其目的在于构建一个有效的压缩感知图像重构框架，充分发掘图像的结构特性，并根据不同的结构特性采用针对性的重构模型，获得更加有效的重构结果，并且减少计算资源的浪费。

实现本发明目的的技术方案是：将图像块分为光滑和非光滑两类，分别对其建模重构。根据光滑块的统计特性，直接采用伪逆解对其直流分量和变化分量进行快速重构。对于非光滑块，利用遗传算法在PCA字典中选择一组原子来对其进行最优重构。其具体步骤包括如下：

(1)输入大小为像素的图像块x的观测向量y，估计该图像块的直流分量和变化分量

(2)根据变化分量的能量，将图像块x进行光滑和非光滑的分类，若该图像块属于光滑块，则执行步骤(3)操作；否则，执行步骤(4)操作；

(3)对于光滑块，直接用直流分量加上变化分量进行重构，得到重构结果为：

$\tilde{x}={\tilde{x}}_{\mathrm{dc}}+{\tilde{x}}_{\mathrm{var}};$

(4)对于非光滑块，利用遗传算法在PCA字典中选择合适的原子组合对该图像块进行重构：

(4.1)构造一个包含18个方向子字典的PCA字典B；

(4.2)初始化一个父代种群种群规模S＝20，个体的长度为16，置计数器t＝0；

(4.3)计算父代种群中每个个体的适应度值，记录适应度值最大的个体为最优个体z_best；

(4.4)对父代种群进行概率为0.8的交叉操作，得到子代种群

(4.5)对子代种群进行概率为0.2的变异操作；

(4.6)计算子代种群中每个个体的适应度值，如果中个体的最大适应度值大于当前最优个体z_best的适应度值，则更新最优个体z_best；

(4.7)在父代种群和子代种群的集合中，选择20个适应度值最大的个体作为新的父代种群并且记录它们的适应度值，置计数器t＝t+1；

(4.8)设最大迭代次数T_max为100，若t<T_max，则返回步骤(4.4)；否则，输出最佳个体z_best的重构结果。

本发明由于集成了遗传算法，稀疏表示，概率统计这些技术手段，因而具有如下优点：

(A)对于光滑块的重构时间极短。

(B)PCA字典和遗传算法的运用充分抓住了非光滑块的结构和方向特性。

(C)重构效果优于传统的正交匹配追踪方法OMP和统计压缩感知方法SCS。

实验证明，与传统压缩感知重构算法OMP和统计压缩感知方法SCS相比，本发明重构的图像具有更好的区域和边界一致性。对于图像中复杂信息的重构，本发明具有更加清晰的视觉效果和更小的人工块效应。

附图说明

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

图2是本发明与现有方法对于4个具有明显方向信息的图像块的重构效果对比；

图3是本发明与现有方法对于4个具有复杂结构信息的图像块的重构效果对比；

图4是本发明与现有方法对于2个具有明显边界和区域信息的局部图像的重构效果对比；

图5是本发明与其它方法对于2个具有复杂纹理信息的局部图像的重构效果对比。

具体实施方式

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

步骤1，输入图像块x的观测向量y，估计该图像块的直流分量和变化分量

1.1)输入大小为像素的图像块x，得到观测向量y：

y＝Φx

其中，Φ为M×N的观测矩阵，M是观测维数，N是信号维数；

1.2)根据观测向量y和观测矩阵Φ，估计直流分量

其中，1是元素全为1的N维向量，上标表示求矩阵的Moore-Penrose伪逆；

1.3)根据观测向量y、观测矩阵Φ和直流分量估计变化分量

步骤2,对图像块x进行光滑和非光滑分类。

2.1)按如下公式计算阈值τ：

$\mathrm{\&tau;}=15+(M/N-0.05)\&times;\frac{30-15}{0.5-0.05},$

其中，M是观测维数，N是信号维数，M/N是采样率；

2.2)将变化分量的能量||x_var||₂与阈值τ进行比较：如果则判定该图像块为光滑块；否则，判定该图像块为非光滑块。

步骤3，根据图像块类别，选择重构方法：若该图像块属于光滑块，则将直流分量和变化分量相加，得到光滑块的估计：若该图像块属于非光滑快，则执行如下步骤。

步骤4，初始化PCA字典B。

4.1)构造人工黑白图：生成18幅大小为65×65像素的人工黑白图{I_i}_1≤i≤18，其中黑色像素点值为0，白色像素点值为1，每幅黑白图I_i的黑色区域与白色区域的分界线均经过图像的中心坐标(33,33)，这18幅人工黑白图分界线的角度分别从0到180度均匀采样；

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

4.3)利用矩阵对角化对图像块样本集合进行主成分(PCA)分解，得到PCA字典；

4.3.1)求每一个图像块样本集合中样本的协方差矩阵：

其中，是的转置，E[·]表示求解数学期望；

4.3.2)对协方差矩阵W_i进行对角化：

$W_i=B_i{S}_i{B}_i^T,$

其中，B_i是第i个方向的PCA子字典，其包含64个原子，S_i是对角矩阵，其对角元素为B_i中原子对应的特征值

4.3.3)将所有子字典进行级联，组合成完整的PCA字典：

B＝[B₁,B₂,...,B₁₈]。

步骤5，设置计数器t＝0，初始化父代种群

5.1)设置种群规模S＝20，即父代种群中含有20个个体：其中每个个体的长度为L＝16，即每个个体有16个元素：其中表示个体z_i的第j个元素；

5.2)将前18个个体{z_i}_1≤i≤18的每个个体按如下方法初始化：设个体的前8个元素取值为[(i-1)×64+1,…,(i-1)×64+8]，后8个元素的每个元素取值为(i-1)×64+9到i×64之间的随机整数；

5.3)将后两个个体{z₁₉,z₂₀}的所有元素分别取1到64×18之间的随机整数。

步骤6，计算父代种群中每个个体z_i的适应度值，记录适应度值最大的个体为最优个体z_best。

6.1)对于父代种群中的每个个体z_i，生成一个子矩阵其中的列向量是以个体z_i中元素为索引从PCA字典B中抽取的原子；

6.2)记录个体z_i中所有元素对应的子字典索引Q，生成一个合成协方差矩阵

$\hat{W}=\underset{q\&Element;Q}{\mathrm{\&Sigma;}}\frac{1}{\left|Q\right|}{W}_q,$

其中，W_q是生成第q个PCA子字典时产生的协方差矩阵，|Q|表示索引Q中元素的数量；

6.3)求解个体z_i对应的系数向量

求解系数向量的方法有两种：基于能量的求解方法和基于线性多高斯的求解方法方法；

6.3.1)若采用基于能量的方法求解系数向量则按如下公式计算：

$\widehat{\mathrm{\&alpha;}}={({\left(\mathrm{\&Phi;}\hat{B}\right)}^T\mathrm{\&Phi;}\hat{B}+{\mathrm{\&sigma;}}^2{\hat{S}}^{-1})}^{-1}{\left(\mathrm{\&Phi;}\hat{B}\right)}^{T}y,$

其中，Φ是观测矩阵，y是观测向量，是PCA字典B的子矩阵，σ是图像块的噪声能量，取值为3，是一对角矩阵，对角元素是中原子对应的特征值；

6.3.2)若采用基于线性多高斯的方法求解系数向量则按如下公式计算：

$\widehat{\mathrm{\&alpha;}}={({\left(\mathrm{\&Phi;}\hat{B}\right)}^T\mathrm{\&Phi;}\hat{B}+{\mathrm{\&sigma;}}^2{\hat{B}}^T{\hat{W}}^{-1}\hat{B})}^{-1}{\left(\mathrm{\&Phi;}\hat{B}\right)}^{T}y,$

其中，Φ是观测矩阵，y是观测向量，是PCA字典B的子矩阵，σ是图像块的噪声能量，取值为3，是一对角矩阵，对角元素是中原子对应的特征值，是合成协方差矩阵；

6.4)计算个体z_i的适应度函数：

$f\left(z_i\right)=\frac{1}{{\left|\right|\mathrm{\&Phi;}\hat{B}\widehat{\mathrm{\&alpha;}}-y\left|\right|}_2^2+{\mathrm{\&sigma;}}^2{\widehat{\mathrm{\&alpha;}}}^T{\hat{B}}^T{\hat{W}}^{-1}\hat{B}\widehat{\mathrm{\&alpha;}}},$

其中，Φ是观测矩阵，y是观测向量，是PCA字典B的子矩阵，是合成协方差矩阵，是系数向量；

6.5)将适应度最大的个体记录为最优个体z_best。

步骤7，对父代种群进行概率为0.8交叉操作，得到子代种群

7.1)新建一个子代种群并设为空集，置i＝1；

7.2)对于个体z_i，生成一个0到1之间的随机数p_c，如果p_c<0.8，执行步骤7.3；否则，执行步骤7.10；

7.3)将个体z_i中的元素按其原子对应的特征值由大到小排列，得到排序后的个体 $z_i=[z_1^i,z_2^i,...z_{16}^i];$

7.4)从父代种群中随机选择另外一个个体z_j，并将个体z_j中的元素按其原子对应的特征值由小到大排列，得到排序后的个体

7.5)生成一个1到16之间的随机整数c,作为个体z_i与z_j的交叉点；

7.6)将个体z_i和个体z_j在交叉点c处断开，得到4个个体片段：

$\begin{array}{cc}{z}_1^i=[z_1^i,...,z_c^i]& z_2^i=[z_{c+1}^i,...,z_{16}^i]\\ z_1^j=[z_1^j,...,z_c^j]& z_2^j=[z_{c+1}^j,...,z_{16}^j]\end{array};$

7.7)将个体片段与组合得到第一新个体z′_i，将个体片段与组合得到第二新个体z′_j：

$\begin{array}{c}{z}_i^{\&prime;}=[z_1^i,...,z_c^i,z_{c+1}^j,...,z_{16}^j]\\ z_j^{\&prime;}=[z_1^j,...,z_c^j,z_{c+1}^i,...,z_{16}^i]\end{array},$

7.8)将这两个新个体z′_i和z′_j添加到子代种群中；

7.9)置i＝i+1，判断父代种群中的所有个体是否完成交叉操作：若i>20，则所有个体已完成交叉操作，输出子代种群否则，返回步骤7.2)。

步骤8，对子代种群进行概率为0.2变异操作。

8.1)对于子代种群中的每个个体z_i，生成一个0到1之间的随机数p_c，如果p_c<0.2，对个体z_i执行下面的变异操作；否则，对个体z_i不执行操作；

8.2)随机生成1到L之间的整数c作为变异位置；

8.3)找到变异位置的元素所属的PCA子字典j：

其中，为第i个个体z_i的第c个元素，表示上取整函数；

8.4)随机生成(j-1)×64+1到j×64之间的整数d，用d替换个体z_i的第c个元素，即 $z_c^i=d.$

步骤9，按步骤6方法计算子代种群中每个个体的适应度值，并将子代个体中的最大适应度值与当前最优个体z_best的适应度值进行比较。如果子代个体中的最大适应度值大于当前最优个体z_best的适应度值，则用子代种群中适应度值最大的个体替换当前最优个体z_best；否则，不更新最优个体。

步骤10，在父代种群和子代种群的集合中，选择20个适应度值最大的个体作为新的父代种群置计数器t＝t+1。

步骤11，设最大迭代次数T_max为100，若t<T_max，则返回步骤7；否则，输出最佳个体z_best的重构结果。

本发明的效果可以通过仿真实验具体说明：

1.实验条件

实验所用计算机处理器为Intel Core(TM)2Duo CPU2.67GHz内存4GB，编程平台是Matlab2013。实验中采用的图像数据为标准图像测试库中的Lena和Barbara图像，大小为512×512。

2.实验内容

实验1：对4个具有明显方向信息的图像块进行重构

分别用本发明方法、正交匹配追踪方法OMP和统计压缩感知方法SCS在采样率0.3的条件下，对4个大小为8×8像素的图像块进行重构，并计算重构结果的峰值信噪比PSNR，重构结果如图2所示，其中：

图2(a)是4个具有明显方向信息的原始图像块，

图2(b)是OMP算法对图2(a)中图像块的重构结果，

图2(c)是SCS算法对图2(a)中图像块的重构结果，

图2(d)是本发明采用基于能量的系数求解方法对图2(a)中图像块的重构结果，

图2(e)是本发明采用线性多高斯的系数求解方法对图2(a)中图像块的重构结果。

实验2：对4个具有复杂结构信息的图像块进行重构

分别用本发明方法、正交匹配追踪方法OMP和统计压缩感知方法SCS在采样率0.3的条件下，对4个大小为8×8像素的图像块进行重构，并计算重构结果的峰值信噪比PSNR，重构结果如图3所示，其中：

图3(a)是4个具有复杂结构信息的原始图像块，

图3(b)是OMP算法对图3(a)中图像块的重构结果，

图3(c)是SCS算法对图3(a)中图像块的重构结果，

图3(d)是本发明采用基于能量的系数求解方法对图3(a)中图像块的重构结果，

图3(e)是本发明采用线性多高斯的系数求解方法对图3(a)中图像块的重构结果。

从图2和图3中可以看出，对于方向和结构比较复杂的图像块，本发明的重构结果相比于其它方法更加接近原始图像块，并且能够获得更高的峰值信噪比PSNR。

实验3：对2个具有明显边界和区域信息的局部图像进行重构

分别用本发明方法、正交匹配追踪方法OMP和统计压缩感知方法SCS在采样率0.25的条件下，重构图像Lena和Barbara中的2个具有明显边界和区域信息的局部图像。每个局部图像被分成8×8的图像块，对每个图像块分别用以上四种方法进行重构，再将图像块的重构结果组合成完整图像，并计算峰值信噪比PSNR。重构结果如图4所示，其中：

图4(a)是Lena局部原始图像，

图4(b)是OMP对图4(a)的重构结果，

图4(c)是SCS对图4(a)的重构结果，

图4(d)是本发明采用基于能量的系数求解方法对图4(a)的重构结果，

图4(e)是本发明采用线性多高斯的系数求解方法对图4(a)的重构结果，

图4(f)是Barbara局部原始图像，

图4(g)是OMP对图4(f)的重构结果，

图4(h)是SCS对图4(f)的重构结果，

图4(i)是本发明采用基于能量的系数求解方法对图4(f)的重构结果，

图4(j)是本发明采用线性多高斯的系数求解方法对B图4(f)的重构结果。

实验4：对2个具有复杂纹理信息的局部图像进行重构

分别用本发明方法、正交匹配追踪方法OMP和统计压缩感知方法SCS在采样率0.25的条件下，重构图像Lena和Barbara中的2个具有复杂纹理信息的局部图像。每个局部图像被分成8×8的图像块，对每个图像块分别用以上四种方法进行重构，再将图像块的重构结果组合成完整图像，并计算峰值信噪比PSNR。重构结果如图5所示，其中：

图5(a)是Lena局部原始图像，

图5(b)是OMP对图5(a)的重构结果，

图5(c)是SCS对图5(a)的重构结果，

图5(d)是本发明采用基于能量的系数求解方法对图5(a)的重构结果，

图5(e)是本发明采用线性多高斯的系数求解方法对图5(a)的重构结果，

图5(f)是Barbara局部原始图像，

图5(g)是OMP对图5(f)的重构结果，

图5(h)是SCS对图5(f)的重构结果，

图5(i)是本发明采用基于能量的系数求解方法对图5(f)的重构结果，

图5(j)是本发明采用线性多高斯的系数求解方法对图5(f)的重构结果。

从图4和图5中可以看出，本发明的方法比OMP和SCS有更好的区域一致性和边界连续性，对于细节的重构也优于其它两种方法，而且人工块效应明显减弱。

综上，无论从视觉效果还是PSNR值上来看，本发明相比较其它两种方法都有明显的提升。

Claims (6)

1.一种基于遗传稀疏优化和贝叶斯估计模型的图像重构方法，包括如下步骤：

(1)输入大小为像素的图像块x的观测向量y，估计该图像块的直流分量和变化分量

(2)根据变化分量的能量，将图像块x进行光滑和非光滑的分类，若该图像块属于光滑块，则执行步骤(3)操作；否则，执行步骤(4)操作；

(3)对于光滑块，直接用直流分量加上变化分量进行重构，得到重构结果为：

$\tilde{x}={\tilde{x}}_{\mathrm{dc}}+{\tilde{x}}_{\mathrm{var}};$

(4)对于非光滑块，利用遗传算法在PCA字典中选择合适的原子组合对图像块x进行重构，包括如下步骤：

(4.1)构造一个包含18个方向子字典的PCA字典B；

(4.2)初始化一个父代种群种群规模S＝20，个体长度为16，置计数器t＝0；

(4.3)计算父代种群中每个个体的适应度值，记录适应度值最大的个体为最优个体z_best；

(4.4)对父代种群进行概率为0.8的交叉操作，得到子代种群

(4.5)对子代种群进行概率为0.2的变异操作；

(4.6)计算子代种群中每个个体的适应度值，如果中个体的最大适应度值大于当前最优个体z_best的适应度值，则更新最优个体z_best；

(4.7)在父代种群和子代种群的集合中，选择20个适应度值最大的个体作为新的父代种群并且记录它们的适应度值，置计数器t＝t+1；

(4.8)设最大迭代次数T_max为100，若t<T_max，则返回步骤(4.4)；否则，输出最佳个体z_best的重构结果。

2.根据权利要求1所述的基于遗传稀疏优化和贝叶斯估计模型的图像重构方法，其中，步骤(1)所述的直流分量和变化分量按如下公式求解：

其中，1是元素全为1的N维向量，Φ是M×N的观测矩阵，y是M维的观测向量，上标表示求矩阵的Moore-Penrose伪逆。

3.根据权利要求1所述的基于遗传稀疏优化和贝叶斯估计模型的图像重构方法，其中，步骤(2)所述的将图像块x进行光滑和非光滑的分类，按如下步骤进行：

(2.1)按如下公式计算计算阈值τ：

$\mathrm{\&tau;}=15+(M/N-0.05)\&times;\frac{30-15}{0.5-0.05},$

其中，M是观测维数，N是信号维数，M/N是采样率。

(2.2)将图像块x的变化分量的能量与阈值τ进行比较：如果则判定该图像块为光滑块；否则，判断该图像块为非光滑块，其中，符号||·||₂表示求向量的2范数。

4.根据权利要求1所述的基于遗传稀疏优化和贝叶斯估计模型的图像重构方法，其中，步骤(4.3)所述的计算父代种群中每个个体的适应度值，其步骤如下：

(4.3.1)对于个体z，生成一个子矩阵其中的列向量是以个体z中元素为索引从PCA字典B中抽取的原子；

(4.3.2)记录个体z中元素对应的子字典索引Q，生成一个合成协方差矩阵

$\hat{W}=\underset{q\&Element;Q}{\mathrm{\&Sigma;}}\frac{1}{\left|Q\right|}{W}_q$

其中，W_q是生成第q个PCA子字典时产生的协方差矩阵，|Q|表示索引Q中元素的数量；

(4.3.3)根据子矩阵和合成协方差矩阵求解个体z对应的系数向量

(4.3.4)按如下公式计算个体z的适应度函数：

$f\left(z\right)=\frac{1}{{\left|\right|\mathrm{\&Phi;}\hat{B}\widehat{\mathrm{\&alpha;}}-y\left|\right|}_2^2+{\mathrm{\&sigma;}}^2{\widehat{\mathrm{\&alpha;}}}^T{\hat{B}}^T{\hat{W}}^{-1}\hat{B}\widehat{\mathrm{\&alpha;}}}$

其中，Φ是观测矩阵，y是观测向量，是PCA字典B的子矩阵，是合成协方差矩阵，是系数向量。

5.根据权利要求4所述的计算个体的适应度值，其步骤(4.3.3)所述的根据子矩阵和合成协方差矩阵求解个体z对应的系数向量采用基于能量的求解方法，其计算公式如下：

$\widehat{\mathrm{\&alpha;}}={({\left(\mathrm{\&Phi;}\hat{B}\right)}^T\mathrm{\&Phi;}\hat{B}+{\mathrm{\&sigma;}}^2{\hat{S}}^{-1})}^{-1}{\left(\mathrm{\&Phi;}\hat{B}\right)}^{T}y,$

其中，Φ是观测矩阵，y是观测向量，是PCA字典B的子矩阵，σ是图像块的噪声能量，取值为3，是一对角矩阵，对角元素是中原子对应的特征值。

6.根据权利要求4所述的计算个体的适应度值，其步骤(4.3.3)所述的根据子矩阵和合成协方差矩阵求解个体z对应的系数向量采用基于线性多高斯的求解方法，其计算公式如下：

$\widehat{\mathrm{\&alpha;}}={({\left(\mathrm{\&Phi;}\hat{B}\right)}^T\mathrm{\&Phi;}\hat{B}+{\mathrm{\&sigma;}}^2{\hat{B}}^T{\hat{W}}^{-1}\hat{B})}^{-1}{\left(\mathrm{\&Phi;}\hat{B}\right)}^{T}y,$

其中，Φ是观测矩阵，y是观测向量，是PCA字典B的子矩阵，σ是图像块的噪声能量，取值为3，是合成协方差矩阵。