第4章醫學圖像小波去噪的MATLAB實現4.1小波基確實定不同的小波基具有不同的時頻特征，用不同的小波基分析同一個問題會產生不同的結果，故小波分析在應用中便存在一個小波基或小波函數的選取和優化問題。我們在應用中要把握小波函數的特征，根據應用需要，選擇適宜的小波基。在小波分析應用中要考查小波函數或小波基的連續性、正交性、對稱性、消失矩、線性相位、時頻窗口的中心和半徑以及時頻窗的面積等，這些特征關系到如何選擇適宜的小波基。本節選取了一些常見的小波基，首先固定小波分解層數和閾值，然后改變小波基，運行結果。通過計算峰值信噪比(PSNR)來判定哪個小波基對醫學圖像去噪效果好。下表為不同小波基去噪前帶噪圖像的峰值信噪比〔PSNR〕和去噪后圖像的峰值信噪比〔PSNR〕，通過峰值信噪比對不同小波基的去噪效果進行評價，從而選出對圖像去噪效果較好的小波基。表4-1不同小波基去噪后圖像的峰值信噪比小波基sym3sym5coif2coif5db2db6帶噪圖像的PSNR17.867817.871717.872017.918417.956117.9500去噪圖像的PSNR22.781324.304218.189621.562720.4057通過去噪效果圖4-1和表4-1以及圖像評價原那么我們可以很容易選出對圖像去噪效果好，而又很好的保持圖像細節的小波基。從圖4-1中我們可以看出選用sym3小波基去噪后噪聲得到了明顯的抑制，但是圖像的細節被弱化了，讀圖有所影響。選用sym5小波基去噪后，噪聲沒有得到很好的抑制，而且圖像細節已明顯消損，對讀圖有所影響。選用coif2小波基對圖像進行去噪后，噪聲得到一定的抑制，圖像的細節保持的也很好。選用coif5小波基對圖像去噪后，圖像細節明顯消損，對讀圖有所影響。選用db2小波基對圖像去噪后圖像的噪聲雖然得到抑制但細節變得模糊，很難區分。選用db6小波基對圖像進行去噪后，圖像失真比擬明顯。從表4-1中可以看出去噪后圖像的PSNR,其中使用coif2小波基去噪后圖像的PSNR最大，通常峰值信噪比PSNR愈大愈好。實驗結果如圖4-1所示：原始圖像加噪圖像sym3小波基sym5小波基coif2小波基sym3小波基sym5小波基coif2小波基coif5coif5小波基db2小波基db6小波基圖4-1不同小波基去噪效果圖綜上所述，coif2小波基去噪效果很好，所以本次課程設計中我選擇coif2小波基進行醫學圖像小波去噪方法研究。相應小波變換語句：[c,s]=wavedec2(A,n,wname);其中：wname為小波基變量，通過改變它來變換程序中運用的小波基。n為小波分解層數A為處理圖像4.2分解層數的選擇在實際的小波圖像去噪過程中，不同信號、不同信噪比下都存在一個去噪效果最好或接近最好的分解層數，分解層數對于去噪效果的影響很大。通常分解層數過多，并且對所有的各層小波空間的系數都進行閾值處理會造成信號的信息喪失嚴重，去噪后的信噪比反而下降，同時導致運算量增大，使處理速度變慢。分解層數過少那么去噪效果不理想，信噪比提高不多。因此分解層數確實定作為一個核心的問題需要解決。11層小波分解2層小波分解3層小波分解4層小波分解5層小波分解6層小波分解圖圖4-21-6層小波分解去噪后的效果圖本論文中通過對醫學圖像進行n層小波分解，從大量的實驗數據可以看出層數超過十層后圖像以很難區分，我通過仿真實驗對醫學圖像進行了1-6層小波分解，確定小波基為coif2。仿真實驗結果如圖4-2所示：不同分解層數圖像去噪效果不同，從上面圖中我們就可以看出，為了更有說服力和可信度，我們將不同分解層數去噪后的峰值信噪比〔PSNR〕列出，見下表：表4-2不同小波分解層數去噪后的峰值信噪比分解層數NN=1N=2N=3N=4N=5N=6PSNR22.408724.723923.624520.283317.575115.7194根據上表繪制不同層次去噪后圖像峰值信噪比走勢圖如下：圖4-3不同層次去噪后圖像峰值信噪比走勢圖由圖4-3可知由本文研究方法確定的最正確分解層數分別為2和3，在沒到達此分解層數之前，分解層數對去噪結果的影響是很大的，去噪后信號的信噪比隨著分解層數的增加而迅速增加，在超過這個層數之后，去噪后信號的信噪比要么沒有太大的變化，要么就會有較大的下滑，由此需要確定一個適宜的分解層數才能到達較好的去噪效果。綜上所述：通過圖4-2、表4-2、圖4-3的綜合評價，我們選取分解層數為2層進行我們后續的研究，力求使我們研究的方法去噪效果更好。以下是相應的仿真程序clear;%清理工作空間I=imread('mm、bmp');%裝載原始圖像X1=rgb2gray(I);X=double(X1);[m,n]=size(X);x=imnoise(X1,'gaussian',0、02);%參加高斯噪聲%分層閾值處理fori=1:6[c,s]=wavedec2(x,i,'coif2');%i層小波分解[thr3,nkeep3]=wdcbm2(c,s,2);%提取i層閾值Duo(:,:,i)=wdencmp('lvd',c,s,'coif2',i,thr3,'h');%i層小波去噪psnr(i)=10*log10(m*n*max(max(X、^2))/sum(sum((double(X)-double(Duo(:,:,i)))、^2)));%i層小波去噪后圖像的峰值信噪比end%顯示各層去噪后圖像fori=1:6figure,imshow(Duo(:,:,i),[])endfigure，plot(psnr)axis([06525])title('不同層次去噪后圖像峰值信噪比走勢圖‘)4.3閾值函數的選擇現在有很多關于閾值函數選擇方面的研究，人們對閾值函數的優缺點大體上有所了解，小波閾值的選取和閾值函數的選擇直接影響著最后的結果。本文選取的方有更好的法克服了硬閾值在閾值點的不連續性和軟閾值過度平滑的缺點，實驗說明，該方法具去噪效果，能在去噪的同時保存圖像的邊緣等重要細節。閾值去噪中，閾值函數表達了對超過和低于閾值的小波系數運用不同處理策略，是閾值去噪中關鍵的一步。經小波分解后，信號的小波系數幅值要大于噪聲的系數幅值。可以認為，幅值比擬大的小波系數一般以信號為主，而幅值比擬小的系數在很大程度上是噪聲。于是，采用閾值的方法可以把信號系數保存，而使大局部噪聲系數減小至零。小波閾值收縮法去噪的具體處理過程為：將含噪信號在各尺度上進行小波分解，設定一個閾值，幅值低于該閾值的小波系數置為0，高于該閾值的小波系數或者完全保存，或者做相應的“收縮(shrinkage)〞最后將處理后獲得的小波系數用逆小波變換進行重構，得到去噪后的圖像。從表4-3中我們可以看出，自適應閾值去噪后圖像的峰值信噪比〔PSNR〕比擬大，去噪效果較好，所以后續研究中我們選用自適應閾值去噪進行研究。表4-3不同閾值函數去噪的峰值信噪比閾值和閾值函數Donoho軟閾值Donoho硬閾值Birge-Massart軟閾值Birge-Massart硬閾值半軟閾值自適應閾值加噪圖像處理前PSNR72.082072.899872.709473.195473.740874.160868.2931圖4-4為不同閾值函數的去噪效果圖，從圖中我們可以看出：硬閾值函數可以很好地保存圖像邊緣等局部特征，但圖像會出現振鈴、偽吉布斯效應等視覺失真。而軟閾值函數處理相對平滑，但可能會造成邊緣模糊等失真現象。相對于以上兩種閾值函數，半軟閾值函數能表現出很好的去噪效果，只不過半軟閾值大都對圖像有針對性，沒有在各個方面都優秀的閾值函數，并且隨著優化程度的提高，算法復雜度也越來越高，實現起來比擬困難。但相對而言改善后的閾值函數都會從某種程度上提高信噪比以及保真度。DonohoDonoho軟閾值Donoho硬閾值Birge-Massart軟閾值Birge-MassartBirge-Massart硬閾值半軟閾值自適應閾值圖4-4不同閾值函數去噪效果圖以下是相應程序代碼%Donoho全局閾值軟閾值公式---------------------------------------------x_soft_r=wdenoise(xr,'gbl','s',thr_r,'coif2',2);x_soft_g=wdenoise(xg,'gbl','s',thr_g,'coif2',2);x_soft_b=wdenoise(xb,'gbl','s',thr_b,'coif2',2);%------------------------------------------------------------------%Donoho全局閾值硬閾值公式---------------------------------------------x_hard_r=wdenoise(xr,'gbl','h',thr_r,'coif2',2);x_hard_g=wdenoise(xg,'gbl','h',thr_g,'coif2',2);x_hard_b=wdenoise(xb,'gbl','h',thr_b,'coif2',2);%--------------------------------------------------------------------%Bige-Massa策略軟閾值公式--------------------------------------------x_soft_lvd_r=wdenoise(xr,'lvd','s',thr_lvd_r,'coif2',2);x_soft_lvd_g=wdenoise(xg,'lvd','s',thr_lvd_g,'coif2',2);x_soft_lvd_b=wdenoise(xb,'lvd','s',thr_lvd_b,'coif2',2);%--------------------------------------------------------%Bige-Massa策略硬閾值公式--------------------------------------------x_hard_lvd_r=wdenoise(xr,'lvd','h',thr_lvd_r,'coif2',2);x_hard_lvd_g=wdenoise(xg,'lvd','h',thr_lvd_g,'coif2',2);x_hard_lvd_b=wdenoise(xb,'lvd','h',thr_lvd_b,'coif2',2);%----------------------------------------------------------------------%半軟閾值--------------------------------------------------------------x1_r=den1(xr,'coif2',2,thr_r);x1_g=den1(xg,'coif2',2,thr_g);x1_b=den1(xb,'coif2',2,thr_b);%-------------------------------------------------------------------%自適應閾值---------------------------------------------------------x4_r=den4(xr,'coif2',2);x4_g=den4(xg,'coif2',2);x4_b=den4(xb,'coif2',2);%-------------------------------------------------------------------%恢復去噪后的圖像========================================================x_soft=cat(3,x_soft_r,x_soft_g,x_soft_b);%Donoho軟閾值x_hard=cat(3,x_hard_r,x_hard_g,x_hard_b);%Donoho硬閾值x_soft_lvd=cat(3,x_soft_lvd_r,x_soft_lvd_g,x_soft_lvd_b);%Birge-Massart軟閾值x_hard_lvd=cat(3,x_hard_lvd_r,x_hard_lvd_g,x_hard_lvd_b);%Birge-Massart硬閾值x1=cat(3,x1_r,x1_g,x1_b);%半軟閾值x1_5=cat(3,x1_5_r,x1_5_g,x1_5_b);%半軟閾值+均值濾波x4=cat(3,x4_r,x4_g,x4_b);%自適應閾值=========================================================================psnr_soft=PSNR_color(x_soft,X);psnr_hard=PSNR_color(x_hard,X);psnr_soft_lvd=PSNR_color(x_soft_lvd,X);psnr_hard_lvd=PSNR_color(x_hard_lvd,X);psnr1=PSNR_color(x1,X);psnr1_5=PSNR_color(x1_5,X);x4=rgb2gray(x4);psnr4=PSNR(x4,X);==================================================================================================================================================functionpsnr=psnr_color(I,J)%計算消噪前后圖像的峰值信噪比psnr=10*log10(255^2/MSE_color(I,J));4.4與傳統去噪方法比擬通過前面的研究確定一個完整去噪算法中影響去噪性能的各種因素。在實際的圖像處理中，為小波閾值去噪法的選擇和改善最終確定了一個比擬好的去噪效果。然后通過與傳統去噪方法進行比擬來證實我們最終確定的方法對醫學圖像去噪有較好的優勢。本文中最后確定的各變量為：coif2小波基、2層小波分解層數、自適應閾值、軟閾值函數等等，通過這些變量確實定，使我們所研究的小波自適應閾值去噪方法的去噪效果更優。由圖4-5我們可以看出自適應閾值去噪后的圖像比傳統方法去噪后圖像效果明顯好。通過軟件仿真我們可以清晰的看出不同方法的優缺點，更加明確了本文研究的小波閾值去噪的優勢。原始圖像加噪圖像3*3的鄰域窗的中值濾波圖像均值濾波圖像均值濾波圖像自適應閾值5*5的鄰域窗的中值濾波圖像維納處理維納處理3維納處理5維納處理7圖4-5小波閾值去噪與傳統去噪方法效果圖以下是相應的程序代碼%對彩色圖像進行去噪I=imread('mm、bmp');%讀入圖像X=im2double(I);%轉換成雙精度類型x_noise=imnoise(X,'gaussian',0、02);%參加高斯噪聲%提取三個通道信息xr=x_noise(:,:,1);%R通道xg=x_noise(:,:,2);%G通道xb=x_noise(:,:,3);%B通道%估計三個通道的閾值[Cr,Sr]=wavedec2(xr,2,'coif2');[Cg,Sg]=wavedec2(xg,2,'coif2');[Cb,Sb]=wavedec2(xb,2,'coif2');thr_r=Donoho(xr);%R通道全局閾值thr_g=Donoho(xg);%G通道全局閾值thr_b=Donoho(xb);%B通道全局閾值thr_lvd_r=Birge_Massart(Cr,Sr);%R通道局部閾值thr_lvd_g=Birge_Massart(Cg,Sg);%G通道局部閾值thr_lvd_b=Birge_Massart(Cb,Sb);%B通道局部閾值%-------------------------------------------------------------------%自適應閾值---------------------------------------------------------x4_r=den4(xr,'coif2',2);x4_g=den4(xg,'coif2',2);x4_b=den4(xb,'coif2',2);%-------------------------------------------------------------------%恢復去噪后的圖像========================================================x4=cat(3,x4_r,x4_g,x4_b);%自適應閾值=========================================================================psnr4=PSNR(x4,X);figure;imshow(x4);title('自適應閾值');%========================================================================functionX=den4(x,wname,n)%"FeatureAdaptiveWaveletShrinkageforImageDenoising"%初始化參數值R=5;%窗口大小alpha=0.1;%控制小波系數收縮減的程度beta=0.3;delta=DELTA(x);%噪方差lambda2=4*delta^2*log(R);%局部閾值[C,S]=wavedec2(x,n,wname);%對圖像進行小波分解%提取每層系數并進行處理fori=n:-1:1cH=detcoef2('h',C,S,i);%水平細節系數cV=detcoef2('v',C,S,i);%垂直細節系數cD=detcoef2('d',C,S,i);%對角線細節系數dim=size(cH);%分別處理三個方向的系數forj=1:dim(1)fork=1:dim(2)S_jk2=energy(cH,j,k,R);cH(j,k)=shrink(cH(j,k),S_jk2,alpha,beta,lambda2);S_jk2=energy(cV,j,k,R);cV(j,k)=shrink(cV(j,k),S_jk2,alpha,beta,lambda2);S_jk2=energy(cD,j,k,R);cD(j,k)=shrink(cD(j,k),S_jk2,alpha,beta,lambda2);endend%再把系數放回去k=size(S,1)-i;first=S(1,1)*S(1,2)+3*sum(S(2:k-1,1)、*S(2:k-1,2))+1;%起始位置add=S(k,1)*S(k,2);%系數長度C(first:first+add-1)=reshape(cH,1,add);C(first+add:first+2*add-1)=reshape(cV,1,add);C(first+2*add:first+3*add-1)=reshape(cD,1,add);endX=waverec2(C,S,wname);%重構圖像%%%%%%delta%%%functiondelta=DELTA(x)%估計噪聲方差%[C,S]=wavedec2(x,1,'db1');%小波分解d=C(pro