1、使用MATLAB設計小波變換程序中的若干問題    使用MATLAB設計小波變換程序中的若干問題在使用MATLAB完成小波變換程序和通過閾值來壓縮圖像的過程中，我和許多同學都是邊學邊用，是從一個接一個的問題中逐步理解小波和MATLAB編寫程序的。因此我愿意就個人遇到和解決問題的經驗與大家討論，希望能夠對遇到同樣問題的人有所幫助。在清華大學林福宗老師倡導的網上互動的學習方式中，老師同學的開誠布公的討論，尤其是林老師啟發大家對出現問題采取的態度和做法，對我今后成長為一名合格的清華的研究生意義重大，謹以此文表示對他指導關心的敬意！1. 內容簡介本文分

2、為三部分:如何使用MATLAB設計小波標準與標準分解；如何完成使用小波變換壓縮圖像；仍需探討的問題。每一部分主要以問題和例子的形式講述，為了便于參照，附錄部分給出部分源代碼供大家參考指正。我個人在完成作業的時候，走了許多彎路，最后，反復讀老師的第3章講義的第3-5節，明白了小波變換的一些非常基礎性的知識，正是如此，我主要是按照第三章的例子，做出使用MATLAB進行Haar變換的實例，這至少對完成非標準Haar小波的3級分解與合成沒有任何問題，而且計算速度遠遠比使用一維變換快，甚至超過dwt2和wavedec2。再者，是通過這些簡短的例子的學習，也能從實踐的角度理解小波變換的基礎知識和道理，我覺

3、得掌握小波的知識要比簡單使用dwt2直接分析出結果更重要。為了讓使用一維變換或打算采用卷積完成分解重構程序的同學有所參考，也給出了我以前寫的代碼和設想供有興趣者參考。在圖像壓縮的任務中，我個人建議采用wdencmp，原因是簡單可靠，能夠生成老師要求小波壓縮與重構的演示及PNG文件。在最后部分著重對使用小波壓縮之后重構圖像的PNG文件會變大問題進行了簡單的猜測性的解釋和討論，希望感興趣的同學能深入研究。2. 如何使用MATLAB設計小波標準與標準分解1. 使用Haar和Db9小波編寫圖像文件的標準和非標準分解重構程序實質上，任務書中規定首先要編寫種程序之一：哈爾（haar）小

4、波的標準分解重構程序；Daubechies9小波的標準分解重構程序；哈爾（haar）小波的非標準分解重構程序；Daubechies9小波的非標準分解重構程序。如果采用標準方法，需要生成圖的系列圖像。如果采用非標準方法則需要生成圖的系列圖像。我個人編寫了非標準的分解重構程序，并總結網上同學們公開的方法，認為可以有種方法實現。但前提是必須仔細閱讀林老師推薦的補充教材第三章“小波與小波變換”的3.3、3.4和3.5節（頁）。如果能使用MATLAB簡單地實踐一下教材的例子，至少完成哈爾小波的標準與非標準程序相當簡單。為此我主要介紹如何實現老師教材給我們的例子是如何在MATLAB中實現的。然后，再引申到

5、使用MATLAB的函數實現方法。為了描述問題簡便，使用黑體表示需在MATLAB的命令窗口（commandwindow）的輸入部分。如果不熟悉下面出現的函數功能和使用方法，請再命令窗口中使用help函數名，如helpdwt2，能夠得到英文的功能和用法說明。2. 哈爾小波變換的MATLAB實例先介紹MATLAB的矩陣(Matrix)和向量(Vector)的賦值方法：在commandwindow的>>符號下輸入：一維向量：=9735=2,5,8,9,7,4,-1,1二維矩陣：A=64236160675795554121351501617474620214342244026273

6、736303133323435292838392541232244451918484915145253111056858595462631和現在就是例3.1例3.2一維向量,A就是3.5.1中的圖像矩陣。在進行Haar小波的一維與二維變換（分解）之前，首先介紹一下MATLAB的矩陣運算。MATLAB是MatrixLaboratory的合寫，它對矩陣運算之功能堪稱一流。由于使用矩陣描述問題更象數學表達式，所以編寫的程序不僅高效，更易讀。所以MATLAB程序應該盡量使用矩陣直接描述。如M是一4X4的矩陣，則B=I*M就完成了使用兩個for循環編寫乘法的程序。分析老師的例3.1，那么可以得到如下結果

7、：（9+7）*1/2（3+5）*1/2（9-7）*1/2（3-5）*1/2如果把其看作矩陣方式的乘法，那么令M為其表述求取平均和差值的系數矩陣，則輸入：M=1/201/201/20-1/2001/201/201/20-1/2把上述的M輸入到MATLAB中,然后使用：C=I*M得到的結果就是841-1。這就是非規范化Haar小波的第一級分解系數。注意C的前兩項84就是近似系數（ApproximationCoefficients）,后兩項1_1就是細節系數（DetailCoefficients）.M就是Haar小波的非規范化（non-normalization）系數矩陣。如果我們執行：M*M_(M

8、_代表的是M的轉置矩陣)，你會發現得到對角線為0.5 為0的4X4矩陣，所以如果讓它變成單位矩陣（對角線為1），必須把原來的1/2增大，變為1/sqrt(2),所以執行：N=M*sqrt(2)后得到新的系數矩陣0.707100.707100.70710-0.7071000.707100.707100.70710-0.7071再執行N*N_就會得到單位矩陣。N就是Haar小波的規范化（normalization）系數矩陣。現在執行:C*M_將得到4.53.51.52.5，顯然是因為M*M_的對角矩陣的值為0.5。現在使用Cn=I*N得到11.31375.65691.4142-1.4142規范化的

9、一維Haar小波一級分解系數。I1=Cn*N_得到9377規范化的一維Haar小波逆變換（重構）結果。可見得到這樣規范化矩陣的好處是為了逆變換時不需要額外的運算，只是需要分解系數矩陣乘以規范化系數矩陣的轉置。如果需要對分解系數矩陣進行第2級分解，那么此時只對C或Cn的前2項的低頻系數進行，此時的Haar的非規范化與規范化系數矩陣M1和N1分別為：M1=1/21/21/2�1/2N1=1/21/21/2�1/2那么第2級分解可以用以下方式進行：C1=C(1:2)*M1Cn1=Cn(1:2)*N1結果是C1=62,Cn1=12.04.0。注意,C1(1:2)表示取出C中

10、第到第個向量，符號“：”的含義是從到的意思，這是MATLAB截取數據的方法。逆變換：C1*M1*2C(3:4)*M*2Cn1*N1Cn(3:4)*N這里采用了中括號的賦值方式直接將第級的重構合并到第級的系數中再重構原始向量I。從上面的實驗可以看出，除了輸入哈爾小波系數矩陣M和N比較麻煩，整個變換和重構過程異常簡單,這就啟示我們只要編寫一個Haar系數生成矩陣的函數，整個的分解問題即將化簡。1. 問題1：如何生成Haar小波系數矩陣？MATLAB的函數可以如以下形勢構造：在MATLAB中使用File?New?M-file,輸入functionM=haarmatrix(rows,cols

11、,flag)%矩陣的行數rows,列數cols,非規范化flag=否則就是規范化矩陣%M為函數返回的Haar系數矩陣M=zeros(rows,cols);%首先制作rowsXcols的矩陣ifflag=0sv=0.5;%non-normalizedelsesv=1.0/sqrt(2.0);%normalizedendhalf=rows/(2);矩陣從中間分開，左半部為近似系數部分，右半部為細節系數部分fori=1:2:half*2M(i,ceil(i/2)=sv;M(i+1,ceil(i/2)=sv;M(i,ceil(i/2)+half)=sv;M(i+1,ceil(i/2)+half)=-s

12、v;endM=sparse(M)%生成稀疏矩陣，提高運算速度注意分號的作用是不在command窗口上顯示M的值。在輸入后，請保存與函數名同名的文件名haarmatrix.m。現在來做例3.2：=2,5,8,9,7,4,-1,1N=haarmatrix(8,8,1)C1=F*NN1=haarmatrix(4,4,1)C2=C1(1:4)*N1N2=haarmatrix(2,2,1)C3=C2(1:2)*N2按圖3-21小波分解的層次結構合成最終系數C=C3,C2(3:4),C1(5:8)得到的C就是老師讓我們使用第三章17頁偽代碼編寫的一維小波變換分解的結果。其逆變換（重構）代碼如下：Fr=C3

13、*N2",C2(3:4)*N1",C1(5:8)*N"注意理解最內部的中擴號實質上是使用C3逆變換的結果（低頻）與C2的后位細節系數一起逆變換得倒C1的前位的低頻近似系數，然后與C1的后位細節系數構造出原始向量Fr。為了驗證我們的一維Haar變換程序，我們使用MATLAB提供的一維離散小波變換函數dwt來進行同樣的計算：第一級分解ca1,cd1=dwt(F,"haar")第級分解ca2,cd2=dwt(ca1,"haar")第級分解ca3,cd3=dwt(ca2,"haar")按圖3-21小波分解的層次結

14、構合成最終系數c=ca3cd3cd2cd1或直接使用MATLAB提供的維多級分解程序c=wavedec(F,3,"haar")比較c與C(注意MATLAB中大小寫變量的區別)，可知我們成功構造了Haar小波分解合成程序。也許這是你可能覺得剛才的系數矩陣比起MATLAB提供的函數來說還是復雜了，但如果把上述過程函數化，做一個matrixDWT和matrixIDWT函數，使用一樣簡單，而且關鍵在于，我們使用矩陣方式解決維圖像的小波非標準變換的演示過程來說會帶來更大便利。2. 問題2如何進行圖像的1級非標準分解和重構？如果我們利用2.2節開始時已經輸入的18頁的圖像矩陣

15、A，那么它的第一級二維變換的非標準分解過程只需步：(為了方便核對MATLAB的dwt2分解結果，使用規范化系數矩陣)M1=haarmatrix(8,8,1);AR1=A*M1;%這就是進逐行分解的圖像AC1=M1"*AR1;%這就是進逐列分解的圖像，也是一次維變換的系數矩陣使用MATLAB的dwt2來驗證cachcvcd=dwt2(A,"haar")c=cacv;chcd%將個系數合并在一起，注意cv與ch的位置與第四章圖4-04區別。對比ca3,ch3,cv3,cd3與AC3，可知矩陣算法正確。細心的同學也許能發現我們的矩陣位置和第章圖-(a)本文圖1解釋的不同

16、，這是因為dwt2采用的算法與我們不同，但這并不影響小波分解系數的使用。詳細解釋見附錄1。圖1、多級分解的各個系數矩陣其重構過程只需要兩步：（參照第三章24頁的公式，但是我們只采用單級重構）AR1=M1*AC1;A=AR1*M1"3. 問題3如何進行圖像的1級非標準分解和重構？如果是進行級非標準分解，那么只需要對前一次系數矩陣的的左上角的近似系數圖-04（a）中的a進行，如對是級分解：M1=haarmatrix(8,8,1);AR1=A*M1;%這就是第1次逐行分解的圖像AC1=M1"*AR1;%這就是第次逐列分解的圖像，M2=haarmatrix(4,4,);A

17、R2=AC1(1:4,1:4)*M2;%第級行分解，只對左上角低頻近似系數部分AC2=M2"*AR2;%第級列分解M3=haarmatrix(2,2,);AR3=AC2(1:2,1:2)*M3;%第3級行分解，只對左上角低頻近似系數部分AC3=M3"*AR3;%第級列分解執行MATLAB提供的維多級小波分解函數wavedec2C,L=wavedec2(A,3,"haar");ca3=appcoef2(C,L,"haar",);ch3,cv3,cd3=detcoef2("all",C,L,3)AC3相當于上圖中(b)

18、左上角的最小的a3v3h3d3，AC2是c2v2h2d2，AC1是c1c2c3c4。那么如果合成圖4-04(c)這樣的圖像只需要把AC3替代AC2的a2,在把AC2替代AC1的a1就可以，相應的代碼：AC2(1:2,1:2)=AC3;AC1(1:4,1:4)=AC2;如何由3級分解系數矩陣重構出原始矩陣？參照1級非標準重構,注意此時從第3級向第1級重構，即首先由A3中的ahvd四個系數重構出A2的低頻部分a,然后重復這種方法：AR3=M3*AC3;A3=AR3*M3"；AC2(1:2,1:2)=A3；%重構的A3就是A2的低頻近似系數，如圖4-04(b)的ahvdAR2=M2*AC2

19、;A2=AR2*M2"；AC1(1:4,1:4)=A2；%重構的A2就是A1的低頻近似系數，如圖4-04(a)的aAR1=M1*AC1;A1=AR1*M1"；%A1就是A注意：上述分解與重構過程與老師第3章中的分解重構過程在一級分解的情況下相同，在第2級分解時就發生了變化。上述分解過程使用了與MATLAB的wavedec2近似的多分辨率分析(multi-resolutionanalysis),產生的矩陣系數除了h與v的位置對調外，二者結果完全相同。然而采用第3章21-24頁介紹的方法（即使是使用23頁的規范化系數矩陣）不能生成與wavedec2相同的矩陣，正是因為如此，采用

20、設置閾值來壓縮圖像時會與使用wavedec2的結果不同。原因分析及對應的24頁的算法的演示程序請參見附錄2。（請林老師務必核對）。3. 三種實現非標準分解的方法：第一種：2.2.3介紹的矩陣A變成256X256的照片矩陣，1:2的矩陣下標變成1:64,1:4變成1:128, 不變，保存系數AR1,AC1,AR2,AC2,AR3,AC3為PNG格式，就是圖3-26所需要的照片，只不過下1級沒有顯示上一級的細節系數（hvd）部分。使用合成技術類似這樣的語句c=cacv;chcd應該很方便得出。當然采用附錄2中的程序可以更直接，更快速地實現老師的范例。第二種：按照第3章偽代碼形式使用MAT

21、LAB一維離散小波變換函數dwt編寫。dwt的用法:cacd=dwt(I,wname);I是你要分解的一維向量，wname是字符串，使用單引號擴起來，如db9或haar。返回系數ca是低頻近似系數，cd是高頻細節系數。由于db9的濾波系數長度為18，dwt需要擴展以減少邊界誤差，重構函數idwt能夠保證恢復到原來的尺寸。所以可以采用不同的擴展邊界模式來使用，缺省值是sym,它在分解后會使ca和cd的長度變為(LI+LF-1)/2,LI是向量I的長度，LF是小波濾波系數的長度。也可以改變邊界擴展模式，如使用per,具體方法是cacd=dwt(I,wname,mode,per).它的ca和cd的長

22、度為LI/2,但是它可能帶來邊界誤差。我編寫了2維非標準分解nsdwt2和重構nsidwt2函數，可參見附錄3。第三種：直接使用MATLAB卷積函數conv2，即仿照dwt2的函數編寫。此處不再詳細討論，有興趣者參見附錄1這三種方式在MATLAB下計算速度最快的是第一種，最慢的是第二種，第三種速度與第一種速度比相差不大。在筆者PIII733上運行分解wbarb圖像過程所用時間（使用TIC,TOC）分別為0.10,0.28,7.23秒。4. 常見問題：如何讀取、轉換和保存圖像文件？X,map=imread(文件名);%文件名必須帶后綴。使用size(X)來判斷圖像矩陣的維數，如果是25

23、6256就是2維的，那么map就是它的色板文件，可以使用imshow(X,map)來顯示它。如果2562563就是3維真彩色圖像（TruecolorImage），使用imshow(X)來顯示。一般我們把2維的彩色圖片稱為偽彩色圖像（Indexedcolorimage）。一般來說，不管圖像文件是不是PNG格式，最好首先使用MATLAB的imwrite方式重新保存：對于真彩色文件，imwrite（X,新文件名）;對于偽彩色文件，imwrite（X,map,新文件名）;由imread讀出的圖像文件采用uint8的類型，必須轉換成double類型之后方可進行運算。轉換方法：X=doubl(X);經過小

24、波變換和重構的X是double類型，在保存圖像前需要轉換成uint8類型，如果處理不好，會產生截斷誤差，許多同學反映閾值為0時的文件大小和原始文件大小不一致可能就是出現在這里，轉換方法：Y=uint8(round(X);這里的round是就近取整，這樣7.99999999就會變成8而不會被截斷變為7。如何處理Indexedcolor文件？原則說，只要可以轉換成double類型我們就可以開始小波變換。但是由于偽彩色圖片的色板map保存了每個象素值對應的rgb的顏色值，因此，如果map中的顏色值不連續，就不能保證小波轉換的效果。可以使用imshow(X,map);colorbar;來顯示圖片的色板

25、是否連續。如果不連續，可以使用MATLAB提供的方法轉換成灰度圖像處理。這是林老師在“教師答疑”中貼出的MATLAB幫助文件ImageConversion，我把它轉換成一個小程序，可以修改使用。Anna指出ind2gray和ind2rgb也能完成同樣的任務。X,map=imread("你的偽彩色圖片文件");X=double(X);ifmin(min(X)=0X=X+1;endimage(X)title("OriginalColorIndexedImage")colormap(map);colorbarR=map(X,1);R=reshape(R,siz

26、e(X);G=map(X,2);G=reshape(G,size(X);B=map(X,3);B=reshape(B,size(X);Xrgb=0.2990*R+0.5870*G+0.1140*B;n=255;%NumberofshadesinnewindexedimageX=round(Xrgb*(n-1)+1;map2=gray(n);image(X),title("ProcessedGrayScaleIndexedImage")colormap(map2),colorbarbaboon=X;map=map2;imwrite(X,map,"gray.png&q

27、uot;);%保存為灰度圖像，文件名為gray.png如何處理真彩色文件？有幾種方法，在MATLAB的ImageConversion中有敘述，但許多同學采用了單獨處理它的3個顏色的2維矩陣。設X為三維矩陣(256,256,3),X(:,:,1)代表紅顏色的2維矩陣X(:,:,2)代表綠2維矩陣,X(:,:,3)代表綠2維矩陣。X,map=imread(真彩色文件);r=double(X(:,:,1);%r是256x256的紅色信息矩陣g=double(X(:,:,2);%g是256x256的綠色信息矩陣b=double(X(:,:,3);%b是256x256的蘭色信息矩陣%開始對rgb分別作小

28、波變換，并分別形成各自的小波系數矩陣%把3個變換的系數矩陣合并成圖像文件Y(:,:,1)=uint8(round(r);Y(:,:,2)=uint8(round(g);Y(:,:,1)=uint8(round(b);可以參照3.1的simplecmp.m程序來使用。3. 如何完成使用小波變換壓縮圖像的任務由于壓縮的目的在于比較使用不同小波壓縮后重構圖像的失真度視覺效果和使用PNG文件保存時的文件大小，如果我們自己編寫的小波變換程序存在小的遺漏，可能對壓縮結果判斷錯誤，所以我認為至少應當使用MATLAB為我們提供的壓縮函數記錄結果，以便與自己設計的閾值處理程序進行比較。一般來說，兩種處

29、理結果應該差別很小，甚至無差別。1. 直接使用wdencmp:這是我使用wdencmp編寫的簡單壓縮處理程序simplecmp，可以直接再運行時輸出圖像、0的個數、PNG文件的大小以及計算時間。如果你的圖片為a.png,不管是真彩色、偽彩色，使用：simplecmp(a.png,haar,3)%進行haar小波的三級分解壓縮合成simplecmp(a.png,db9,3)%進行db9小波的三級分解壓縮合成下面是具體程序清單。我只對高頻細節系數部分進行硬閾值設置。functionsimplecmp(fname,wname,level)rgb,map=imread(fname);fig=

30、figure;colorbar;axison;axisequal;set(fig,"position",1020790580,"name","壓縮程序演示")iflength(size(rgb)=3rgbcmp(rgb,wname,level);elseindcmp(rgb,map,wname,level);end%-functionrgbcmp(rgb,wname,level)%這是壓縮真彩色圖像rgb=double(rgb);THR=051020;otxt=sprintf("%s_zeros.txt",wnam

31、e);fid=fopen(otxt,"w");fprintf(fid,"%20s%15s%15s%15s","文件名","大小","閾值","0數");forT=THRTIC;rgb(:,:,1),cxc,lxc,perf0,perfl2=wdencmp("gbl",rgb(:,:,1),wname,level,T,"h",1);num0=length(find(abs(cxc)<0.0000001);rgb(:,:,2),cxc

32、,lxc,perf0,perfl2=wdencmp("gbl",rgb(:,:,2),wname,level,T,"h",1);num0=num0+length(find(abs(cxc)<0.0000001);rgb(:,:,3),cxc,lxc,perf0,perfl2=wdencmp("gbl",rgb(:,:,3),wname,level,T,"h",1);num0=num0+length(find(abs(cxc)<0.0000001);x=uint8(round(rgb);oname=spr

33、intf("%s_%d.png",wname,T);imwrite(x,oname);tmp=imfinfo(oname);fs=tmp.FileSize;fprintf(fid,"%20s%15d%15d%15d",oname,fs,T,num0);e_t=TOC;sTitle=sprintf("%s閾值%d的文件大小%d,0數%u,用時%f,任意鍵繼續.",wname,T,fs,num0,e_t);image(x)title(sTitle);pauseendfclose(fid);%-functionindcmp(x,map,wn

34、ame,level)%這是壓縮偽彩色圖像THR=051020;x=double(x);otxt=sprintf("%s_zeros.txt",wname);fid=fopen(otxt,"w");fprintf(fid,"%20s%15s%15s%15s","文件名","大小","閾值","0數");forT=THRTIC;y,cxc,lxc,perf0,perfl2=wdencmp("gbl",x,wname,level,T,&quo

35、t;h",1);num0=length(find(abs(cxc)<0.0000001);y=uint8(round(y);oname=sprintf("%s_%d.png",wname,T);imwrite(y,map,oname);tmp=imfinfo(oname);fs=tmp.FileSize;fprintf(fid,"%20s%15d%15d%15d",oname,fs,T,num0);e_t=TOC;sTitle=sprintf("%s閾值%d的文件大小%d,0數%u,用時%f,任意鍵繼續.",wname

36、,T,fs,num0,e_t);imshow(y,map)title(sTitle);pauseendfclose(fid);2. 自己編寫閾值設置函數：可以使用wthresh(X,h,T)，X是要進行閾值設置的矩陣，h表示使用硬閾值處理方式，閾值T=5，如對X=6,-652-2進行閾值設置的結果是66000，如果設置軟s，結果是11000,大于T的，也被縮小。注意:低頻系數是對圖像重構質量最重要的系數，一般不需要設置。4. 仍需探討的問題：為什么使用PNG存儲經小波變換后的重構圖像變大？我曾在清華大學的多媒體課程的教師答疑中寫了“老師：尊重事實:DB9閾值10的PNG文件

37、就是比原文件大”和“續一：尊重事實:DB9閾值10的PNG文件就是比原文件大”，在林老師的鼓勵和指導下，我進行了繼續試驗、分析，與劉趙璧（Anna）同學進行了探討，并得到了Lily（姓名還不知道）同學的幫助，同時同學們也做了各自不同的實驗，現在的實驗結果可以說基本上比較明確，那就是有些圖像就是會變大，這與圖像的種類、紋理等密切相關。林老師曾經鼓勵我去研究一下PNG的壓縮方法，無奈我資質不夠，至今在這方面的進展不大。由于臨近期末考試，作業也要抓緊，所以我暫且將沒有搞明白的內容擱置，待寒假期間再進行，希望對這些問題有各種看法也有興趣研究的同學對此發表意見。以下是我最近試驗、分析和閱讀到的一些相關信

38、息。1. 試驗結果我首先根據老師第三章的Haar矩陣算法推演出DB9的系數矩陣，并實現了分解重構及閾值處理程序，對幾種照片進行了比較，然后使用3.1節的simplecmp進行了相同照片的實驗，結果相當一致。細小差別是因為我的程序對邊界的擴展與MATLAB不一樣，在設置閾值后引起了邊界上小部分不一致造成的。表一：真彩色圖像百合花的處理結果閾值 PngHaar(Mat/Mine) 0數Haat(Mat/Mine) PNGDb9(MAT/Mine) 0數Db9(MAT/Mine) 95973/95973   0

39、 95973/95973 27524/24268 95973/95973 27/95 74552/74292 135838/136063 101882/101992 167412/16566210 51976/51504 163423/163741 98411/98861 199200/19573020 32474/32346 180167/180267 92295/93660 220629/217214從對比表中我們能夠看到2個程序的

40、結果相當一致，因此，我不再給出兩種程序的對比，而是使用simplecmp直接處理的結果說明。將百合花圖像使用I,map=rgb2ind(x,255);轉換成為彩色圖像處理，在將偽彩色圖像轉換為連續變換的灰度圖像（如2.4常見問題中討論的方法）進行處理：表二：百合花的偽彩色圖像和處理后的灰度（gray）圖像的處理結果閾值 PngHaar(Index/Gray) 0數Haar(Index/Gray) PNDb9(Index/Gray) 0數Db9(Index/Gray) 48535/43235   0 485

41、35/43235 6096/7430 48535/43235 18/225 53207/36450 9473/43626 60362/49927 7009/5285210 58025/23602 13362/54344 64916/47813 13202/6588120 60193/14347 21948/60039 66020/46014 24468/73494其他偽彩色與進行加工的灰度圖的結果與此完全一致，這也就說明了如果偽彩色文件的色板不是

42、單調性遞增就不適合小波分解。“Thecolorbartotherightoftheimageisnotsmoothanddoesnotmonotonicallyprogressfromdarktolight.Thistypeofindexedimageisnotsuitablefordirectwaveletdecompositionwiththetoolboxandneedstobepreprocessed.”。我對Facets進行同樣的實驗，結果與此一致。這種處理的結果可以從圖像象素值的連續性來理解。這是處理與不處理的圖像的中間一行的數據圖。另外，不連續的圖像質量在壓縮后會被極大地破壞圖2

43、偽彩色文件變化前后的第128行數據的連續性情況對比2. 分析多種試驗圖片基本能夠反映類似的結果，雖然IndexedColorimage有時令Haar小波的分解重構圖像出現增大現象，單經過處理之后，這種現象就會消失。然而對于DB9可以看到無論真彩色還是處理后的灰度圖像都在閾值510處超過原始圖像的大小，能不能因此得出DB9不適合進行圖像壓縮的結論呢？有一些同學確實這樣認為，但我認為這種觀點因為忽略了如何利用小波進行壓縮和還原的過程，這也正是第四章老師為我們講述的那些編碼算法而造成的。在推薦材料1中也有類似的說明。圖3、JPEG2000的基本結構看一下上圖就可以明白為什么PNG不能衡量小

44、波壓縮的效率問題。上圖的圖像原始數據首先經過正變換（ForwardTransform）就是小波變換的得到小波系數，變換的小波系數經過閾值處理后進行量化，編碼后得到壓縮的圖像文件JPEG2000，如果你沒有JPEG2000的顯示程序，那么你就不能看到它。它的顯示程序就是由解碼器從壓縮數據中解出編碼，進行反量化，得到小波系數，再實施逆變換（InverseTransform）就是小波系數重構。最終得到圖像的原始數據。因此衡量小波變換的效率是應該看你選擇的小波能不能分解出適合“編碼器”壓縮的小波系數，這種編碼器不是PNG的LZ77，因為LZ77壓縮小波分解系數的效率不是最好的。這種高效編碼器在第四章可

45、以找到。那么我們存儲PNG文件的目的是什么呢？我認為壓縮與去噪(de-noising)是同一種方法的兩種提法。他們都使用了設置閾值的方法。我們可以仔細分析經過重構的PNG圖片的質量來體會這種消除噪音的效果，也可以評定小波壓縮后的圖片的視覺質量，同時PNG的文件大小也可以讓我們從LZ77算法的本質來理解小波變換壓縮后的重構圖像的內容變化情況。比如，我們可以從表2中的灰度圖像在haar變換取閾值20時出現塊狀象素，文件大小變為14347，而db9卻為46014，超過原始的PNG大小，但并不出現塊狀而是具有波狀的特征。這本身說明了采用Harr小波壓縮或去噪后重構的圖像中相同的串增多，便于PNG方式壓

46、縮，而db9則在相同閾值的情況下不會象Haar那樣制造馬賽克，說明了它的平滑性，這也能幫助我們理解小波的特性。當然，當閾值繼續增加后，超過某一界限，即使DB9也仍然會使PNG文件大小減小。這本身也就是由雙尺度（Dyadic）小波變換的兩種濾波器決定的。低通濾波結果相當于平均值，高通濾波結果相當于差值，差值能夠保證重構圖像的細節部分丟失最小，如果差值部分被閾值略去的過多，細節就會越來越少，平均意義的值約來越多，直到多到某一個臨界值時（該圖像的閾值取到40），重構的圖像也可能出現較多的相同數字串，這就會提高PNG的壓縮結果。下圖是我對Haar（藍色）小波取閾值為20，db9（紅色）小波閾值取40時

47、第128行1：32列的數據曲線與原始數據（黑色）曲線的對比。可以看出也db9在閾值=40時出現了較多的平均值，但比haar在閾值=20時的曲線要少的多。圖4、haar（藍色）和db9（紅色）壓縮后重構圖像的第128行，1：32列的數據曲線不過，MATLAB給我們提供了量化的方法來決定如何選取閾值。在HELPWaveletToolbox:AdvanceConcepts:ChoosingtheOptimalDecomposition中提到了幾種利用“熵”的概念來衡量如何選取合適的分解級。感興趣的同學還可以參看wentropy,wdcbm2,wpdec的幫助。文獻1中也提到了衡量壓縮質量的客觀化方法

48、MSE,PNSR并指出小波的重構濾波器的長度越長，形狀越規則越能夠提供良好的壓縮性能。上面對PNG的討論因為沒有足夠的算法分析和程序解讀，同時也沒有準確的試驗數據，因此只能作為猜測。但衡量小波壓縮效率的方法我堅持認為不能以PNG文件大小來解說，如果采用圖像文件大小來衡量，應該以JPEG2000來衡量。致謝：感謝林老師的多次指導，感謝廣東中山站的劉趙璧同學的幫助和建議并校對修改了許多錯誤，Lily（我還不知道她的姓名）為我提供的寶貴資料，以及曹紅軍、李澄鈺、鄭南、施吉鴻等在網絡上共享他們程序的同學。感謝張子亮以及網上許多同學的討論對我的啟發。附錄：1、dwt2的分解系數與第3章使用矩陣方法的分解

49、系數不一致的原因第3章小波與小波變換procedureNonstandardDecomposition(C:array.,.hh11ofreals)C=C/hwhileh>1doforrow1tohdoDecompositionStep(Crow,1.h)endfor%CischangedtoAverageandDifference%A-Lefthalf,D-righthalfforcol1tohdoDecompositionStep(C1.h,col)endfor%Cischangedtocach,cv,cd,ca-TopLeft,%ch-bottomleft,cv-topright,

50、cd-bottomrighth=h/2endwhileendprocedure1）兩個DecompositionStep都處理(C)，那么C在row處理后，原來的C已經改變為行處理后的。應當聲明。2）我仔細閱讀了教材，并使用MATLAB編寫上述算法，發現老師沒有說明產生了幾個系數。這是我分析MATLAB的dwt,dwt2,idwt2之后并親自編寫自己的nsdwt2和nsidwt2時發現的。您的原文“首先對圖像的每一行計算像素對的均值和差值，然后對每一列計算像素對的均值和差值。”設平均產生結果為A,差為D,試想每一行的平均值是h/2個A,差值是h/2個D,那么仍然存儲在C中（兩者都假設下采樣），

51、此時仍然使用DecompositionStep(C1.h,col),就會出現h/4個AA,h/4個AD(A的均值與差值）,h/4個DA,h/4個DD（A的均值與差值），這實際上就是說，2維與1維只有A與D兩個系數相比，產生了4個系數，這也就是使用cachcvcd=dwt2(C,wname)的得到的低頻系數，水平、垂直、對角系數。對于最后生成的矩陣C,左上ca,左下ch,右上cv,右下cd。雖然您在前面章節介紹了左上角矩陣AA為低頻，其余均為高頻細節系數，但如果使用dwt編寫上述算法時，很容易只考慮低頻部分而丟掉cv,cd。-dwt2源程序:前面的源代碼略去首先擴展邊界，sym->size

52、EXT=(18-1),per->18/2y=wextend("2D",dwtEXTM,x,sizeEXT,sizeEXT);在進行低頻Lo_D系數的每行的過濾得到z，這相當于行變換，只不過沒有下采樣，就是上面的A部分z=wconv("row",y,Lo_D);然后再用行變換得低頻矩陣z再一次進行列方向低頻濾波卷積下采樣，得到的低頻部分就是ca=a,對z得高頻濾波就是ch=h,再對比上面對AA,ADa=convdown(z,Lo_D,sizeKEPT,shift);h=convdown(z,Hi_D,sizeKEPT,shift);在進行高頻Hi_D

53、系數的每行的過濾得到z，這相當于行變換，只不過沒有下采樣，就是上面的D部分z=wconv("row",y,Hi_D);然后再用行變換得低頻矩陣z再一次進行列方向Lo_D濾波卷積下采樣，得到的高頻中的低頻部分就是cv=v,對z得高頻濾波就是cd=d對角部分,再對比上面對DA,DDv=convdown(z,Lo_D,sizeKEPT,shift);d=convdown(z,Hi_D,sizeKEPT,shift);如果你想使用dwt2確得到圖3-26那樣的圖像，比使用dwt要快的dwt2中的關鍵代碼的解釋：typedwt2.m與下面的源代碼解釋相對照，稍加修改，然后換一個函數名

54、稱，作業就完成了。可行的方法：對第一個z進行下采樣a1=dyaddown(wkeep(z,sizeKEPT)就是圖3-26的第一次行分解的低頻部分，對第2個的z進行上述過程得到d1就是圖3-26的高頻部分。然后合成圖像：img3_26=uint8(round(a1d1);這個速度絕對與dwt2一致，因為它就是dwt2.2、采用第3章21-24頁介紹的方法不能生成與wavedec2相同的矩陣。原來認為老師第三章的haar的2級，3級小波矩陣就只對低頻（ca)部分分解,實際上不是，ch部分也參加了計算。但是因為矩陣是正交的，重構沒有問題，但壓縮可能就會出現問題。舉例：4X4A與2級系數MA=a11

55、a12a13a14M=1/21/200a21a22a23a241/2-1/200a31a32a33a340010a41a42a43a440001此時a11a12a21a22為ca;a31a32a41a42是ch;A*M=a11*1/2+a12*1/2a11*1/2-a12*1/2a13a14a21*1/2+a22*1/2a21*1/2-a22*1/2a23a24-a31*1/2+a32*1/2a31*1/2-a31*1/2a33a34a41*1/2+a42*1/2a41*1/2-a41*1/2a43a44由此，ch參加了運算。下面是根據老師第3章21-24頁編寫的標準和非標準分解重構程序。可以

56、看出它分解的高頻部分越來越亮0的數目不是越來越多。但如果將X的賦值語句的注釋符號%去掉，可以看到與老師3.5節一樣的結果。functionshowhaardecrec%ThisprogramiscompletelyaccordingProfessorLinfuzong"sadditionalmaterials%chapter3,section3.52-Dhaarwavelettransformandsynthesis.%Hereweusethe.matfilewbarbwhichshouldbeinthedirectorywavedemo%Itusefunctionmakehaarm

57、atrixmakethehaarwavelettransformmatrixat3differentlevel.%andthen1)showingstandarddecomposition/reconstructionprocessanddisplaytheresult%inafigure2)showingnon-standarddecompositionprocessanddisplaytheresultina%figure3)showingdirectwavelettransform/synthesisusingtheproductmatrixof3-levels%matrixanddis

58、playtheimage.loadwbarb;colormap(map);%X=642361606757%955541213515016%1747462021434224%4026273736303133%3234352928383925%4123224445191848%4915145253111056%858595462631X=double(X);maxX=max(max(X);X=wcodemat(X,maxX+1);fig=figure;set(fig,"position",00800600);image(X);TICrows,cols=size(X);M1=makehaarmatrix(rows,cols,1,1);M2=makehaarmatrix(rows,cols,2,1);M3=makehaarmatrix(rows,cols,3,1);S1=X*M1;%standarddecrowS2=S1*M2;%standarddecrowS3=S2*M3;%standarddecrowS4=M1"*S3;%standarddeccolS5=M2"*S4;%standarddeccolS6=M3"*S5;%standarddeccoluiwai