[發明專利]基于小波變換和引導濾波器的醫學超聲圖像去噪方法

在審

• 申請號: 201410561594.3
• 申請日: 2014-10-21
• 公開（公告）號: CN104318527A
• 公開（公告）日: 2015-01-28
• 發明（設計）人: 張聚;林廣闊;胡標標;吳麗麗;楊立;林忠款;胡立偉
• 申請（專利權）人: 浙江工業大學
• 主分類號: G06T5/00
• 分類號: G06T5/00
• 代理公司: 杭州天正專利事務所有限公司 33201
• 代理人: 王兵;黃美娟

1.一種基于小波變換和引導濾波器的醫學超聲圖像去噪方法,包括以下步驟：

步驟1)醫學超聲圖像模型的建立；

超聲成像系統采集的包絡信號由兩部分組成，一是有意義的體內組織的反射信號，另一部分是噪聲信號；其中噪聲信號可分為相乘噪聲與相加噪聲；相乘噪聲與超聲信號成像的原理有關，主要來源于隨機的散射信號；相加噪聲認為是系統噪聲，如傳感器的噪聲等；超聲成像系統初步得到的包絡信號為f^pre，它的一般模型如下

f^pre＝g^pren^pre+w^pre????(1)

這里，上標^pre表示系統初步得到的信號；函數g^pre表示無噪聲信號，n^pre和w^pre分別表示相乘噪聲和相加噪聲，式中n^pre是噪聲的主要成分；

和相乘噪聲n^pre相比，相加噪聲w^pre所占比重很小，因此將w^pre忽略后的模型為

f^pre＝g^pren^pre???(2)

為了適應超聲成像系統顯示屏幕的動態顯示范圍，對超聲成像系統采集到的包絡信號進行對數壓縮處理；此時相乘的式(2)模型將變為相加的模型，如下

log(f^pre)＝log(g^pre)+log(n^pre)???(3)

此時，得到的信號log(f^pre)即是通常看到的醫學超聲圖像；

步驟2)對第一步得到的對數變換后的圖像進行小波分解，得到四個頻域(LL¹、LH¹、HL¹和HH¹)；對低頻域LL¹繼續進行小波分解，再得到四個頻域(LL²、LH²、HL²和HH²)；然后重復這個步驟，直到分解最大層數J；

由于小波變換是線性變換，因此式(3)模型經過二維離散小波變換后得到下面模型：

Wl,kj(log(fpre))=Wl,kj(log(gpre))+Wl,kj(log(npre))---(4)]]>

其中和分別表示含有噪聲圖像的小波系數、無噪聲圖像的小波系數和斑點噪聲的小波系數；其中上標j為小波變換的分解層數，下標(l,k)為小波域內的坐標；為了方便表示，將式(4)簡化為

Fl,kj=Gl,kj+Nl,kj---(5)]]>

對于離散的二維圖像f(n,m)，對其進行二維小波分解的步驟為：首先對圖像的每一行像素進行一維離散小波分解，然后對再圖像的每一列進行一維離散小波分解，這樣便將一幅圖像分解為四個子頻帶信號；同理，對二維小波進行重構與上述過程非常相似，即按照相反的順序進行處理，即可得到，因此這里不再贅述；接下來將對小波分解后的一些子頻帶分量，即相應的小波分解系數作簡單的分析；LL0為原始信號，圖像的信息都集中在這里；每次小波分解都會得到四個子頻帶，對LL0進行一級小波分解后得到LL1、LH1、HL1和HH1四個子頻帶；

LL1分量是對原始信號LL0的列和行進行小波分解后得到的低頻分量，即一級小波分解后近似部分，它包含了原始圖像最多的低頻信息；

LH1是一次小波分解后的垂直方向上的高頻分量，即它包含了圖像水平方向上的近似信息和垂直方向上的邊緣等高頻信息；

HL1是一次小波分解后的水平方向上的高頻分量，即它包含了圖像垂直方向上的近似信息和水平方向上的邊緣等高頻信息；

HH1是一次小波分解后對角方向上的高頻分量，即它包含了圖像水平和垂直方向上的邊緣等高頻信息；

經過小波分解后的無噪信號的小波系數服從廣義拉普拉斯分布，其概率分布如下

pG(g)=v2sΓ(1/v)exp(-|gs|v),s,v>0,---(6)]]>

式中，是伽馬函數，v為形狀參數，s為尺度參數，u為位置參數；當v＝1，u＝0時，式(6)將變為拉普拉斯分布，它是廣義拉普拉斯分布的特殊模型；

同時斑點噪聲的小波系數服從零均值高斯分布

pN(n)=12πσNexp(-n22σN2)---(7)]]>

式中σ_N為小波域內噪聲的標準差；

步驟3)利用引導濾波器對最后一層的低頻部分(LL^J)中的小波系數做濾波處理；

引導濾波器對低頻域內的小波系數作濾波處理，基本原理如下式所示

qi=ΣjWij(I)pj---(8)]]>

式(8)中，I為引導圖像，p為輸入圖像，q為輸出圖像，W_ij為關于引導圖像I的函數，i和j為像素點的位置，I由具體問題確定，可以令I＝p；

假設在窗口w_k內，中心點為后k，q為I的線性變換，如式(9)所示

qi=akIi+bk,∀i∈wk---(9)]]>

在圖像濾波中，希望能在達到濾波效果的前提下最小化輸入圖像和輸出圖像的差異，減小原始圖像細節的損失，故通過最小化p和q的差來確定系數a_k和b_k，即使式(10)最小

E(ak,bk)=Σi∈wk[(akIi+bk-pi)2+ϵak2]---(10)]]>

式(10)中，ε是正則化參數，目的是為了防止a_k過大；求解式(10)，得

ak=1|w|Σi∈wkIipi-μkp‾kσk2+ϵ]]>

bk=p‾k-akμk---(11)]]>

p‾k=1|w|Σi∈wkpi]]>

式中，μ_k和分別為I在w_k中的均值和方差；|w|為w_k中的像素個數，是輸入圖像p在w_k中的均值；求出該線性模型后，帶入整幅圖像，由于每一個像素點會有多個包含該像素點的窗口w_k，所以當在不同窗口w_k計算時，q_i的值會不同；故需要對其進行平均處理

qi=1|w|Σk,i∈wk(akIi+bk)=a‾iIi+b‾i---(12)]]>

式中，a‾i=1|w|Σak,]]>b‾i=1|w|Σbk;]]>

綜上所述，核函數W_ij可以定義如下

Wij=1|w|2Σk,(i,j)∈wk(1+(Ii-μk)(Ij-μk)σk2+ϵ)---(13)]]>

由以上原理可知，引導濾波器去噪的過程如下：

(3.1)輸入圖像p；

(3.2)輸入濾波窗口w_k的尺寸和正則化參數ε；

(3.3)計算I,p和I*p的均值；

(3.4)計算(I,p)的協方差；

(3.5)計算(I*I)的均值并計算I的方差；

(3.6)計算系數a,b；

(3.7)分別計算a和b的均值；

(3.8)得到輸出圖像q；

步驟4)對每一層的高頻部分(LH^j、HL^j和HH^j，j＝1,2,...,J)的小波系數進行閾值法收縮處理；

無噪信號的小波系數服從廣義拉普拉斯分布，小波域內的斑點噪聲部分服從高斯分布；選擇v＝1，則式(6)變為拉普拉斯分布

pG(g)=12sexp(-|g|s),s>0---(17)]]>

為了得到小波域內的信號估計值，使用貝葉斯最大后驗估計的方法；在后驗概率的計算過程中，使用貝葉斯公式如下

pG|F(g|f)=1pF(f)pF|G(f|g)·pG(g)=1pF(f)pN(f-g)·pG(g)---(18)]]>

將式(7)、式(17)帶入上式(18)，得到

pG|F(g|f)=1pF(f)·122πsσN×exp(2σN2|g|-s(f-g)22sσN2)---(19)]]>

為了得到最大后驗概率，將ln(p_G|F(g|f))對g求一次導數的方程置零，最后得到

g^=sign(f)·max(|f|-σN2s,0)---(20)]]>

為g的估計，并且假設f和無噪信號g同號；這樣就得到新的收縮方法

g^=0,f<=Tjsign(f)·max(|f|-σn2s,0),f>Tj---(21)]]>

式中只有尺度s是未知的，可由下式確定

s=[0.5(σF,j2-σN2)]2---(22)]]>

其中σ_F,j為噪聲圖像小波系數在j層的標準差；

步驟5)作小波逆變換處理，得到去噪后的醫學超聲圖像；對第5步得到的超聲圖像作指數變換，得到去噪后的超聲包絡信號。