一種在擴(kuò)散加權(quán)圖像降噪中的算法

易三莉 賀建峰 邵黨國 劉正剛

（1.云南昆明理工大學(xué)信息工程與自動化學(xué)院，昆明，650500；2.云南昆明理工大學(xué)外國語言文化學(xué)院，昆明，650500）

引 言

擴(kuò)散張量成像（Diffusion tensor imaging，DTI）技術(shù)的出現(xiàn)是醫(yī)學(xué)影像技術(shù)領(lǐng)域中一個非常大的進(jìn)展，它實(shí)現(xiàn)了對腦白質(zhì)纖維的無創(chuàng)成像。DTI是根據(jù)Stejskal-Tanner方程對擴(kuò)散加權(quán)圖像（Diffusion weighted imaging，DWI）的計算 得來［1，2］。由于大腦內(nèi)部具有豐富的白質(zhì)結(jié)構(gòu)，因而決定了擴(kuò)散張量圖像具有多邊界的特點(diǎn)，并且DWI圖像中的邊界信息對于大腦白質(zhì)結(jié)構(gòu)的反映尤其重要。因而要求對擴(kuò)散張量圖像進(jìn)行降噪的同時，能較好地保存圖像的邊界信息。針對DWI圖像多邊界，數(shù)據(jù)量大，低信噪比且噪聲為萊斯分布的特點(diǎn)，國內(nèi)外研究人員提出了許多對DWI進(jìn)行降噪的算法。如McGraw［3］等應(yīng)用全變差（Total variation，TV）濾波方法進(jìn)行擴(kuò)散張量圖像的降噪，Basu［4］等將基于貝葉斯理論的最大后驗(yàn)估計濾波器對擴(kuò)散張量圖像進(jìn)降噪，以及Fernandez［5］等人采用了各向異性維納濾波的方法，國內(nèi)張相芬［6］等人應(yīng)用向量復(fù)擴(kuò)散模型對DWI圖像進(jìn)行降噪等。

自適應(yīng)維納濾波器是一種經(jīng)典的線性降噪濾波器，其抗噪性能優(yōu)良并且計算簡單。然而，由于DWI圖像具有豐富的邊界信息，且其噪聲分布為萊斯噪聲，直接將自適應(yīng)維納濾波方法應(yīng)用于對DWI圖像的降噪并不能取得較好的效果。這是因?yàn)椋菏紫龋珼WI圖像的噪聲主要集中于圖像的高頻部分，圖像的低頻部分含有的噪聲較高頻部分而言則要小很多。自適應(yīng)維納濾波算法在對圖像進(jìn)行降噪的過程中，并不能分辨圖像中的不同頻率成份，因而只能對所有頻率成份都不加區(qū)分地進(jìn)行濾波降噪處理。二維經(jīng)驗(yàn)?zāi)B(tài)分解算法則可將圖像分解為表示不同頻率成份的細(xì)節(jié)圖像以及輪廓圖像，它是一種將Hilbert-Huang變換用于圖像等二維信號處理的方法［7］。其次，自適應(yīng)維納濾波器主要針對高斯噪聲具有較好的抗噪性能，對于萊斯噪聲則并不十分理想。通過使用改進(jìn)的維納濾波算法能夠使其更適用于對DWI圖像的降噪。本文提出了一種將二維經(jīng)驗(yàn)?zāi)B(tài)分解（Bidimensional empivical mode decomposition，BEMD）和改進(jìn)的維納濾波器相結(jié)合的算法，并將該算法應(yīng)用于DWI圖像的降噪。

1 BEMD與改進(jìn)的維納濾波法算法

1.1 二維經(jīng)驗(yàn)?zāi)B(tài)分解

二維經(jīng)驗(yàn)?zāi)B(tài)分解是在一維經(jīng)驗(yàn)?zāi)B(tài)分解的基礎(chǔ)上發(fā)展而來的。一維經(jīng)驗(yàn)?zāi)B(tài)分解算法（Empirical mode decomposition，EMD），最 初 是 由Huang等人于1998年提出的用于對復(fù)雜信號進(jìn)行平穩(wěn)化處理的信號處理方法［8］。它首先對一個時間序列進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，將信號分解為有限個內(nèi)稟模式函數(shù)（Intrinsic mode function，IMF），然后對各個分量IMF做希爾伯特變換從而得到其瞬時頻率與振幅。從本質(zhì)上來看，該方法的處理結(jié)果是將信號中不同尺度的波動或趨勢逐級分解開來。由于這種分解是基于局部特征尺度，因而它具有良好的局部自適應(yīng)性及多尺度的優(yōu)點(diǎn)。由于EMD在一維信號的處理尤其是非平穩(wěn)信號的處理上具有較好的效果并得到廣泛的應(yīng)用，近年來一些國內(nèi)外研究者將其推廣到二維圖像分析并且取得了一定的成果。如瑞典的Linderhed［9］在2005年提出了基于樣條插值的二維EMD方法，并將該方法用在圖像壓縮中；Damerval［10］提出了基于三角剖分和立方樣條插值的二維EMD，將它用于圖像降噪并取得較好的效果；Nunes［11］等人基于數(shù)學(xué)形態(tài)學(xué)重構(gòu)進(jìn)行圖像的BEMD分解，并將其用于紋理圖像及醫(yī)學(xué)圖像的分析。國內(nèi)葛光濤［12］等通過對BEMD算法引入新的篩分停止準(zhǔn)則實(shí)現(xiàn)更好的保留圖像細(xì)節(jié)。

對于二維圖像，采用BEMD進(jìn)行分解必須基于以下假設(shè)［13］：（1）二維數(shù)據(jù)平面至少包含一個極大值點(diǎn)和一個極小值點(diǎn)，或整個二維平面沒有極值點(diǎn)，但在進(jìn)行一階或多階求導(dǎo)運(yùn)算后能出現(xiàn)一個極大值點(diǎn)和一個極小值點(diǎn)；（2）特征尺度用極值點(diǎn)之間距離定義。

對一幅圖像f（x，y）進(jìn)行二維經(jīng)驗(yàn)?zāi)B(tài)分解，將它分解為有限個二維IMF以及最終的余項(xiàng)函數(shù)。其分解方法與一維處理基本相同，所得到的二維IMF具有如下特征：IMF的極值點(diǎn)數(shù)目等于過零點(diǎn)的數(shù)目，或者最多相差為1；IMF的局部極大值所構(gòu)成的上包絡(luò)面與局部極小值所構(gòu)成的下包絡(luò)面的均值曲面為零。BEMD算法的具體實(shí)現(xiàn)如下：

（1）賦初值，令i＝1，k＝1，分別表示第i個IMF以及為了計算IMFi所進(jìn)行的第k次循環(huán)，余項(xiàng)的初值為：rik（x，y）＝f（x，y）。

（2）計算rik（x，y）的極大值包絡(luò)hup（x，y）以及它的極小值包絡(luò)hlow（x，y）。

（3）計算rik（x，y）的均值曲面：mean（x，y）＝ （hup（x，y）＋hlow（x，y））/2。

（4） 根 據(jù)rik（x，y） 以 及 它 的 均 值 曲 面mean（x，y） 計 算 其 差 值 函 數(shù)：Dik（x，y） ＝rik（x，y）-mean（x，y）。

（5）對差值函數(shù)Dik（x，y）進(jìn)行判斷，看它是否滿足上文所述IMF函數(shù)的特征：當(dāng)它滿足時，就可得到內(nèi)稟模式函數(shù)：IMFi（x，y）＝Dik（x，y），令r（i＋1）1（x，y）＝ri1（x，y）-IMFi（x，y），并且令i＝i＋1，k＝1，重復(fù)（2）～（4）計算下一個內(nèi)稟模式函數(shù)；當(dāng)不滿足時，令ri（k＋1）（x，y）＝Dik（x，y），并且令k＝k＋1，重復(fù)（2）～（4）直到差值函數(shù)滿足內(nèi)稟模式函數(shù)為止。

綜上所述，對于圖像f（x，y）通過BEMD分解為n個內(nèi)稟模式函數(shù)IMFi及其余項(xiàng)函數(shù)r（x，y）可以用下式來表示

通過BEMD算法可將圖像分解為一系列細(xì)節(jié)信息和趨勢信息。上式中的內(nèi)稟模式函數(shù)IMFi即表示圖像的細(xì)節(jié)信息，對應(yīng)圖像中的邊緣、噪聲等高頻信息，i的值越小表示越早分解出來的內(nèi)稟模式函數(shù)，它所含有的頻率也就越高，余項(xiàng)函數(shù)r（x，y）則表示了圖像的基本結(jié)構(gòu)、基本變化的趨勢信息［8］。

1.2 改進(jìn)的自適應(yīng)維納濾波

由于擴(kuò)散加權(quán)圖像中的噪聲為萊斯噪聲，而維納濾波器是針對高斯噪聲進(jìn)行降噪的，直接用維納濾波器對擴(kuò)散加權(quán)圖像進(jìn)行降噪并不能得到理想的降噪效果［14］。萊斯校正技術(shù)在核磁共振圖像中應(yīng)用較為廣泛，對于含有萊斯噪聲的信號，當(dāng)信噪比較高時它趨向于高斯分布，因此可以在對擴(kuò)散加權(quán)圖像進(jìn)行降噪之前先通過萊斯校正技術(shù)對擴(kuò)散加權(quán)圖像中的數(shù)據(jù)進(jìn)行校正處理，從而使處理后的圖像可運(yùn)用維納濾波算法對其進(jìn)行降噪處理。

通過式（2）進(jìn)行萊斯校正之后所得到的圖像f則可通過自適應(yīng)維納濾波器對其進(jìn)行降噪處理。對于圖像f（x，y）（其中（x，y）表示像素點(diǎn)f（x，y）在圖像中的位置），自適應(yīng)維納濾波器就是尋找圖像的估計值（x，y）滿足它與圖像f（x，y）之間的均方誤差為最小。因而自適應(yīng)維納濾波器對該圖像的降噪模型表示如下

式中：g（x，y）是濾波結(jié)果，μ（x，y）為圖像鄰域均值，w是自適應(yīng)維納濾波器，其表達(dá)式為

1.3 BEMD與維納濾波的結(jié)合算法

由于DWI圖像是由不同的頻率成份組成的，圖像的高頻成份代表著圖像的細(xì)節(jié)信息，而圖像的低頻成份代表圖像的趨勢信息。對于一幅含噪的DWI圖像，其噪聲主要集中于它的高頻部分，而圖像的低頻部分含有的噪聲較高頻部分而言則要小很多。事實(shí)上，由于自適應(yīng)維納濾波算法對圖像中所有頻率成份都不加區(qū)分地進(jìn)行濾波降噪處理并不合理。對于圖像中不含噪或者含噪非常少的圖像成份而言，不僅不能改善圖像質(zhì)量，反而由于濾波器的過濾波作用使圖像質(zhì)量變差。本文通過將二維經(jīng)驗(yàn)?zāi)B(tài)分解和改進(jìn)的自適應(yīng)維納濾波相結(jié)合的方法對DWI圖像進(jìn)行降噪處理。

首先，使用上文所述的BEMD算法，將含噪圖像分解為四個不同頻率成份的子圖像IMF1，IMF2，IMF3及REF。它們分別對應(yīng)的頻率成份為從高到低：IMF1體現(xiàn)了原始含噪圖像的高頻信號部分，主要反映圖像的邊界等細(xì)節(jié)；第2個分量IMF2所包含高頻成份較分量IMF1次之，但仍然具有豐富的邊界信息；第3個分量IMF3則包含較多的低頻成份；圖像的余項(xiàng)部分即REF則體現(xiàn)了圖像中的低頻成份，代表圖像趨勢信息。

然后采用改進(jìn)的自適應(yīng)維納濾波器對不同子圖像進(jìn)行處理。事實(shí)上，由于噪聲主要集中于圖像中的高頻部分，因而對IMF1和IMF2采用改進(jìn)的自適應(yīng)維納濾波算法進(jìn)行降噪，從而得到更準(zhǔn)確的邊界信息。而對于IMF3和REF，它們主要體現(xiàn)的是圖像的趨勢信息，因而保持原始圖像不變。

最后，將處理后的對應(yīng)各頻率成份子圖像進(jìn)行相加得到降噪后的圖像。

2 實(shí)驗(yàn)結(jié)果與分析

將本文所提出的BEMD與改進(jìn)的維納濾波器相結(jié)合的算法用于DWI圖像的降噪，來對該算法的降噪性能進(jìn)行分析比較。實(shí)驗(yàn)中采用GE1.5T磁共振系統(tǒng)對人腦進(jìn)行DWI成像，掃描參數(shù)為：TR＝12 000ms，TE＝107ms，圖像大小為256×256，體素大小為1mm×1mm×4mm，梯度編碼方向數(shù)為13，b＝1 000s/mm2，其參考DWI圖像為DWI0每一層DWI數(shù)據(jù)集中包括一幅未加權(quán)的參考圖像DWI0，以及13幅擴(kuò)散加權(quán)的圖像DWI1～13，共計14幅DWI圖像。本文通過對編碼方向?yàn)?的圖像DWI0進(jìn)行相關(guān)的實(shí)驗(yàn)分析。

2.1 DWI圖像的BEMD分析

使用BEMD方法將DWI圖像分解為3個IMF圖像及一個余項(xiàng)函數(shù)圖像，共計4個子圖像。分解后的結(jié)果如圖1所示，其中圖1（a～e）分別表示DWI0，IMF1，IMF2，IMF3以及余項(xiàng)函數(shù)。從圖1可以看到，IMF1幾乎包含了圖像的的所有細(xì)節(jié)信息，而IMF2所含的細(xì)節(jié)信息較IMF1要少，IMF3則包含更少的細(xì)節(jié)信息，而在余項(xiàng)圖像中幾乎看不到圖像的細(xì)節(jié)信息而只是圖像的趨勢信息。

圖1 DWI0及BEMD算法對其分解結(jié)果Fig.1 Decomposed results of DWI0of BEMD method

2.2 降噪效果比較

通過兩個實(shí)驗(yàn)來對本文所提出的降噪算法進(jìn)行分析比較。

首先，在第一個實(shí)驗(yàn)中，將萊斯噪聲加入到原始DWI圖像中，萊斯噪聲的大小設(shè)置為：σ2Ra＝0.05：0.05：0.25（其中σ2Ra為萊斯噪聲的方差），然后將自適應(yīng)維納濾波方法、改進(jìn)的維納濾波法以及本文所提出的基于BEMD算法的濾波方法分別應(yīng)用于含噪的DWI圖像的降噪中，其結(jié)果如表1所示。從表1中可以看到與BEMD相結(jié)合的算法相較于其他算法，隨著信噪比的增加，具有更高的峰值信噪比（Peak signal to noise ratio，PSNR），從而其降噪性能更好。

表1 降噪結(jié)果的PSNR比較Table 1 PSNR of denoising results of different methods

其次，在第二個實(shí)驗(yàn)中，將方差σ2Ra＝0.1的萊斯噪聲加入到原始的DWI0圖像中，并將上述三種算法應(yīng)用于對該含噪DWI圖像的降噪。降噪結(jié)果如圖2所示，其中圖2（a）為含噪的DWI圖像，圖2（b）為自應(yīng)用維納降噪圖像；圖2（c）為使用改進(jìn)的維納濾波進(jìn)行降噪的圖像；圖2（d）為使用本文所提出的BEMD與改進(jìn)的維納濾波相結(jié)合的算法所得到的降噪圖像。

圖2 不同降噪算法對含噪DWI0圖像的降噪Fig.2 Denoised results of DWI0of different denoising methods

從圖2中，可以看出圖2（d）的降噪效果最為理想，即BEMD與改進(jìn)的維納濾波器相結(jié)合的算法降噪后所得到的圖像最接近原始圖像。因此，可得出與上個實(shí)驗(yàn)相同的結(jié)論，即本文所提出的算法能夠更好地對DWI圖像進(jìn)行降處理。

3 結(jié)束語

本文針對DWI圖像的特點(diǎn)，提出了BEMD和改進(jìn)的維納濾波器相結(jié)合的降噪算法，并將該算法應(yīng)用于DWI圖像的降噪中。通過采用BEMD算法將具有豐富邊界信息的DWI圖像分解為含有不同頻率成份的內(nèi)稟模式函數(shù)IMF，得到分別代表圖像高頻成份以及代表圖像低頻成份的子圖像；根據(jù)噪聲主要集中于圖像高頻部分的特點(diǎn)，采用改進(jìn)的維納濾波算法對不同子圖像進(jìn)行不同的降噪處理；最后將處理后的子圖像相加得到最終的降噪圖像。在本文的實(shí)驗(yàn)中，將該算法應(yīng)用于DWI圖像的降噪中，并將其與維納濾波器的降噪結(jié)果進(jìn)行比較。通過比較，可以看出本文算法具有更峰值信噪比，即本文所提方法較維納濾波算法具有更好的降噪效果。

［1］Basser P J，Mattiello J，LeBihan D.MR diffusion tensor spectroscopy and imaging［J］.Biological Physics，1994，66：259-267.

［2］Stejskal E O，Tanner J E.Spin diffusion measurements：spin echoes in the presence of a time-dependent field gradient［J］.Journal of Chemical Physics，1965，42：288-292.

［3］McGraw T，Vermuri B C，Chen Y，et al.DT-MRI denoising and neuronal fiber tracking ［J］.Medical Image Analysis，2004，8：95-111.

［4］Saurav B，Thomas F，Ross W，Rician noise removal in diffusion tensor MR［J］.MICCAI，LNCS4190，2006：117-125.

［5］Martin-Fernandez M，Munoz-Moreno E，Cammoun L，et al.Sequential anisotropic multichannel Wiener filtering with Rician bias correction applied to 3Dregularization of DWI data［J］.Medical Image Analysis，2009，13：19-35.

［6］張相芬，田蔚風(fēng)，葉宏.DTI圖像恢復(fù)的向量復(fù)擴(kuò)散模型［J］.計算機(jī)工程與設(shè)計，2009，30（3）：634-635.

Zhang Xiangfen，Tian Weifeng，Ye Hong.Restoring DTI images by vector-valued complex diffusion model［J］.Computer Engineering and Design，2009，30（3）：634-635.

［7］Nunes J C，Bouaoune Y，Delechelle E，et al.Image analysis by bidimensional empirical mode decomposi-tion［J］.Image and Vision computing，2003，21：1019-1026.

［8］Huang N E，Shen Z，Long S，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proc R Soc Lond Ser A，1998，454：903-905.

［9］Linderhed A.Compression by image empirical mode decomposition［C］∥Image Processing，IEEE International Conference on Image Processing.Piscataway，USA：IEEE，2005：553-556.

［10］Damerval C，Sylvain M，Valerie P.A fast algorithm for bidimensional EMD ［J］.Signal Processing Letters，IEEE，2005，12（10）：701-704.

［11］Nunes J C，Bouaoune Y，Delechelle E.Image analysis by bidi-mensional empirical mode decomposition［J］.Image and Vision Computing，2003，21（12）：1019-1026.

［12］葛光濤，桑恩方，劉卓夫，等.一種新的BEMD篩分停止準(zhǔn)則［J］.數(shù)據(jù)采集與處理，2010，25（2）：195-200.

Ge Guangtao，Sang Enfang，Liu Zhoufu，et al.Novel BEMD criterion for stopping sifting process［J］.Journal of Data Acquisition and Processing，2010，25（2）：195-200.

［13］周欣，李衷怡.基于二維經(jīng)驗(yàn)?zāi)Ｊ椒纸獾膱D像去噪［J］.計算機(jī)與數(shù)字工程，2007，35（11）：93-94.

Zhou Xin，Li Zhongyi.Image denoise based on BEMD method［J］.Computer and Digital Engineering，2007，35（11）：93-94.

［14］Martin M F，Westin C F，Alberola C L.3DBayesian regularization of diffusion tensor MRI using multivariate Gaussian Markov random fields［J］.Lecture Notes in Computer Science，2004，3216：351-359.

［15］Matzner R.An SNR estimation algorithm for comoplex baseband signals using higher order statistics［J］.Facta Universitatis（Nis），1993，6：41-52.

［16］Koay C G，Basser P J.Analytically exact correction scheme for signal extraction from noisy magnitude MR signals ［J］.Journal of Magnetic Resonance，2006，179：317-322.