楚春雨

渤海大學(xué)工學(xué)院，遼寧錦州121013

前言

磁共振擴散成像是目前唯一能夠在活體上測量組織內(nèi)水分子擴散運動與成像的無創(chuàng)方法，它通過測量和量化組織中水分子的擴散信息來探測組織的微觀結(jié)構(gòu)。水分子沿不同方向的擴散信息包含在一組不同擴散梯度方向的擴散加權(quán)圖像（Diffusion Weighted Image,DWI）中，通過對擴散函數(shù)進行建模可以解析出每個體素內(nèi)的纖維結(jié)構(gòu)信息（主要是纖維的走向）。根據(jù)每個體素內(nèi)的纖維走向，利用纖維束追蹤技術(shù)可重建出組織纖維束的整體三維結(jié)構(gòu)，并可從中提取有效的臨床統(tǒng)計特征，從而用于醫(yī)療診斷及相關(guān)研究等。

擴散張量成像（Diffusion Tensor Imaging,DTI）［1］是最早提出的一種擴散模型，目前已被廣泛用于臨床和醫(yī)學(xué)研究中。該方法將擴散函數(shù)建模為D(g)=gTDg，其中g(shù)是擴散加權(quán)梯度方向，D是二階對稱正定張量。這種建模方法雖然簡單、穩(wěn)定，但所采用的二階對稱張量僅能描述體素內(nèi)的單一平均纖維走向，而無法用于描述存在纖維交叉、分叉等體素內(nèi)復(fù)雜纖維結(jié)構(gòu)。

目前，已有一些方法能夠通過更復(fù)雜的模型來估計體素內(nèi)的多個纖維走向，從而解決DTI方法的限制。這些方法主要包括：多張量模型［2-4］、高階擴散張量［5］、Q-Ball 成像［6-9］、擴散譜成像［10-12］、球面反卷積（Spherical Deconvolution,SD）［13-17］、獨立分量分析［18］、混合擴散成像［19］等。其中，基于SD的方法由于具備不需指定纖維分布的數(shù)量、計算效率高且能夠在低角分辨率成像條件下估計體素內(nèi)纖維走向分布等優(yōu)點，得到較廣泛的應(yīng)用，但是基于SD 的方法對噪聲非常敏感。為了解決SD 方法對圖像噪聲敏感的問題，通常有兩種改進方案：一是在SD 模型中引入正則化技術(shù)；二是在使用SD方法進行體素內(nèi)纖維結(jié)構(gòu)估計前，先采用有效的降噪方法對DWI 圖像進行預(yù)處理。非局部均值（Non-Local Mean, NLM）作為一種經(jīng)典的圖像平滑方法已被應(yīng)用于DWI 降噪［20］，在使用NLM方法對DWI降噪時，只是簡單地對一組不同擴散加權(quán)梯度方向的DWI 圖像分別進行NLM 平滑，往往難以得到較好的結(jié)果。

本研究針對NLM 方法直接用于DWI 降噪存在的不足，提出一種面向DWI 的NLM 方法，該方法在對DWI平滑過程中既考慮到不同位置體素間的關(guān)聯(lián)性，同時也考慮到不同擴散加權(quán)梯度方向數(shù)據(jù)間的關(guān)聯(lián)性。原始DWI 數(shù)據(jù)經(jīng)本研究提出的方法處理后，再采用SD 方法進行體素內(nèi)纖維結(jié)構(gòu)估計，可有效提高體素內(nèi)纖維結(jié)構(gòu)估計結(jié)果的準(zhǔn)確性。

1 方法

1.1 基于SD的體素內(nèi)纖維結(jié)構(gòu)估計

SD方法［13］是目前應(yīng)用較廣的體素內(nèi)纖維結(jié)構(gòu)估計方法，該方法將磁共振擴散信號建模為：

其中，S為擴散加權(quán)信號，R為響應(yīng)函數(shù)，F(xiàn)為纖維走向分布函數(shù)，?為卷積操作。通過構(gòu)造合適的響應(yīng)函數(shù)R，就可以通過反卷積方法求得纖維走向分布函數(shù)F。在該模型下，采用纖維走向分布來表征體素內(nèi)纖維結(jié)構(gòu)。更多關(guān)于SD 的細(xì)節(jié)信息參見文獻(xiàn)［13］。

1.2 經(jīng)典的非局部均值算法

非局部均值算法由Buades 等于2005年首次提出，該算法基于圖像信息的自然冗余進行加權(quán)平均降噪。在NLM算法的理論公式中，體素xi的非局部均值NLM(xi)是圖像I上所有體素的加權(quán)平均，即：

其中，v(xj)為體素xj的灰度值，ω(xi,xj)為體素xj對體 素xi的權(quán)重系數(shù) ，ω(xi,xj)∈[0,1] 并且∑xj∈Iω(xi,xj)=1。

在NLM算法的實際應(yīng)用過程中，對與xi相似體素的搜索通常限定在其鄰域Vi中，因此式（2）可寫為：

對于在Vi中的每個體素xi，權(quán)重ω(xi,xj) 定義為：

其中，Ni和Nj分別為體素xi和xj的鄰域，Z(i)=∑jω(xi,xj)為歸一化常數(shù)，為圖像噪聲標(biāo)準(zhǔn)差的估計值，h為濾波器參數(shù)，距離d定義為：

其中，N為鄰域Ni中體素的數(shù)量，其值等于鄰域Nj中體素的數(shù)量，yk和zk分別為鄰域Ni和Nj中的第k個體素。對于灰度圖像，

1.3 面向磁共振DWI的NLM算法

對于一組DWI 圖像，其灰度值（擴散加權(quán)信號值）可表示為S(p,g)，p為體素的位置，g為擴散加權(quán)梯度方向?？紤]角度信息的DWI 的NLM 平滑算法可表示為：

其中，pk為pi的鄰域Vpi中體素的位置，G為所有擴散加權(quán)梯度方向的集合，權(quán)重ωDWI定義為：

其中，Npi和Npk分別為體素pi和pk的鄰域，ZD為歸一化常數(shù)，σ?D為DWI 圖像噪聲標(biāo)準(zhǔn)差的估計值，hD為濾波器參數(shù)，距離Dis 定義為：

其中，ND為鄰域Npi中體素的數(shù)量，其值等于鄰域Npk中體素的數(shù)量，qm和rm分別為鄰域Npi和Npk中第m個體素的位置，Δ(S(qm,gj),S(rm,gl))=‖qm-rm‖2lg2(|gjTgl|)。

2 實驗與分析

本研究提出的基于NLM平滑的體素內(nèi)纖維結(jié)構(gòu)估計方法分別在數(shù)值仿真數(shù)據(jù)、仿真實體數(shù)據(jù)上進行了實驗測試，并與不對DWI 圖像進行降噪直接采用SD方法［13］以及先利用經(jīng)典的NLM方法［20］對DWI圖像進行降噪后再使用SD 方法（NLM+SD）進行了對比。在數(shù)值仿真數(shù)據(jù)上，通過定量比較不同方法得到的體素內(nèi)纖維結(jié)構(gòu)的平均角度誤差來評價不同方法的準(zhǔn)確性，本研究中平均角度誤差定義為估計得到的纖維走向與實際纖維走向之間夾角的平均值；在仿真實體數(shù)據(jù)上，根據(jù)已知的纖維路徑的走向來推斷體素內(nèi)纖維結(jié)構(gòu)的實際情況，并將不同方法得到的結(jié)果與之比較，從而評價不同方法的準(zhǔn)確性。

2.1 數(shù)值仿真數(shù)據(jù)實驗

實驗采用多張量模型［2］的簡單形式二張量模型生成一組數(shù)值仿真數(shù)據(jù)，二張量模型是將一個體素內(nèi)沿擴散加權(quán)梯度方向g的擴散加權(quán)信號S(g) 建模為兩項擴散張量項的和的形式：

其中，b為擴散敏感因子，f為體積分?jǐn)?shù)，D1和D2為兩擴散張量，擴散張量可以根據(jù)指定的特征向量和特征值進行構(gòu)造。在數(shù)值仿真中，為了簡化操作，通常設(shè)定S0=1，f=0.5，并假定擴散張量的第二和第三特征值相同，同時令D1和D2的特征值相同。因此只要給定b值，一定數(shù)量的擴散加權(quán)梯度方向，兩個擴散張量的主特征方向，以及張量的特征值，即可生成仿真的擴散加權(quán)信號。在本仿真中，給定b值為1 000 s/mm2，30個擴散梯度方向，兩個擴散張量的主特征方向如圖1a所示，張量的特征值為（1.2,0.3,0.3）×10-3mm2/s。由于數(shù)值仿真得到的是無噪聲數(shù)據(jù)，為了測試算法的抗噪性能，向數(shù)值仿真數(shù)據(jù)中加入不同強度的萊斯噪聲，從而得到不同信噪比的數(shù)值仿真數(shù)據(jù)。

圖1b～圖1d給出不同方法在一組信噪比為25的數(shù)值仿真數(shù)據(jù)上的估計結(jié)果。從圖1b可以看出直接采用SD 方法估計得到的體素內(nèi)纖維結(jié)構(gòu)呈現(xiàn)出較為混亂的分布，這說明SD 方法極易受到噪聲的影響，在噪聲強度不高的情況下，也會導(dǎo)致估計結(jié)果的誤差大大增加。從圖1c 和圖1d 的結(jié)果可以看出，在進行體素內(nèi)纖維結(jié)構(gòu)估計前對DWI數(shù)據(jù)進行平滑能夠有效提升纖維結(jié)構(gòu)估計結(jié)果的質(zhì)量，減小估計誤差。通過對比圖1c 和圖1d 的結(jié)果可以看出，NLM+SD方法一些邊緣區(qū)域存在模糊現(xiàn)在，例如圖1c中紅色方框所標(biāo)注的體素，對于一些實際只有單一纖維走向的體素，NLM+SD方法卻得到了兩個纖維走向；而對于一些實際上有兩個體積分?jǐn)?shù)相同的纖維走向的體素，NLM+SD 方法得到的結(jié)果卻顯示其中一個纖維走向的體積分?jǐn)?shù)較小。而本研究方法在這組仿真數(shù)據(jù)上并沒有顯示出這些問題，可以看出本研究方法明顯優(yōu)于NLM+SD方法。

圖1 不同方法估計得到的體素內(nèi)纖維結(jié)構(gòu)Fig.1 Intravoxel fiber structures obtained by different methods

為了定量比較不同方法的估計結(jié)果，通過計算圖1b～圖1d所示估計結(jié)果的平均角度誤差，得到SD方法的平均角度誤差為14.8°，NLM+SD方法的平均角度誤差為7.5°，而本研究方法的平均角度誤差為4.8°，顯然本研究方法得到的結(jié)果誤差低于其它方法。為了比較不同噪聲強度對3種方法（SD方法、NLM+SD方法和本研究方法）的估計結(jié)果的影響，在5種信噪比條件下，不同方法估計結(jié)果的平均角度誤差如圖2所示。與SD方法和NLM+SD方法相比，本研究方法的估計結(jié)果的平均角度誤差約降低25%～30%。

2.2 仿真實體數(shù)據(jù)實驗

實驗采用的仿真實體數(shù)據(jù)來源于MICCAI 2009國際會議，該數(shù)據(jù)最初是用于人腦白質(zhì)纖維追蹤算法競賽，關(guān)于該組數(shù)據(jù)的更多細(xì)節(jié)信息參見文獻(xiàn)［21］。該比賽的數(shù)據(jù)集中包含多組不同參數(shù)的DWI數(shù)據(jù)，本實驗中采用其中一組b值為1 500 s/mm2、64個擴散加權(quán)梯度方向的數(shù)據(jù)。

圖2 不同信噪比條件下3種方法的平均角度誤差Fig.2 Mean angular errors of 3 methods at different signal-tonoise ratios

為了比較經(jīng)典NLM方法與本研究提出的面向DWI數(shù)據(jù)的NLM方法，圖3a～圖3c分別給出未經(jīng)降噪處理的DWI數(shù)據(jù)的各向異性分?jǐn)?shù)（Fractional Anisotropy,FA）圖、采用經(jīng)典NLM方法處理后的DWI數(shù)據(jù)對應(yīng)的FA圖以及采用本研究提出的面向DWI數(shù)據(jù)的NLM方法處理后的DWI數(shù)據(jù)對應(yīng)的FA圖。從圖中可以明顯看出，采用經(jīng)典NLM方法或本研究方法處理后，相應(yīng)的FA圖的噪聲明顯降低。與經(jīng)典NLM方法的結(jié)果（圖3b）相比，本研究方法的結(jié)果（圖3c）在存在纖維束的區(qū)域內(nèi)的均勻性更好，在纖維束區(qū)域與背景邊界更清晰，沒有明顯的模糊現(xiàn)象。

圖3 不同方法得到的FA圖和體素內(nèi)纖維結(jié)構(gòu)Fig.3 Fractional anisotropy maps and intravoxel fiber structures obtained by different methods

為了更好地比較不同方法的效果，選取如圖3a～圖3c 中紅色方框所標(biāo)識的感興趣區(qū)域，該區(qū)域包含兩束纖維交叉的情況，分別采用3 種方法（SD 方法、NLM+SD方法和本研究方法）對該區(qū)域內(nèi)的體素進行纖維結(jié)構(gòu)估計，得到的估計結(jié)果如圖3d～圖3f 所示?？梢钥闯?，在纖維交叉處SD方法的結(jié)果比較混亂，說明SD 方法容易受噪聲影響，估計結(jié)果的誤差較大，而NLM+SD 方法和本研究方法的結(jié)果均具有較好的局部一致性。對比圖3d～圖3f 中紅色橢圓標(biāo)記的區(qū)域，該區(qū)域的體素內(nèi)纖維結(jié)構(gòu)應(yīng)為單一纖維走向，而NLM+SD 方法卻得到存在體素內(nèi)纖維交叉的結(jié)果，這說明NLM+SD 方法在存在纖維交叉的區(qū)域和不存在纖維交叉的區(qū)域的邊界處出現(xiàn)了一定的模糊現(xiàn)象，而本研究方法的結(jié)果則較少存在這種情況。

3 結(jié)論

本研究提出一種基于NLM平滑的體素內(nèi)纖維結(jié)構(gòu)估計方法，該方法通過建立面向DWI數(shù)據(jù)的NLM平滑方法，在考慮DWI 數(shù)據(jù)不同位置的體素間的關(guān)聯(lián)性同時也考慮到不同擴散加權(quán)梯度方向的數(shù)據(jù)間的關(guān)聯(lián)性，對DWI 數(shù)據(jù)進行NLM 平滑后再采用SD方法估計體素內(nèi)纖維結(jié)構(gòu)。數(shù)值仿真數(shù)據(jù)和仿真實體數(shù)據(jù)實驗結(jié)果表明，與SD 和NLM+SD 方法相比，采用本研究方法得到的體素內(nèi)纖維結(jié)構(gòu)的平均角度誤差更小，且較少存在邊緣模糊現(xiàn)象。