吳 錫周激流畢務忠謝明元羅代升

1(成都信息工程學院電子工程系，成都 610225)2(四川大學電子信息學院，成都 610065)

采用偏微分方程的磁共振彌散張量各向異性擴散濾波去噪

吳 錫1，2*周激流2畢務忠1謝明元1羅代升2

1(成都信息工程學院電子工程系，成都 610225)2(四川大學電子信息學院，成都 610065)

彌散張量磁共振成像(DTI)是無創(chuàng)研究大腦白質結構及其他人體纖維狀組織結構的主要工具。由于合成DTI數(shù)據(jù)的彌散加權成像數(shù)據(jù)(DWI)易受噪聲干擾，需要有效去噪以保證DTI數(shù)據(jù)精度和后續(xù)應用的實現(xiàn)。使用各向異性擴散濾波理論，綜合考慮各方向通道DWI數(shù)據(jù)的幾何形態(tài)和結構特點，重構其特征向量和特征值，獲得統(tǒng)一平滑的結構張量，以期在有效去噪的基礎上最大程度地保持DTI數(shù)據(jù)幾何結構和特征。利用所提出方法在合成彌散張量數(shù)據(jù)上進行仿真，并在真實腦部DTI數(shù)據(jù)上進行實驗。仿真和實驗結果表明，該方法能有效減少噪聲對DTI數(shù)據(jù)的影響，較之現(xiàn)有的時頻分析去噪方法，可更準確地恢復DTI數(shù)據(jù)，減少主分量方向的偏差和部分各向異性值的誤差。

彌散張量磁共振成像;各向異性擴散濾波;偏微分方程

Abstract:Diffusion tensor magnetic resonance imaging(DTI)is the main non-invasive utility to reveal the information of local diffusivity of white matter and other fibrous human tissues.Because the diffusion weighted image(DWI)which is used to acquire DTI is sensitive to noise，effective noise removal is required to improve the accuracy of the DTI data and its subsequent applications.Instead of denoising the DWI in different direction separately，an anistropic filtering method synthetically considering the structure information and characteristic of the DTI was proposed in this paper.The eigenvalue and eigenvector of the smoothing structure tensor were reconstructed to remove the noise and keep the structure characteristic simultaneously.Simulations using a synthetic DTI dataset and experiments using an in vivo brain DTI dataset were performed.The results demonstrated that the noise could be removed significantly and the direction bias of the main eigenvector and the error of fractional anisotropy were reduced effectively compared with the common methods.

Key words:diffusion tensor magnetic resonance imaging;anisotropic smoothing;partial differential equation

引言

彌散張量磁共振成像(diffusion tensor magnetic resonance imaging，DT-MRI或者 DTI)，是對大腦白質結構進行無創(chuàng)研究的主要工具［1］。利用腦白質中水分子彌散效應，對腦白質三維結構的每一體素使用一個3×3對稱正定矩陣描述，該矩陣即為體素彌散張量。該矩陣通過對每一體素，采集6個非共面的擴散敏感梯度磁場方向的回波衰減信號測量值構成的彌散加權圖像(diffusion weighted image，DWI)和一個不施加擴散敏感梯度磁場的MR信號參考測量值解得。在實際掃描過程中，一般采集多于6個的擴散敏感方向信號，在每個方向上采用更多的編碼幅度，以此來抑制成像噪聲。在得到多于6個的擴散編碼方向的測量結果后，可采用最小二乘擬合方法求解彌散矩陣，獲得DTI數(shù)據(jù)。DTI也同時大量應用在腦損傷、腦發(fā)育［2］及人體其他纖維狀組織(如骨骼?。?］、筋腱［4］等)研究。

由于DWI在快速成像時極易受系統(tǒng)噪聲干擾，同時在掃描成像過程中由于呼吸、心跳等會產(chǎn)生生理性噪聲，通過DWI獲得的DTI數(shù)據(jù)受噪聲干擾極大［5］，另外DTI的主要應用如纖維追蹤成像等對于噪聲極為敏感，少量隨機噪聲即會影響全部成像結果，因此DTI數(shù)據(jù)應用的首要工作就是去噪，以獲得正確的彌散張量。

現(xiàn)有的DTI去噪方法一般分為3類:第一類是對DWI數(shù)據(jù)使用傳統(tǒng)時域或者頻域手段進行去噪［6］，但由于是分別對不同方向通道的DWI數(shù)據(jù)進行濾波，沒有綜合考慮DTI數(shù)據(jù)的形態(tài)和結構等信息，容易丟失信息;第二類是在DWI估計 DTI的過程中，對該過程進行正則化去噪［7］，但由于 DWI數(shù)據(jù)本身的復雜性，使該正則過程的適應受到很多約束;第三類是對 DTI數(shù)據(jù)進行去噪［8］，由于 DTI數(shù)據(jù)是二階張量，將極大地增加去噪的計算量和計算復雜度。

本研究提出采用偏微分方程的DTI各向異性擴散去噪方法，綜合考慮所有不同通道DWI數(shù)據(jù)所包含的DTI結構和形態(tài)信息，重構特征值和特征向量;對每個方向通道的DWI數(shù)據(jù)，使用統(tǒng)一平滑結構張量進行去噪，既能保持DTI的結構和邊沿信息，同時又能有效去噪。下面，首先介紹采用偏微分方程的DTI各向異性擴散去噪的基本原理，然后闡述采用這種方法在人工合成數(shù)據(jù)集中的仿真試驗和結果，最后介紹采用這種方法對真實腦部DTI數(shù)據(jù)進行去噪的實驗和結果。

1 采用偏微分方程的DTI各向異性濾波去噪

采用偏微分方程的各向異性濾波，可用下述模型表示［9］為

式中，ρ為卷積核的尺度參數(shù)。

G的特征向量 v1，v2，v3及其對應的特征值λ1，λ2，λ3表征 ρ內灰度梯度的變化，最大特征值及其對應的特征向量代表灰度梯度變化最大的方向，反之亦然。當λ1≈λ2≈λ3≈0時，為各向同性結構(即平坦區(qū)域);當 λ1＞λ2≈λ3≈0時，則為各向異性結構(邊緣區(qū)域)。結構張量T由灰度的梯度張量G按照一定比例重構特征值獲得，有

式中，A為梯度張量G的各向異性系數(shù)，a為歸一化參數(shù)，C為閾值參數(shù)。

在平坦區(qū)域，A?C，λ3→a，使用各向同性濾波;在邊緣區(qū)域，A?C，λ3→1，使用各向異性濾波，在梯度小的方向進行大尺度濾波，即在邊沿的切向方向增強濾波。

不同于傳統(tǒng)DWI去噪僅考慮各獨立方向通道的數(shù)據(jù)特點，對于不同方向通道的 DWI數(shù)據(jù)，綜合考慮DTI結構信息，使用統(tǒng)一的平滑結構張量，使在一個DWI方向通道上失去的結構信息可以從其他DWI數(shù)據(jù)彌補，其梯度張量 G′可由各 DWI方向通道的梯度和與高斯核卷積獲得，即

平滑結構張量T′由梯度張量G′重構獲得。不同于傳統(tǒng)各向異性濾波主要針對線狀邊沿，DTI由于其生理特點，其結構主要以帶狀分布為主，特別在腦白質纖維束內部，彌散信息變化較大，容易出現(xiàn)除傳統(tǒng)各向同性結構(λ1≈λ2≈λ3≈0)和各向異性結構(λ1≈λ2＞λ3≈0)以外的特殊結構，其特征值滿足 λ1≥λ2≥λ3＞0。針對腦白質 DTI數(shù)據(jù)的特殊性，在各向同性結構區(qū)域，希望在各個方向都有較強擴散去噪;在邊緣區(qū)域，希望沿邊沿切向方向有較強擴散去噪，垂直于邊沿盡量保持不要擴散;最后對于腦白質帶狀纖維束內部進行相應小尺度各向異性擴散，以保持數(shù)據(jù)精度。

根據(jù)上述分析，對梯度張量G′進行重構，T′的特征向量為

T′的特征值重構為:

其中，λ1，λ2由下式獲得，有式中 c(x，y，z)為權函數(shù)，且為的減函數(shù)，滿足 0≤c(x，y，z)≤1，令其為c(x，y，z)

由上述特征值和特征向量，即可重構平滑結構張量T′。當位于各向同性結構時，各方向特征值均較大，保證各方向均進行較大的擴散濾波;當位于各向異性結構時，沿邊緣方向特征值大，垂直邊緣特征值小，保證邊緣增強去噪;當位于腦白質纖維束內部區(qū)域時，特征值均較小，使用較小的擴散保持數(shù)據(jù)精度。

圖1 合成DTI數(shù)據(jù)用本方法和Parker方法的去噪結果比較。(a)合成DTI三維數(shù)據(jù);(b)合成DTI二維放大結構;(c)加入5%高斯噪聲的結果;(d)加入10%高斯噪聲的結果;(e)和(f)分別加入5%和10%高斯噪聲的本方法去噪結果;(g)和(f)分別加入5%和10%高斯噪聲Parker方法去噪結果Fig.1 Denoising results of the proposed method and Parker’s method using the synthetic phantom.(a)A 3D view of one slice of the synthetic dataset;(b)Enlarged 2D view of the boxed region in(a);(c)Result of adding 5%Gaussian noise;(d)Result of adding 10%Gaussian noise;(e)and(f)Result of the proposed method to 5%and 10%Gaussian noise;(g)and(h)Result of the Parker’s method to 5%and 10%Gaussian noise

2 仿真與實驗結果

在合成數(shù)據(jù)集和腦部真實DTI數(shù)據(jù)集中，分別使用本方法和Parker的方法［6］進行去噪，并對結果進行比較。

首先，采用所提出的算法在合成張量數(shù)據(jù)上進行了仿真，實驗平臺為Matlab 7.0.4，Windows Vista操作系統(tǒng)，處理器為Intel Core2 Duo CPU(1.80 GHz)，系統(tǒng)內存4 GB。圖1(a)為模擬數(shù)據(jù)集，掃描矩陣分辨率為64像素×64像素，層數(shù)為20層，合成纖維束張量的 FA為0.8，彌散度為2.1×10-5cm2/s。如圖1(a)所示，合成彌散張量數(shù)據(jù)分為橫向和縱向兩組，其余部分為非彌散區(qū)域。圖1(b)為圖1(a)框中圖像的放大結果，其中的直線段表征彌散張量，其長短和方向與該體素彌散張量主分量相同，非彌散區(qū)域用圓點表示。圖1(c)和(d)分別為上述DTI數(shù)據(jù)集加入5%和10%高斯噪聲的結果，由于噪聲干擾DTI數(shù)據(jù)的彌散大小和方向都發(fā)生了變化，隨著噪聲增加，出現(xiàn)了明顯的錯誤 DTI數(shù)據(jù)。圖1(e)和(f)是使用本方法分別對受5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集去噪結果，圖1(g)和(h)是使用Parker方法分別對受5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集去噪結果。

由圖1(e)和(g)可知，對5%高斯噪聲分別使用本方法和Parker方法，均可在一定程度上抑制噪聲的影響。相比而言，使用本方法的去噪效果更明顯，除少數(shù)體素外，DTI數(shù)據(jù)的彌散大小和方向基本與未受噪聲干擾的DTI數(shù)據(jù)一致，特別是在非彌散區(qū)域，基本抑制了噪聲影響，見圖1(e)。由圖1(f)和(g)可知，在噪聲較大的情況下，兩種方法均可有效去噪。相比而言，本方法的去噪結果其彌散大小和方向與未加噪數(shù)據(jù)基本吻合，見圖1(f)。

DTI數(shù)據(jù)的主要應用之一為纖維追蹤成像，主要使用DTI數(shù)據(jù)彌散方向信息進行追蹤成像。因此，對DTI數(shù)據(jù)方向信息進行復原和保護是DTI去噪的主要目的。使用本方法，分別對加入5%和10%高斯噪聲的DTI數(shù)據(jù)去噪，對方向信息的復原效果見圖2，數(shù)據(jù)參數(shù)同圖1。橫軸為迭代次數(shù)，縱軸為所有體素DTI數(shù)據(jù)主分量方向與標準DTI數(shù)據(jù)主分量方向角度偏差的均值。

圖2 使用本方法的DTI數(shù)據(jù)方向信息復原效果分析Fig.2 Variation of mean angular difference of the proposed method in different noise

圖3 真實腦部DTI數(shù)據(jù)用本方法和Parker方法的去噪效果比較。(a)原始FA圖;(b)加入5%高斯噪聲的FA圖;(c)加入10%高斯噪聲的FA圖;(d)和(e)用本方法加入5%和10%高斯噪聲的去噪結果;(f)和(g)用Parker方法加入5%和10%高斯噪聲的去噪結果Fig.3 Filtering results using the proposed method and Parker’s method in real human in vivo dataset.(a)One slice FA map of the in vivo dataset;(b)FA map of adding 5%Gaussian noise;(c)FA map of adding 10%Gaussian noise;(d)and(e)Result of the proposed method to 5%and 10%Gaussian noise;(f)and(g)Result of the Parker’s method to 5%and 10%Gaussian noise

由圖2可知，較之標準的DTI數(shù)據(jù)集，被5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集其平均彌散方向分別有大約4.5°和9.5°的偏差。使用本方法進行去噪后，隨著迭代次數(shù)增加，平均角度偏差逐漸減小，大約在第10個迭代周期后基本達到穩(wěn)定，平均角度偏差分別下降至2°和5°。

此外，采用所提出的算法在人腦部DTI數(shù)據(jù)上進行了成像實驗。人腦部 DTI數(shù)據(jù)使用 Philips Intera Achieva 3T MRI掃描儀獲得，每次掃描矩陣分辨率為128像素×128像素，掃描層數(shù)為53層，層厚2 mm。彌散加權數(shù)據(jù)集由31個不同梯度方向數(shù)據(jù)(彌散加權值1 000 s/mm2)和1個未加權基本數(shù)據(jù)(彌散加權值 0 s/mm2)構成。使用 Philips Diffusion Registration Tool(Release 0.4)進行圖像配準，消除偏移等影響。張量彌散模型使用傳統(tǒng)的最小二乘模型［1］，并據(jù)此計算部分各向異性(fractional anisotropy，F(xiàn)A)。

FA是衡量DTI數(shù)據(jù)各向異性大小的主要指標，可客觀描述體素彌散的大小程度。圖3(a)為上述真實人體DTI數(shù)據(jù)中第25層的 FA圖，圖3(b)和(c)為分別在原始DWI數(shù)據(jù)中加入5%和10%高斯噪聲后獲得的DTI數(shù)據(jù)第25層FA圖。由圖可知，由于噪聲影響，F(xiàn)A圖的準確性下降。

使用本方法，分別對5%和10%的高斯噪聲干擾數(shù)據(jù)去噪，獲得的FA圖見圖3(d)和(e)。分別使用Parker方法，對5%和10%的高斯噪聲干擾數(shù)據(jù)去噪，結果FA圖。見圖3(f)和(g)。由圖可知，兩種方法在不同噪聲情況下，都可以在一定程度上去噪復原圖像，相比而言，本方法可以在更大程度上去除噪聲影響。圖3(b)和(c)在胼胝體壓部周圍的非彌散區(qū)域，由于噪聲影響產(chǎn)生一些錯誤的彌散信息，使用本方法根據(jù)平滑結構張量進行不同尺度的濾波，去噪效果明顯，這些錯誤的彌散區(qū)域基本得到修正。同時，本方法在彌散區(qū)域根據(jù)白質結構自適應地進行各向異性濾波，在去除噪聲的情況下盡可能地考慮結構對其的影響，由圖3(d)和(e)可見，胼胝體壓部的結構保持良好。

較之DTI數(shù)據(jù)的大小信息，去除噪聲對DTI方向信息的影響更加重要。使用本方法，分別對加入5%和10%高斯噪聲的DTI數(shù)據(jù)去噪，表1為對方向信息的復原效果(參數(shù)設置同上)，為各數(shù)據(jù)集體素彌散張量主分量方向與標準數(shù)據(jù)集彌散張量主分量方向偏離值的均值。

表1 本方法和Parker方法真實DTI數(shù)據(jù)集濾波去噪方向信息的復原效果比較Tab.1 Variation of mean angular differencein the proposed method and Parker’s

由表1可知，較之標準 DTI數(shù)據(jù)集，被5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集其主軸平均彌散方向分別有平均5.271°和9.198°的偏差。本方法和Parker方法均能降低噪聲對彌散張量數(shù)據(jù)方向信息的影響，相比而言，本方法對方向信息復原效果更好，平均角度偏離值分別下降至2.087°和4.671°，較之Parker方法，復原結果方向信息更接近于標準數(shù)據(jù)集。

3 結論

使用偏微分方程的方法，對磁共振彌散張量成像的DWI數(shù)據(jù)進行各向異性擴散去噪。不同于現(xiàn)有方法對多方向通道的DWI數(shù)據(jù)獨立去噪，本研究綜合DTI的結構和數(shù)據(jù)特征，對各方向通道的DWI數(shù)據(jù)，按照各向異性擴散原理，重構統(tǒng)一的平滑結構特征張量，重新設計相應的特征值和特征向量。實現(xiàn)在非彌散的各向同性結構區(qū)域進行大尺度的各向同性擴散;在邊緣區(qū)域，平行于邊沿方向進行大尺度擴散去噪，垂直于邊沿方向進行小尺度擴散;最后對于腦白質帶狀纖維束內部進行相應小尺度擴散，以保持數(shù)據(jù)精度。在合成數(shù)據(jù)和真實人腦部DTI數(shù)據(jù)集上進行去噪仿真和實驗，較之現(xiàn)有方法效果更好。

附錄

R方向張量擴散公式

張量矩陣D

部分各向異性FA

角度偏差均值Ψ

FA變化率均值ζ

［1］Basser P，Mattiello J，LeBihan D.MR diffusion tensor spectroscopy and imaging［J］.Biophysical Journal，1994，66(1):259-267.

［2］Susumu M，Zhang Jiangyang.Principles of diffusion tensor primer imaging and its applications to basic neuroscience research［J］.Neuron，2006，51(5):527-539.

［3］Heemskerk AM，Sinha TK，Wilson KJ，et al.Quantitative assessment of DTI-based muscle fiber tracking and optimal tracking parameters［J］.Magn Reson Med，2009，61(2):467-72.

［4］Fechete R，Demco DE，Eliav U，et al.Self-diffusion anisotropy of water in sheep Achilles tendon［J］.NMR in Biomedicine，2005，18(8):577-586.

［5］Anderson AW.Theoretical analysis of the effects of noise on diffusion tensor imaging［J］.Magnetic Resonance in Medicine，2001，46(6):1174-1188.

［6］Parker GJ，Schnabel JA，Symms MR，et al.Nonlinear smoothing for reduction of systematic and random errors in diffusion tensor imaging［J］.J.Magn Reson Imag，2000，11:702-710.

［7］白衡，高玉蕊，王世杰，等.DTI擴散張量的一種穩(wěn)健估計方法［J］.計算機研究與發(fā)展，2008，45(7):1232-1238.

［8］Pennec X，F(xiàn)illard P，Ayache N.A Riemannian framework for tensor computing［J］.Int J Comput Vis，2006，66:41-66.

［9］Weickert J.Coherence-enhancing diffusion filtering［J］.Int J Comput Vis，1999，31:111-127.

Diffusion Tensor Magnetic Resonance Imaging Denoising Using Anisotropic Partial Differential Function Filter

WU Xi1，2*ZHOU Ji-Liu2BI Wu-Zhong1XIE Ming-Yuan1LUO Dai-Sheng2
1(Department of Electronic Engineering，Chengdu University of Information Technology，Chengdu 610225，China)2(School of Electronics and Information Engineering，Sichuan University，Chengdu 610065，China)

TP391

A

0258-8021(2010)04-0481-05

10.3969/j.issn.0258-8021.2010.04.001

2010-01-12，

2010-05-12

四川省教育廳科技創(chuàng)新重大培育項目(09ZZ004)

*通訊作者。 E-mail:xiey@cuit.edu.cn