在多变量观测下基于边缘约束的图像重构方法

摘要

本发明公开了一种在多变量观测下基于边缘约束的图像重构方法，主要解决现有技术对压缩感知图像重构不准确和低鲁棒的问题。其实现过程为：1)接收观测矩阵、低频小波分解系数和多变量测量矩阵；2)通过边缘检测和相关性指导得到非零系数组支撑集合；3)在多变量高斯模型的基础上，根据观测矩阵、多变量测量矩阵、基础协方差和残差协方差矩阵，运用吉布斯采样的方法重构非零系数组支撑集合中的高频小波系数；4)将低频小波分解系数和重构的高频小波系数进行小波逆变换得到重构图像。本发明与OMP和BEPA方法相比，具有重构图像质量高、鲁棒性好的优点，可用于自然图像和医学图像的重构。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及统计压缩感知图像重构方法，可用于对自然图像进行重构。

背景技术

近几年，在信号处理领域出现了一种新的数据理论压缩感知CS，该理论在数据采集的同时实现压缩，突破了传统奈奎采集斯特采样定理的限制，为数据采集技术带来了革命性的变化，使得该理论在压缩成像系统、军事密码学、无线传感等领域有着广阔的应用前景。压缩感知理论主要包括信号的稀疏表示、信号的观测和信号的重构等三个方面。其中设计快速有效的重构算法是将CS理论成功推广并应用于实际数据模型和采集系统的重要环节。

在压缩采样领域，小波基是一组很好的稀疏基。图像经过小波分解后得到的分解系数，分为低频部分和高频部分，低频部分包含原始图像的低频稀疏，通常认为是非稀疏的，而高频部分包含图像的水平、垂直、对角信息，具有良好的稀疏性。目前，经常采用在小波域下对低频全部保留，对高频进行压缩观测的采样方法。该采样方法的优点是可有效提高重构图像质量。

Lihan He等人在文献“Exploiting Structure in Wavelet-Based Bayesian Compressive Sensing”中提出基于小波树结构的贝叶斯压缩感知图像重构方法。该方法对多尺度小波系数分层建立单高斯模型，并通过吉布斯采样重构图像。但该方法将图像展开成列向量进行观测重构，不仅没有结合原始图像数据的先验，并且对计算机内存要求很高，限制了处理图像的大小。

Jiao Wu等人在文献“Multivariate Compressive Sensing for Image Reconstruction in the Wavelet Domain:Using Scale Mixture Models”中提出基于混合尺度模型的多变量压缩感知图像重构MPA。该方法对小波系数构造多变量分布模型，抓住小波系数具有聚集性这一特点，对其统计相关性进行建模，但该方法忽略了保留下来的小波低频系数对图像重构的指导作用，从而导致其不具有鲁棒性，且重构出的图像不够准确。

发明内容

本发明的目的在于针上述已有技术的不足，提出一种在多变量观测下基于边缘约束的图像重构方法，以充分利用保留的低频小波系数对图像重构的指导作用，提高重构图像的准确 性。

现本发明目的技术思路是：通过对多变量测量矩阵建立多变量高斯模型，抓住小波的聚集性；通过边缘检测和相关性的联合，指导确定小波系数的非零支撑；通过对非零支撑系数利用吉布斯采样方法依次迭代更新，实现高质量的压缩感知图像重构。

根据上述思路，本发明的技术方案包括如下步骤：

1.一种在多变量观测下基于边缘约束的图像重构方法，包括如下步骤：

(1)接收方接收图像发送方发送的正交随机高斯观测矩阵Φ、低频小波分解系数L、水平高频子带多变量测量矩阵Y₁、垂直高频子带多变量测量矩阵Y₂和对角高频子带多变量测量矩阵Y₃，将三个高频子带多变量测量矩阵统一用Y表示；

(2)根据接收的观测矩阵Φ、低频小波分解系数L和高频子带多变量测量矩阵Y，通过边缘检测和相关性的指导得到非零系数组索引集合：u＝{s₁,s₂,...,s_i,...,s_c}，其中s_i代表第i个非零系数组的索引，i＝1,2,...,c，c为小于Φ的列数：

(2.1)将接收的低频小波分解系数L和三个全部为零的高频子带进行小波逆变换，得到边缘模糊图像；

(2.2)对边缘模糊图像进行边缘检测，得到边缘位置；

(2.3)提取边缘模糊图像中对应边缘位置的像素得到模糊边缘；

(2.4)对模糊边缘进行一层小波变换，得到模糊边缘小波高频系数和模糊边缘小波低频系数；将模糊边缘小波高频系数绝对值大于阈值h的位置设为1，小于阈值h的位置设为0得到初始模糊位置矩阵，将初始模糊位置矩阵按照多变量矩阵的形式排列成M×Q维的多变量模糊位置矩阵E，其中M为Φ的列数，Q为Y的列数，阈值h＝0.2；

(2.5)根据观测矩阵Φ的转置和高频子带多变量测量矩阵Y相乘得到的相关性矩阵Φ^T*Y,将该相关性矩阵的绝对值|Φ^T*Y|和多变量模糊位置矩阵E进行加权求和，得到系数重要性矩阵V＝|Φ^T*Y|+w*E，其中w为加权系数；

(2.6)将系数重要性矩阵V的每一行相加，得到系数重要性向量，将系数重要性向量中最大的c个元素的索引取出构成非零系数组索引集合u＝{s₁,s₂,...,s_i,...,s_c}，其中 i＝1,2,...,c，s_i代表第i个非零系数组的索引，c为小于M的整数；

(3)初始化外部迭代次数t＝1，模型求解迭代总次数为N₁，系数求解迭代总次数为N₂；初始化待重构小波高频系数矩阵和总叠加结果O均为M×Q维的零矩阵，其中M为Φ的列数，Q为Y的列数；初始化待重构小波高频系数矩阵的每一行均服从基础协方差矩阵Ω的多变量高斯分布；初始化残差的每一行均服从残差协方差矩阵Π的多变量高斯分布；

(4)根据观测矩阵Φ、高频子带多变量测量矩阵Y、基础协方差矩阵Ω、残差协方差矩阵Π和非零系数组索引集合u，通过吉布斯采样方法计算待重构小波高频系数矩阵

(4.1)初始化迭代次数i＝1；

(4.2)根据基础协方差矩阵Ω、残差协方差矩阵Π、观测矩阵Φ得到待重构小波高频系数矩阵的第s_i行的协方差矩阵：

其中为观测矩阵Φ的第s_i列，为的转置；

(4.3)根据残差协方差矩阵Π、观测矩阵Φ和待重构小波高频系数矩阵第s_i行的协方差矩阵得到待重构小波高频系数矩阵第s_i行的均值向量：

其中为观测矩阵Φ的第k列，x_k为待重构小波高频系数矩阵的第k行，其中k为大于等于1小于等于M且不等于s_i的整数；

(4.4)根据待重构小波高频系数矩阵第s_i行的均值向量和协方差矩阵建立对应的多变量高斯模型，计算待重构小波高频系数矩阵的第s_i行系数 $x_{s_i}=\mathrm{Gaussian}({\mu }_{s_i},{\Lambda }_{s_i}),$其中，表示随机生成一个服从均值向量为协方差矩阵为的多变量高斯分布的向量；

(4.5)将迭代次数i与非零系数组索引集合u元素个数c进行比较：如果迭代次数i<c， 则迭代次数i自加1，返回步骤(4.2)，否则，输出待重构小波高频系数矩阵

(5)根据待重构小波高频系数矩阵计算基础协方差矩阵Ω和残差协方差矩阵Π：

$\Omega ={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+c,a_0+\frac{1}{2}\mathrm{diag}({\tilde{X}}^T*\tilde{X}))\right)\right)}^{-1},$

$\Pi ={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+\frac{M}{2},a_0+\frac{1}{2}\mathrm{diag}\left({(Y-\Phi *\tilde{X})}^T(Y-\Phi *\tilde{X})\right))\right)\right)}^{-1},$

其中参数a₀为Q维常数向量，向量每个元素的值均为0.000001，*为相乘操作，·^T为转置操作，·^-1为求逆操作，diag(·)操作表示将括号内矩阵的对角线元素取出组成向量，Gamma(·,·)操作表示生成一个以括号内第一项为形状参数，括号内第二项为尺度参数的伽马分布的向量，其中括号内的两项与产生的向量均为Q维，diag0(·)操作表示产生一个方阵，方阵的对角线元素为括号内的向量元素，非对角线元素为0；

(6)将外部迭代次数t与模型求解迭代总次数N₁进行比较：如果外部迭代次数t≤N₁，则t＝t+1，返回步骤(4)，否则执行步骤(7)；

(7)设上代叠加结果等于总叠加结果O，将外部迭代次数t分别与模型求解迭代总次数N₁和系数求解迭代总次数N₂进行比较：如果外部迭代次数N₁<t<N₁+N₂，则计算总叠加结果外部迭代次数t＝t+1，返回步骤(4)，否则计算总叠加结果计算最终重构小波高频系数矩阵执行步骤(8)；

(8)根据保留的低频小波分解系数L和最终重构小波高频系数矩阵进行小波逆变换，得到原图的重构图。

本发明通过建立多变量模型高斯抓住小波的聚集性，在求解时使用基于吉布斯采样方法对多变量高斯模型的参数自适应修正，并且用保留的小波低频系数通过边缘检测得到的模糊位置对小波系数的非零支撑进行指导，使图像的重构质量和鲁棒性得到了显著提高。

附图说明

图1是本发明的总流程图；

图2是本发明中将小波高频子带系数转化为多变量系数矩阵的示意图；

图3是本发明中确定模糊边缘位置的示意图；

图4是本发明与现有技术在采样率为40％时对Boat图像的重构结果图；

图5是用本发明与现有技术重构的Lena图像的峰值信噪比随采样率变化曲线图。

具体实施方式

以下结合附图和具体实施方式对本发明进行详细说明。

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

步骤一，发送方发送观测矩阵、低频小波分解系数和高频多变量测量矩阵。

(1a)图像发送方在小波域观测图像，对低频小波分解系数全部保留作为低频小波分解系数的观测，用正交随机高斯观测矩阵Φ对水平高频子带系数A₁、垂直高频子带系数A₂和对角高频子带系数A₃分别进行多变量观测，得到水平高频子带多变量测量矩阵Y₁、垂直高频子带多变量测量矩阵Y₂和对角高频子带多变量测量矩阵Y₃，其中多变量观测的过程为：

(1a1)，将水平高频子带系数A₁、垂直高频子带系数A₂、对角水平高频子带系数A₃都分为M个大小为个像素块，每个像素块拉成行向量，转化成的M×Q维水平多变量系数矩阵X₁、垂直多变量系数矩阵X₂、对角多变量系数矩阵X₃，如图2所示，将256×256个像素的高频子带系数A按照3×3个像素分块，由于不能整分高频子带系数A的最后两行和最后两列用零添补，多变量系数矩阵X的每一行对应每个3×3的像素块拉成行；

(1a2)，对水平多变量系数矩阵X₁、垂直多变量系数矩阵X₂、对角多变量系数矩阵X₃进行观测，得到水平高频子带多变量测量矩阵Y₁＝Φ*X₁、垂直高频子带多变量测量矩阵Y₂＝Φ*X₂、对角高频子带多变量测量矩阵Y₃＝Φ*X₃；

(1b)图像发送方发送正交随机高斯观测矩阵Φ、低频小波分解系数L、水平高频子带多变量测量矩阵Y₁、垂直高频子带多变量测量矩阵Y₂和对角高频子带多变量测量矩阵Y₃给接收方。

步骤二，根据接收的观测矩阵Φ、低频小波分解系数L和高频子带多变量测量矩阵Y，通过边缘检测和相关性的指导得到非零系数组索引集合u。

(2a)接收方接收图像发送方发送的观测矩阵Φ、低频小波分解系数L、水平高频子带多变量测量矩阵Y₁、垂直高频子带多变量测量矩阵Y₂和对角高频子带多变量测量矩阵Y₃，由于三个高频子带系数的重构方式一样，因此三个高频子带多变量测量矩阵统一用Y表示；

(2b)将接收的低频小波分解系数L和三个全部为零的高频子带进行小波逆变换，得到边缘模糊图像，如图3中的小波逆变换过程所示；

(2c)对边缘模糊图像用Canny边缘检测方法进行边缘检测，得到边缘位置，提取边缘模糊图像中对应边缘位置的像素得到模糊边缘，如图3中的Canny边缘检测过程所示；

(2d)对模糊边缘进行一层小波变换，得到模糊边缘小波高频系数和模糊边缘小波低频系数，如图3中小波分解所示；将模糊边缘小波高频系数绝对值大于阈值h的位置设为1，小于阈值h的位置设为0得到初始模糊位置矩阵，如图3中的得到模糊位置矩阵所示,将初始模糊位置矩阵按照多变量矩阵的形式排列成M×Q维的多变量模糊位置矩阵E，其中M为Φ的列数，Q为Y的列数，阈值h＝0.2；

(2e)根据观测矩阵Φ的转置和高频子带多变量测量矩阵Y相乘得到的相关性矩阵Φ^T*Y,将该相关性矩阵的绝对值|Φ^T*Y|和多变量模糊位置矩阵E进行加权求和，得到系数重要性矩阵V＝|Φ^T*Y|+w*E，其中w为加权系数；

(2f)将系数重要性矩阵V的每一行相加，得到系数重要性向量，将系数重要性向量中最大的c个元素的索引取出构成非零系数组索引集合u＝{s₁,s₂,...,s_i,...,s_c}，其中i＝1,2,...,c，s_i代表第i个非零系数组的索引，c为小于M的整数。

步骤三，初始化迭代次数和多变量高斯模型。

(3a)初始化外部迭代次数t＝1，模型求解迭代总次数为N₁，系数求解迭代总次数为N₂；

(3b)初始化待重构小波高频系数矩阵和总叠加结果O均为M×Q维的零矩阵，其中M为Φ的列数，Q为Y的列数；

(3c)初始化待重构小波高频系数矩阵的每一行均服从基础协方差矩阵Ω的多变量高斯分布，其中Ω为Q×Q维矩阵，其对角线元素为1，非对线元素为0；

(3d)初始化残差的每一行均服从残差协方差矩阵Π的多变量高斯分布，其中Π为Q×Q维矩阵，其对角线元素为0.01，非对角线元素为0。

步骤四，根据观测矩阵Φ、高频子带多变量测量矩阵Y、基础协方差矩阵Ω、残差协方差矩阵Π和非零系数组索引集合u，通过吉布斯采样方法计算待重构小波高频系数矩阵 

(4a)初始化迭代次数i＝1；

(4b)根据基础协方差矩阵Ω、残差协方差矩阵Π、观测矩阵Φ得到待重构小波高频系数矩阵的第s_i行的协方差矩阵：

其中为观测矩阵Φ的第s_i列，为的转置，·^-1为求逆操作；

(4c)根据残差协方差矩阵Π、观测矩阵Φ和待重构小波高频系数矩阵第s_i行的协方差矩阵得到待重构小波高频系数矩阵第s_i行的均值向量：

其中为观测矩阵Φ的第k列，x_k为待重构小波高频系数矩阵的第k行，其中k为大于等于1小于等于M且不等于s_i的整数；

(4d)根据待重构小波高频系数矩阵第s_i行的均值向量和协方差矩阵建立对应的多变量高斯模型，计算待重构小波高频系数矩阵的第s_i行系数 $x_{s_i}=\mathrm{Gaussian}({\mu }_{s_i},{\Lambda }_{s_i}),$其中，表示随机生成一个服从均值向量为协方差矩阵为的多变量高斯分布的向量；

(4e)将迭代次数i与非零系数组索引集合u元素个数c进行比较：如果迭代次数i<c，则迭代次数i自加1，返回步骤(4.2)，否则，输出待重构小波高频系数矩阵

步骤5，根据待重构小波高频系数矩阵计算基础协方差矩阵Ω和残差协方差矩阵Π：

$\Omega ={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+c,a_0+\frac{1}{2}\mathrm{diag}({\tilde{X}}^T*\tilde{X}))\right)\right)}^{-1},$

$\Pi ={\left(\mathrm{diag}0\left(\mathrm{Gamma}(a_0+\frac{M}{2},a_0+\frac{1}{2}\mathrm{diag}\left({(Y-\Phi *\tilde{X})}^T(Y-\Phi *\tilde{X})\right))\right)\right)}^{-1},$

其中，参数a₀为Q维常数向量，向量每个元素的值均为0.000001，*为相乘操作，·^T为转置操作，·^-1为求逆操作，diag(·)操作表示将括号内矩阵的对角线元素取出组成向量， Gamma(·,·)操作表示生成一个以括号内第一项为形状参数，括号内第二项为尺度参数的伽马分布的向量，其中括号内的两项与产生的向量均为Q维，diag0(·)操作表示产生一个方阵，方阵的对角线元素为括号内的向量元素，非对角线元素为0。

步骤6，将外部迭代次数t与模型求解迭代总次数N₁进行比较：如果外部迭代次数t≤N₁，则t＝t+1，返回步骤4，否则，执行步骤7。

步骤7，设上代叠加结果等于总叠加结果O，将外部迭代次数t分别与模型求解迭代总次数N₁和系数求解迭代总次数N₂进行比较：如果外部迭代次数N₁<t<N₁+N₂，则计算总叠加结果并设外部迭代次数t＝t+1，返回步骤4；否则，计算总叠加结果 计算最终重构小波高频系数矩阵执行步骤8。

步骤8，根据保留的低频小波分解系数L和最终重构小波高频系数矩阵经过小波逆变换，得到重构图像。

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

1、仿真条件：本发明的仿真在windows 7，SPI，CPU Intel(R)Core(TM)i5-3470，基本频率3.20GHz，软件平台为Matlab R2011b上运行，仿真选用的是512×512的四幅标准测试自然图像Lena、Peppers、Boat、Barbara，对于OMP算法分块大小32×32,对于MPA方法和本发明方法使用阈值为0.15的Canny边缘检测，c＝1800，相关性矩阵的绝对值和模糊位置矩阵的加权值ω＝1.5，迭代次数N₁＝100，N₂＝200。

2、仿真内容与结果：

仿真1：在固定采样率40％下，用本发明与现有的OMP、MPA方法对标准测试自然图像在小波域下进行重构，图像Boat的重构视觉效果如图4所示，其中图4(a)为Boat原图，图4(b)是图4(a)的局部放大图，图4(c)、图4(e)和图4(g)分别是OMP、MPA和本发明方法的重构图，图4(d)、图4(f)和图4(h)分别是图4(c)、图4(e)和图4(g)的局部放大图。

从重构图和局部放大图可以看出，本发明的重构图像的边缘部分保持的较好，平滑部分的噪声也比OMP、MPA的重构图像的少得多。

在固定采样率40％下，用本发明与现有的OMP、MPA方法对四幅大小为512*512的图像Lena、Peppers、Barbara、Boat分别进行五次重构，五次重构重构结果的峰值信噪比PSNR 的平均值如表1所示。

表1 自然图像大小为512*512，采样率40％的重构结果

从表1可以看出，本发明重构图像的平均PSNR值比OMP、MPA方法均高，表明重构图像的质量好。

仿真2：在采样率分别为30％、35％、40％、45％、50％的情况下，用本发明与现有的OMP、BEPA方法在小波域下对大小为512×512的Lena图像的进行重构，5次重构结果的峰值信噪比PSNR的平均值如表2所示。

表2 Lena图像在不同采样率下用OMP、MPA和本发明方法的重构结果

从表2数据可以看出，本发明方法在采样率为30％、35％、40％、45％、50％下得到的结果图的峰值信噪比PSNR都要高于OMP和MPA方法得到的PSNR，即本发明的方法的重构图像质量要比OMP和MPA方法高。

根据表2数据得到OMP，MPA和本发明方法重构出的Lena图像的PSNR随采样率变化的趋势图如图5所示，图5中的横坐标表示采样率，纵坐标表示峰值信噪比PSNR(dB)值。

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

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