基于塊變換與哈爾小波編碼的

1、醫(yī)學(xué)圖像處理實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)十一：基于塊變換和哈爾小波編碼的圖像壓縮與重建摘要 本次實(shí)驗(yàn)的實(shí)驗(yàn)?zāi)康募爸饕獌?nèi)容是： 基于塊變換的圖像壓縮重建 基于小波變換的圖像壓縮重建一、技術(shù)討論1.1實(shí)驗(yàn)原理圖像信息熵：它表征圖像灰度分布的聚集特性，卻不能反映圖像灰度分布的空間特征，為了表征這種空間特征，可以在一維熵的基礎(chǔ)上引入能夠反映灰度分布空間特征的特征量來(lái)組成圖像的二維熵。一階熵：相鄰兩個(gè)像素的像素值和作為統(tǒng)計(jì)標(biāo)準(zhǔn)；二階熵：相鄰三個(gè)像素的像素值和作為統(tǒng)計(jì)標(biāo)準(zhǔn)。保真度準(zhǔn)則：表征電子設(shè)備輸出再現(xiàn)輸入信號(hào)的相似程度。保真度越高，無(wú)線電接收機(jī)輸出的聲音或電視機(jī)輸出的影像越逼真。塊變換編碼:映射變換：函數(shù)變換，常用的

2、稱為正交變換-量化-編碼小波編碼：小波變換（不損失數(shù)據(jù)，它是EZW算法具有漸進(jìn)性的基礎(chǔ)）-量化（對(duì)數(shù)據(jù)產(chǎn)生損失，損失程度取決于量化閾值的大?。?熵編碼（對(duì)每個(gè)輸入字精確的確定它的概率，并根據(jù)概率生成一個(gè)和合適的代碼，使輸出的碼流小雨輸入碼流）。1.2實(shí)驗(yàn)方法1.void perform( const Mat &InputArray, Mat &OutputArray,const int reserve,const int flag);method-perform(originalimage,dst,1,DCT); /depend on dct or dft2. 調(diào)用對(duì)象的成員函數(shù) Imagen

3、tropy來(lái)求信息熵，ZEROORDER為零階熵，F(xiàn)IRST_ORDER為一階熵，method-Imagentropy(originalimage,FIRST_ORDER);3. 使用對(duì)象method的成員函數(shù)FidelityCriteria求保真度；二、結(jié)果與討論2.1實(shí)驗(yàn)結(jié)果1. 改變塊變換編碼中所使用變換的種類（DFT與DCT），結(jié)果圖如下：DFT： 圖一（a） （b）DCT： 圖二（a） （b）2. 改變哈爾小波編碼中分級(jí)參數(shù)level與置零系數(shù)閾值res，結(jié)果圖如下：Level=1 res=1 圖三（a） （b） Level=1 res=0.01 圖四（a） （b）Level=2 r

4、es=0.01 圖五（a） （b）Level=3 res=0.01 圖六（a） （b）3.了解：可利用塊變換編碼的程序計(jì)算圖像的0、1、2階信息熵，通過(guò)改變主程序中計(jì)算信息熵的函數(shù)最后一個(gè)輸入?yún)?shù)ZERO_ORDER為FIRST_ORDER或SECOND_ORDER實(shí)現(xiàn)。2.2實(shí)驗(yàn)討論1. 塊變換編碼中，比較使用不同變換運(yùn)行程序后出現(xiàn)的均方根誤差ems與信噪比SNR的結(jié)果值，簡(jiǎn)單解釋出現(xiàn)差別的原因，并對(duì)兩種變換的壓縮效果做出評(píng)價(jià)；由圖一和圖二可以看出，DCT變換結(jié)果優(yōu)于DFT。離散傅里葉變換DFT是傅里葉變換在時(shí)域和頻域上都呈離散的形式。實(shí)際應(yīng)用中通常采用快速傅里葉變換計(jì)算DFT。離散余弦變換

5、DCT是與傅里葉變換相關(guān)的一種變換，DCT相當(dāng)于一個(gè)長(zhǎng)度大概是它兩倍的DFT，但DCT是對(duì)一個(gè)實(shí)偶函數(shù)進(jìn)行變換的。從形式上來(lái)看，DCT是線性可逆變換。當(dāng)圖像壓縮比增大時(shí)，壓縮效率減小，圖像質(zhì)量降低。但是，當(dāng)壓縮比較高時(shí)，這種塊的獨(dú)立處理帶來(lái)了壓縮效率上的限制和方塊效應(yīng)問(wèn)題，方塊效應(yīng)成為限制DCT編碼質(zhì)量的主要因素。2.試闡釋哈爾小波編碼中分級(jí)參數(shù)level與置零系數(shù)閾值res對(duì)結(jié)果的影響；根據(jù)圖四、圖五、圖六的對(duì)照試驗(yàn)可以看出：分級(jí)參數(shù)對(duì)于圖像變換后的清晰度有很大影響，level值越大，則圖像越模糊。這和理論分析結(jié)果相互對(duì)應(yīng)；根據(jù)圖三和圖四的對(duì)照試驗(yàn)可以看出：res的值對(duì)圖像的壓縮效果不明顯（

6、由于實(shí)驗(yàn)代碼問(wèn)題），理論分析出置零系數(shù)閾值是對(duì)圖像壓縮影響最大的系數(shù)。小于res的值置零，大于res的值不變，因此res值越大，壓縮效果越明顯。附錄（實(shí)驗(yàn)代碼）.pro文件程序：#-# Project created by QtCreator 2014-04-04T14:35:39#-QT += coreQT -= guiTARGET = 2dfourieCONFIG += consoleCONFIG -= app_bundleTEMPLATE = appSOURCES += main.cppINCLUDEPATH+=d:Qtopencv2.2includeopencvd:Qtopencv2.

7、2includeopencv2d:Qtopencv2.2includeLIBS+=d:Qtopencv2.2liblibopencv_calib3d220.dll.ad:Qtopencv2.2liblibopencv_contrib220.dll.ad:Qtopencv2.2liblibopencv_core220.dll.ad:Qtopencv2.2liblibopencv_features2d220.dll.ad:Qtopencv2.2liblibopencv_flann220.dll.ad:Qtopencv2.2liblibopencv_gpu220.dll.ad:Qtopencv2.2

8、liblibopencv_highgui220.dll.ad:Qtopencv2.2liblibopencv_imgproc220.dll.ad:Qtopencv2.2liblibopencv_legacy220.dll.ad:Qtopencv2.2liblibopencv_ml220.dll.ad:Qtopencv2.2liblibopencv_objdetect220.dll.ad:Qtopencv2.2liblibopencv_video220.dll.aWavelet文件程序如下：頭文件code.h:#ifndef CODE_H#define CODE_H#include #inclu

9、de #include #include using namespace cv;Mat WDT( const Mat &_src, const String _wname, const int _level );void wavelet( const String _wname, Mat &_lowFilter, Mat &_highFilter );Mat IWDT( const Mat &_src, const String _wname, const int _level );Mat waveletDecompose( const Mat &_src, const Mat &_lowFi

10、lter, const Mat &_highFilter );Mat waveletReconstruct( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter );#endif/ CODE_HCode.cpp:#includecode.h/ 小波變換Mat WDT( const Mat &_src, const String _wname, const int _level ) /int reValue = THID_ERR_NONE; Mat src = Mat_(_src); Mat dst = Mat:zeros

11、( src.rows, src.cols, src.type() ); int N = src.rows; int D = src.cols; / 高通低通濾波器 Mat lowFilter; Mat highFilter; wavelet( _wname, lowFilter, highFilter ); / 小波變換 int t=1; int row = N; int col = D; while( t=_level ) /先進(jìn)行行小波變換 for( int i=0; irow; i+ ) / 取出src中要處理的數(shù)據(jù)的一行 Mat oneRow = Mat:zeros( 1,col, s

12、rc.type() ); for ( int j=0; jcol; j+ ) oneRow.at(0,j) = src.at(i,j); oneRow = waveletDecompose( oneRow, lowFilter,highFilter ); / 將src這一行置為oneRow中的數(shù)據(jù) for ( int j=0; jcol; j+ ) dst.at(i,j) = oneRow.at(0,j); / 小波列變換 for ( intj=0; jcol; j+) / 取出src數(shù)據(jù)的一列輸入 Mat oneCol = Mat:zeros( row, 1,src.type() ); fo

13、r ( int i=0; irow; i+ ) oneCol.at(i,0)= dst.at(i,j); oneCol=( waveletDecompose( oneCol.t(),lowFilter,highFilter ).t(); for ( inti=0; irow; i+ ) dst.at(i,j) = oneCol.at(i,0); / 更新 row /= 2; col /=2; t+; src=dst; normalize(dst,dst,0,1,CV_MINMAX); return dst;/ 小波逆變換Mat IWDT( const Mat &_src, const Stri

14、ng _wname, const int _level ) /int reValue = THID_ERR_NONE; Mat src= Mat_(_src); Mat dst=Mat:zeros(src.rows,src.cols,src.type() ); int N =src.rows; int D = src.cols; /高通低通濾波器 Mat lowFilter; Mat highFilter; wavelet( _wname,lowFilter, highFilter ); /小波變換 intt=1; int row =N/std:pow( 2., _level-1); int

15、col= D/std:pow(2., _level-1); while ( row=N&col=D ) /小波列逆變換 for ( int j=0;jcol;j+) /取出src數(shù)據(jù)的一行輸入 Mat oneCol =Mat:zeros( row,1,src.type() ); for ( int i=0; irow;i+) oneCol.at(i,0)=src.at(i,j); oneCol = ( waveletReconstruct( oneCol.t(),lowFilter,highFilter ) ).t(); for(int i=0;irow;i+) dst.at(i,j) = o

16、neCol.at(i,0); /行小波逆變換 for( int i=0;irow; i+ ) /取出src中要處理的數(shù)據(jù)的一行 Mat oneRow = Mat:zeros( 1,col,src.type(); for ( intj=0; jcol; j+ ) oneRow.at(0,j) =dst.at(i,j); oneRow=waveletReconstruct(oneRow, lowFilter, highFilter ); /將src這一行置為oneRow中的數(shù)據(jù) for( int j=0;jcol; j+ ) dst.at(i,j) =oneRow.at(0,j); row*=2;

17、 col *= 2; src=dst; return dst;/*.*/ 調(diào)用函數(shù)/ 生成不同類型的小波，現(xiàn)在只有haar，sym2void wavelet(constString _wname, Mat &_lowFilter, Mat &_highFilter) if( _wname=haar|_wname=db1 ) intN=2; _lowFilter = Mat:zeros(1, N,CV_32F ); _highFilter =Mat:zeros(1, N, CV_32F ); _lowFilter.at(0,0)=1/sqrtf(N); _lowFilter.at(0,1)=1/

18、sqrtf(N); _highFilter.at(0,0)=-1/sqrtf(N); _highFilter.at(0, 1)= 1/sqrtf(N); if (_wname =sym2 ) int N=4; floath =-0.483,0.836, -0.224, -0.129; floatl=-0.129,0.224, 0.837,0.483; _lowFilter =Mat:zeros( 1,N,CV_32F ); _highFilter=Mat:zeros(1, N,CV_32F ); for (int i=0;iN;i+) _lowFilter.at(0,i)= li; _high

19、Filter.at(0,i)=hi; /小波分解Mat waveletDecompose(const Mat &_src, const Mat &_lowFilter, constMat &_highFilter ) assert(_src.rows=1 &_lowFilter.rows=1&_highFilter.rows=1 ); assert(_src.cols=_lowFilter.cols &_src.cols=_highFilter.cols ); Mat src = Mat_(_src); intD= src.cols; Mat lowFilter=Mat_(_lowFilter

20、); Mat highFilter=Mat_(_highFilter); / 頻域?yàn)V波，或時(shí)域卷積；ifft(fft(x) *fft(filter)=cov(x,filter) Mat dst1=Mat:zeros( 1, D, src.type(); Mat dst2 =Mat:zeros( 1,D, src.type() ); filter2D(src,dst1,-1,lowFilter ); filter2D(src,dst2,-1, highFilter); / 下采樣 Mat downDst1 =Mat:zeros(1,D/2,src.type(); Mat downDst2 =Ma

21、t:zeros(1,D/2, src.type() ); resize(dst1, downDst1,downDst1.size(); resize(dst2, downDst2, downDst2.size() ); / 數(shù)據(jù)拼接 for ( inti=0;iD/2; i+) src.at(0, i)=downDst1.at(0, i); src.at(0, i+D/2)=downDst2.at(0, i); return src;/小波重建Mat waveletReconstruct(const Mat &_src,constMat &_lowFilter, constMat &_high

22、Filter) assert( _src.rows=1 &_lowFilter.rows=1& _highFilter.rows=1 ); assert(_src.cols=_lowFilter.cols & _src.cols=_highFilter.cols ); Mat src=Mat_(_src); int D =src.cols; Mat lowFilter = Mat_(_lowFilter); Mat highFilter =Mat_(_highFilter); /插值; Mat Up1 = Mat:zeros(1,D,src.type() ); Mat Up2=Mat:zero

23、s(1, D,src.type();/ /插值為0/ for( inti=0,cnt=1;iD/2; i+,cnt+=2)/ / Up1.at(0, cnt)=src.at(0,i); /前一半/ Up2.at(0, cnt) =src.at(0, i+D/2 );/ 后一半/ /線性插值 Mat roi1( src, Rect(0,0,D/2,1); Mat roi2( src, Rect(D/2,0,D/2, 1); resize(roi1,Up1, Up1.size(),0,0, INTER_CUBIC ); resize(roi2,Up2,Up2.size(), 0,0,INTER_C

24、UBIC ); / 前一半低通，后一半高通 Mat dst1=Mat:zeros(1, D,src.type() ); Mat dst2= Mat:zeros(1,D, src.type() ); filter2D(Up1,dst1,-1, lowFilter); filter2D( Up2, dst2, -1, highFilter ); /結(jié)果相加 dst1 =dst1+ dst2; normalize(dst1,dst1,0,1,CV_MINMAX); returndst1;Main.cpp:#includecode.hintmain(int argc, char *argv) int

25、level=1; Mat src=imread(d:/lena.tif,CV_LOAD_IMAGE_GRAYSCALE); Mat dst=WDT(src,haar,level); float res=0.01; Mat te = Mat_(dst); for (int i=round(dst.rows/(pow(2,level);idst.rows;i+) for (int j=round(dst.cols/(pow(2,level);jdst.cols;j+) te.at(i,j)=(abs(te.at(i,j)res)*te.at(i,j); Mat th=IWDT(te,haar,le

26、vel); imshow(dst,dst); imshow(th,th); waitKey(0); return 0;Transform-coding文件程序如下：頭文件transformcoding.h#ifndefTRANFORM_CODING_H#define TRANFORM_CODING_H#include #include #include #include using namespace cv;using namespace std;enum DCT=0,DFT=1,ZERO_ORDER=2,FIRST_ORDER=3,SECOND_ORDER=4;class Tranformc

27、odingpublic: double ems,SNR,H; /構(gòu)造函數(shù)，初始化成員變量 Tranformcoding() this-ems=0; this-SNR=0; this-H=0; this-index=0; Tranformcoding() delete this; void FidelityCriteria(const Mat &compressed, const Mat &decompressed); void perform( const Mat &InputArray, Mat &OutputArray,const int reserve,const int flag);

28、void Imagentropy(const Mat &InputArray,const int flag);private: void Cosinecoding(Mat &InputArray,Mat &Z,const int reserve); void Fouriercoding( Mat &InputArray, Mat &Z,const int reserve); int index; Point maxidex;#endif/TRANFORM_CODING_H源文件main.cpp:#include tranformcoding.hint main(int argc, char *

29、argv) Mat originalimage = imread(/home/timlentse/Qtproject/lena.tif,CV_LOAD_IMAGE_GRAYSCALE); / Tranformcoding是我定義的一個(gè)類，使用這個(gè)類來(lái)實(shí)例化一個(gè)對(duì)象,名字叫method Tranformcoding *method=new Tranformcoding; Mat dst; /前面我們已經(jīng)有了method這個(gè)對(duì)象了，當(dāng)然可以使用它來(lái)調(diào)用它本身的成員函數(shù)，來(lái)執(zhí)行我們預(yù)先定義好的行為 /例如method-Imagentropy 執(zhí)行信息熵的計(jì)算 method-Imagentropy(o

30、riginalimage,ZERO_ORDER); / i want to output the second_order entroy cout-the entropy =Hperform 執(zhí)行塊轉(zhuǎn)換編碼 method-perform(originalimage,dst,8,DFT); /depend on dct or dft imshow(the original image,originalimage); imshow (the effect, dst); waitKey(0); /method-FidelityCriteria 執(zhí)行保真度計(jì)算 method-FidelityCrite

31、ria(originalimage,dst); cout-the ems =emsendl-and the SNR =SNRendl; return 0;Transformcoding.cpp：#include tranformcoding.h/這是計(jì)算圖像保真度（包括ems,SNR)的子程序，你需要輸入兩幅圖像，即原圖和壓縮的圖/第一個(gè)參數(shù)是原圖，第二個(gè)參數(shù)是壓縮后的圖void Tranformcoding:FidelityCriteria(const Mat &compressed, const Mat &decompressed) if (compressed.size != decom

32、pressed.size| compressed.type()!= decompressed.type() cout-endl; couterror:nothing done because the two images are not in the same size or type endl It is recommended to check your two imagesendl; cout-endl; return ; SNR = norm(decompressed, NORM_L2); ems = norm(compressed, decompressed, NORM_L2); S

33、NR = pow(SNR/ems), 2); ems = sqrt(ems*ems/(compressed.rows* compressed.cols);/以下是求圖像信息熵的子函數(shù) ,flag標(biāo)志位：/-ZERO ORDER ：求零階熵/-FIRST ORDER ：求一階熵/-SECOND ORDER ：求二階熵void Tranformcoding:Imagentropy(const Mat &InputArray, const int flag) int array766= 0; /像素值統(tǒng)計(jì)部分，放在array數(shù)組里面 int row= InputArray.rows; int col

34、= InputArray.cols; if (InputArray.isContinuous() col= col*row; row= 1; for (int i(0); irow; +i) for(intj(0); jcol; +j) const uchar*p =InputArray.ptr(i); if (flag=2) index=pj; else if(flag= 3) index=pj+(j!=col-1)*pj+1; else if(flag= 4) index=(j!=0)*pj-1+pj+(j!=col-1)*pj+1; arrayindex+; /計(jì)算信息熵 double

35、pr(0); for(intk(0); k766; k+) pr= arrayk; if (pr!=0) pr/=(InputArray.cols* InputArray.rows); H=- (log(pr)/log(2)*pr + H; /定義perform函數(shù)來(lái)響應(yīng)用戶的需求，perform函數(shù)是用來(lái)實(shí)施塊轉(zhuǎn)換編碼的/根據(jù)標(biāo)志位flag和保留系數(shù)個(gè)數(shù)來(lái)執(zhí)行塊轉(zhuǎn)換編碼voidTranformcoding:perform( constMat &InputArray, Mat &OutputArray, const int reserve, const int flag) OutputArra

36、y= Mat_(InputArray); float t88= 16,11,10,16,24,40,51,61, 12,12,14,19,26,58,60,55, 14,13,16,24,40,57,69,56, 14,17,22,29,51,87,80,62, 18,22,37,56,68,109,103,77, 24,35,55,64,81,104,113,92, 49,64,78,87,103,121,120,101, 72,92,95,98,112,100,103,99; Mat Z = Mat(8,8,CV_32F,t); for(intu(0); uOutputArray.rows

37、/8;+u) for (int v(0);vOutputArray.cols/8; +v) Mat temp(OutputArray,Rect(v*8,u*8,8,8); if (flag=1) Fouriercoding(temp,Z, reserve); elseif(flag=0) Cosinecoding(temp, Z, reserve); else std:cout-no matching fuctioncall,maybethe flagiswrong-std:endl; OutputArray.convertTo(OutputArray,CV_8U); return; OutputArray.convertTo(OutputArray, CV_8U);/這是余弦變換編碼子函數(shù)voidTranfor