圖像的分解與合成

實驗八　　圖像的分解與合成

　　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 = \sum\limits_{i = 1}^m {p(d_i) \log_2 {p(d_i)}} \]

單位是bit/灰度級，圖像的熵表示像素的各個灰度級為數的統計平均值，它給出了對此輸入灰度級集合進行編碼時所需的平均位數的下限。

　　（2）平均字碼長度：設 \( \beta_i \) 為數字圖像中灰度級 \( d_i \) 對應的碼字長度（二進制代碼的位數）。其相應出現的概率為 \( p(d_i) \) ，則該數字圖像賦予的平均碼字長度為 \[  R = \sum\limits_{i = 1}^m {\beta_i p(d_i)} \]

　　（3）編碼效率：當H接近R時，編碼效果好，當R遠大於H時，則編碼效果差 \[ \eta = \frac{H}{R} \times  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

　　可計算信源的熵、平均碼長、編碼效率及冗余度如下：

\begin{aligned}
　　\text { 熵: } 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\\
　　\text { 平均碼長: } R(X) &= \sum_{k=1}^{6} \beta_{k} p_{k} \\
　　&=2 \times 0.25 +2 \times 0.25 +2 \times 0.20 +3 \times 0.15 +4 \times 0.10+4 \times 0.05 \\
　　&=2.45 \\
　\text {編碼效率: } \eta &= \frac{H(x)}{R(x)}=\frac{2.42}{2.45}\times 100\%= 98.8\% \\
　\text {冗余度: } 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 \times 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 \right \}^{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\times\lg\frac{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 \times 8 \) 或  \( 16 \times 16 \) 。

　　3. 變換系數的編碼：對圖像子塊進行變換后，得到的是變換系。可采用“區域編碼”或者“閾值編碼”等。

　　（1）區域編碼：選出能量集中的區域，舍棄其他區域（對其置零）。可采用不同的量化器。

　　（2）閾值編碼：記變換系數 \( F(u,v) \) ，量化和去門限后的結果 \( \hat{F}(u,v) \) ，量化矩陣中的相應元素 \( S(u,v) \)，則 \[ \hat{F}(u,v)=INT\left [ \frac{F(u,v)}{S(u,v)} \pm 0.5 \right ]   \] 閾值編碼中保留系數的位置是可變的，因此需要對這些系數的位置信息和幅度信息一起編碼和傳輸，才能在接受端正確恢復圖像。常用的對位置編碼方法是基於對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;
}