本發(fā)明屬于圖像處理技術(shù)領(lǐng)域。具體來說，涉及一種以去除MRI萊斯噪聲，提高信噪比為目的的基于卡方無偏估計的鄰域收縮MRI去噪方法。

背景技術(shù)：

現(xiàn)代醫(yī)學(xué)中，磁共振成像(MRI)作為一項基本的醫(yī)學(xué)成像技術(shù)，已經(jīng)成為了醫(yī)生診斷和治療的重要輔助手段，磁共振圖像在獲取過程中由于硬件電路和人體的因素會引入噪聲，噪聲降低了MRI的質(zhì)量，使得一些組織的邊界變得模糊，細(xì)微結(jié)構(gòu)難以辨別，增加了對圖像細(xì)節(jié)的識別分析和處理的難度，影響醫(yī)學(xué)診斷。因此，針對噪聲MRI用軟件去噪是一件實(shí)用方便且性價比高的事情。

常用的MRI獲取方法是對頻率域(k域)內(nèi)感興趣的目標(biāo)進(jìn)行采樣，然后對頻率域數(shù)據(jù)進(jìn)行離散傅里葉反變換，得到包含實(shí)部和虛部的復(fù)數(shù)圖像數(shù)據(jù)，再轉(zhuǎn)換為幅值數(shù)據(jù)，在實(shí)部和虛部引入的加性高斯白噪聲，轉(zhuǎn)化為幅值圖像后，呈現(xiàn)萊斯(Rician)分布。當(dāng)信噪比較高時，萊斯分布傾向于高斯分布，當(dāng)信噪比低時，傾向于瑞利分布，因此針對低信噪比的MRI，用傳統(tǒng)的去除高斯噪聲的方法進(jìn)行噪聲估計和去噪會產(chǎn)生較大的偏差。

目前，MRI去噪中一個重要問題是保留有用的邊緣、細(xì)節(jié)信息，因為醫(yī)學(xué)圖像事關(guān)病患安全，任何微小的失誤都是不允許的。MRI去噪通常有三種方式，一是在幅值域針對萊斯分布的模型去噪，二是在幅值的平方域去噪，也有對三維的MRI進(jìn)行去噪的。1999年Nowak提出，平方使得噪聲圖像產(chǎn)生了一個與信號無關(guān)的、加性的偏置，經(jīng)小波變換后該偏置在尺度系數(shù)中保留，直接減去偏置即可提升去噪效果，而且平方幅值圖像服從兩個自由度的非中心卡方分布。Nowak將平方圖像的尺度系數(shù)去除偏置，對高頻小波系數(shù)用卡方分布的性質(zhì)計算信號方差，再應(yīng)用類似維納濾波的逐點(diǎn)比例收縮方法，得到優(yōu)于傳統(tǒng)小波收縮的去噪效果。Anand等人(2010年)在Nowak的基礎(chǔ)上，將平方幅值圖像經(jīng)平穩(wěn)小波變換后，高低頻分別采用NeighShrinkSURE和雙邊濾波器進(jìn)行去噪，對高頻系數(shù)的處理不再是逐點(diǎn)處理，而是考慮了像素的鄰域。在低頻系數(shù)中，含有大部分的圖像能量，信噪比高，而雙邊濾波器良好的保持邊緣的能力，可以去模糊，該方法較好的濾除了噪聲。Luisier等人(2012年)在論文“CURE for noisy Magnetic Resonance Imgaes：Chi-square unbiased risk estimation”針對平方MRI的非中心卡方分布模型，推導(dǎo)出針對卡方分布的均方誤差的無偏估計表達(dá)式(CURE)，應(yīng)用到兩種線性閾值擴(kuò)展去噪方法(LET)中，其一為圖像域的CURE估計確定小波域UWT變換下的LET參數(shù)，其二為用小波域的CURE估計確定小波域的LET的參數(shù)，由于非中心卡方分布的特殊性質(zhì)，小波域的CURE估計只能應(yīng)用在未歸一化的哈爾(Haar)小波變換。

技術(shù)實(shí)現(xiàn)要素：

為了去除噪聲，提高M(jìn)RI的信噪比，本發(fā)明推導(dǎo)了在小波域僅取決于小波系數(shù)的閾值收縮方法的CURE估計表達(dá)式，使其可以適用于鄰域收縮去噪方法上，結(jié)合低頻部分的雙邊濾波器處理和快速的循環(huán)移位技術(shù)，取得了較好的MRI去噪效果。

基于卡方無偏估計的鄰域收縮MRI去噪方法，包括如下步驟：

(1)在MRI的背景區(qū)域估計噪聲標(biāo)準(zhǔn)差σ；

(2)對噪聲圖像m進(jìn)行平方，然后除以噪聲標(biāo)準(zhǔn)差σ的平方，得圖像y

(3)對y進(jìn)行未歸一化的平穩(wěn)哈爾小波變換(SWT)，分解L層，得高頻系數(shù)和低頻系數(shù)，高頻系數(shù)有L*3個子帶；

(4)對低頻系數(shù)用雙邊濾波器去模糊，對高頻系數(shù)用基于卡方無偏估計的鄰域收縮法(NeighShrinkCURE)去噪；

(5)對去噪后的小波系數(shù)連同低頻系數(shù)進(jìn)行循環(huán)移位、小波反變換(ISWT)，得到若干循環(huán)移位的去噪結(jié)果；

(6)將上述若干循環(huán)移位的去噪結(jié)果通過非線性映射、求平均得到最終的去噪圖像。

所述的NeighShrinkCURE包括兩個部分，其一是小波域的取決于小波系數(shù)的閾值收縮方法的CURE估計表達(dá)式的推導(dǎo)，其二是鄰域收縮方法的CURE估計值的計算。

1、小波域關(guān)于小波系數(shù)的閾值收縮方法的CURE估計表達(dá)式的推導(dǎo)過程如下：

如果用表示一個N_x×N_y大小的圖像的索引集合，用μ＝[μ_n]_n∈Ω∈C^N表示理想(無噪)的大小為N＝N_x×N_y的MR復(fù)數(shù)圖像，m＝[m_n]_n∈Ω∈C^N表示含噪聲的MR圖像，去噪問題轉(zhuǎn)換為從m中估計μ的問題。在MR數(shù)據(jù)獲取過程中，實(shí)部和虛部會分別引入標(biāo)準(zhǔn)差為σ的高斯噪聲，其統(tǒng)計特性可表述為：

表示取復(fù)數(shù)的實(shí)部，表示取復(fù)數(shù)的虛部。我們的目標(biāo)是估計μ_n的幅值，定義x和y如下：

則y服從非中心的卡方分布，非中心參數(shù)為x，自由度K＝2，去噪問題又可以轉(zhuǎn)換為從非中心卡方分布的數(shù)據(jù)y中估計非中心參數(shù)x的問題。在設(shè)計去噪函數(shù)f(y)來估計x時，一個自然的評價標(biāo)準(zhǔn)是均方誤差(MSE)：

從式中可以看出，計算MSE需要無噪信號x，實(shí)際情況中這是無法預(yù)知的，當(dāng)即y服從非中心參數(shù)為x，自由度為K時，Luisier等人根據(jù)非中心卡方分布的統(tǒng)計特性，推導(dǎo)出了MSE的無偏估計一卡方無偏估計(CURE)，只要f(y)滿足連續(xù)或是平滑的限制，就可以由f(y)計算MSE的無偏估計CURE，并將其推廣到小波域，提出未歸一化的哈爾小波變換下的CURE估計。本方法只考慮未歸一化哈爾小波變換下的CURE估計，未歸一化的哈爾小波的意思是基本哈爾小波變換的尺度函數(shù)小波函數(shù)而未歸一化的哈爾小波尺度函數(shù)[1，1]，小波函數(shù)[1，-1]，這是為了使小波系數(shù)仍然滿足非中心卡方分布的性質(zhì)，在小波反變換時，重建濾波器也要做相應(yīng)的改變。在Luisier論文中的相關(guān)定理是這樣表述的，對圖像y做未歸一化哈爾小波變換，得到第j層的尺度系數(shù)s^j和小波系數(shù)w^j，用θ(w，s)＝θ^j(w^j，s^j)表示對無噪信號x的小波變換系數(shù)ω＝ω^j的估計，則第j層小波系數(shù)的MSE可以由下式估計：

式中K_j＝2^jK，K是非中心卡方分布的自由度，在磁共振儀器用單線圈采集成像時，K＝2。去噪函數(shù)取決于噪聲小波系數(shù)w、尺度系數(shù)s，而很多常用的去噪算法諸如軟閾值收縮、二變量收縮都是不考慮尺度系數(shù)的，在不考慮尺度系數(shù)的情況下，本發(fā)明推導(dǎo)出僅考慮小波系數(shù)w的閾值收縮方法的CURE估計，即用θ(w)代替θ(w，s)估計無噪信號的小波系數(shù)。

為簡化問題，我們把圖像展成一維的信號，噪聲信號y、理想信號x通過第j層小波分解所得尺度系數(shù)分別為s^j、高頻系數(shù)分別為則根據(jù)未歸一化哈爾小波變換的性質(zhì)可得：

式中K_j＝2^jK＝2^j+1。

第j層小波系數(shù)的MSE可表示為：

考慮第一層分解j＝1，噪聲信號y＝s^j-1，理想信號K＝K_j/2，則term2可表述為：

E{ω_nθ_n(w)}＝E{x_2nθ_n(w)}-E{x_2n-1θ_n(w)} (7)

Luisier已經(jīng)證明，當(dāng)f(y)連續(xù)可導(dǎo)時：

套用公式(6)求(5)：

代入term2：

套用公式(10)求term3：

最后將term1、term2、term3代入(8)可得僅取決于小波系數(shù)w的閾值收縮方法的小波域CURE估計：

其余層的分解亦可以得到類似的證明，而二維圖像的小波變換是對水平和垂直方向分別做小波變換，上述公式亦可直接推廣到二維圖像，所做的改動只是K_j＝2^2·jK。

2、如上述在推導(dǎo)出僅考慮小波系數(shù)w的閾值方法的小波域CURE估計之后，應(yīng)用到鄰域收縮去噪方法中，主要問題是閾值函數(shù)的導(dǎo)數(shù)計算以及閾值函數(shù)的平滑。

鄰域收縮法考慮當(dāng)前待處理的小波系數(shù)w_ij為中心的平方鄰域B_ij，給定平方鄰域系數(shù)和設(shè)λ為當(dāng)前小波子帶的全局閾值，為了表述方便，用w_n表示w_ij，S_n表示S_ij，則閾值收縮公式可表述為：

式中，θ(w_n)為去噪后的小波系數(shù)，用CURE估計來確定最佳的閾值λ，即λ_opt＝argminCURE(λ)，窗口大小取3*3，CURE估計需要計算θ(w_n)的一階和二階微分，令g＝max(u，0)，通過分層求導(dǎo)步驟如下：

CURE估計對去噪函數(shù)有平滑的要求，g＝max(u，0)須做平滑，利用不完全貝塔函數(shù)做平滑，設(shè)平滑的去噪函數(shù)為g＝f(u)＝u.*I(u)，其中

式中，a、b用于控制平滑的程度，為了保證平滑的曲線中心對稱，a、b須取值相同，a＝b＝80。

3、所述的循環(huán)移位(Cycle Spinning，CS)是為了消除由于平穩(wěn)小波變換缺乏平移不變性而產(chǎn)生的偽吉布斯現(xiàn)象。傳統(tǒng)的循環(huán)平移(Cycle Spinning)，即對含噪信號進(jìn)行“循環(huán)平移-閾值去噪-逆向循環(huán)平移”。由于對每次平移后的信號進(jìn)行閾值去噪會使偽吉布斯現(xiàn)象出現(xiàn)在不同的地方，因此針對行和列方向上的每組平移量都會得到一個不同的去噪結(jié)果μ_i，j，i、j分別為行和列方向上的平移量，對所有去噪結(jié)果進(jìn)行線性平均將得到抑制偽吉布斯現(xiàn)象的去噪結(jié)果。然而，傳統(tǒng)平移方法每次平移都要調(diào)用去噪函數(shù)，耗費(fèi)時間長，為了加快運(yùn)行速度，本發(fā)明考慮了只針對去噪后的小波系數(shù)的循環(huán)平移，這樣可以大大節(jié)省時間，因為平穩(wěn)小波變換每個子帶的尺寸和原圖相同，對圖像平移與對去噪系數(shù)的平移影響不大，尤其在圖像尺寸較大時，僅在邊緣處會有所不同，事實(shí)證明在平穩(wěn)小波變換下，這種快速循環(huán)平移方法和基本方法結(jié)果是一樣的，運(yùn)行時間卻大大減少。

附圖說明

下面結(jié)合附圖對本發(fā)明做進(jìn)一步的說明：

圖1基于卡方無偏估計的鄰域收縮MRI去噪方法流程圖；

圖2鄰域收縮法的MSE和CURE估計對比圖；

圖3傳統(tǒng)循環(huán)平移和快速循環(huán)平移示意圖；

圖4 T2w，7％噪聲第40切片噪聲、無噪和無噪圖像局部放大圖

圖5 T2w，7％噪聲第40切片去噪效果對比圖，包括去噪圖像，殘差圖像(噪聲圖減去去噪圖)和局部放大圖像。

圖6 T2w，7％噪聲第40到第100切片的去噪指標(biāo)對比

圖7 T2w第40切片模擬加不同強(qiáng)度Rician噪聲去噪指標(biāo)對比

具體實(shí)施方式

本發(fā)明提供了一種基于卡方無偏估計的鄰域收縮MRI去噪方法，當(dāng)需要對噪聲MRI進(jìn)行去噪處理時，如圖1所示依次執(zhí)行如下步驟：

步驟1 在MRI的背景區(qū)域估計噪聲標(biāo)準(zhǔn)差σ：

μ是選中的背景區(qū)域像素值的均值。

步驟2 對噪聲圖像m進(jìn)行平方，然后除以噪聲標(biāo)準(zhǔn)差σ的平方，得圖像y；

y＝m²/σ² (2)

步驟3 對y進(jìn)行未歸一化的平穩(wěn)哈爾小波變換，分解L層，得高頻系數(shù)和低頻系數(shù)，高頻系數(shù)有L*3個子帶；

步驟4 對低頻系數(shù)用雙邊濾波器去模糊；對每個子帶在預(yù)定的閾值搜索范圍內(nèi)用NeighShrinkCURE方法去噪，在閾值范圍內(nèi)用CURE估計來確定最佳的閾值λ，即λ_opt＝arg min CURE(λ)，再用進(jìn)行去噪；

步驟5 對去噪后的小波系數(shù)連同低頻系數(shù)進(jìn)行循環(huán)移位，CS＝9，即分別向右移i位，向下移j位，i∈[0，2]，j∈[0，2]，進(jìn)行小波反變換，再反平移相應(yīng)的位數(shù)得

步驟6 由得

步驟7 對得到的求平均的最終的去噪圖像

下面通過三組實(shí)驗來詳細(xì)分析該方法的性能：

從MRI數(shù)據(jù)庫(http：//brainweb.bic.mni.mcgill.ca/brainweb/)下載T1加權(quán)(T1w)、T2加權(quán)(T2w)，和質(zhì)子密度加權(quán)(PD)的不同噪聲等級的8比特MRI，原始數(shù)據(jù)為三維立體的人腦模擬MRI掃描數(shù)據(jù)，分辨率為181*217*181，將本發(fā)明的方法(NeighShrinkCURE+BF/CS＝9)和Anand等人的雙邊濾波器結(jié)合基于Stein無偏估計的鄰域收縮法(NeighShrinkSURE+BF)、Luisier等人的CURE-LET haar/CS＝16方法做對比。

實(shí)驗1：取T1w，5％噪聲、PD，7％噪聲、T2w，7％噪聲的MRI的第40個切片進(jìn)行去噪，該切片含有較完整的人腦邊緣和細(xì)節(jié)信息，用上述三種算法進(jìn)行去噪，得去噪指標(biāo)峰值信噪比(PSNR)和結(jié)構(gòu)相似度(SSIM)對比如表1所示。其中T2w，7％的噪聲圖像如圖3所示，去噪效果對比圖像如圖4所示。

表1不同去噪方法指標(biāo)對比

實(shí)驗2：對T2w，7％噪聲的MRI的第40到第100個切片用上述三種算法分別進(jìn)行去噪，去噪指標(biāo)PSNR和SSIM對比如圖5所示。

實(shí)驗3：對人為加噪聲的T2w，第40切片進(jìn)行試驗，將噪聲等級從3％逐漸增大到15％，間隔2％，得到不同去噪方法的去噪指標(biāo)對比。如圖6所示

通過上述三組實(shí)驗說明本發(fā)明算法能夠很好的實(shí)現(xiàn)圖像去噪，通過對比，證明本發(fā)明的算法基本能得到最高的PSNR和SSIM，視覺效果好，且魯棒性強(qiáng)，適用不同噪聲級的MRI去噪，能清晰地保留MRI的邊緣和細(xì)節(jié)信息。