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

吳 錫周激流畢務(wù)忠謝明元羅代升

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

T′的特征值重構(gòu)為:

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

由上述特征值和特征向量，即可重構(gòu)平滑結(jié)構(gòu)張量T′。當(dāng)位于各向同性結(jié)構(gòu)時(shí)，各方向特征值均較大，保證各方向均進(jìn)行較大的擴(kuò)散濾波;當(dāng)位于各向異性結(jié)構(gòu)時(shí)，沿邊緣方向特征值大，垂直邊緣特征值小，保證邊緣增強(qiáng)去噪;當(dāng)位于腦白質(zhì)纖維束內(nèi)部區(qū)域時(shí)，特征值均較小，使用較小的擴(kuò)散保持?jǐn)?shù)據(jù)精度。

圖1 合成DTI數(shù)據(jù)用本方法和Parker方法的去噪結(jié)果比較。(a)合成DTI三維數(shù)據(jù);(b)合成DTI二維放大結(jié)構(gòu);(c)加入5%高斯噪聲的結(jié)果;(d)加入10%高斯噪聲的結(jié)果;(e)和(f)分別加入5%和10%高斯噪聲的本方法去噪結(jié)果;(g)和(f)分別加入5%和10%高斯噪聲Parker方法去噪結(jié)果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í)驗(yàn)結(jié)果

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

首先，采用所提出的算法在合成張量數(shù)據(jù)上進(jìn)行了仿真，實(shí)驗(yàn)平臺(tái)為Matlab 7.0.4，Windows Vista操作系統(tǒng)，處理器為Intel Core2 Duo CPU(1.80 GHz)，系統(tǒng)內(nèi)存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)框中圖像的放大結(jié)果，其中的直線段表征彌散張量，其長(zhǎng)短和方向與該體素彌散張量主分量相同，非彌散區(qū)域用圓點(diǎn)表示。圖1(c)和(d)分別為上述DTI數(shù)據(jù)集加入5%和10%高斯噪聲的結(jié)果，由于噪聲干擾DTI數(shù)據(jù)的彌散大小和方向都發(fā)生了變化，隨著噪聲增加，出現(xiàn)了明顯的錯(cuò)誤 DTI數(shù)據(jù)。圖1(e)和(f)是使用本方法分別對(duì)受5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集去噪結(jié)果，圖1(g)和(h)是使用Parker方法分別對(duì)受5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集去噪結(jié)果。

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

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

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

圖3 真實(shí)腦部DTI數(shù)據(jù)用本方法和Parker方法的去噪效果比較。(a)原始FA圖;(b)加入5%高斯噪聲的FA圖;(c)加入10%高斯噪聲的FA圖;(d)和(e)用本方法加入5%和10%高斯噪聲的去噪結(jié)果;(f)和(g)用Parker方法加入5%和10%高斯噪聲的去噪結(jié)果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可知，較之標(biāo)準(zhǔn)的DTI數(shù)據(jù)集，被5%和10%高斯噪聲干擾的DTI數(shù)據(jù)集其平均彌散方向分別有大約4.5°和9.5°的偏差。使用本方法進(jìn)行去噪后，隨著迭代次數(shù)增加，平均角度偏差逐漸減小，大約在第10個(gè)迭代周期后基本達(dá)到穩(wěn)定，平均角度偏差分別下降至2°和5°。

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

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

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

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

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

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

3 結(jié)論

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

附錄

R方向張量擴(kuò)散公式

張量矩陣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擴(kuò)散張量的一種穩(wěn)健估計(jì)方法［J］.計(jì)算機(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)新重大培育項(xiàng)目(09ZZ004)

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