小波去噪MATLAB实现

第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