实验八 图像的分解与合成
1）对一图像按块进行离散余弦变换（即分成8×8的小块，对每块进行DCT），利用JPEG建议的两化矩阵对DCT系数进行量化，观察8×8小块交换系数。
2）对变换系数进行区域选择，然后进行逆量化和逆变换（IDCT），重建原图像，计算重建图像的PSNR及图像压缩的压缩比。
3）选择一个小波函数，对一幅图像进行二级分解，显示各变换系数，并进行小波反变换，重建原图像。
图像的压缩编码的概述
一副640×480分辨率的彩色图像（24bit/像素），其数据量为900KB，如果以每秒30帧的速度播放，则1秒钟播放的数据量为：640×480×24bit×30帧=210Mbit=26.4MB，需要210Mbit/s的通信回路；如果存放在650MB的光盘中，在不考虑音频信号的情况下，每张光盘也只能存储24s的视频。对于电话线传送二值图像的传真，如果以每秒200dpi(点/英寸)的分辨率传输，一张A4稿纸内容的数据量为（200×210/25.4）×（200×297/25.4）bit=3866948bit，按14.kbit/s的电话线传输速率，需要传送的时间是263s(4.4min)。可见，研究图像压缩是有必要的。图像压缩的理论基础是信息论。从信息论的角度来看，压缩是去掉信息中冗余成分，即保留不确定信息，去掉确定信息。
（1）空域冗余：如果用水平方向的任何一行像素预测垂直方向的其他行像素，都能准确预测出其他行数据，或者其他行数据完全能够用一行数据复制得到。
（2）时域冗余：两帧图像越接近，说明图像序列携带的动态信息越少，换言之，第2帧相对第1帧而言，存在大量的时域冗余。对于视频压缩而言，时域冗余是最主要的冗余。
（3）频域冗余：大多数图像的频谱具有低通特性，低频部分的系数能够提供绝大部分的图像信息。保留低频部分的系数，去掉高频部分的系数，可保持大部分图像能量。
（4）信息熵冗余：图像中像素灰度分布的不均匀性导致信息熵冗余。若将出现概率大的灰度级用长度较短的码表表示，将出现概率小的灰度级用长度较长的码表示，有可能使编码的平均长度下降。
（5）接收端图像设备分辨率较低，则可降低图像分辨率。
（6）用户关心的图像区域有限，可对其余部分图像采用空间和灰度级上的粗化。
（7）根据人眼的视觉特性，对不敏感区域进行降分辨率编码。
20世纪80年代以前，图像压缩编码主要是根据传统的信息源编码方法，研究的内容是有关信息熵、编码方法及数据压缩比；20世纪20年代以后，它突破了信息源编码理论，结合分形、模型基、神经网络、小波变换等数学工具，充分利用视觉系统和图像信源的各种特性。对图像压缩编码的分类，从解码结果对图像的保真程度，可将图像的压缩编码分为两大类：
无损压缩编码（冗余度压缩、可逆压缩）——是一种在解码时精确第恢复原图像，没有任何损失的编码方法，但压缩比不大，通常只能获得1~5倍的压缩率。
有损压缩编码（熵压缩、不可逆压缩）——解码时只能近似原图像，不能无失真地恢复原图像，压缩比大，但有信息损失。
从具体编码技术出发：预测编码、变换编码、统计编码、轮廓编码、模型编码等。
压缩编码系统评价
（1）图像熵：设数字图像像素的灰度级集合为有 ({d_1,d_2,...,d_m}) ，( p(d_1),p(d_2),...,p(d_m) )为其对应的概率。按信息论中信息源信息熵的定义，图像的熵定义为 [ H = sumlimits_{i = 1}^m {p(d_i) log_2 {p(d_i)}} ]
单位是bit/灰度级，图像的熵表示像素的各个灰度级为数的统计平均值，它给出了对此输入灰度级集合进行编码时所需的平均位数的下限。
（2）平均字码长度：设 ( eta_i ) 为数字图像中灰度级 ( d_i ) 对应的码字长度（二进制代码的位数）。其相应出现的概率为 ( p(d_i) ) ，则该数字图像赋予的平均码字长度为 [ R = sumlimits_{i = 1}^m {eta_i p(d_i)} ]
（3）编码效率：当H接近R时，编码效果好，当R远大于H时，则编码效果差 [ eta = frac{H}{R} imes 100 \% ]
（4）压缩比：编码前后平均码长之比，记n为编码前每个灰度级所用的平均码长，通常用自然二进制码表示时的比特数，如灰度级为256时，n为8bit，二进制码长度为3。即 [ r = frac{n}{R} ]
基于压缩编码参数的基本评价：最好的编码结果应使R等于或很接近于H，这样的编码方法，称为最佳编码，它既不丢失信息，又占用很少的比特数，比如后面将要介绍的霍夫曼编码。若要求编码结果 ( R<H ) 则必然会丢失信息而引起图像失真。一般而言，压缩比大，则说明被压缩掉的数据量大。一个编码系统要研究的问题是设法减小编码平均长度R，使编码效率 ( eta ) 尽量趋于1，而冗余度尽量趋于0.举例来说，一个符号的信源X，其霍夫曼编码为：
| 符号 | ( u_1 ) | ( u_2 ) | ( u_3 ) | ( u_4 ) | ( u_5 ) | ( u_6 ) |
|
概率 |
0.25 |
0.25 | 0.20 | 0.15 | 0.10 | 0.05 |
| 码字 | 01 | 10 | 11 | 000 | 0010 | 0011 |
可计算信源的熵、平均码长、编码效率及冗余度如下：
egin{aligned}
ext { 熵: } H(X) &= - sum_{k=1}^{6} log_{2} p_{k} \
&=-0.25 log_2 {0.25} -0.25 log_2 {0.25} -0.20 log_2 {0.20} - 0.15 log_2 {0.15} - \
& qquad 0.10 log_2 {0.10} - 0.05 log_2 {0.05} \
&=2.42\
ext { 平均码长: } R(X) &= sum_{k=1}^{6} eta_{k} p_{k} \
&=2 imes 0.25 +2 imes 0.25 +2 imes 0.20 +3 imes 0.15 +4 imes 0.10+4 imes 0.05 \
&=2.45 \
ext {编码效率: } eta &= frac{H(x)}{R(x)}=frac{2.42}{2.45} imes 100\%= 98.8\% \
ext {冗余度: } r &=1-eta =1-98.8\%= 1.2\%
end{aligned}
可见，对于上述信源X的霍夫曼编码，其编码效率已达 ( 98.8\% ) ，只有 ( 1.2\% ) 的冗余度。
基于保真度准则的评价：可分为客观保真度准则和主观保真度准则（这里主要介绍客观保真度准则）。其中常用的客观保真度准则有输入图像与输出图像的方均根误差、输入图像与输出图像的方均根信噪比和峰值信噪比三种。
方均根误差的定义：令 ( f(x,y) ) 代表输入图像，( hat{f} (x,y) ) 代表对 ( f(x,y) ) 先压缩又解压后得到的图像，对于任意的x和y， ( f(x,y) ) 和 ( hat{f} (x,y) ) 的误差定义为 [ e(x,y)=hat{f}(x,y)-f(x,y) ] 若两幅图像大小均为 ( M imes N ) ，则它们之间的总误差为 [ sum_{x=0}^{M-1}sum_{y=0}^{N-1}[hat{f}(x,y)-f(x,y)] ] 这样 ( f(x,y) ) 和 ( hat{f} (x,y) ) 之间的方均根误差 ( e_{rms} ) 为 [ e_{rms}=left { frac{1}{MN} sum_{x=0}^{M-1}sum_{y=0}^{N-1}[hat{f}(x,y)-f(x,y)]^2 ight }^{1/2} ]
方均根信噪比的定义：如果将( hat{f} (x,y) ) 看作原始图像 ( f(x,y) ) 和噪声图像 ( e(x,y) ) 之和，那么输出图像的方均信噪比 ( SNR_{ms} ) 为 [SNR_{ms}=sum_{x=0}^{M-1}sum_{y=0}^{N-1}hat{f}(x,y)^2 / sum_{x=0}^{M-1}sum_{y=0}^{N-1}[hat{f}(x,y)-f(x,y)]^2 ]如果对上式求平方根，就得到方均根信噪比 ( SNR_{ms} ) ，实际使用时常将 SNR 归一化并用分贝（dB）表示。
峰值信噪比的定义：如果令 ( f_{max}=max{ f(x,y),x=0,1,cdots,M-1, y=0,1,cdots,N-1} ) ，即 ( f_{max} ) 为图像灰度的最大值，则峰值信噪比PSNR定义为 [ PSNR=10 imeslgfrac{f_{max}^2}{frac{1}{MN}sum_{x=0}^{M-1}sum_{y=0}^{N-1}[hat{f}(x,y)-f(x,y)]^2 } ]
变换编码
变换编码的方法很多，如离散傅里叶变换（Discrete Fourier Transform，DFT）、离散余弦变换（Discrete Cosine Transform，DCT）、离散哈达玛变换（Discrete Hadamard Transform，DHT）等。变换编码的基本原理是将原来在空域描述的图像信号，变换到另外一些正交空间中去，用变换系数来表示原始图像，并对变换系数进行编码。变换本身并不产生数据压缩，只是提供用于编码压缩的变换系数矩阵，只有对变换系数采用量化和熵编码才能产生压缩作用。统计表明，在变换域中，图像信号的绝大部分能量集中在低频部分，编码中如果略去那些能量很小的高频分量，或者给这写高频分量分配合适的比特数，就可以明显减小图像传输或存储的数据量。
残余相关准则，变换域内变换系数具有相关性称为残余相关性，它代表经过正交变换后图像相关性被消弱的程度。显然，如果变换方法及有关参数选择恰当，则变换系数矩阵所对应的相关系数会很小，原图像的相关性可得到充分的消除，冗余度可得到很好的压缩。
方均误差准则，是一种将解码后的重建图像与未经压缩的原始图像之间的方均误差作为衡量各种正交变换效果的准则。
主观评价标准，人眼辨别。
上述几种评价标准中，三种准则都与正交变换本身的特性有关，但是后两种准则还与变换系数样本的选择及量化特性有关，因此可将方均误差准则和主观评价准则看做整个正交变换编码系统特性的评价标准。
1. 选择变换方法：理论上，K-L变换是最优的正交变换，但是由于K-L基的不确定性，使用起来会很不方便，只能作为理论上的比较标准。实际上离散余弦变换是图像压缩中最常用的，它的性能最接近K-L变换，并且具有快速算法。
2. 确定子块图像的大小：实际应用中，图像子块一般取 ( 8 imes 8 ) 或 ( 16 imes 16 ) 。
3. 变换系数的编码：对图像子块进行变换后，得到的是变换系。可采用“区域编码”或者“阈值编码”等。
（1）区域编码：选出能量集中的区域，舍弃其他区域（对其置零）。可采用不同的量化器。
（2）阈值编码：记变换系数 ( F(u,v) ) ，量化和去门限后的结果 ( hat{F}(u,v) ) ，量化矩阵中的相应元素 ( S(u,v) )，则 [ hat{F}(u,v)=INTleft [ frac{F(u,v)}{S(u,v)} pm 0.5 ight ] ] 阈值编码中保留系数的位置是可变的，因此需要对这些系数的位置信息和幅度信息一起编码和传输，才能在接受端正确恢复图像。常用的对位置编码方法是基于对0值（非保留系数）的游程长度的编码，也就是通过指明一个保留系数之前有多少个0，来确保该系数的位置。
栈溢出的代码：
#include <opencv2/core.hpp> #include <opencv2/highgui.hpp> #include <iostream> #include <math.h> #include <complex> using namespace std; using namespace cv; const int height = 128, width = 128, channel = 3; const double M_PI = 3.1415926535; // DCT hyper-parameter int T = 8; int K = 8; // DCT coefficient struct dct_str { double coef[height][width][channel]; // 栈溢出的罪魁祸首 }; // Discrete Cosine transformation dct_str dct(Mat img, dct_str dct_s) { double I; double F; double Cu, Cv; for (int ys = 0; ys < height; ys += T) { for (int xs = 0; xs < width; xs += T) { for (int c = 0; c < channel; c++) { for (int v = 0; v < T; v++) { for (int u = 0; u < T; u++) { F = 0; if (u == 0) { Cu = 1. / sqrt(2); } else { Cu = 1; } if (v == 0) { Cv = 1. / sqrt(2); } else { Cv = 1; } for (int y = 0; y < T; y++) { for (int x = 0; x < T; x++) { I = (double)img.at<Vec3b>(ys + y, xs + x)[c]; F += 2. / T * Cu * Cv * I * cos((2. * x + 1) * u * M_PI / 2. / T) * cos((2. * y + 1) * v * M_PI / 2. / T); } } dct_s.coef[ys + v][xs + u][c] = F; } } } } } return dct_s; } // Inverse Discrete Cosine transformation Mat idct(Mat out, dct_str dct_s) { double f; double Cu, Cv; for (int ys = 0; ys < height; ys += T) { for (int xs = 0; xs < width; xs += T) { for (int c = 0; c < channel; c++) { for (int y = 0; y < T; y++) { for (int x = 0; x < T; x++) { f = 0; for (int v = 0; v < K; v++) { for (int u = 0; u < K; u++) { if (u == 0) { Cu = 1. / sqrt(2); } else { Cu = 1; } if (v == 0) { Cv = 1. / sqrt(2); } else { Cv = 1; } f += 2. / T * Cu * Cv * dct_s.coef[ys + v][xs + u][c] * cos((2. * x + 1) * u * M_PI / 2. / T) * cos((2. * y + 1) * v * M_PI / 2. / T); } } f = fmin(fmax(f, 0), 255); out.at<Vec3b>(ys + y, xs + x)[c] = (uchar)f; } } } } } return out; } // Main int main(int argc, const char* argv[]) { // read original image Mat img = imread("1.jpg", IMREAD_COLOR); if (img.empty()) { cout << "Image empty!" << endl; return -1; } // DCT coefficient dct_str dct_s; // output image Mat out = Mat::zeros(height, width, CV_8UC3); // DCT dct_s = dct(img, dct_s); // IDCT out = idct(out, dct_s); imwrite("out.jpg", out); //imshow("answer", out); //waitKey(0); destroyAllWindows(); return 0; }