彌散加權磁共振圖像(DWI)是使用特殊自旋平面回波序列進行快速成像,它易被噪聲干擾,需要有效去噪以保證后續應用。目前去噪方法一般為常用圖像去噪方法擴展,缺乏針對DWI多不同梯度磁場方向數據構成特點的針對性應用。本文提出一種DWI Rician噪聲的線性最小均方誤差(LMMSE)復原方法,使用局部信息的統計特征,對DWI的Rician噪聲進行有效估計,并綜合多梯度磁場方向改進使用LMMSE進行復原。在合成模擬DWI數據和真實人體腦部DWI數據上進行的仿真和實驗表明,本文方法較之目前使用較多的逐梯度方向去噪方法能夠更好去除DWI中Rician噪聲,有效改善計算獲得的DTI大小和方向信息的有效性和準確性。
引用本文: 吳錫, 何晉, 朱明. 綜合多梯度磁場方向彌散加權磁共振圖像線性最小均方誤差去噪. 生物醫學工程學雜志, 2014, 31(1): 7-12. doi: 10.7507/1001-5515.20140002 復制
引言
彌散磁共振成像(diffusion magnetic resonance imaging,dMRI)是對活體腦部神經組織結構進行研究的非介入成像技術,通過測量組織中水分子的三維彌散對其生理機能進行定量分析和數字成像,并可利用纖維成像技術重建腦白質神經纖維的三維結構[1]。常用的彌散張量磁共振圖像(diffusion tensor magnetic resonance image,DT-MRI或者DTI)對三維組織的每一體素使用一個3×3對稱正定矩陣構成張量模型進行描述,該模型需要至少6個非共面擴散敏感梯度磁場方向的回波衰減信號測量值構成的彌散加權磁共振圖像(diffusion weighted magnetic resonance image,DW-MRI或者DWI)和一個不施加擴散敏感梯度磁場的磁共振信號參考測量值解得[2]。由于DWI成像使用特殊的自旋平面回波序列易受噪聲干擾[3-4],掃描獲得的DWI數據需要進行有效去噪處理以保證DTI數據計算的準確性和后續應用的可靠性[5]。
現有DWI去噪方法一般由經典灰度圖像去噪推廣而來,使用各向異性濾波[6-7]、極大似然估計[8]等逐個對不同梯度磁場方向DWI進行去噪。現有DWI去噪方法的缺陷主要有兩方面:首先,經典灰度圖像去噪一般針對噪聲模型為高斯噪聲,而DWI的噪聲一般滿足Rician分布,直接使用高斯噪聲去噪方法會產生誤差,需要進行補償[9];另一方面,由于DWI掃描特殊性,僅使用單一梯度磁場DWI去噪并不能獲得最優復原結果[10-11]。
本文使用線性最小均方誤差(linear minimum mean square error,LMMSE)復原去噪方法直接基于Rician噪聲進行去噪,克服常用基于高斯噪聲去噪算法缺陷,同時綜合多梯度方向信息獲得最優去噪結果。下面,首先對本文方法的基本原理進行分析,然后分別使用合成和真實DWI數據對本文方法的有效性進行驗證分析。
1 綜合多梯度磁場方向LMMSE去噪
不同于普通灰度圖像,DWI成像結果為復數,當該信號的實部和虛部同時被方差相等的非相關零均值高斯噪聲時,該信號的模服從Rician分布[3]。令M為測量數據,A為原始數據,則DWI數據集根據掃描矩陣的大小可視為P×Q×R的三維矩陣(P、Q、R分別為其行、列和層數),每一體素存儲一個包含7個(或以上)參數的列向量,滿足Rician分布的DWI數據及其概率密度函數表示為
${{M}_{i}}=\sqrt{{{\left( {{A}_{i}}+{{n}_{i.c}} \right)}^{2}}+n_{i.c}^{2}},$ |
${{P}_{M}}\left( {{M}_{i}}/{{A}_{i}},\sigma \right)=\frac{{{M}_{i}}}{{{\sigma }^{2}}}{{e}^{-\frac{M_{i}^{2}+A_{i}^{2}}{2{{\sigma }^{2}}}}}{{I}_{0}}\left( \frac{{{A}_{i}}{{M}_{i}}}{{{\sigma }^{2}}} \right)u\left( {{M}_{i}} \right),$ |
其中i為梯度磁場編號,ni,c和ni,s為非相關零均值高斯噪聲,方差均為σ2,其二階矩和四階矩為:
${{{{\mu }''}}_{2}}=E\left\{ M_{i}^{2} \right\}=2{{\sigma }^{2}}+A_{i}^{2}$ |
${{{{\mu }''}}_{4}}=E\left\{ M_{i}^{4} \right\}=8{{\sigma }^{4}}+8{{\sigma }^{2}}A_{i}^{2}+A_{i}^{4}$ |
LMMSE的原理即使用測量數據M獲得原始數據A的最優估計值A′,其方法由[12]可得:
${A}'=E\left\{ A \right\}+{{C}_{AM}}C_{AM}^{-1}\left( M-E\{M\} \right)$ |
由于DWI數據中的噪聲呈Rician分布,無法直接通過其一階矩計算,因此最優估計值A′只能通過其平方間接獲得,如下式
$A_{i}^{'2}=E\left\{ A_{i}^{2} \right\}+{{C}_{A_{i}^{2}M_{i}^{2}}}C_{M_{i}^{2}M_{i}^{2}}^{-1}\left( M_{i}^{2}-E\left\{ M_{i}^{2} \right\} \right),$ |
其中CAi2Mi2和CMi2Mi2分別為互協方差和自協方差,則DWI數據集中體素p的濾波器輸出為:
$A_{i}^{'2}\left( p \right)=E\left\{ A_{i}^{2}\left( p \right) \right\}+{{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)C_{M_{i}^{2}M_{i}^{2}}^{-1}\left( p \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right)$ |
CAi2Mi2和CMi2Mi2分別由測量值的矩通過下式獲得:
$\begin{align} & {{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ \left( A_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right) \right\} \\ & =E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}-4{{\sigma }^{2}}\left( {{E}^{2}}\left\{ {{M}_{i}}\left( p \right) \right\}-{{\sigma }^{2}} \right) \\ \end{align}$ |
$\begin{align} & {{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ {{\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right)}^{2}} \right\} \\ & =E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\} \\ \end{align}$ |
將式(8)、(9)帶入式(7),得
$A_{i}^{'2}\left( p \right)=E\left\{ M_{i}^{2}\left( p \right) \right\}-2{{\sigma }^{2}}+K\left( p \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right),$ |
其中,K(p)為:
$K\left( p \right)=1-\frac{4{{\sigma }^{2}}\left( E\left\{ M_{i}^{2}\left( p \right) \right\}-{{\sigma }^{2}} \right)}{E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}}$ |
由于測量數據M的矩無法直接獲得,因此需要對當前體素p的鄰域采樣進行估計:
$E\left\{ M_{i}^{r}\left( p \right) \right\}=\frac{1}{S}\sum\limits_{q\in S}{M_{i}^{r}\left( q \right)},$ |
其中S為進行采樣選擇的鄰域區域。
令Mi2和Ai2分別為[Mi-n2,…Mi2,…Mi+n2]T和[Ai-n2,…Ai2,…Ai+n2]T,其中n為常數,Mi2和Ai2表征以第i個梯度磁場為中心的多梯度磁場信息,其互協方差和自協方差分別為:
$\begin{align} & {{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ \left( A_{i}^{2}\left( p \right)-E\left\{ A_{i}^{2}\left( p \right) \right\} \right) \right.\left. \left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right) \right\} \\ & =E\left\{ M_{i}^{2}\left( p \right) \right\}=E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}- \\ & 4{{\sigma }^{2}}\left( {{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}-{{\sigma }^{2}} \right) \\ \end{align}$ |
$\begin{align} & {{C}_{M_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ M_{i}^{2}\left( p \right)-E{{\left\{ M_{i}^{2}\left( p \right) \right\}}^{2}} \right\} \\ & =E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\} \\ \end{align}$ |
則式(10)改寫為
$A_{i}^{'2}\left( p \right)=E\left\{ M_{i}^{2}\left( p \right) \right\}-2{{\sigma }^{2}}+K\left( p \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right),$ |
其中噪聲方差σ2根據Rician噪聲特點,對原始數據為0區域(圖像背景)估計獲得。由式(3) 可知,當A2=0時,有μ2″=E{Mi2}=2σ2,則
${{\sigma }^{2}}=\frac{1}{2N}\sum\limits_{i=1}^{N}{M_{i}^{2}},$ |
其中 N為估計所用背景區域大小。
2 仿真與實驗結果
在合成和真實腦部DWI數據集中,分別使用綜合多梯度磁場方向LMMSE、逐梯度磁場方向LMMSE和各向異性濾波[6]三種方法進行去噪復原比較。
首先,在合成DWI數據集上進行仿真,仿真平臺為MATLAB 7.0.4,Windows 7 Professional操作系統,Intel Core(TM) i7 CPU處理器和4G內存。合成DWI數據集大小為64×64×1,由6個不同方向的DWI數據集和1個標準方向DWI數據集使用最小二乘模型計算獲得,合成纖維結構的FA為0.8,彌散度為2.1×10-5 cm2/s,幾何結構為一條正弦曲線,寬度為10體素。
由圖 1可知,三種方法均能在一定程度上去除DWI噪聲,相比而言,基于Rician噪聲模型的LMMSE方法綜合數據的統計信息,較之局部灰度差異的各向異性濾波方法可以獲得更好的去噪效果,特別綜合多梯度磁場方向提供更多信息獲得最優去噪效果。

(a) 未加噪DWI;(b) 加入5% Rician噪聲;(c) 綜合多梯度磁場方向LMMSE;(d) 逐梯度磁場方向LMMSE;(e) 各向異性濾波
Figure1. Denoising results of DWI in the third direction(a) original DWI; (b) add 5% Rician noise; (c) result using LMMSE with multiple magnitude directions; (d) result using LMMSE with single magnitude directions; (e) result of anisotropic method
除上述定性分析以外,對三種去噪方法的效果進行定量分析,分別計算圖像信噪比(signal noise ratio,SNR)和結構相似度(structure similarity index,SSIM)[13],前者描述圖像的總體噪聲情況,后者強調結構信息的復原效果,計算方法如下式
$SNR=10lg\left[ \frac{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{g{{\left( i,j \right)}^{2}}}}}{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{\left\lceil g\left( i,j \right)-f{{\left( i,j \right)}^{2}} \right\rceil }}} \right],$ |
$SSIM\left( x,y \right)=\frac{\left( 2{{\mu }_{x}}{{\mu }_{y}} \right)\left( 2{{\sigma }_{xy}}+{{c}_{2}} \right)}{\left( \mu _{x}^{2}+\mu _{x}^{2}+{{c}_{1}} \right)\left( \sigma _{x}^{2}+\sigma _{x}^{2}+{{c}_{2}} \right)},$ |
其中μx和μy分別為圖像x,y的均值,σx和σy分別為x,y的標準差,σxy為x,y的協方差。
圖 2為DWI數據集中第3個梯度磁場方向數據在不同噪聲情況下使用三種方法去噪的SNR和SSIM,其中橫軸為噪聲比例,縱軸為SNR和SSIM。

由圖 2可知,隨著噪聲的逐漸增加,DWI的SNR和SSIM都逐漸下降,三種方法均在不同程度上去除噪聲影響,綜合多梯度磁場信息的本文方法在不同噪聲等級下均獲得最高SNR和SSIM。
圖 3(a)為上述DTI數據集的三維結構,圖 3(b)為(a)中部分區域放大圖像,橢圓表征每個體素的彌散張量,橢圓主軸大小和方向與彌散張量主分量的大小和方向相同。圖 3(c)為加入5% Rician噪聲后的結果,如圖所示,由于噪聲影響,本來排列整齊大小均勻的DTI數據均發生明顯變化,其表征彌散大小的橢圓形狀變化明顯,表征彌散方向的橢圓主軸方向也發生不同程度的偏差。

(a) 原始DTI數據;(b) 原始DTI數據放大部分;(c) 加入5% Rician噪聲后DTI數據;(d) 綜合多梯度磁場方向LMMSE;(e) 逐梯度磁場方向LMMSE;(f) 各向異性濾波
Figure3. Denoising result of the synthetic DTI data(a) original DTI data; (b) enlarged area of the original DTI data,(c) add 5% Rician noise; (d) result using LMMSE with multiple magnitude directions; (e) result using LMMSE with single magnitude directions; (f) result of anisotropic method
由圖 3可知,經過三種方法去噪后,計算獲得DTI其張量結構和方向均有明顯改善,張量結構與未加噪標準DTI基本相同,方向基本規律。仔細觀察可見,本文綜合多磁場梯度獲得的去噪后DTI較之其他結果效果更好,張量橢圓無明顯偏差。
除上述對去噪DTI質量影響的定性分析以外,還定量地對上述三種方法DWI去噪對DTI的大小和方向信息復原能力進行了比較,DTI的彌散大小信息由局部各向異性(fractional anisotropy,FA)描述,DTI方向信息保持效果由其主分量的角度偏移確定。
由表 1可知,使用經過三種方法去噪后的DWI計算獲得的DTI的大小和方向信息均有較大提高,較之標準未加噪DTI其角度偏移值和FA的平均變化值均明顯下降,相比而言本文方法效果更為明顯,角度偏移值和FA的平均變化值均較逐梯度磁場方向LMMSE方法進一步下降。

其次,采用本文提出的算法在人腦部DTI數據上進行了成像實驗。人腦部DTI數據使用Philips Intera Achieva 3T MRI 掃描儀獲得,每次掃描矩陣分辨率為128×128,掃描層數為53層,層厚2 mm。DWI由6個不同梯度方向數據(彌散加權值1 000 s/mm2)和1個未加權基本數據(彌散加權值0 s/mm2)構成。
圖 4為真實腦部DWI數據三種方法去噪效果比較。其中圖 4第一行為真實人體DWI數據計算出的DTI數據第28層FA圖,描述了各體素的彌散度大小。第二行為上述對應結果的DTI圖,且僅進行部分放大顯示,背景為FA圖,DTI數據每個體素張量使用橢圓表征。
如圖 4所示,經過三種方法復原的DTI數據 FA值噪聲有明顯下降,灰度級大小變化均勻,沒有異常FA突變;具體分析放大區域,圖 4(g)中橢圓大小和方向均較為凌亂的部分經過三種方法復原后排列均較為整齊,特別是使用本文綜合多梯度磁場方法的圖 4(h)去噪效果更為明顯,與未加噪區域張量結構基本類似。

(a) 原始數據FA圖;(b) 加入10% Rician噪聲FA圖;(c) 綜合多梯度磁場方向LMMSE;(d) 逐梯度磁場方向LMMSE;(e) 各向異性濾波;(f) 原始DTI數據放大圖;(g) 加入10%Rician噪聲DTI數據;(h) 綜合多梯度磁場方向LMMSE;(i) 逐梯度磁場方向LMMSE;(j) 各向異性濾波
Figure4. Denoising results of the in vivo DWI data(a) original FA map; (b) add 10% Rician data; (c) FA map using LMMSE with multiple magnitude directions; (d) FA map using LMMSE with single magnitude directions; (e) FA map using anisotropic method.; (f) enlarged DTI data; (g) add 10% Rician noise; (h) result using LMMSE with multiple magnitude directions; (i) result using LMMSE with single magnitude directions; (j) result of anisotropic method
3 結論
本文針對目前DWI逐磁場梯度方向數據的不足,針對其噪聲呈Rician分布的特點,基于Rician分布改進LMMSE,并綜合多梯度磁場方向DWI數據進行去噪。模擬和真實DWI數據集中的仿真實驗結果均說明,本方法可以有效去除DWI中的噪聲,提高計算獲得的DTI大小和方向信息的準確性和可靠性。
引言
彌散磁共振成像(diffusion magnetic resonance imaging,dMRI)是對活體腦部神經組織結構進行研究的非介入成像技術,通過測量組織中水分子的三維彌散對其生理機能進行定量分析和數字成像,并可利用纖維成像技術重建腦白質神經纖維的三維結構[1]。常用的彌散張量磁共振圖像(diffusion tensor magnetic resonance image,DT-MRI或者DTI)對三維組織的每一體素使用一個3×3對稱正定矩陣構成張量模型進行描述,該模型需要至少6個非共面擴散敏感梯度磁場方向的回波衰減信號測量值構成的彌散加權磁共振圖像(diffusion weighted magnetic resonance image,DW-MRI或者DWI)和一個不施加擴散敏感梯度磁場的磁共振信號參考測量值解得[2]。由于DWI成像使用特殊的自旋平面回波序列易受噪聲干擾[3-4],掃描獲得的DWI數據需要進行有效去噪處理以保證DTI數據計算的準確性和后續應用的可靠性[5]。
現有DWI去噪方法一般由經典灰度圖像去噪推廣而來,使用各向異性濾波[6-7]、極大似然估計[8]等逐個對不同梯度磁場方向DWI進行去噪。現有DWI去噪方法的缺陷主要有兩方面:首先,經典灰度圖像去噪一般針對噪聲模型為高斯噪聲,而DWI的噪聲一般滿足Rician分布,直接使用高斯噪聲去噪方法會產生誤差,需要進行補償[9];另一方面,由于DWI掃描特殊性,僅使用單一梯度磁場DWI去噪并不能獲得最優復原結果[10-11]。
本文使用線性最小均方誤差(linear minimum mean square error,LMMSE)復原去噪方法直接基于Rician噪聲進行去噪,克服常用基于高斯噪聲去噪算法缺陷,同時綜合多梯度方向信息獲得最優去噪結果。下面,首先對本文方法的基本原理進行分析,然后分別使用合成和真實DWI數據對本文方法的有效性進行驗證分析。
1 綜合多梯度磁場方向LMMSE去噪
不同于普通灰度圖像,DWI成像結果為復數,當該信號的實部和虛部同時被方差相等的非相關零均值高斯噪聲時,該信號的模服從Rician分布[3]。令M為測量數據,A為原始數據,則DWI數據集根據掃描矩陣的大小可視為P×Q×R的三維矩陣(P、Q、R分別為其行、列和層數),每一體素存儲一個包含7個(或以上)參數的列向量,滿足Rician分布的DWI數據及其概率密度函數表示為
${{M}_{i}}=\sqrt{{{\left( {{A}_{i}}+{{n}_{i.c}} \right)}^{2}}+n_{i.c}^{2}},$ |
${{P}_{M}}\left( {{M}_{i}}/{{A}_{i}},\sigma \right)=\frac{{{M}_{i}}}{{{\sigma }^{2}}}{{e}^{-\frac{M_{i}^{2}+A_{i}^{2}}{2{{\sigma }^{2}}}}}{{I}_{0}}\left( \frac{{{A}_{i}}{{M}_{i}}}{{{\sigma }^{2}}} \right)u\left( {{M}_{i}} \right),$ |
其中i為梯度磁場編號,ni,c和ni,s為非相關零均值高斯噪聲,方差均為σ2,其二階矩和四階矩為:
${{{{\mu }''}}_{2}}=E\left\{ M_{i}^{2} \right\}=2{{\sigma }^{2}}+A_{i}^{2}$ |
${{{{\mu }''}}_{4}}=E\left\{ M_{i}^{4} \right\}=8{{\sigma }^{4}}+8{{\sigma }^{2}}A_{i}^{2}+A_{i}^{4}$ |
LMMSE的原理即使用測量數據M獲得原始數據A的最優估計值A′,其方法由[12]可得:
${A}'=E\left\{ A \right\}+{{C}_{AM}}C_{AM}^{-1}\left( M-E\{M\} \right)$ |
由于DWI數據中的噪聲呈Rician分布,無法直接通過其一階矩計算,因此最優估計值A′只能通過其平方間接獲得,如下式
$A_{i}^{'2}=E\left\{ A_{i}^{2} \right\}+{{C}_{A_{i}^{2}M_{i}^{2}}}C_{M_{i}^{2}M_{i}^{2}}^{-1}\left( M_{i}^{2}-E\left\{ M_{i}^{2} \right\} \right),$ |
其中CAi2Mi2和CMi2Mi2分別為互協方差和自協方差,則DWI數據集中體素p的濾波器輸出為:
$A_{i}^{'2}\left( p \right)=E\left\{ A_{i}^{2}\left( p \right) \right\}+{{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)C_{M_{i}^{2}M_{i}^{2}}^{-1}\left( p \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right)$ |
CAi2Mi2和CMi2Mi2分別由測量值的矩通過下式獲得:
$\begin{align} & {{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ \left( A_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right) \right\} \\ & =E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}-4{{\sigma }^{2}}\left( {{E}^{2}}\left\{ {{M}_{i}}\left( p \right) \right\}-{{\sigma }^{2}} \right) \\ \end{align}$ |
$\begin{align} & {{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ {{\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right)}^{2}} \right\} \\ & =E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\} \\ \end{align}$ |
將式(8)、(9)帶入式(7),得
$A_{i}^{'2}\left( p \right)=E\left\{ M_{i}^{2}\left( p \right) \right\}-2{{\sigma }^{2}}+K\left( p \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right),$ |
其中,K(p)為:
$K\left( p \right)=1-\frac{4{{\sigma }^{2}}\left( E\left\{ M_{i}^{2}\left( p \right) \right\}-{{\sigma }^{2}} \right)}{E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}}$ |
由于測量數據M的矩無法直接獲得,因此需要對當前體素p的鄰域采樣進行估計:
$E\left\{ M_{i}^{r}\left( p \right) \right\}=\frac{1}{S}\sum\limits_{q\in S}{M_{i}^{r}\left( q \right)},$ |
其中S為進行采樣選擇的鄰域區域。
令Mi2和Ai2分別為[Mi-n2,…Mi2,…Mi+n2]T和[Ai-n2,…Ai2,…Ai+n2]T,其中n為常數,Mi2和Ai2表征以第i個梯度磁場為中心的多梯度磁場信息,其互協方差和自協方差分別為:
$\begin{align} & {{C}_{A_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ \left( A_{i}^{2}\left( p \right)-E\left\{ A_{i}^{2}\left( p \right) \right\} \right) \right.\left. \left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right) \right\} \\ & =E\left\{ M_{i}^{2}\left( p \right) \right\}=E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}- \\ & 4{{\sigma }^{2}}\left( {{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\}-{{\sigma }^{2}} \right) \\ \end{align}$ |
$\begin{align} & {{C}_{M_{i}^{2}M_{i}^{2}}}\left( p \right)=E\left\{ M_{i}^{2}\left( p \right)-E{{\left\{ M_{i}^{2}\left( p \right) \right\}}^{2}} \right\} \\ & =E\left\{ M_{i}^{4}\left( p \right) \right\}-{{E}^{2}}\left\{ M_{i}^{2}\left( p \right) \right\} \\ \end{align}$ |
則式(10)改寫為
$A_{i}^{'2}\left( p \right)=E\left\{ M_{i}^{2}\left( p \right) \right\}-2{{\sigma }^{2}}+K\left( p \right)\left( M_{i}^{2}\left( p \right)-E\left\{ M_{i}^{2}\left( p \right) \right\} \right),$ |
其中噪聲方差σ2根據Rician噪聲特點,對原始數據為0區域(圖像背景)估計獲得。由式(3) 可知,當A2=0時,有μ2″=E{Mi2}=2σ2,則
${{\sigma }^{2}}=\frac{1}{2N}\sum\limits_{i=1}^{N}{M_{i}^{2}},$ |
其中 N為估計所用背景區域大小。
2 仿真與實驗結果
在合成和真實腦部DWI數據集中,分別使用綜合多梯度磁場方向LMMSE、逐梯度磁場方向LMMSE和各向異性濾波[6]三種方法進行去噪復原比較。
首先,在合成DWI數據集上進行仿真,仿真平臺為MATLAB 7.0.4,Windows 7 Professional操作系統,Intel Core(TM) i7 CPU處理器和4G內存。合成DWI數據集大小為64×64×1,由6個不同方向的DWI數據集和1個標準方向DWI數據集使用最小二乘模型計算獲得,合成纖維結構的FA為0.8,彌散度為2.1×10-5 cm2/s,幾何結構為一條正弦曲線,寬度為10體素。
由圖 1可知,三種方法均能在一定程度上去除DWI噪聲,相比而言,基于Rician噪聲模型的LMMSE方法綜合數據的統計信息,較之局部灰度差異的各向異性濾波方法可以獲得更好的去噪效果,特別綜合多梯度磁場方向提供更多信息獲得最優去噪效果。

(a) 未加噪DWI;(b) 加入5% Rician噪聲;(c) 綜合多梯度磁場方向LMMSE;(d) 逐梯度磁場方向LMMSE;(e) 各向異性濾波
Figure1. Denoising results of DWI in the third direction(a) original DWI; (b) add 5% Rician noise; (c) result using LMMSE with multiple magnitude directions; (d) result using LMMSE with single magnitude directions; (e) result of anisotropic method
除上述定性分析以外,對三種去噪方法的效果進行定量分析,分別計算圖像信噪比(signal noise ratio,SNR)和結構相似度(structure similarity index,SSIM)[13],前者描述圖像的總體噪聲情況,后者強調結構信息的復原效果,計算方法如下式
$SNR=10lg\left[ \frac{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{g{{\left( i,j \right)}^{2}}}}}{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{\left\lceil g\left( i,j \right)-f{{\left( i,j \right)}^{2}} \right\rceil }}} \right],$ |
$SSIM\left( x,y \right)=\frac{\left( 2{{\mu }_{x}}{{\mu }_{y}} \right)\left( 2{{\sigma }_{xy}}+{{c}_{2}} \right)}{\left( \mu _{x}^{2}+\mu _{x}^{2}+{{c}_{1}} \right)\left( \sigma _{x}^{2}+\sigma _{x}^{2}+{{c}_{2}} \right)},$ |
其中μx和μy分別為圖像x,y的均值,σx和σy分別為x,y的標準差,σxy為x,y的協方差。
圖 2為DWI數據集中第3個梯度磁場方向數據在不同噪聲情況下使用三種方法去噪的SNR和SSIM,其中橫軸為噪聲比例,縱軸為SNR和SSIM。

由圖 2可知,隨著噪聲的逐漸增加,DWI的SNR和SSIM都逐漸下降,三種方法均在不同程度上去除噪聲影響,綜合多梯度磁場信息的本文方法在不同噪聲等級下均獲得最高SNR和SSIM。
圖 3(a)為上述DTI數據集的三維結構,圖 3(b)為(a)中部分區域放大圖像,橢圓表征每個體素的彌散張量,橢圓主軸大小和方向與彌散張量主分量的大小和方向相同。圖 3(c)為加入5% Rician噪聲后的結果,如圖所示,由于噪聲影響,本來排列整齊大小均勻的DTI數據均發生明顯變化,其表征彌散大小的橢圓形狀變化明顯,表征彌散方向的橢圓主軸方向也發生不同程度的偏差。

(a) 原始DTI數據;(b) 原始DTI數據放大部分;(c) 加入5% Rician噪聲后DTI數據;(d) 綜合多梯度磁場方向LMMSE;(e) 逐梯度磁場方向LMMSE;(f) 各向異性濾波
Figure3. Denoising result of the synthetic DTI data(a) original DTI data; (b) enlarged area of the original DTI data,(c) add 5% Rician noise; (d) result using LMMSE with multiple magnitude directions; (e) result using LMMSE with single magnitude directions; (f) result of anisotropic method
由圖 3可知,經過三種方法去噪后,計算獲得DTI其張量結構和方向均有明顯改善,張量結構與未加噪標準DTI基本相同,方向基本規律。仔細觀察可見,本文綜合多磁場梯度獲得的去噪后DTI較之其他結果效果更好,張量橢圓無明顯偏差。
除上述對去噪DTI質量影響的定性分析以外,還定量地對上述三種方法DWI去噪對DTI的大小和方向信息復原能力進行了比較,DTI的彌散大小信息由局部各向異性(fractional anisotropy,FA)描述,DTI方向信息保持效果由其主分量的角度偏移確定。
由表 1可知,使用經過三種方法去噪后的DWI計算獲得的DTI的大小和方向信息均有較大提高,較之標準未加噪DTI其角度偏移值和FA的平均變化值均明顯下降,相比而言本文方法效果更為明顯,角度偏移值和FA的平均變化值均較逐梯度磁場方向LMMSE方法進一步下降。

其次,采用本文提出的算法在人腦部DTI數據上進行了成像實驗。人腦部DTI數據使用Philips Intera Achieva 3T MRI 掃描儀獲得,每次掃描矩陣分辨率為128×128,掃描層數為53層,層厚2 mm。DWI由6個不同梯度方向數據(彌散加權值1 000 s/mm2)和1個未加權基本數據(彌散加權值0 s/mm2)構成。
圖 4為真實腦部DWI數據三種方法去噪效果比較。其中圖 4第一行為真實人體DWI數據計算出的DTI數據第28層FA圖,描述了各體素的彌散度大小。第二行為上述對應結果的DTI圖,且僅進行部分放大顯示,背景為FA圖,DTI數據每個體素張量使用橢圓表征。
如圖 4所示,經過三種方法復原的DTI數據 FA值噪聲有明顯下降,灰度級大小變化均勻,沒有異常FA突變;具體分析放大區域,圖 4(g)中橢圓大小和方向均較為凌亂的部分經過三種方法復原后排列均較為整齊,特別是使用本文綜合多梯度磁場方法的圖 4(h)去噪效果更為明顯,與未加噪區域張量結構基本類似。

(a) 原始數據FA圖;(b) 加入10% Rician噪聲FA圖;(c) 綜合多梯度磁場方向LMMSE;(d) 逐梯度磁場方向LMMSE;(e) 各向異性濾波;(f) 原始DTI數據放大圖;(g) 加入10%Rician噪聲DTI數據;(h) 綜合多梯度磁場方向LMMSE;(i) 逐梯度磁場方向LMMSE;(j) 各向異性濾波
Figure4. Denoising results of the in vivo DWI data(a) original FA map; (b) add 10% Rician data; (c) FA map using LMMSE with multiple magnitude directions; (d) FA map using LMMSE with single magnitude directions; (e) FA map using anisotropic method.; (f) enlarged DTI data; (g) add 10% Rician noise; (h) result using LMMSE with multiple magnitude directions; (i) result using LMMSE with single magnitude directions; (j) result of anisotropic method
3 結論
本文針對目前DWI逐磁場梯度方向數據的不足,針對其噪聲呈Rician分布的特點,基于Rician分布改進LMMSE,并綜合多梯度磁場方向DWI數據進行去噪。模擬和真實DWI數據集中的仿真實驗結果均說明,本方法可以有效去除DWI中的噪聲,提高計算獲得的DTI大小和方向信息的準確性和可靠性。