CN103037212A - 基于视觉感知的自适应块压缩传感图像编码方法 - Google Patents

Abstract

一种低复杂的基于视觉感知的自适应块压缩传感图像编码方法，其特征是：I.块压缩传感：1.编码：读入图像并分块；块观测；将观测值和观测率发至解码端；2.解码：由观测率得观测矩阵；求初始解；维纳滤波；用PL更新；DCT变换；双变量收缩阈值处理；IDCT；用PL更新；反复迭代到初始阶段结束；II.视觉分析：分析每个重构的DCT系数块；块分类；送入反馈信道；III.视觉感知指导的自适应块压缩传感：1.编码：根据反馈来的结果自适应观测；将观测值和观测率送至解码端；2.解码：组合观测值；上一阶段恢复图像做初值；维纳滤波；用PL更新；DCT变换；双变量收缩阈值处理；IDCT；用PL更新；迭代到解码结束；进行下一阶段分析、自适应观测和重构，到恢复图像满足要求。

Description

基于视觉感知的自适应块压缩传感图像编码方法

技术领域

本发明属于图像编码方法技术领域，具体涉及一种低复杂度的基于视觉感知的自适应块压缩传感图像编码方法。

背景技术

压缩传感(Compressive Sensing，简写为CS)理论是近几年来刚刚兴起的一种全新的信号处理方法。其优点在于CS对信号的采样数据量远远小于传统采样方法所获取的数据量，突破了香农采样定理的瓶颈，使得高分辨率信号的采集成为可能。CS理论要求信号本身是可压缩的或在某种变换域下可以稀疏表示，采用非自适应线性投影来保持信号的原始结构，然后通过数值最优化问题准确重构原始信号。

由于CS中的采样一般通过一个随机矩阵来实现，一旦应用于二维图像，将面临如下问题：一是观测矩阵比较大，将对内存资源要求高，并使计算量增大；二是整幅图像的变换忽略了图像在空间位置和时间域上的局部稀疏性，造成系统性能下降。为此，Lu Gan提出了块压缩传感(Block Compressive Sensing，简写为BCS)，对图像进行分块观测，从而降低了观测矩阵的大小；同时，一些快速有效的恢复算法也接连被提出，如M.A.T.Figueiredo等人提出的GPSR算法、T.T.Do等人提出的稀疏自适应匹配追踪算法等。但是这种块压缩传感算法并没有考虑到图像中不同块的特性及其对恢复质量的影响。专利(申请号：201110074327.X)提出了一种基于自适应块压缩传感的渐进图像编码方法，根据块分类使用不同的观测率对图像块进行采样。但是这种方法中，一方面由于块分类在编码端进行，其复杂性影响了CS固有的编码简单特性；另一方面这种分类方法仅针对图像像素适用，忽略了人眼的感知特性。

近年来，人类视觉特性被大量应用于图像处理中以进一步提高系统的主观、客观恢复性能。如：国内学者路文提出了基于HDW的图像质量评价方法，梁慧提出将视觉分析用到视频隐写中。由于视觉感知信息可以通过分析某种变换域下(如：离散余弦变换DCT，小波变换以及contourlet变换)的系数来获得，以更能够捕捉到图像中的细节和变换信息。而CS重构也都在变换域中进行，这就很方便进行CS域的视觉分析，从而提高图像的恢复质量。然而，目前尚未见到如何将视觉特性应用在CS中。

发明内容

本发明的目的是提供一种低复杂度的、利用人类视觉感知特性来指导自适应块压缩传感(Adaptive Block compressive sensing，简写为ABCS)的图像编码方法，可以克服现有技术中存在的诸多缺点。

本发明的技术方案是：基于视觉感知的自适应块压缩传感图像编码方法，其特征是包括下列步骤：

I.块压缩传感(BCS)，具体包括下列步骤：

1.BCS编码：

(1).读入一副N×N大小图像x，将其分成多个W×W大小不重叠的图像块；

(2).生成一个W²×W²大小的正交高斯随机矩阵作为种子观测阵Φ；

(3).对每一个图像块x_i利用公式(1)进行CS观测：

y_i＝Φ_W.x_i，          (1)

其中，Φ_W是第i个图像块的观测矩阵，其大小为m×W²，由种子观测阵Φ的前m行组成，

MR为预设的观测率；y_i是对应x_i的观测值，所有y_i组成整幅图像x的观测值集合y；

(4).将观测值集合y以及观测率MR发送到解码端；

2.BCS解码，具体步骤如下：

(1).在解码端，预先存储一个与编码端相同的种子观测阵Φ，通过编码端发送过来的观测率MR抽取Φ的前m行，得到观测矩阵Φ_W；

(2).由公式(2)得到整幅图像的初始解：

$x^{\left(0\right)}={\mathrm{\Φ}}_0^{T}y,---\left(2\right)$

其中，Φ₀是由Φ_W组成的块对角矩阵，

是Φ₀的转置矩阵，y是由编码端发送过来的观测值集合；

(3).对x^(j)(j代表迭代次数，第一次迭代时为x⁽⁰⁾)进行维纳滤波，以消除图像的块效应；

(4).对于滤波后的图像中每一个块通过PL算法进行更新，如公式(3)所示：

${\hat{x}}_i^{\left(j\right)}=x_i^{\left(j\right)}+{\mathrm{\Φ}}_W^T(y_i-{\mathrm{\Φ}}_W{x}_i^{\left(j\right)});---\left(3\right)$

(5).利用公式(4)对每个块

进行DCT变换，得到

的稀疏表示，

${{\mathrm{\θ}}_i}^{\left(j\right)}=\mathrm{\ψ}{\hat{x}}_i^{\left(j\right)},---\left(4\right)$

其中，θ_i ^(j)为

在DCT基下的系数，Ψ为DCT变换基；

(6).对每个θ_i ^(j)利用公式(5)进行双变量收缩阈值处理，得到更加稀疏的系数：

${{\widehat{\mathrm{\θ}}}_i}^{\left(j\right)}=\frac{{(\sqrt{{\left({{\mathrm{\θ}}_i}^{\left(j\right)}\right)}^2+{\left({{\mathrm{\θ}}_i}^{(j-1)}\right)}^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(j\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({{\mathrm{\θ}}_i}^{\left(j\right)}\right)}^2+{\left({{\mathrm{\θ}}_i}^{(j-1)}\right)}^2}}\·{{\mathrm{\θ}}_i}^{\left(j\right)},---\left(5\right)$

其中，若g≤0，则(g)+＝0，若g＞0，则(g)+＝g；λ是一个收敛控制常数；

是θ_i ^(j)的中位数估计值；

是θ_i ^(j)的边缘方差；

(7).利用公式(6)对

进行DCT变换，得到

的稀疏表示，

${\mathrm{\θ}}^{\left(j\right)}=\mathrm{\ψ}{\hat{x}}^{\left(j\right)},---\left(6\right)$

其中，θ^(j)为

在DCT基下的系数，Ψ为DCT变换基；

(8).对θ^(j)利用公式(7)进行双变量收缩阈值处理，得到更加稀疏的系数：

${\widehat{\mathrm{\θ}}}^{\left(j\right)}=\frac{{(\sqrt{{\left({\mathrm{\θ}}^{\left(j\right)}\right)}^2+{\left({\mathrm{\θ}}^{(j-1)}\right)}^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(j\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({\mathrm{\θ}}^{\left(j\right)}\right)}^2+{\left({\mathrm{\θ}}^{(j-1)}\right)}^2}}\·{\mathrm{\θ}}^{\left(j\right)},---\left(7\right)$

其中，若g≤0，则(g)+＝0，若g＞0，则(g)+＝g；λ是一个收敛控制常数；

是θ^(j)的中位数估计值；是θ^(j)的边缘方差；

(9).利用公式(8)对

进行反变换(IDCT)，得到本次迭代的近似解：

${\stackrel{\‾}{x}}^{\left(j\right)}={\mathrm{\ψ}}^{-1}{\widehat{\mathrm{\θ}}}^{\left(j\right)};---\left(8\right)$

(10).对于

中的每一个块

利用公式(9)进行更新：

$x_i^{(j+1)}={\stackrel{\‾}{x}}_i^{\left(j\right)}+{\mathrm{\Φ}}_W^T(y_i-{\mathrm{\Φ}}_W{\stackrel{\‾}{x}}_i^{\left(j\right)});---\left(9\right)$

(11).重复进行步骤(3)-(10)，直到得到的解满足终止条件(如公式(10)所示)，停止迭代，得到初始阶段的所有块的最优DCT系数

并恢复图像R₀：

|D^(j+1)-D^(j)|＜10^-4，         (10)

其中， $D^{\left(j\right)}=\frac{1}{N}{\left|\right|x^{\left(j\right)}-{\hat{x}}^{(j-1)}\left|\right|}_2.$

II.视觉特性分析，具体包括下列步骤：

(1).将每个重构的W×W大小的DCT系数块

分成四个部分：直流(DC)、低频(LF)、中频(MF)和高频(HF)，并分别标记四个部分的系数绝对值之和为DC、L、M和H；

(2)根据下列五种情况将DCT系数块分成平滑块、边沿块和纹理块。

情况A：若M+H≤μ₁，那么这个块是平滑的，被划分为平滑块；否则进入情况B的判断，其中μ₁＝125；

情况B：若M+H＞μ₂不成立，那么进入情况C1的判断；否则进入情况C2的判断。其中μ₂＝900；

情况C1：若满足下列条件之一：①L/M≥α₁且(L+M)/H≥β₁；②L/M≥β₁且(L+M)/H≥α₁；③(L+M)/H≥γ.那么被划分为边沿块；否则进入情况D的判断。其中α₁＝2.3，β₁＝1.6，γ＝4；

情况C2：若满足下列条件之一：①L/M≥α₂且(L+M)/H≥β₂；②L/M≥β₂且(L+M)/H≥α₂；③(L+M)/H≥γ.那么被划分为边沿块；否则被划分为纹理块。其中α₂＝1.4，β₂＝1.1，γ＝4；

情况D：若M+H＞κ，那么被划分为纹理块，否则这个块是平滑的，被划分为平滑块。其中κ＝290。

此外，进一步划分边沿块：若L+M≤400，则该块判为边沿1，否则被判为边沿2；

(3).将块分类结果送入反馈信道。

III.视觉感知指导的自适应块压缩传感(ABCS)，具体包括下列步骤：

1.ABCS编码：

(1).根据每个阶段反馈信道发送来的块分类结果，对不同类型的块使用不同的观测率进行自适应观测，如公式(11)所示：

y_ni＝Φ_Wnix_i，              (11)

其中，Φ_Wni是第n个阶段中第i个块x_i的观测矩阵，由抽取种子观测阵Φ中的

行所得，相应的观测率为MR_ni；对于平滑块，MR_ni取较大的值；对纹理块，MR_ni取较小的值；而对边沿块，MR_ni则取中间值，其中边沿1的MR_ni值略大于边沿2的；

(2).将自适应观测值y_ni与相应的观测率MR_ni发送到解码端；

2.ABCS解码，具体步骤如下：

(1).在解码端，将当前接收到的每一个块的所有观测值(即：包含存储器中之前阶段发送过来的块观测值)组合起来，如公式(12)：

$Y_{\mathrm{ni}}=\left(\begin{array}{c}{y}_{\mathrm{ni}}\\ \·\\ \·\\ \·\\ y_{1i}\\ y_i\end{array}\right),---\left(12\right)$

其中，Y_ni表示第i个图像块在n个阶段中的所有观测值，将所有Y_ni组成整幅图像的观测值集合Y_n。同时，利用编码端发送来的观测率MR_ni得到相应的观测矩阵Φ_Wni(Φ_Wni的抽取方法跟公式(9)中相同)，并将其组合起来，得到对应于Y_ni的观测阵Φ_ni，如公式(13)所示：

${\mathrm{\Φ}}_{\mathrm{ni}}=\left(\begin{array}{c}{\mathrm{\Φ}}_{\mathrm{Wni}}\\ \·\\ \·\\ \·\\ {\mathrm{\Φ}}_{W1i}\\ {\mathrm{\Φ}}_W\end{array}\right);---\left(13\right)$

(2).将上一阶段恢复的图像R_n-1作为当前阶段迭代的初始解，以加快重构的速度；

(3).对x^(j)(j代表迭代次数，第一次迭代时x^(j)＝R₀)进行维纳滤波，以消除图像的块效应；

(4).对于滤波后的x^(j)中的每一个块

通过PL算法进行更新精确，如公式(14)所示：

${\hat{x}}_i^{\left(j\right)}=x_i^{\left(j\right)}+{\mathrm{\Φ}}_{\mathrm{ni}}^T(Y_{\mathrm{ni}}-{\mathrm{\Φ}}_{\mathrm{ni}}{x}_i^{\left(j\right)});---\left(14\right)$

(5).利用公式(4)对每个块进行DCT变换，得到的稀疏表示θ_i ^(j)；

(6).利用公式(5)对每个θ_i ^(j)进行双变量收缩阈值处理，得到更加稀疏的系数

(7).利用公式(6)对

进行DCT变换，得到

的稀疏表示θ^(j)；

(8).利用公式(7)对θ^(j)进行双变量收缩阈值处理，得到更加稀疏的系数

(9).通过公式(8)对

进行反变换(IDCT)，得到本次迭代的近似解

(10).对于中的每一个块

通过公式(15)进行更新：

$x_i^{(j+1)}={\stackrel{\‾}{x}}_i^{\left(j\right)}+{\mathrm{\Φ}}_{\mathrm{ni}}^T(Y_{\mathrm{ni}}-{\mathrm{\Φ}}_{\mathrm{ni}}{\stackrel{\‾}{x}}_i^{\left(j\right)});---\left(15\right)$

(11).重复进行步骤(3)-(10)，直到得到的解满足终止条件(如公式(10)所示)，停止迭代，得到本阶段的所有块的最优DCT系数

并恢复图像R_n；

(12).根据解码端对恢复图像质量的要求，决定是否重复进行视觉特性分析和视觉感知指导的自适应块压缩传感(ABCS)，得到质量渐进提高的恢复图像。

本发明的优点及积极效果：

1、本发明的渐进图像编码系统，同时采用了块压缩传感技术、人类视觉特性分析技术、视觉分析指导的自适应块压缩传感技术、以及分阶段编码技术，有效的提高了图像的视觉质量。

2、与现有CS图像编码方法相比，本发明的创新点及其优势体现在以下几点：(1)将人类视觉特性与自适应块压缩传感技术结合，既降低了观测矩阵对存储空间的要求，又考虑到了人眼对不同类图像块的不同敏感度，从而有效地提高了图像的恢复质量；(2)将视觉分析转移到解码端，利用CS重构在某种变换域下的特性，直接对重构的变换域系数进行视觉分析，通过反馈信道将分析结果传送到编码端以指导观测率的分配，保持了CS“低复杂度编码”的特性；(3)采用分阶段的编码技术，考虑到了信道的带宽变化和对图像质量的不同要求，实现了图像的渐进编码，最终得到满足人眼视觉要求的恢复图像。

3、本发明很适合于目前新出现的一些要求“低复杂度编码”的通信设备以及信道带宽较窄或容易发生拥塞和对图像质量有可伸缩要求的应用场合，如无线传感器网络、手机照相、无线视频监控等。

附图说明

图1是本发明的系统实施框图；

图2是本发明中DCT系数块的频率区域划分；

图3是本发明中块分类的流程图；

图4是本发明中观测矩阵的构造示意图；

图5是本发明的编解码算法程序流程图；

图6是本发明的不同等效观测率下恢复图像放大部分的对比图；

图7是本发明的不同等效观测率下测试图像Lena的块分类图。

具体实施方式

本发明提出的基于视觉感知的自适应块压缩传感图像编码系统由软件仿真实现。

在图1中我们给出了本发明所提出系统的一个实施框图。输入图像后，对图像进行分块，对分好的图像块x_i首先采用传统的块压缩传感(BCS)方法进行观测(所有块使用相同的观测率即Φ_Wni＝Φ_W)，观测值y_i传输到解码端进行整体重构和分块重构，其中，整体重构后的系数通过后处理得到初始阶段的恢复图像，而分块重构后得到的DCT系数块进行第一阶段的视觉分析并分类，同时将分类结果由反馈信道传送到编码端；根据分级反馈来的分析结果对图像块进行自适应观测(不同类图像块使用不同的观测率即Φ_Wni不同)，观测值y_ni传输到解码端，与前n-1阶段的观测值共同用于整体与分块重构，以得到整幅图像的DCT系数

和块DCT系数块

根据解码端对图像的恢复质量要求决定是否重复进行视觉分析及自适应观测和重构。随着观测率的增加，每一阶段都可以得到比上一阶段更高质量的恢复图像，进而获得了图像的渐进恢复效果。

图2给出了一个8×8的DCT系数块的频率区域划分。将每个分块重构后得到的DCT系数块划分成直流(DC)、低频(LF)、中频(MF)和高频(HF)四个部分，分别标记四个部分的系数绝对值之和为DC、L、M和H。

图3则给出了块分类的流程图。其中μ₁＝125，μ₂＝900，{α₁，β₁}为{2.3，1.6}，{α₂，β₂}为{1.4，1.1}，γ＝4，κ＝290.根据每个DCT系数块中M+H以及L/M和(L+M)/H与预定阈值之间的关系，通过五种情况的判断，将每个块划分为平滑块、边沿块和纹理块。另外，进一步划分了边沿块：若L+M≤400，则为边沿1，否则为边沿2。

图4给出了观测矩阵Φ_ni的构造示意图。预先构造一个正交的随机高斯矩阵Φ作为种子观测阵。在传统的块压缩传感方法(BCS)中，观测矩阵Φ_W是Φ的前M_W行组成的矩阵；而在自适应块压缩传感(ABCS)中，观测矩阵Φ_Wni是抽取Φ中的

行组成的矩阵，MR_ni为预定的观测率。重构时，由于要用到所有的观测值，相应的观测矩阵Φ_ni则由Φ_W、Φ_W1i、......Φ_Wni组成。

本发明算法的程序流程图如图5所示，主要包括如下所述的“BCS编解码”、“视觉特性分析”和“ABCS编解码”三部分。下面以测试图像Lena为例具体说明：

I.块压缩传感(BCS)

本发明首先采用传统的块压缩传感方法(BCS)，使用一个较小的固定观测率0.1对Lena的原图像进行分块观测，然后在解码端进行分块重构和整体重构，得到初始阶段重构的DCT系数和一个基本清晰的恢复图像(PSNR＝23.35dB)，具体操作步骤为：

1.BCS编码：

第一步：读入256×256大小的Lena图像记为x，将其分成1024个8×8大小不重叠的图像块；

第二步：生成一个8²×8²大小的正交高斯随机矩阵作为种子观测阵Φ；第三步：对每一个图像块利用公式(1)进行CS观测：

y_i＝Φ_W.x_i，           (1)

其中，Φ_W是第i个图像块x_i的观测矩阵，其大小为m×8²，

由种子观测阵Φ的前6行组成，MR为设定的观测率0.1；y_i是对应x_i的观测值，所有y_i组成整幅图像x的观测值集合y；

第四步：将观测值集合y以及观测率MR发送到解码端；

2.BCS解码，利用基于DCT变换的迭代重构算法分别进行分块和整体重构，分块重构后的DCT系数用于第一阶段的视觉分析，而整体重构后得到一个低质量的恢复图像，具体步骤如下：

第一步：在解码端，预先存储一个与编码端相同的种子观测阵Φ，通过编码端发送过来的观测率0.1抽取Φ的前6行，得到观测矩阵Φ_W；

第二步：由公式(2)得到整幅图像的初始解：

$x^{\left(0\right)}={\mathrm{\Φ}}_0^{T}y,---\left(2\right)$

其中，

是由Φ_W组成的块对角矩阵，

是Φ₀的转置矩阵，y是由编码端发送过来的观测值集合；

第三步：对x^(j)(j代表迭代次数，本发明预设最大迭代次数为200，第一次迭代时为x⁽⁰⁾)进行维纳滤波，以消除图像的块效应；

第四步：对于滤波后图像中的每一个块

通过PL算法进行更新精确，如公式(3)所示：

${\hat{x}}_i^{\left(j\right)}=x_i^{\left(j\right)}+{\mathrm{\Φ}}_W^T(y_i-{\mathrm{\Φ}}_W{x}_i^{\left(j\right)});---\left(3\right)$

第五步：利用公式(4)对每个块

进行DCT变换，得到

的稀疏表示，

${{\mathrm{\θ}}_i}^{\left(j\right)}=\mathrm{\ψ}{\hat{x}}_i^{\left(j\right)},---\left(4\right)$

其中，θ_i(j)为

在DCT基下的系数，Ψ为DCT变换基；

第六步：对每个θ_i ^(j)利用公式(5)进行双变量收缩阈值处理，得到更加稀疏的系数：

其中，若g≤0，则(g)+＝0，若g＞0，则(g)+＝g；λ是一个收敛控制常数，本发明中取λ＝6；

是θ_i ^(j)的中位数估计值；

是θ_i ^(j)的边缘方差；

这里特别注意的是：第五、六步属于分块重构，是为了得到迭代后的最优系数

用来进行第一个阶段的视觉特性分析，以便指导编码端的自适应观测，所以不需要对

进行反变换和恢复处理。

第七步：利用公式(6)对

进行DCT变换，得到的稀疏表示，

${\mathrm{\θ}}^{\left(j\right)}=\mathrm{\ψ}{\hat{x}}^{\left(j\right)},---\left(6\right)$

其中，θ^(j)为

在DCT基下的系数，Ψ为DCT变换基；

第八步：对θ^(j)利用公式(7)进行双变量收缩阈值处理，得到更加稀疏的系数：

${\widehat{\mathrm{\θ}}}^{\left(j\right)}=\frac{{(\sqrt{{\left({\mathrm{\θ}}^{\left(j\right)}\right)}^2+{\left({\mathrm{\θ}}^{(j-1)}\right)}^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(j\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({\mathrm{\θ}}^{\left(j\right)}\right)}^2+{\left({\mathrm{\θ}}^{(j-1)}\right)}^2}}\·{\mathrm{\θ}}^{\left(j\right)},---\left(7\right)$

其中，若g≤0，则(g)+＝0，若g＞0，则(g)+＝g；λ是一个收敛控制常数，本发明中取λ＝6；

是θ^(j)的中位数估计值；

是θ^(j)的边缘方差；

第九步：利用公式(8)对进行反变换(IDCT)，得到本次迭代的近似解：

${\stackrel{\‾}{x}}^{\left(j\right)}={\mathrm{\ψ}}^{-1}{\widehat{\mathrm{\θ}}}^{\left(j\right)};---\left(8\right)$

第十步：对于中的每一个块

利用公式(9)进行更新：

$x_i^{(j+1)}={\stackrel{\‾}{x}}_i^{\left(j\right)}+{\mathrm{\Φ}}_W^T(y_i-{\mathrm{\Φ}}_W{\stackrel{\‾}{x}}_i^{\left(j\right)});---\left(9\right)$

这里第七到十步为整体重构，是为了得到一个基本清晰的恢复图像。

第十一步：重复进行第三到第十步，直到得到的解满足终止条件(如公式(10)所示)，停止迭代，如果到达预设的最大迭代次数时还未满足终止条件，也将停止迭代，得到初始阶段的所有块的最优DCT系数

并恢复图像R₀：

|D^(j+1)-D^(j)|＜10^-4，    (10)

其中， $D^{\left(j\right)}=\frac{1}{256}{\left|\right|x^{\left(j\right)}-{\hat{x}}^{(j-1)}\left|\right|}_2.$

II.视觉特性分析

人类视觉系统(HVS)是一个高度复杂的系统，但并不是一个完美的系统，其特性体现在：人眼对各种能量区域的敏感程度是不同的。一般而言，纹理能量高的区域不易引起视觉注意，可以少发送一些观测值以减少码率；反之，平滑区域较易引起视觉注意，可以多发送一些观测值以精确恢复。此外，由于人眼往往会刻意捕捉目标的边沿，以刻画目标形状，所以对边沿的敏感度处于平滑和纹理之间。本发明主要考虑这些特性，试图对不同类型的图像区域采用不同数量的观测值，将稀疏系数块划分为平滑块、边沿块和纹理块。

另外，为了保持CS“低复杂度编码”的特性，本发明将视觉特性分析移到解码端，利用CS可以重构在DCT域下的性质，直接对每个阶段分块重构得到的DCT系数块进行视觉分析，并利用反馈信道逐级地将分析结果传送到编码端来指导观测率的分配，从而渐进地提高图像的恢复质量。具体包括下列步骤：

第一步：在解码端，对每个重构的8×8大小的DCT系数块

进行分析，分成四个部分：直流(DC)、低频(LF)、中频(MF)和高频(HF)，如图2所示，并分别标记上述四个部分的系数绝对值之和为DC、L、M和H，其中，中、高频系数能量和(即M+H)反映了纹理特性；

第二步：根据下列五种情况将DCT系数块分成平滑块、边沿块和纹理块，如图3所示。

情况A：若M+H≤μ₁，那么这个块是平滑的，被划分为平滑块；否则进入情况B的判断，其中μ₁＝125；

情况B：若M+H＞μ₂不成立，那么进入情况C1的判断；否则进入情况C2的判断。其中μ₂＝900；

情况C1：若满足下列条件之一：①L/M≥α₁且(L+M)/H≥β₁；②L/M≥β₁且(L+M)/H≥α₁；③(L+M)/H≥γ.那么被划分为边沿块；否则进入情况D的判断。其中α₁＝2.3，β₁＝1.6，γ＝4；

情况C2：若满足下列条件之一：①L/M≥α₂且(L+M)/H≥β₂；②L/M≥β₂且(L+M)/H≥α₂；③(L+M)/H≥γ.那么被划分为边沿块；否则被划分为纹理块。其中α₂＝1.4，β₂＝1.1，γ＝4；

情况D：若M+H＞κ，那么被划分为纹理块，否则这个块是平滑的，被划分为平滑块。其中κ＝290。

此外，进一步划分边沿块：若L+M≤400，则该块判为边沿1，否则被判为边沿2。这样做的主要目的是使块分类更精确，从而准确指导编码端的自适应观测；

第三步：将块分类结果送入反馈信道；

III.视觉感知指导的自适应块压缩传感(ABCS)

本发明利用反馈信道分级传送来的块分类信息，分阶段地来指导自适应观测和重构-平滑块分配较大的观测率，其次是边沿1，再次是边沿2，最后纹理块分配较小的观测率，渐进地提高图像的恢复质量，本发明做了三个阶段的实验：第一阶段为利用观测率为0.1使用BCS方法得到的块分类结果，自适应地增加观测率和重构，等效观测率为0.336时，得到PSNR＝28.12dB的恢复图像；第二阶段为利用等效观测率为0.336使用ABCS方法得到的块分类结果，自适应增加观测率和重构，等效观测率为0.475时，得到PSNR＝30.57dB的恢复图像；第三阶段是利用等效观测率为0.475得到的块分类结果，自适应增加观测率和重构，等效观测率为0.507时，得到PSNR＝30.93dB的恢复图像。具体包括下列步骤：

1.ABCS编码：

第一步：分别根据反馈信道发送来的三个阶段的块分类结果(依次为利用观测率为0.1使用BCS方法得到的，利用等效观测率为0.336以及0.475使用ABCS方法得到的)，对不同类型的块使用不同的观测率进行自适应观测，如公式(11)所示：

y_ni＝Φ_Wnix_i，         (11)

其中，Φ_Wni是第n个阶段中第i个块x_i的观测矩阵，这里n取1，2，3，由抽取种子观测阵Φ中的

行所得，相应的观测率为MR_ni，如图4所示；对于平滑块，MR_ni取较大的值；对纹理块，MR_ni取较小的值；而对边沿块，MR_ni则取中间值，其中边沿1的MR_ni值略大于边沿2的；

第二步：计算等效观测率EMR：用所有块的全部观测率之和(包括之前阶段的观测率)除以块数1024，即求平均(三个阶段依次得到0.336，0.475，0.507)；

这里求EMR主要是作为传统BCS算法的观测率，以便于与本发明的性能作比较；

第三步：将自适应观测值y_ni与相应的观测率MR_ni发送到解码端；

2.ABCS解码，利用基于DCT变换的迭代重构算法分三阶段进行分块和整体重构，分块重构后的DCT系数用于视觉分析，而整体重构后得到三幅质量渐进提高的恢复图像，具体步骤如下：

第一步：在解码端，将当前接收到的每一个块的所有观测值(即：包含存储器中之前阶段发送过来的块观测值)组合起来，如公式(12)：

$Y_{\mathrm{ni}}=\left(\begin{array}{c}{y}_{\mathrm{ni}}\\ \·\\ \·\\ \·\\ y_{1i}\\ y_i\end{array}\right),---\left(12\right)$

其中，Y_ni表示第i个图像块在n个阶段中的所有观测值，将所有Y_ni组成整幅图像的观测值集合Y_n。同时，利用编码端发送来的观测率MR_ni得到相应的观测矩阵Φ_Wni(Φ_Wni的抽取方法跟公式(9)中相同)，并将其组合起来，得到对应于Y_ni的观测阵Φ_ni，如公式(13)所示：

${\mathrm{\Φ}}_{\mathrm{ni}}=\left(\begin{array}{c}{\mathrm{\Φ}}_{\mathrm{Wni}}\\ \·\\ \·\\ \·\\ {\mathrm{\Φ}}_{W1i}\\ {\mathrm{\Φ}}_W\end{array}\right);---\left(13\right)$

第二步：将上一阶段恢复的图像R_n-1(即PSNR依次为23.35dB，28.12dB，30.57dB的恢复图像R₀，R₁，R₂)作为当前阶段迭代的初始解，以大大加快重构的速度；

第三步：对x^(j)(j代表迭代次数，本发明预设最大迭代次数为200，第一次迭代时x^(j)＝R₀)进行维纳滤波，以消除图像的块效应；

第四步：对于滤波后的x^(j)中的每一个块

通过PL算法进行更新精确，如公式(14)所示：

${\hat{x}}_i^{\left(j\right)}=x_i^{\left(j\right)}+{\mathrm{\Φ}}_{\mathrm{ni}}^T(Y_{\mathrm{ni}}-{\mathrm{\Φ}}_{\mathrm{ni}}{x}_i^{\left(j\right)});---\left(14\right)$

第五步：利用公式(4)对每个块

进行DCT变换，得到

的稀疏表示θ_i ^(j)；

第六步：利用公式(5)对每个θ_i ^(j)进行双变量收缩阈值处理，得到更加稀疏的系数

特别注意的是：第五、六步属于分块重构，是为了得到迭代后的最优系数

如果图像恢复质量仍不满足要求，将继续用来进行视觉特性分析，以便指导编码端下一阶段的自适应观测和重构，也不需要对

进行反变换和恢复处理。

第七步：利用公式(6)对

进行DCT变换，得到

的稀疏表示θ^(j)；

第八步：利用公式(7)对θ^(j)进行双变量收缩阈值处理，得到更加稀疏的系数

第九步：通过公式(8)对

进行反变换(IDCT)，得到本次迭代的近似解

第十步：对于

中的每一个块

通过公式(15)进行更新：

$x_i^{(j+1)}={\stackrel{\‾}{x}}_i^{\left(j\right)}+{\mathrm{\Φ}}_{\mathrm{ni}}^T(Y_{\mathrm{ni}}-{\mathrm{\Φ}}_{\mathrm{ni}}{\stackrel{\‾}{x}}_i^{\left(j\right)});---\left(15\right)$

这里第七到十步为整体重构，是为了得到质量比上一阶段提高的恢复图像。

第十一步：重复进行第三到第十步，直到得到的解满足终止条件(如公式(10)所示)，停止迭代，如果到达预设的最大迭代次数时还未满足终止条件，也将停止迭代，得到本阶段的所有块的最优DCT系数并恢复图像R_n；

第十二步：根据解码端对恢复图像质量的要求，决定是否重复进行视觉特性分析和视觉感知指导的自适应块压缩传感(ABCS)，得到质量渐进提高的恢复图像。由于第一阶段和第二阶段得到的Lena恢复图像(PSNR分别为28.12dB和30.57dB)仍不满足质量要求，所以继续进行视觉分析，自适应地增加观测率，完成第三阶段的重构，恢复的图像(PSNR＝30.93dB)满足了视觉质量的要求，采样结束。

我们对本发明提出的渐进图像编码方案做了初步测试实验，我们采用三幅数字图像处理的标准测试图像(Lena、Cameraman和Barbara)做输入图像。假设一无损传输信道。使用联想ThinkPadE40笔记本电脑作算法处理，笔记本参数为：Intel(R)，Core(TM)i3 CPU，M 350，2.27GHz，2.27GHz，2.00GB内存。软件平台为MatlabR2008a，用Matlab语言编程实现了渐进图像编码方案。

在实验中，我们使用了峰值信噪比(PSNR)和频域结构相似度(FSSM)共同作为实验结果的质量评价测度。频域结构相似度(FSSM)源于近年来著名的结构相似度评价方法(SSIM)，将频域信息作为图像的主要结构信息，根据人眼对不同频率分量的敏感程度不同，对DCT变换后的各频率分量加权后得到图像的频域函数，然后由频域函数、亮度函数和对比度函数计算得到结构相似度。其性能更符合人眼视觉特性。具体计算方法如下：

FSSM(x，y)＝[l(x，y)]^α[c(x，y)]^β[s_F(x，y)]^γ，   (14)

$l(x,y)=\frac{2{\mathrm{\μ}}_x{\mathrm{\μ}}_y+C_1}{{{\mathrm{\μ}}_x}^2+{{\mathrm{\μ}}_y}^2+C_1},---\left(15\right)$

$C(x,y)=\frac{2{\mathrm{\σ}}_x{\mathrm{\σ}}_y+C_2}{{{\mathrm{\σ}}_x}^2+{{\mathrm{\σ}}_y}^2+C_2},---\left(16\right)$

$s_F(x,y)=\frac{2{\mathrm{\σ}}_{\mathrm{Fxy}}+C_3}{{{\mathrm{\σ}}_{\mathrm{Fx}}}^2+{{\mathrm{\σ}}_{\mathrm{Fy}}}^2+C_3},---\left(17\right)$

其中，x和y分别为参考图像和测试图像中的对应位置的块；l(x，y)、c(x，y)和s_F(x，y)分别为x和y的亮度函数、对比度函数和频域结构度函数；α＞0、β＞0和γ＞0是控制3部分权重的参数，本发明中设α＝β＝γ＝1。μ_x和μ_y分别为x和y的均值；σ_x和σ_y分别为x和y的标准差；σ_Fx和σ_Fy分别为x和y的频域加权标准差；σ_Fxy为x和y的频域加权协方差。计算方法如下：

μ_F＝∑_u，v w_uvF_AC(u，v)，                            (18)

σ_F ²＝∑_u，v w_uv(F_AC(u，v)-μ_F)²，                   (19)

σ_Fxy＝∑_u，v w_uv(F_ACx(u，v)-μ_Fx)(F_ACy(u，v)-μ_Fy)，(20)

其中，F_AC(u，v)是对图像块进行DCT变换后的交流系数，w_uv是对频率F_AC(u，v)的感知加权系数，∑_u，v w_uv＝1.

最后，对各图像块的FSSM利用公式(21)求均值可得到整幅图像的基于频域的结构相似度：

$\mathrm{FSSM}(X,Y)=\frac{1}{M}{\mathrm{\Σ}}_{i=1}^M\mathrm{FSSM}(x_i,y_i),---\left(21\right)$

其中，M为整幅图像中图像块的个数。

实验中，我们将图像分成了8×8大小的块，测试了在相同的观测率(等效观测率)下，本发明方案和传统BCS算法的对比，包括图像PSNR、FSSM以及平滑块和边沿块的平均PSNR、FSSM等性能参数的对比。

图6给出了三幅测试图像在不同等效观测率下，BCS算法和本发明方案恢复图像放大区域的对比图，放大区域主要包括平滑块和边沿块，即原图像中红色方框内的区域。可以明显看出，在相同的等效观测率下，本发明方案恢复的图像有更好的视觉质量。

另外，图7给出了测试图像Lena在不同等效观测率下的块分类图。黑色代表平滑块，灰色代表纹理块，白色代表边沿块。可以看到，随着观测率的增加，分类图越来越接近于原图像的分类图，也就是说，块分类越来越准确，恢复图像质量也得到了渐进的提高。

表1给出了上述三幅图像的PSNR、FSSM以及平滑块和边沿块的平均PSNR、FSSM的测试比较结果。可以看出，本发明方案比BCS算法均有明显的提高，尤其是对于视觉上比较敏感的平滑块和边沿块。特别注意的是，对Lena图像，在等效观测率为0.507时，本发明方案得到的整幅恢复图像的PSNR比使用BCS算法的下降了约0.1dB，但是FSSM仍有着一定的提高，仍然证明了本发明方案恢复的图像具有更好的主观和客观质量。

表1

Claims (1)

1.一种基于视觉感知的自适应块压缩传感图像编码方法，是将视觉特性应用在压缩传感(CS)中，利用CS重构在某种变换域下的特性，在解码端直接对每次重构的变换域系数进行视觉分析，并通过反馈信道逐级地将分析结果传送到编码端来指导观测率分配，其特征在于具体包括下列步骤：

I.块压缩传感(BCS)，具体包括下列步骤：

1.BCS编码：

(1).读入一副N×N大小图像x，将其分成多个W×W大小不重叠的图像块；

(2).生成一个W²×W²大小的正交高斯随机矩阵作为种子观测阵Φ；

(3).对每一个图像块x_i利用公式(1)进行CS观测：

y_i＝Φ_W.x_i，        (1)

其中，Φ_W是第i个图像块的观测矩阵，其大小为m×W²，由种子观测阵Φ的前m行组成，

MR为预设的观测率；y_i是对应x_i的观测值，所有y_i组成整幅图像x的观测值集合y；

(4).将观测值集合y以及观测率MR发送到解码端；

2.BCS解码，具体步骤如下：

(1).在解码端，预先存储一个与编码端相同的种子观测阵Φ，通过编码端发送过来的观测率MR抽取Φ的前m行，得到观测矩阵Φ_W；

(2).由公式(2)得到整幅图像的初始解：

$x^{\left(0\right)}={\mathrm{\Φ}}_0^{T}y,---\left(2\right)$

其中，Φ₀是由Φ_W组成的块对角矩阵，

是Φ₀的转置矩阵，y是由编码端发送过来的观测值集合；

(3).对x^(j)(j代表迭代次数，第一次迭代时为x⁽⁰⁾)进行维纳滤波，以消除图像的块效应；

(4).对于滤波后的图像中每一个块

通过PL算法进行更新，如公式(3)所示：

${\hat{x}}_i^{\left(j\right)}=x_i^{\left(j\right)}+{\mathrm{\Φ}}_W^T(y_i-{\mathrm{\Φ}}_W{x}_i^{\left(j\right)});---\left(3\right)$

(5).利用公式(4)对每个块进行DCT变换，得到

的稀疏表示，

${{\mathrm{\θ}}_i}^{\left(j\right)}=\mathrm{\ψ}{\hat{x}}_i^{\left(j\right)},---\left(4\right)$

其中，θ_i ^(j)为

在DCT基下的系数，Ψ为DCT变换基；

(6).对每个θ_i ^(j)利用公式(5)进行双变量收缩阈值处理，得到更加稀疏的系数：

${{\widehat{\mathrm{\θ}}}_i}^{\left(j\right)}=\frac{{(\sqrt{{\left({{\mathrm{\θ}}_i}^{\left(j\right)}\right)}^2+{\left({{\mathrm{\θ}}_i}^{(j-1)}\right)}^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(j\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({{\mathrm{\θ}}_i}^{\left(j\right)}\right)}^2+{\left({{\mathrm{\θ}}_i}^{(j-1)}\right)}^2}}\·{{\mathrm{\θ}}_i}^{\left(j\right)},---\left(5\right)$

其中，若g≤0，则(g)+＝0，若g＞0，则(g)+＝g；λ是一个收敛控制常数；

是θ_i ^(j)的中位数估计值；

是θ_i ^(j)的边缘方差；

(7).利用公式(6)对

进行DCT变换，得到

的稀疏表示，

${\mathrm{\θ}}^{\left(j\right)}=\mathrm{\ψ}{\hat{x}}^{\left(j\right)},---\left(6\right)$

其中，θ^(j)为

在DCT基下的系数，Ψ为DCT变换基；

(8).对θ^(j)利用公式(7)进行双变量收缩阈值处理，得到更加稀疏的系数：

${\widehat{\mathrm{\θ}}}^{\left(j\right)}=\frac{{(\sqrt{{\left({\mathrm{\θ}}^{\left(j\right)}\right)}^2+{\left({\mathrm{\θ}}^{(j-1)}\right)}^2}-\mathrm{\λ}\frac{\sqrt{3{\mathrm{\σ}}^{\left(j\right)}}}{{\mathrm{\σ}}_{\mathrm{\ξ}}})}_{+}}{\sqrt{{\left({\mathrm{\θ}}^{\left(j\right)}\right)}^2+{\left({\mathrm{\θ}}^{(j-1)}\right)}^2}}\·{\mathrm{\θ}}^{\left(j\right)},---\left(7\right)$

其中，若g≤0，则(g)+＝0，若g＞0，则(g)+＝g；λ是一个收敛控制常数；是θ^(j)的中位数估计值；

是θ^(j)的边缘方差；

(9).利用公式(8)对

进行反变换(IDCT)，得到本次迭代的近似解：

${\stackrel{\‾}{x}}^{\left(j\right)}={\mathrm{\ψ}}^{-1}{\widehat{\mathrm{\θ}}}^{\left(j\right)};---\left(8\right)$

(10).对于

中的每一个块

利用公式(9)进行更新：

$x_i^{(j+1)}={\stackrel{\‾}{x}}_i^{\left(j\right)}+{\mathrm{\Φ}}_W^T(y_i-{\mathrm{\Φ}}_W{\stackrel{\‾}{x}}_i^{\left(j\right)});---\left(9\right)$

(11).重复进行步骤(3)-(10)，直到得到的解满足终止条件(如公式(10)所示)，停止迭代，得到初始阶段的所有块的最优DCT系数

并恢复图像R₀：

|D^(j+1)-D^(j)|＜10^-4，      (10)

其中， $D^{\left(j\right)}=\frac{1}{N}{\left|\right|x^{\left(j\right)}-{\hat{x}}^{(j-1)}\left|\right|}_2.$

II.视觉特性分析，具体包括下列步骤：

(1).将每个重构的W×W大小的DCT系数块

分成四个部分：直流(DC)、低频(LF)、中频(MF)和高频(HF)，并分别标记四个部分的系数绝对值之和为DC、L、M和H；

(2)根据下列五种情况将DCT系数块分成平滑块、边沿块和纹理块。

情况A：若M+H≤μ₁，那么这个块是平滑的，被划分为平滑块；否则进入情况B的判断，其中μ₁＝125；

情况B：若M+H＞μ₂不成立，那么进入情况C1的判断；否则进入情况C2的判断。其中μ₂＝900；

情况C1：若满足下列条件之一：①L/M≥α₁且(L+M)/H≥β₁；②L/M≥β₁且(L+M)/H≥α₁；③(L+M)/H≥γ.那么被划分为边沿块；否则进入情况D的判断。其中α₁＝2.3，β₁＝1.6，γ＝4；

情况C2：若满足下列条件之一：①L/M≥α₂且(L+M)/H≥β₂；②L/M≥β₂且(L+M)/H≥α₂；③(L+M)/H≥γ.那么被划分为边沿块；否则被划分为纹理块。其中α₂＝1.4，β₂＝1.1，γ＝4；

情况D：若M+H＞κ，那么被划分为纹理块，否则这个块是平滑的，被划分为平滑块。其中κ＝290。

此外，进一步划分边沿块：若L+M≤400，则该块判为边沿1，否则被判为边沿2；

(3).将块分类结果送入反馈信道。

III.视觉感知指导的自适应块压缩传感(ABCS)，具体包括下列步骤：

1.ABCS编码：

(1).根据每个阶段反馈信道发送来的块分类结果，对不同类型的块使用不同的观测率进行自适应观测，如公式(11)所示：

y_ni＝Φ_Wnix_i，           (11)

其中，Φ_Wni是第n个阶段中第i个块x_i的观测矩阵，由抽取种子观测阵Φ中的

行所得，相应的观测率为MR_ni；对于平滑块，MR_ni取较大的值；对纹理块，MR_ni取较小的值；而对边沿块，MR_ni则取中间值，其中边沿1的MR_ni值略大于边沿2的；

(2).将自适应观测值y_ni与相应的观测率MR_ni发送到解码端；

2.ABCS解码，具体步骤如下：

(1).在解码端，将当前接收到的每一个块的所有观测值(即：包含存储器中之前阶段发送过来的块观测值)组合起来，如公式(12)：

$Y_{\mathrm{ni}}=\left(\begin{array}{c}{y}_{\mathrm{ni}}\\ \·\\ \·\\ \·\\ y_{1i}\\ y_i\end{array}\right),---\left(12\right)$

其中，Y_ni表示第i个图像块在n个阶段中的所有观测值，将所有Y_ni组成整幅图像的观测值集合Y_n。同时，利用编码端发送来的观测率MR_ni得到相应的观测矩阵Φ_Wni(Φ_Wni的抽取方法跟公式(9)中相同)，并将其组合起来，得到对应于Y_ni的观测阵Φ_ni，如公式(13)所示：

${\mathrm{\Φ}}_{\mathrm{ni}}=\left(\begin{array}{c}{\mathrm{\Φ}}_{\mathrm{Wni}}\\ \·\\ \·\\ \·\\ {\mathrm{\Φ}}_{W1i}\\ {\mathrm{\Φ}}_W\end{array}\right);---\left(13\right)$

(2).将上一阶段恢复的图像R_n-1作为当前阶段迭代的初始解，以加快重构的速度；

(3).对x^(j)(j代表迭代次数，第一次迭代时x^(j)＝R₀)进行维纳滤波，以消除图像的块效应；

(4).对于滤波后的x^(j)中的每一个块

通过PL算法进行更新精确，如公式(14)所示：

${\hat{x}}_i^{\left(j\right)}=x_i^{\left(j\right)}+{\mathrm{\Φ}}_{\mathrm{ni}}^T(Y_{\mathrm{ni}}-{\mathrm{\Φ}}_{\mathrm{ni}}{x}_i^{\left(j\right)});---\left(14\right)$

(5).利用公式(4)对每个块

进行DCT变换，得到

的稀疏表示θ_i ^(j)；

(6).利用公式(5)对每个θ_i ^(j)进行双变量收缩阈值处理，得到更加稀疏的系数

(7).利用公式(6)对

进行DCT变换，得到

的稀疏表示θ^(j)；

(8).利用公式(7)对θ^(j)进行双变量收缩阈值处理，得到更加稀疏的系数

(9).通过公式(8)对

进行反变换(IDCT)，得到本次迭代的近似解

(10).对于中的每一个块

通过公式(15)进行更新：

$x_i^{(j+1)}={\stackrel{\‾}{x}}_i^{\left(j\right)}+{\mathrm{\Φ}}_{\mathrm{ni}}^T(Y_{\mathrm{ni}}-{\mathrm{\Φ}}_{\mathrm{ni}}{\stackrel{\‾}{x}}_i^{\left(j\right)});---\left(15\right)$

(11).重复进行步骤(3)-(10)，直到得到的解满足终止条件(如公式(10)所示)，停止迭代，得到本阶段的所有块的最优DCT系数

并恢复图像R_n；

(12).根据解码端对恢复图像质量的要求，决定是否重复进行视觉特性分析和视觉感知指导的自适应块压缩传感(ABCS)，得到质量渐进提高的恢复图像。

Images