Rician噪聲水平場(chǎng)的估計(jì)及其在MR圖像去噪中的應(yīng)用

余麗玲 陽(yáng) 維* 馮衍秋 劉 閩 馮前進(jìn) 陳武凡

1(南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，廣州 510515)

2(深圳出入境檢驗(yàn)檢疫局工業(yè)品檢測(cè)技術(shù)中心，深圳 518067)

引言

MRI技術(shù)提供了豐富的、有價(jià)值的人體結(jié)構(gòu)信息。然而，MR圖像在采集的過(guò)程中經(jīng)常被隨機(jī)噪聲干擾，MR圖像中的噪聲一般認(rèn)為服從Rician分布，是一種與信號(hào)相關(guān)的非加性噪聲。此外，隨著多線圈并行采集技術(shù)(敏感編碼、廣義自動(dòng)校準(zhǔn)部分并行采集)的采用［1］，MR圖像中的Rician噪聲水平在空間上可能不均勻。噪聲不僅會(huì)影響醫(yī)生的觀察，還會(huì)影響后續(xù)處理和分析的效果，如圖像分割、配準(zhǔn)、可視化等，抑制MR圖像中的Rician噪聲是后續(xù)處理和分析的重要步驟。

現(xiàn)有的Rician噪聲抑制方法，通常假設(shè)噪聲水平在空間上為常數(shù)。經(jīng)過(guò)適當(dāng)?shù)母脑?，非局部均值?］和 BM3D［3］算法，已被研究人員應(yīng)用于 Rician噪聲的抑制。其中，F(xiàn)oi等提出了利用方差穩(wěn)定變換(variance－stabilization transformations，VST)進(jìn)行Rician噪聲抑制的方法［4］，利用估計(jì)的全局噪聲水平，對(duì)MR圖像幅值進(jìn)行變換，使得變換后圖像中的噪聲與信號(hào)是獨(dú)立的，進(jìn)而采用BM3D算法進(jìn)行噪聲抑制，然而該方法僅適用于噪聲水平為常數(shù)的Rician噪聲圖像。Manjón等提出了依據(jù)局部Rician噪聲水平進(jìn)行自適應(yīng)調(diào)整的非局部均值［5］(adaptive non－local means，ANLM)算法，用于抑制 MR圖像中空間變化的Rician噪聲，取得了較好的去噪效果，但該方法較為耗時(shí)，且需調(diào)整的參數(shù)較多。

針對(duì)空間變化的Rician噪聲，本研究的思路是:對(duì)空間變化的噪聲水平建模并進(jìn)行估計(jì)，然后對(duì)圖像各處的幅值依據(jù)不同的局部噪聲水平進(jìn)行方差穩(wěn)定變換，進(jìn)而結(jié)合有效的去噪算法抑制這種空間變化的Rician噪聲。對(duì)于研究目的而言，其關(guān)鍵是如何有效估計(jì)MR圖像中不同位置的噪聲水平。

MR圖像中Rician噪聲水平的估計(jì)，較為常用的是基于背景區(qū)域的方法［6］。假設(shè)背景區(qū)域的信號(hào)值為零，采用最大似然或者幅值圖像的二階矩進(jìn)行估計(jì)Rician噪聲的水平，但不適用于空間變化噪聲水平的估計(jì)。本研究提出了一種基于稀疏性約束的Rician噪聲水平場(chǎng)的估計(jì)方法和空間自適應(yīng)方差穩(wěn)定變換方法，結(jié)合BM3D算法，有效地實(shí)現(xiàn)了對(duì)空間變化Rician噪聲的抑制。

1 方差穩(wěn)定變換

1.1 Rician噪聲分布

MR成像過(guò)程中，原始K空間數(shù)據(jù)被復(fù)數(shù)高斯噪聲干擾。由于Fourier變換的正交性和線性，經(jīng)Fourier逆變換后圖像實(shí)部和虛部的噪聲仍呈高斯分布。這樣，MR幅值圖像為兩個(gè)獨(dú)立高斯隨機(jī)變量的平方根，其噪聲不再服從高斯分布，而是服從Rician分布。設(shè)M為MR幅值，則M服從分布［2］為

式中，S為無(wú)噪聲時(shí)的幅值，σ為實(shí)圖像和虛圖像中高斯噪聲的標(biāo)準(zhǔn)偏差，I0為零階第一類(lèi)貝塞爾函數(shù)。Rician噪聲是信號(hào)相關(guān)的，既不是加性的也不是乘性的，針對(duì)加性高斯白噪聲的去噪方法不能直接應(yīng)用于MR圖像的去噪。

Rician噪聲抑制的一種可行方案，是利用Rician分布的最大似然在最大后驗(yàn)(maximum a posteriori，MAP)框架下對(duì) MR圖像進(jìn)行恢復(fù)，但MAP框架一般僅可融入有限的先驗(yàn)，而且其優(yōu)化過(guò)程相對(duì)耗時(shí)。另一種可行方案，是先將圖像進(jìn)行變換，使得在變換域中噪聲與信號(hào)獨(dú)立，在變換域中對(duì)去噪后圖像進(jìn)行逆變換，如小波變換、方差穩(wěn)定變換。本研究采用方差穩(wěn)定變換。方差穩(wěn)定變換可使圖像中的噪聲與信號(hào)近似不相關(guān)［4］，并且可通過(guò)相應(yīng)的逆變換得到信號(hào)的無(wú)偏估計(jì)。

1.2 方差穩(wěn)定變換

方差穩(wěn)定變換［4，7］是為了使異方差數(shù)據(jù)轉(zhuǎn)換為同方差的，更容易地運(yùn)用標(biāo)準(zhǔn)方法來(lái)解決問(wèn)題，去除Rician噪聲方差和無(wú)噪聲圖像幅值的依賴性。

圖1為所采用的方差穩(wěn)定變換及其對(duì)應(yīng)的逆變換，其中f為方差穩(wěn)定變換，D為去噪后逆變換之前的圖像，Vf為逆變換。對(duì)MR幅值圖像進(jìn)行方差穩(wěn)定變換處理后，使得σf(θ)=std{f(M)|S，σ}=1，圖像各處噪聲的方差基本一致，這樣可將VST變換后的Rician噪聲作為加性噪聲進(jìn)行抑制處理。對(duì)方差穩(wěn)定變換圖像進(jìn)行噪聲抑制后，需經(jīng)過(guò)對(duì)應(yīng)的方差穩(wěn)定逆變換，才能得到無(wú)偏的去噪圖像。

圖1 方差穩(wěn)定變換及其逆變換。(a)方差穩(wěn)定變換;(b)方差穩(wěn)定逆變換Fig.1 VST and inverse VST.(a)VST;(b)Inverse VST

2 噪聲水平場(chǎng)估計(jì)

對(duì)于水平空間變化的Rician噪聲，依據(jù)局部的噪聲水平逐像素進(jìn)行VST，即實(shí)現(xiàn)了空間自適應(yīng)VST?？臻g各處不同的噪聲水平在空間上構(gòu)成一個(gè)場(chǎng)，稱(chēng)之為噪聲水平場(chǎng)(noise level field，NLF)。顯然，進(jìn)行空間自適應(yīng)VST的關(guān)鍵在于NLF的準(zhǔn)確估計(jì)。利用NLF進(jìn)行空間自適應(yīng)VST和去噪的流程如圖2所示。主要包括4個(gè)步驟:估計(jì)NLF;利用NLF進(jìn)行空間自適應(yīng) VST;對(duì)變換的圖像使用BM3D算法進(jìn)行去噪;對(duì)去噪后圖像進(jìn)行VST逆變換。

圖2 空間自適應(yīng)VST和去噪的流程圖Fig.2 The flow chart of spatially adaptive VST and denosing

2.1 噪聲水平場(chǎng)模型

采用多項(xiàng)式描述噪聲水平場(chǎng)［8］

式中，{at，l}為多項(xiàng)式系數(shù)，xi，yi為像素點(diǎn) i的坐標(biāo)，K為多項(xiàng)式階數(shù)。為了方便計(jì)算，{at，l}用列向量A表示，{}用列向量C表示，B表示估計(jì)的噪聲水平場(chǎng)。式(2)可以擴(kuò)展到三維的NLF。

式中，λ為設(shè)定的正則化系數(shù)，B'為噪聲水平的局部估計(jì)值，其估計(jì)方法將在下一小節(jié)做詳細(xì)介紹。式(3)為l1正則化約束優(yōu)化問(wèn)題，解決此類(lèi)問(wèn)題的優(yōu)化方法很多［9］，采用 L1_LS 算法［10］求解多項(xiàng)式系數(shù)A，進(jìn)而由式(2)得到NLF。L1_LS算法結(jié)合了截?cái)嗯ｎD內(nèi)點(diǎn)法、預(yù)條件共軛梯度算法，比普通的內(nèi)點(diǎn)法使用方向或共軛梯度法計(jì)算效率要高，求解時(shí)間要短。

2.2 局部Rician噪聲水平的估計(jì)方法

由式(3)可見(jiàn)，為了估計(jì)參數(shù)A，需要得到噪聲水平的局部估計(jì)。傳統(tǒng)的噪聲水平估計(jì)方法一般假設(shè)噪聲水平是空間不變的常數(shù)。對(duì)于MR圖像中的Rician噪聲，較為常用的是基于背景區(qū)域的估計(jì)方法，但這些方法對(duì)于空間變化的噪聲水平不再適用。假設(shè)圖像局部的噪聲水平為常數(shù)，采用修正的中值絕對(duì)偏差(median absolute deviation，MAD)估計(jì)方法估計(jì)圖像局部的噪聲水平［11］。

假設(shè)噪聲服從高斯分布，對(duì)局部噪聲水平進(jìn)行MAD 估計(jì)［12－13］

式中，i，j表示圖像中的像素點(diǎn)，Ni，Nj表示點(diǎn) i，點(diǎn) j的鄰域，S(i)，S(j)是像素點(diǎn)的信號(hào)幅值，median(S)為中值濾波，為估計(jì)的噪聲標(biāo)準(zhǔn)偏差。與其它噪聲方差估計(jì)方法相比，MAD估計(jì)不但具有較強(qiáng)的魯棒性，而且運(yùn)算量較小。當(dāng)信噪比足夠大(＞5)時(shí)，Rician分布的標(biāo)準(zhǔn)偏差估計(jì)與高斯分布的趨于一致，無(wú)須進(jìn)行修正;當(dāng)信噪比較小時(shí)，式(4)的估計(jì)將會(huì)產(chǎn)生一定的偏差，須對(duì)估計(jì)的進(jìn)行修正。這種修正方法是基于對(duì)Rician噪聲信噪比(signal noise ratio，SNR)的迭代估計(jì)［14］。用 MAD估計(jì)出的標(biāo)準(zhǔn)偏差初始化修正過(guò)程為

式中，θ為SNR值，ζ(θ)為修正因子，定義如下

式中，I0，I1分別是零階、一階修正的貝塞爾函數(shù)。修正因子ζ(θ)通過(guò)迭代估計(jì)直到收斂或者達(dá)到給定的迭代次數(shù)。使用|θt－θt－1|作為收斂條件，迭代過(guò)程可表示為

式中，〈S〉是所給像素信號(hào)幅值的平均值，t為迭代次數(shù)。在梯度大的像素點(diǎn)上估計(jì)的局部噪聲水平與實(shí)際噪聲水平相比，存在著較大偏差。因此，在式(3)估計(jì)噪聲水平場(chǎng)時(shí)，不使用梯度大和修正后信噪比仍然較小的局部噪聲水平估計(jì)。

2.3 空間自適應(yīng)方差穩(wěn)定變換去噪

利用估計(jì)的NLF對(duì)圖像中各像素進(jìn)行方差穩(wěn)定變換

式中，β為可調(diào)節(jié)參數(shù)，B為估計(jì)的噪聲水平場(chǎng)，i，j表示圖像的位置。式(8)實(shí)際上相當(dāng)于對(duì)MR圖像進(jìn)行同態(tài)化和歸一化處理，使得整幅圖像的噪聲水平處處近似為1，噪聲與MR幅值和空間位置近似獨(dú)立。這樣，可采用針對(duì)加性高斯白噪聲的去噪方法進(jìn)行處理，然后進(jìn)行方差穩(wěn)定逆變換，最終得到無(wú)偏的去噪圖像。采用BM3D(BM4D)算法對(duì)方差穩(wěn)定變換后的圖像進(jìn)行去噪。BM3D是當(dāng)前公認(rèn)對(duì)加性高斯白噪聲去噪性能良好的算法，它包含了非局部去噪的思想，同時(shí)又用到了變換域?yàn)V波的方法。通過(guò)空間自適應(yīng)VST，可有效利用BM3D的去噪能力，實(shí)現(xiàn)空間變化Rican噪聲的抑制。

3 實(shí)驗(yàn)

為驗(yàn)證所提出方法的有效性，進(jìn)行仿真MR圖像實(shí)驗(yàn)和真實(shí)MR圖像實(shí)驗(yàn)，并與VST－BM3D、VSTBM4D和ANLM算法進(jìn)行比較。為方便起見(jiàn)，本方法簡(jiǎn)記為SAVST (spatially adaptivevariancestabilization transformations)。

3.1 仿真實(shí)驗(yàn)

仿真實(shí)驗(yàn)中，多項(xiàng)式的階數(shù)K設(shè)定為5，正則化系數(shù)λ為1，β為1.2。仿真的噪聲水平場(chǎng)由離散點(diǎn)三次方插值生成［5］，仿真的MR噪聲圖像的峰值信噪比(Peak signal to noise ratio，PSNR)為24 dB。

式中，RMSE為原始無(wú)噪聲圖像與去噪后圖像的均方差根誤差。

采用本方法估計(jì)噪聲水平場(chǎng)如圖3所示。由圖3可以看到，估計(jì)的噪聲水平場(chǎng)較好地逼近實(shí)際的噪聲水平場(chǎng)。噪聲水平場(chǎng)估計(jì)的精度采用均方根誤差(root mean square error，RMSE)和平均相對(duì)誤差(mean relative error，MRE)進(jìn)行度量。表1列出對(duì)于3種不同的仿真噪聲水平場(chǎng)的估計(jì)精度，平均相對(duì)誤差小于0.2%。

圖3 仿真噪聲水平場(chǎng)的估計(jì)結(jié)果。(a)噪聲圖像;(b)仿真的噪聲水平場(chǎng);(c)估計(jì)的噪聲水平場(chǎng)Fig.3 The estimated results of simulated NLF.(a)The noisy image;(b)The simulated NLF;(c)The estimated NLF

表1 噪聲水平場(chǎng)估計(jì)的精度Tab.1 The precision of NLF estimation

仿真實(shí)驗(yàn)圖像數(shù)據(jù)來(lái)自BrainWeb。對(duì)T1加權(quán)圖像、T2加權(quán)圖像，添加仿真 Rician噪聲水平從1%到15%進(jìn)行實(shí)驗(yàn)，仿真的NLF如圖4(b)所示。采用PSNR評(píng)價(jià)去噪效果。

圖4 不同方法噪聲估計(jì)結(jié)果。(a)原始圖像;(b)仿真噪聲水平場(chǎng);(c)噪聲圖像;(d)仿真的噪聲;(e)VSTBM3D的去噪效果;(f)VST-BM3D估計(jì)的噪聲;(g)本方法的去噪效果;(h)本方法估計(jì)的噪聲Fig.4 The estimated results of different methods.(a)Original image;(b)Simulated NLF;(c)Noisy image;(d)Simulated noise;(e)Denoised result of VST-BM3D;(f)Estimated noise of VST-BM3D;(g)Denoised result of the proposed method;(h)Estimated noise of the proposed method

圖4顯示了仿真噪聲水平為3%時(shí)，兩組圖像使用不同方法得到的去噪結(jié)果。其中，VST－BM3D使用噪聲水平為常數(shù)的方差穩(wěn)定變換去噪方法，圖4第一行為T(mén)1加權(quán)圖像的仿真實(shí)驗(yàn)結(jié)果，第二行為T(mén)2加權(quán)圖像的仿真實(shí)驗(yàn)結(jié)果。由圖4可見(jiàn)，本方法對(duì)于噪聲水平在空間上變化的Rician噪聲抑制效果較好，產(chǎn)生的模糊較小，并能保持原圖像的細(xì)節(jié)特征和邊緣信息。而VST－BM3D在信噪比較低的區(qū)域(如圖4(e)腦室區(qū)域)，噪聲未被有效抑制。對(duì)T1加權(quán)圖像、T2加權(quán)圖像采用本文方法(SAVST)和VST－BM3D方法進(jìn)行去噪實(shí)驗(yàn)，兩種去噪方法在不同噪聲水平下，對(duì)應(yīng)的去噪后PSNR值如圖5所示。其中，橫軸表示加噪水平，從1% ～15%，縱軸表示去噪后圖像的PSNR值。由圖可知，與VST－BM3D方法相比，本方法對(duì)應(yīng)的PSNR約提高了2 dB。

圖5 二維圖像去噪性能比較。(a)T1圖像;(b)T2圖像Fig.5 Comparison of denoising performance using different methods.(a)T1 image;(b)T2 image

本方法可應(yīng)用在三維MR圖像上，對(duì)同一噪聲水平場(chǎng)的三維腦部 MR圖像，采用 VST－BM4D、ANLM和本方法進(jìn)行去噪實(shí)驗(yàn)。圖6為在不同噪聲水平下，幾種方法去噪后PSNR的比較。從圖5和圖6可看出，本方法對(duì)噪聲水平空間變化的二維和三維MR圖像進(jìn)行去噪，對(duì)應(yīng)的PSNR明顯高于VST－BM3D方法。ANLM結(jié)合了噪聲水平的局部估計(jì)，其去噪性能也略高于VST－BM3D方法，但相對(duì)耗時(shí)。對(duì)于圖像大小為181像素×217像素×181像素的腦部數(shù)據(jù)，本方法所需的運(yùn)行時(shí)間約為340 s，而 VST－BM3D 方法為 378 s。

圖6 三維圖像去噪性能比較。(a)T1圖像;(b)T2圖像Fig.6 Comparison of denoising performance using different method.(a)T1 image;(b)T2 image

3.2 真實(shí)MR圖像實(shí)驗(yàn)

所用真實(shí)的乳腺M(fèi)R圖像，圖像大小為384像素×384像素，像素大小為1 mm×1 mm。真實(shí)乳腺M(fèi)R圖像的噪聲水平場(chǎng)估計(jì)見(jiàn)圖7。

圖7 乳腺M(fèi)R圖像噪聲水平場(chǎng)估計(jì)。(a)噪聲圖像;(b)估計(jì)的噪聲水平場(chǎng)Fig.7 The experimental results of breast MR/NLF estimation.(a)Noisy image;(b)Estimated NLF

由于真實(shí)噪聲圖像無(wú)對(duì)應(yīng)的參考圖像，為了評(píng)價(jià)去噪效果，采用Q度量［15］評(píng)價(jià)去噪圖像的圖像質(zhì)量。當(dāng)參考圖像不能獲得的時(shí)候，Q度量可較好地評(píng)價(jià)圖像質(zhì)量。Q度量值越大，圖像質(zhì)量越好。為了使去噪效果達(dá)到最優(yōu)，利用Q度量調(diào)節(jié)式(8)中的參數(shù)β。圖8顯示了對(duì)兩幅乳腺M(fèi)R圖像使用不同方法進(jìn)行去噪的結(jié)果，每一行代表一幅含噪聲乳腺M(fèi)R圖像的實(shí)驗(yàn)結(jié)果，相應(yīng)的Q值如圖8(f)所示，可看出當(dāng)參數(shù)β在1.2左右時(shí)，Q值達(dá)到最大值，去噪效果趨于最優(yōu)。本去噪方法對(duì)應(yīng)的Q值高于VST－BM3D，而且對(duì)于細(xì)節(jié)信息和邊沿結(jié)構(gòu)有更好的保存能力。

4 討論

本研究要解決的問(wèn)題是抑制MR圖像中空間變化的Rician噪聲。提出的方法是通過(guò)估計(jì)Rician噪聲水平場(chǎng)和方差穩(wěn)定變換，將Rician噪聲轉(zhuǎn)化為常數(shù)水平的加性噪聲，然后使用一般的去噪算法實(shí)現(xiàn)對(duì)噪聲的抑制。本文選取的去噪算法為BM3D算法(因其性能良好)，其他針對(duì)加性高斯白噪聲的去噪算法，如非局部均值算法、高斯混合尺度模型(scale mixtures of Gaussians，GSM)算法［16］也可以替代BM3D算法。

圖8 兩幅真實(shí)MR圖像的處理結(jié)果(上行為一幅，下行為另一幅)。(a)噪聲圖像;(b)VST-BM3D的去噪效果;(c)(a)與(b)之間的殘差圖像;(d)本方法的去噪效果，(e)(a)與(d)之間的殘差圖像;(f)調(diào)節(jié)β和Q值的變化Fig.8 The experimental results of two real MR images(The top row concerns one image，the bottom row concerns the other one).(a)Noisy image;(b)Result of VST-BM3D;(c)Residual image between(a)and(b);(d)Result of our proposed method;(e)Residual image between(a)and(d);(f)The corresponding Q values with varying β

噪聲水平場(chǎng)的準(zhǔn)確估計(jì)對(duì)于去噪效果有很大影響。文中，局部噪聲水平使用了修正的MAD方法進(jìn)行估計(jì)。MAD估計(jì)對(duì)服從高斯分布的噪聲具有較高的精度，但當(dāng)信噪比較低時(shí)，Rician分布不服從高斯分布，需要對(duì)估計(jì)的噪聲標(biāo)準(zhǔn)方差進(jìn)行修正，才能得到無(wú)偏的噪聲標(biāo)準(zhǔn)方差，從而提高了噪聲水平場(chǎng)的估計(jì)精度。估計(jì)噪聲水平場(chǎng)時(shí)，相關(guān)參數(shù)的設(shè)置對(duì)估計(jì)精度有較大影響，當(dāng)多項(xiàng)式階數(shù)K較小時(shí)，會(huì)得到偏低的估計(jì)，當(dāng)K較大時(shí)，會(huì)得到偏高的估計(jì);而正則化系數(shù)λ大小決定了噪聲水平場(chǎng)的平滑程度，需要依據(jù)經(jīng)驗(yàn)或先驗(yàn)知識(shí)進(jìn)行調(diào)節(jié)，對(duì)于實(shí)際問(wèn)題，可通過(guò)優(yōu)化Q度量確定最優(yōu)的λ。

5 結(jié)論

針對(duì)MR圖像中噪聲水平空間變化的Rician噪聲，本研究提出了一種噪聲水平場(chǎng)的估計(jì)方法，并用于抑制MR圖像中空間變化的Rician噪聲。通過(guò)Rician噪聲水平的局部估計(jì)和稀疏性約束模型估計(jì)噪聲水平場(chǎng)，然后對(duì)圖像各處的幅值依據(jù)不同的局部噪聲水平進(jìn)行方差穩(wěn)定變換，使得噪聲與信號(hào)幅值和空間位置無(wú)關(guān)，然后利用BM3D算法的強(qiáng)大去噪能力，實(shí)現(xiàn)了對(duì)空間變化Rician噪聲的有效抑制。實(shí)驗(yàn)結(jié)果表明，本方法能有效估計(jì)MR圖像中的Rician噪聲水平場(chǎng)，與使用全局噪聲水平估計(jì)的方差穩(wěn)定變換去噪方法、自適應(yīng)非局部均值去噪算法相比，本方法對(duì)應(yīng)的峰值信噪比和 Q度量均較高。

［1］陳武凡.并行磁共振成像的回顧、現(xiàn)狀與發(fā)展前景［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2005，24(6):649－654.

［2］Buades A，Coll B，Morel JM.A review of image denoising algorithm，with a new one［J］.SIAM Multiscale Modeling and Simulation，2005，4(2):490－530.

［3］Dabov K，F(xiàn)oi A，Katkovnik V，at al.Image denoising by sparse3D transform－domain collaborative filtering ［J］.IEEE Trans Image Process，2007，16(8):2080－2095.

［4］Foi A.Noise estimation and removal in MR imaging:the variance－stabilization approach ［C］// IEEE International Symposium on Biomedical Imaging.Chicago:IEEE，2011:1809－1814.

［5］Manjón JV，Coupé P，Martí－ Bonmatí L，et al.Adaptive nonlocal means denoising of MR images with spatially varying noise levels［J］.Magn Reson Imaging，2010，31(1):192 －203.

［6］Sijbers J，Poot D，Dekker AJ，et al.Automatic estimation of the noise variance from the histogram of a magnetic resonance image［J］.Phys Med Biol，2007，52(5):1335 －1348.

［7］FoiA. Optimization ofvariance－stabilizing transformations［OL］.http://www.cs.tut.fi/～foi/.

［8］Zheng Yuanjie.Estimation of image bias field with sparsity constraints［C］//Computer Vision and Pattern Recognition(CVPR).San Francisco:IEEE Computer Society，2010:255 －262.

［9］王浩.求解正則化問(wèn)題的算法研究［D］.北京:北京航空航天大學(xué)，2010.

［10］Kim SJ，Koh K，Lustig M，at al.An interior－point method for large－scale l1－regularized leastsquares ［J］.IEEE Signal Processing，2007，1(4):606 －617.

［11］Rousseeuw PJ，Verboven S.Robust estimation in very small samples［J］.Comp Stat Data Anal，2002，40(4):741 －758.

［12］Maximov II，F(xiàn)arrher E，Grinberg F，at al.Spatially variable Rician noise in magnetic resonance imaging［J］.Med Image Anal，2012，16(2):536 －548.

［13］Coupé P，Manjón JV，Gedamu E，at al.Robust Rician noise estimation for MR images［J］.Med Image Anal，2010，14(4):483－493.

［14］Koay CG，Basser PJ.Analytically exact correction scheme for signal extraction from noisy magnitude MR signals［J］.J Magn Reson，2006，179(2):317－322.

［15］Zhu Xiang，Milanfar P. Automatic parameter selection fordenoising algorithms using a no－reference measure of image content［J］.IEEE Trans Image Process，2010，19(12):3116 －3132.

［16］Portilla J，Strela V，Wainwright M，et al.Image senoising using scale mixtures of gaussians in the wavelet domain ［J］.IEEE Trans on Image Process，2003，12(11):1338－1351.