Matlab實現小波去噪
小波去噪MATLAB實現
第二章 影象小波去噪理論
第4章 醫學影象小波去噪的MATLAB實現
4。1 小波基的確定
不同的小波基具有不同的時頻特徵，用不同的小波基分析同一個問題會產生不同的結果，故小波分析在應用中便存在一個小波基或小波函式的選取和最佳化問題。我們在應用中要把握小波函式的特徵，根據應用需要，選擇合適的小波基。在小波分析應用中要考查小波函式或小波基的連續性、正交性、對稱性、消失矩、線性相位、時頻視窗的中心和半徑以及時頻窗的面積等，這些特徵關係到如何選擇合適的小波基。本節選取了一些常見的小波基，首先固定小波分解層數和閾值，然後改變小波基，執行結果。透過計算峰值信噪比（PSNR）來判定哪個小波基對醫學影象去噪效果好。
下表為不同小波基去噪前帶噪影象的峰值信噪比（PSNR）和去噪後圖像的峰值信噪比（PSNR），透過峰值信噪比對不同小波基的去噪效果進行評價，從而選出對影象去噪效果較好的小波基。
透過去噪效果圖4-1和表4-1以及影象評價原則我們可以很容易選出對影象去噪效果好，而又很好的保持影象細節的小波基。從圖4-1中我們可以看出選用sym3小波基去噪後噪聲得到了明顯的抑制，但是影象的細節被弱化了，讀圖有所影響。選用sym5小波基去噪後，噪聲沒有得到很好的抑制，而且影象細節已明顯消損，對讀圖有所影響。選用coif2小波基對影象進行去噪後，噪聲得到一定的抑制，影象的細節保持的也很好。選用coif5小波基對影象去噪後，影象細節明顯消損，對讀圖有所影響。選用db2小波基對影象去噪後圖像的噪聲雖然得到抑制但細節變得模糊，很難辨別。選用db6小波基對影象進行去噪後，影象失真比較明顯。綜上所述，coif2小波基去噪效果很好，所以本次課程設計中我選擇coif2小波基進行醫學影象小波去噪方法研究。
相應小波變換語句：［c，s］=wavedec2（A，n，wname）；
其中：wname為小波基變數，透過改變它來變換程式中運用的小波基。
n為小波分解層數
A為處理影象
4。2 分解層數的選擇
在實際的小波影象去噪過程中，不同訊號、不同信噪比下都存在一個去噪效果最好或接近最好的分解層數，分解層數對於去噪效果的影響很大。 通常分解層數過多，並且對所有的各層小波空間的係數都進行閾值處理會造成訊號的資訊丟失嚴重，去噪後的信噪比反而下降，同時導致運算量增大，使處理速度變慢。
分解層數過少則去噪效果不理想，信噪比提高不多。 因此分解層數的確定作為一個核心的問題需要解決。
本論文中透過對醫學影象進行n層小波分解，從大量的實驗資料可以看出層數超過十層後圖像以很難辨別，我透過模擬實驗對醫學影象進行了1-6層小波分解，確定小波基為coif2。模擬實驗結果如圖4-2所示：
不同分解層數影象去噪效果不同，從上面圖中我們就可以看出，為了更有說服力
由圖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）； %加入高斯噪聲
%分層閾值處理
for i=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
% 顯示各層去噪後圖像
for i=1：6
figure，
imshow（Duo（：，：，i），［］）
end
figure，
plot（psnr）
axis（［0 6 5 25］）
title（‘不同層次去噪後圖像峰值信噪比走勢圖‘）
4。3 閾值函式的選擇
現在有很多關於閾值函式選擇方面的研究，人們對閾值函式的優缺點大體上有所瞭解，小波閾值的選取和閾值函式的選擇直接影響著最後的結果。本文選取的方有更好的法克服了硬閾值在閾值點的不連續性和軟閾值過度平滑的缺點，實驗表明，該方法具去噪效果，能在去噪的同時保留影象的邊緣等重要細節。閾值去噪中，閾值函式體現了對超過和低於閾值的小波係數運用不同處理策略，是閾值去噪中關鍵的一步。經小波分解後，訊號的小波係數幅值要大於噪聲的係數幅值。可以認為，幅值比較大的小波係數一般以訊號為主，而幅值比較小的係數在很大程度上是噪聲。於是，採用閾值的辦法可以把訊號係數保留，而使大部分噪聲係數減小至零。小波閾值收縮法去噪的具體處理過程為：將含噪訊號在各尺度上進行小波分解，設定一個閾值，幅值低於該閾值的小波係數置為0，高於該閾值的小波係數或者完全保留，或者做相應的“收縮（shrinkage）”最後將處理後獲得的小波係數用逆小波變換進行重構，得到去噪後的影象。
從表4-3中我們可以看出，自適應閾值去噪後圖像的峰值信噪比（PSNR）比較大，去噪效果較好，所以後續研究中我們選用自適應閾值去噪進行研究。
圖4-4為不同閾值函式的去噪效果圖，從圖中我們可以看出：硬閾值函式可以很好地保留影象邊緣等區域性特徵，但影象會出現振鈴、偽吉布斯效應等視覺失真。而軟閾值函式處理相對平滑，但可能會造成邊緣模糊等失真現象。相對於以上兩種閾值函式，半軟閾值函式能表現出很好的去噪效果，只不過半軟閾值大都對影象有針對性，沒有在各個方面都優秀的閾值函式，並且隨著最佳化程度的提高，演算法複雜度也越來越高，實現起來比較困難。但相對而言改善後的閾值函式都會從某種程度上提高信噪比以及保真度。
硬閾值 Birge-Massart Donoho軟閾值 Donoho 軟閾值
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）；
========================================================================= ========================================================================= function psnr = 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（’自適應閾值‘）；
% =====================================================================
function X = den4（x， wname， n）
% “Feature Adaptive Wavelet Shrinkage for Image Denoising”
% 初始化引數值
R = 5； % 視窗大小
alpha = 0。1； %控制小波係數收縮減的程度 beta = 0。3；
delta = DELTA（x）； % 噪方差
lambda2 = 4 * delta^2 * log（R）； %區域性閾值
［C， S］ = wavedec2（x， n， wname）； %對影象進行小波分解 %提取每層係數並進行處理
for i = n ： -1 ： 1
cH = detcoef2（’h‘， C， S， i）； %水平細節係數
cV = detcoef2（’v‘， C， S， i）； %垂直細節係數
cD = detcoef2（’d‘， C， S， i）； % 對角線細節係數
dim = size（cH）；
%分別處理三個方向的係數
for j = 1 ： dim（1）
for k = 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
function delta = 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
function S_jk2 = energy（cM， j， k， R）
% 計算小波係數附近的能量
dim = size（cM）；
%邊界判斷
row_min = （j-1 < fix（R/2）） * （1-j） + （j-1 >= fix（R/2）） * fix（-R/2）；
row_max = （dim（1）-j < fix（R/2）） * （dim（1）-j） + （dim（1）-j >= fix（R/2）） * fix（R/2）； col_min = （k-1 < fix（R/2）） * （1-k） + （k-1 >= fix（R/2）） * fix（-R/2）；
col_max = （dim（2）-k < fix（R/2）） * （dim（2）-k） + （dim（2）-k >= fix（R/2）） * fix（R/2）； s = 0；
for m = row_min ： row_max
for n = col_min ： col_max
s = cM（j + m， k + n）^2 + s；
end
end
S_jk2 = s / R^2；
%%% shrink
function d_jk = shrink（d， S_jk2， alpha， beta， lambda2）
% 處理小波係數
if S_jk2 >= beta * lambda2
d_jk = d * （1 - alpha * lambda2 / S_jk2）；
else
d_jk = 0；
End
4。5 本章小結
本章透過MATLAB實驗模擬對前面相應的理論知識進行實踐，根據要確定的變數編寫相應的程式對醫學影象進行處理，透過處理後的影象效果和去噪後圖像的峰值信噪比確定 去噪效果較好的變數，最終確定本文中所有的變數，已達到對影象去噪效果最優的選擇。透過本章的總結我們對小波閾值去噪的細節更加清晰，能夠自己除錯程式，並改變相應的引數使影象去噪效果更好。而且我們還學會了判斷醫學影象去噪好壞的評價標準，能夠更好的評價一副影象的質量。實際應用中最常用的閾值函式是軟閾值函式。在閾值及閾值函式的選取中，沒有最好的，只有最合適的，在運用過程中只有透過不斷的試驗，才能找到合適的處理方法。
