小波去噪MATLAB实现.docx
《小波去噪MATLAB实现.docx》由会员分享，可在线阅读，更多相关《小波去噪MATLAB实现.docx（19页珍藏版）》请在冰豆网上搜索。
小波去噪MATLAB实现
第4章医学图像小波去噪的MATLAB实现
4.1小波基的确定
不同的小波基具有不同的时频特征，用不同的小波基分析同一个问题会产生不同的结果，故小波分析在应用中便存在一个小波基或小波函数的选取和优化问题。
我们在应用中要把握小波函数的特征，根据应用需要，选择合适的小波基。
在小波分析应用中要考查小波函数或小波基的连续性、正交性、对称性、消失矩、线性相位、时频窗口的中心和半径以及时频窗的面积等，这些特征关系到如何选择合适的小波基。
本节选取了一些常见的小波基，首先固定小波分解层数和阈值，然后改变小波基，运行结果。
通过计算峰值信噪比（PSNR）来判定哪个小波基对医学图像去噪效果好。
下表为不同小波基去噪前带噪图像的峰值信噪比（PSNR）和去噪后图像的峰值信噪比（PSNR），通过峰值信噪比对不同小波基的去噪效果进行评价，从而选出对图像去噪效果较好的小波基。
表4-1不同小波基去噪后图像的峰值信噪比
小波基
sym3
sym5
coif2
coif5
db2
db6
带噪图像的PSNR
17.8678
17.8717
17.8720
17.9184
17.9561
17.9500
去噪图像的PSNR
19.4783
22.7813
24.3042
18.1896
21.5627
20.4057
通过去噪效果图4-1和表4-1以及图像评价原则我们可以很容易选出对图像去噪效果好，而又很好的保持图像细节的小波基。
从图4-1中我们可以看出选用sym3小波基去噪后噪声得到了明显的抑制，但是图像的细节被弱化了，读图有所影响。
选用sym5小波基去噪后，噪声没有得到很好的抑制，而且图像细节已明显消损，对读图有所影响。
选用coif2小波基对图像进行去噪后，噪声得到一定的抑制，图像的细节保持的也很好。
选用coif5小波基对图像去噪后，图像细节明显消损，对读图有所影响。
选用db2小波基对图像去噪后图像的噪声虽然得到抑制但细节变得模糊，很难辨别。
选用db6小波基对图像进行去噪后，图像失真比较明显。
从表4-1中可以看出去噪后图像的PSNR,其中使用coif2小波基去噪后图像的PSNR最大，通常峰值信噪比PSNR愈大愈好。
实验结果如图4-1所示：
原始图像加噪图像
sym3小波基sym5小波基coif2小波基
coif5小波基db2小波基db6小波基
图4-1不同小波基去噪效果图
综上所述，coif2小波基去噪效果很好，所以本次课程设计中我选择coif2小波基进行医学图像小波去噪方法研究。
相应小波变换语句：
[c,s]=wavedec2（A,n,wname）;
其中：
wname为小波基变量，通过改变它来变换程序中运用的小波基。
n为小波分解层数
A为处理图像
4.2分解层数的选择
在实际的小波图像去噪过程中，不同信号、不同信噪比下都存在一个去噪效果最好或接近最好的分解层数，分解层数对于去噪效果的影响很大。
通常分解层数过多，并且对所有的各层小波空间的系数都进行阈值处理会造成信号的信息丢失严重，去噪后的信噪比反而下降，同时导致运算量增大，使处理速度变慢。
分解层数过少则去噪效果不理想，信噪比提高不多。
因此分解层数的确定作为一个核心的问题需要解决。
1层小波分解2层小波分解3层小波分解
4层小波分解5层小波分解6层小波分解
图4-21-6层小波分解去噪后的效果图
本论文中通过对医学图像进行n层小波分解，从大量的实验数据可以看出层数超过十层后图像以很难辨别，我通过仿真实验对医学图像进行了1-6层小波分解，确定小波基为coif2。
仿真实验结果如图4-2所示：
不同分解层数图像去噪效果不同，从上面图中我们就可以看出，为了更有说服力和可信度，我们将不同分解层数去噪后的峰值信噪比（PSNR）列出，见下表：
表4-2不同小波分解层数去噪后的峰值信噪比
分解层数N
N=1
N=2
N=3
N=4
N=5
N=6
PSNR
22.4087
24.7239
23.6245
20.2833
17.5751
15.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:
6
figure,
imshow（Duo（:
:
i）,[]）
end
figure，
plot（psnr）
axis（[06525]）
title（'不同层次去噪后图像峰值信噪比走势图‘）
4.3阈值函数的选择
现在有很多关于阈值函数选择方面的研究，人们对阈值函数的优缺点大体上有所了解，小波阈值的选取和阈值函数的选择直接影响着最后的结果。
本文选取的方有更好的法克服了硬阈值在阈值点的不连续性和软阈值过度平滑的缺点，实验表明，该方法具去噪效果，能在去噪的同时保留图像的边缘等重要细节。
阈值去噪中，阈值函数体现了对超过和低于阈值的小波系数运用不同处理策略，是阈值去噪中关键的一步。
经小波分解后，信号的小波系数幅值要大于噪声的系数幅值。
可以认为，幅值比较大的小波系数一般以信号为主，而幅值比较小的系数在很大程度上是噪声。
于是，采用阈值的办法可以把信号系数保留，而使大部分噪声系数减小至零。
小波阈值收缩法去噪的具体处理过程为：
将含噪信号在各尺度上进行小波分解，设定一个阈值，幅值低于该阈值的小波系数置为0，高于该阈值的小波系数或者完全保留，或者做相应的“收缩（shrinkage）”最后将处理后获得的小波系数用逆小波变换进行重构，得到去噪后的图像。
从表4-3中我们可以看出，自适应阈值去噪后图像的峰值信噪比（PSNR）比较大，去噪效果较好，所以后续研究中我们选用自适应阈值去噪进行研究。
表4-3不同阈值函数去噪的峰值信噪比
阈值和阈值函数
Donoho软阈值
Donoho硬阈值
Birge-Massart软阈值
Birge-Massart硬阈值
半软
阈值
自适
应阈值
加噪图像处理前
PSNR
72.0820
72.8998
72.7094
73.1954
73.7408
74.1608
68.2931
图4-4为不同阈值函数的去噪效果图，从图中我们可以看出：
硬阈值函数可以很好地保留图像边缘等局部特征，但图像会出现振铃、伪吉布斯效应等视觉失真。
而软阈值函数处理相对平滑，但可能会造成边缘模糊等失真现象。
相对于以上两种阈值函数，半软阈值函数能表现出很好的去噪效果，只不过半软阈值大都对图像有针对性，没有在各个方面都优秀的阈值函数，并且随着优化程度的提高，算法复杂度也越来越高，实现起来比较困难。
但相对而言改善后的阈值函数都会从某种程度上提高信噪比以及保真度。
Donoho软阈值Donoho硬阈值Birge-Massart软阈值
Birge-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:
1
cH=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）;
end
end
%再把系数放回去
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）;
end
X=waverec2（C,S,wname）;%重构图像
%%%
%%%delta
%%%
functiondelta=DELTA（x）
%估计噪声方差
%
[C,S]=wavedec2（x,1,'db1'）;%小波分解
d=C（prod（S（1,:
））+2*prod（S（2,:
））+1:
end）;%HH子带系数
delta=median（abs（d））/0、6745;%计算delta
%%%
%%%energy
%%%
functionS_jk2=energy（cM,j,k,R）
%计算小波系数附近的能量
%
dim=size（cM）;
%边界判断
row_min=（j-1=fix（R/2））*fix（-R/2）;
row_max=（dim
（1）-j（1）-j）+（dim
（1）-j>=fix（R/2））*fix（R/2）;
col_min=（k-1=fix（R/2））*fix（-R/2）;
col_max=（dim
（2）-k（2）-k）+（dim
（2）-k>=fix（R/2））*fix（R/2）;
s=0;
form=row_min:
row_max
forn
升级会员 