基于字典基函數(shù)框架的纖維方向分布模型重建

吳 燁 馮遠(yuǎn)靜 李 斐 高成鋒

(浙江工業(yè)大學(xué)信息工程學(xué)院, 杭州 310023)

基于字典基函數(shù)框架的纖維方向分布模型重建

吳 燁 馮遠(yuǎn)靜*李 斐 高成鋒

(浙江工業(yè)大學(xué)信息工程學(xué)院, 杭州 310023)

腦白質(zhì)纖維方向分布模型重建是腦纖維成像的重要過程之一, 為纖維束跟蹤提供精確的纖維方向估計(jì)。傳統(tǒng)方法的約束條件往往依賴于先驗(yàn)的纖維方向信息, 限制了計(jì)算效率與精度的提高。提出一種基于字典基函數(shù)框架的纖維方向分布函數(shù)(fODF) 估計(jì)方法, 在單位球面上建立均勻的字典基函數(shù)分布, 構(gòu)造基于字典基函數(shù)的球面反卷積模型，直接通過非負(fù)最小二乘方法求得基函數(shù)系數(shù), 避免偽峰與復(fù)雜的算法以及更高階的大規(guī)模病態(tài)逆問題的計(jì)算，使得纖維方向分布的稀疏性和非負(fù)性可以通過幾個(gè)基函數(shù)的加權(quán)和很容易地表示。模擬磁共振數(shù)據(jù)與臨床腦部數(shù)據(jù)的實(shí)驗(yàn)結(jié)果顯示，該方法在同等條件下與目前流行的Super-CSD和QBI方法相比，計(jì)算復(fù)雜度大大降低, 計(jì)算時(shí)間僅為Super-CSD的1/3; 角度分辨率擴(kuò)大到20°以內(nèi)(其中, Super-CSD為40°, QBI為60°)。由于該方法簡(jiǎn)單高效, 角度誤差不會(huì)隨角度等參數(shù)的變化而跳變，也因?yàn)榻嵌确直媛屎徒嵌日`差性能的提升, 使得纖維方向的錯(cuò)誤重構(gòu)概率降低到8%以下，可為纖維跟蹤技術(shù)提供精確的局部纖維走向。

字典基函數(shù); 球面反卷積; 纖維方向分布; 方向分布模型

引言

人腦是世界上最復(fù)雜的系統(tǒng)之一,大約由1 000億個(gè)神經(jīng)元和100萬億條神經(jīng)聯(lián)絡(luò)通路組成[1]。盡管相關(guān)的科學(xué)研究不斷深化,但諸如阿爾茲海默癥、帕金森綜合癥、孤獨(dú)癥、癲癇病、精神分裂癥以及創(chuàng)傷性腦損傷等神經(jīng)類疾病的根源仍然無法解釋。2013年初,歐盟、美國(guó)等先后啟動(dòng)了腦研究計(jì)劃,旨在從時(shí)間和空間上展示復(fù)雜的神經(jīng)環(huán)路之間的相互作用和連通關(guān)系,腦科學(xué)的研究達(dá)到了一個(gè)前所未有的高度[2]。在大腦微結(jié)構(gòu)成像中,擴(kuò)散加權(quán)磁共振成像(diffusion weighted magnetic resonance imaging,DW-MRI)基于磁共振成像(magnetic resonance imaging,MRI)技術(shù),是目前唯一有效的腦部結(jié)構(gòu)無創(chuàng)檢測(cè)技術(shù)。DW-MRI在常規(guī)MRI序列掃描的基礎(chǔ)上,在Cartesian坐標(biāo)系的3個(gè)正交方向施加兩組巨大的對(duì)稱梯度脈沖信號(hào),用來反映水分子的擴(kuò)散運(yùn)動(dòng)[3]。腦纖維成像技術(shù)以DW-MRI為基礎(chǔ),旨在準(zhǔn)確地對(duì)大腦功能區(qū)的纖維束分布進(jìn)行重構(gòu)及成像分析。

目前,纖維跟蹤技術(shù)是唯一的一種能在活體無創(chuàng)條件下獲得腦纖維走向的方法,其基本思想是利用大腦體素內(nèi)水分子擴(kuò)散的各向異性對(duì)臨床數(shù)據(jù)進(jìn)行方向分布的建模,精確計(jì)算纖維在體素中的擴(kuò)散方向,利用纖維跟蹤算法對(duì)ROI內(nèi)的纖維方向進(jìn)行連接[4]。其中,擴(kuò)散張量成像(diffusion tensor imaging,DTI)技術(shù)標(biāo)定了腦白質(zhì)水分子的擴(kuò)散運(yùn)動(dòng)狀態(tài),是目前纖維重構(gòu)中最常用的纖維方向估計(jì)模型[5-6],其基本思想是用一個(gè)二階張量來描述纖維方向分布情況。正因?yàn)槿绱?，DTI方法無法描述腦纖維存在的諸如交叉、分叉、扇形和瓶頸型等復(fù)雜結(jié)構(gòu)。高角度分辨率成像(high angular resolution diffusion imaging,HARDI)技術(shù)的出現(xiàn)，成功地彌補(bǔ)了DTI成像的不足。現(xiàn)有的HARDI方法主要包括以擴(kuò)散譜成像(diffusion spectrum imaging,DSI)[7],Q-ball成像(Q-ball imaging,QBI)[8]、高階張量(high order tensor,HOT)[9-10]、球面反卷積(spherical deconvolution,SD)[11-12]為主導(dǎo)的多纖維重建方法。廣義上講,這些方法通過計(jì)算纖維方向分布函數(shù) (fiber orientation distribution function,fODF) 建立擴(kuò)散模型。其中,DSI方法將擴(kuò)散函數(shù)(水分子的位移概率密度函數(shù))表示為擴(kuò)散信號(hào)在Q空間的傅里葉變換,通過擴(kuò)散信號(hào)直接計(jì)算得到擴(kuò)散函數(shù)，但是對(duì)采樣的敏感程度和采樣方向數(shù)要求較高,難以應(yīng)用于實(shí)際臨床中; Q-ball成像是一種相對(duì)于DSI來講比較成熟的多纖維成像方法,通過測(cè)量擴(kuò)散信號(hào)的FRT變換(funk radon transform,FRT)得到水分子的位移概率密度函數(shù)與零階貝塞爾函數(shù)加權(quán)卷積后的徑向積分[13],雖然成功地降低了采樣時(shí)間,但對(duì)纖維方向的估計(jì)存在一定的隨機(jī)性,這是因?yàn)镼BI得到的擴(kuò)散函數(shù)被一個(gè)零階貝塞爾函數(shù)調(diào)整過,擴(kuò)散函數(shù)極值頻譜拓寬,纖維方向不明確,水分子的位移概率密度函數(shù)與纖維方向之間的關(guān)系并沒有得到證實(shí),所以得到的纖維方向并不一定與真實(shí)纖維方向一致; SD模型將擴(kuò)散信號(hào)表示為纖維方向分布函數(shù)與單條纖維響應(yīng)函數(shù)的卷積[11]，不需要任何的先驗(yàn)信息就能夠準(zhǔn)確地重構(gòu)出腦白質(zhì)纖維走向,但該模型直接從原始信號(hào)中通過反卷積得到,必然涉及高度病態(tài)性的高維逆問題計(jì)算,所以對(duì)噪聲較為敏感; CSD技術(shù)[14]在原有的球面卷積方法中引入了迭代非負(fù)約束的優(yōu)化過程,提高了角度纖維識(shí)別能力,抗噪性顯著提高,但容易產(chǎn)生偽峰現(xiàn)象[15],且隨著階數(shù)的提高重建的穩(wěn)定性大大降低,因此無法穩(wěn)定重建小角度的纖維分叉。

近年來，隨著信息技術(shù)進(jìn)入大數(shù)據(jù)時(shí)代,磁共振數(shù)據(jù)也進(jìn)入了高維高密度層次。雖然高精度的磁共振數(shù)據(jù)為精確重構(gòu)復(fù)雜的腦纖維網(wǎng)絡(luò)帶來了可能,但同時(shí)也帶來了數(shù)據(jù)分析上更大的挑戰(zhàn)[16],纖維方向分布的群體建模過程變得相當(dāng)復(fù)雜。傳統(tǒng)的用提高階數(shù)和復(fù)雜的數(shù)值優(yōu)化算法來逼近纖維分叉角度的精確識(shí)別,顯然無法滿足臨床應(yīng)用的需要。筆者提出了一種基于字典基函數(shù)框架的纖維方向分布重建方法,建立了字典基函數(shù)的球面反卷積模型。用球面字典基的簡(jiǎn)單加權(quán)形式，避免了偽峰與復(fù)雜的算法,避免了更高階的大規(guī)模病態(tài)逆問題的計(jì)算。在單位球空間上，建立基函數(shù)字典來搜索纖維方向的最優(yōu)稀疏分布,以更小的計(jì)算量和更精確的角度分辨率為纖維跟蹤技術(shù)提供纖維的精確走向。該方法分別在模擬磁共振數(shù)據(jù)與臨床腦部數(shù)據(jù)中驗(yàn)證了可行性與高效性。

1 字典基函數(shù)框架的纖維方向分析

1.1球面反卷積模型

(1)

由于受數(shù)據(jù)采集噪聲大等因素的影響,纖維方向分布帶有一定的擴(kuò)散性。纖維方向分布的計(jì)算依賴于大規(guī)模的病態(tài)逆問題。為了使重建的纖維方向分布更接近于真實(shí)的纖維方向,減少這種大規(guī)模的逆問題對(duì)重建的影響,筆者提出一種利用字典基函數(shù)框架來實(shí)現(xiàn)纖維方向分布的重建方法。

1.2方向基函數(shù)框架的纖維方向分析

單位時(shí)間內(nèi),水分子沿纖維方向的各向異性擴(kuò)散呈現(xiàn)近紡錘體分布[19],通過基函數(shù)擬合離散fODF，可以直觀地表示纖維方向的擴(kuò)散分布特性。但是，纖維的線性走向?qū)λ肿訑U(kuò)散各向異性存在約束限制,使得真實(shí)纖維方向信息多數(shù)僅集中于fODF極值的鄰域范圍。2012年，在MICCAI會(huì)議中，Merlet等將該現(xiàn)象稱為“纖維方向的稀疏性分布”[20]。目前,盡可能地提高有效纖維方向信息的比例,是構(gòu)建高分辨率纖維方向分布的途徑之一。

假設(shè)式(1)中所有非共線的纖維之間相對(duì)獨(dú)立,那么可以將體素內(nèi)的所有纖維主方向fri表示為

(2)

為了表示纖維在球面每個(gè)方向上的可能性分布,建立了一種預(yù)先定義的基函數(shù)字典(f1,f2,…,fm)近似模型(見圖1)。那么，纖維方向分布模型可以表示為一組基函數(shù)字典(f1,f2,…,fm)的線性加權(quán)組合，即

(3)

式中,ω=[ω1,ω2,…,ωm]T是加權(quán)向量。

字典基函數(shù)是構(gòu)成纖維方向的所有可能性分布的一組基。通過在球面上均勻采樣m個(gè)離散點(diǎn)，可以得到均勻分布的單位向量υ，那么所建立的基函數(shù)字典均勻分布在整個(gè)Q空間上。

圖1 一種字典基函數(shù)在單位球面空間的方向分布Fig.1 The orientation distribution of a dictionary basis funtion on a unit sphere

1.3字典基函數(shù)的建立

在數(shù)學(xué)分析中,正矢函數(shù)通過緯度來精確計(jì)算球面點(diǎn)距離,重新定義了余弦函數(shù)的球坐標(biāo)形式,常用于航海中,尤其在小角度和短距離的計(jì)算上尤為適用。正矢函數(shù)可定義為

vercosin(θ)=1+cosωθ(θ∈[0,2π])

(4)

式中，ω∈+是周期參數(shù),vercosin(θ)∈[0,2]在極坐標(biāo)下是一個(gè)原點(diǎn)對(duì)稱函數(shù)。

?(γ,υ))

(α∈+,?∈[0,π]

(5)

式中，?(γi,υj)∈[0,π]是單位重建向量γi∈S2與單位基準(zhǔn)向量υj∈S2的相對(duì)夾角α∈+是一個(gè)微調(diào)參數(shù)。

為了提高纖維方向分布的極值比例,使其具有更強(qiáng)的各向異性表現(xiàn)能力,拓展了函數(shù)的表示形式，有

Γμ(γ,υ)=Ζ·(1+αcos2?(γ,υ))μ

(α,μ∈+?∈[0,π])

(6)

式中，Ζ是歸一化參數(shù)。

Γμ(γ,υ)中小于1的方向?qū)⑼ㄟ^冪運(yùn)算而衰減,這提高了模型的徑向緊湊程度,使得各重建點(diǎn)的fODF值比例更加鮮明(見圖2)。

圖2 正矢函數(shù)在球坐標(biāo)下的分布。(a)μ=1；(b)μ=5；(c)μ=20Fig.2 The distribution of versine function on spherical coordinate. (a)μ=1；(b)μ=5；(c)μ=20

2 算法實(shí)現(xiàn)

2.1字典基函數(shù)框架的球面反卷積模型

在單位半球面上均勻采樣m個(gè)重建向量υj。令γ←υ， 那么方向分布函數(shù)可以表示為

(7)

沿梯度方向,g的擴(kuò)散衰減信號(hào)在方向υ上可以表示為

(8)

式中,系數(shù)ω是唯一的未知量,那么式(8)可以表示為線性方程,即

S=Aω+η

(9)

因?yàn)橹亟ㄏ蛄康牟蓸訑?shù)n遠(yuǎn)高于磁共振數(shù)據(jù)的梯度方向的采樣數(shù)m，假定問題式(9)有無窮多解,然而所尋求的是在可行范圍內(nèi)符合纖維方向分布的最優(yōu)解。此外,為了減小噪聲η對(duì)重建的影響,可以將式(9)轉(zhuǎn)化為能量方程來求解，即

(10)

2.2數(shù)學(xué)模型求解及特性分析

纖維沿?cái)U(kuò)散方向的稀疏分布和纖維方向的非負(fù)擴(kuò)散是兩大不可否認(rèn)的物理性質(zhì)。因此,有理由相信: 盡可能地平衡纖維方向的稀疏性與非負(fù)性的相對(duì)約束，是實(shí)現(xiàn)高角度分辨率纖維交叉的條件。在球坐標(biāo)下，建立了纖維方向的字典空間,通過少量的加權(quán)系數(shù)來重建纖維方向的擴(kuò)散模型(見圖3)，避免求解時(shí)的強(qiáng)制稀疏約束?？紤]到纖維方向的非負(fù)性約束,通過非負(fù)最小二乘法，得到與纖維主方向最接近的字典基函數(shù)，有為了提高fODF計(jì)算的精確性,通常會(huì)遍歷整個(gè)球面進(jìn)行球面向量的積分運(yùn)算。但是，對(duì)連續(xù)球面充分采樣是不可能的，且連續(xù)球面的積分過于復(fù)雜。筆者采用一種離散近似的方法,對(duì)單位正20面體進(jìn)行?次細(xì)分，獲得(10×4?+2)個(gè)重建向量(見圖4)。細(xì)分次數(shù)越高，離散誤差越小，但運(yùn)算時(shí)間將延長(zhǎng)。離散誤差雖無法忽略，但通常限制在一定范圍內(nèi)的角度誤差是允許的。例如，四次分形獲得2 562個(gè)單位重建向量κi，令γ←κ,有

(11)

圖3 模擬50°纖維分叉的fODF重構(gòu),在每個(gè)采樣方向上的字典基函數(shù)所對(duì)應(yīng)的加權(quán)系數(shù)ω的分布情況,由于研究的均勻采樣是無序采樣,因此橫坐標(biāo)僅表示每個(gè)采樣方向?qū)?yīng)于采樣方向組中的位置。 (a)μ=50; (b)μ=100Fig.3 The weighted coefficients ω of corresponding basis function along each sampling direction for simulated 50° crossing fiber. Abscissa indicates the position of each sampling direction in the whole sampling group for the unordered uniform sampling. (a)μ=50; (b)μ=100

(12)

式中，ωi是式(11)的解。

通過靈活的基函數(shù)來提高纖維方向信息的比例,從而消除不必要的纖維存在,是一種間接但高效的方法,稱之為字典稀疏性。利用字典基函數(shù)框架建立纖維方向分布模型,可以避免對(duì)復(fù)雜數(shù)值優(yōu)化算法的需求,從而大幅降低計(jì)算時(shí)間。

圖4 模擬50°纖維分叉重構(gòu)。(a)2 562個(gè)均勻采樣的重建方向中重構(gòu)fODF的結(jié)果,白線為真實(shí)纖維方向; (b)重構(gòu)的fODF的極值方向Fig.4 Simulate the reconstruction of fiber cross of 50°. (a) The fODF reconstruction from 2 562 uniform sampling directions,the white line is the real fiber direction; (b)The extremum direction of fODF

2.3實(shí)驗(yàn)方法

根據(jù)S?derman提出的方法[21]產(chǎn)生模擬加權(quán)磁共振信號(hào)，有

(13)

在模擬數(shù)據(jù)實(shí)驗(yàn)中,信號(hào)S通過以下方式添加Rician的隨機(jī)噪聲，即

(14)

式中，ζ1,ζ2～N(0,δ2)，δ=S0/SNR。在實(shí)驗(yàn)中,每次計(jì)算100次取平均結(jié)果。

臨床數(shù)據(jù)實(shí)驗(yàn)利用3-T GE系統(tǒng)從51個(gè)梯度方向采集真實(shí)人腦數(shù)據(jù),其中包括50個(gè)軸向切片,每個(gè)切片層的體素分辨率為256像素×256像素,b=2 000 s/mm2。在臨床數(shù)據(jù)的實(shí)驗(yàn)中,主要關(guān)注的是重構(gòu)的角度分辨率與重構(gòu)的效率。

方法的實(shí)現(xiàn)主要分為6個(gè)步驟。

步驟1:在單位半球上均勻采樣321個(gè)單位向量υ∈321×3；

步驟2:構(gòu)建單位球上均勻分布的字典基函數(shù)；

步驟3:按本文2.1節(jié)和2.2節(jié)求解體素的基函數(shù)系數(shù)；

步驟4:4次分型單位正20邊形,采樣2 562個(gè)頂點(diǎn)向量κ∈2562×3；

步驟5:按本文2.2節(jié)計(jì)算纖維方向分布F∈2562；

步驟6:按每個(gè)fODF對(duì)應(yīng)的向量κ進(jìn)行擬合,得到每個(gè)體素的纖維方向分布模型。

2.4評(píng)價(jià)標(biāo)準(zhǔn)

本研究的目標(biāo)之一是提高纖維方向的重建正確性，文獻(xiàn)[22-23]中廣泛使用了以下度量指標(biāo)來評(píng)估纖維方向重建的效率與纖維方向的準(zhǔn)確性。

1)錯(cuò)誤纖維方向的概率: 量化一個(gè)體素內(nèi)實(shí)際的纖維束數(shù)量M的先對(duì)重構(gòu)準(zhǔn)確度，有

(15)

為了真實(shí)反映算法的有效性，通常對(duì)同一模擬數(shù)據(jù)多次計(jì)算求平均值，或在某一真實(shí)區(qū)域內(nèi)求得平均體素的Pd。

2)角度誤差: 量化體素中重構(gòu)的纖維方向角的精度，角度誤差可計(jì)算為

(16)

最終的值是每一個(gè)重構(gòu)的纖維方向與標(biāo)準(zhǔn)方向?qū)Ρ群蟮玫降恼`差平均值。其中，極值方向的檢測(cè)是通過一個(gè)局部最大值搜索過程，僅考慮纖維方向附近在一個(gè)15°圓錐范圍內(nèi)的方向,因?yàn)檫@種評(píng)價(jià)標(biāo)準(zhǔn)對(duì)于小的交叉角度并不敏感，通常在30°～90°的范圍內(nèi)進(jìn)行評(píng)價(jià)。另外，fODF值小于最大值得10%的識(shí)別為偽方向,要將其過濾。

3 結(jié)果

3.1模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果

纖維方向分布函數(shù)值可以看作是每一個(gè)重建方向上纖維方向分布的可能性,在真實(shí)纖維方向分布的主方向上fODF值最大。圖5顯示，當(dāng)纖維交叉角度很小時(shí)，fODF的極值方向通過指數(shù)μ的提高將逐漸分離;當(dāng)纖維交叉角度較大時(shí)，fODF的極值將通過μ的提高而逼近纖維的真實(shí)方向，這能夠有效地減小角度誤差。從圖5(a)中可以看出,當(dāng)交叉纖維無法重構(gòu)時(shí),fODF將集中于兩個(gè)方向的中間,這在纖維跟蹤時(shí)很容易被誤判為單纖維,因此對(duì)基函數(shù)加冪將使得基函數(shù)更細(xì)致化,使得在加權(quán)重構(gòu)時(shí)能夠更精確地獲得纖維方向分布模型。

圖5 fODF在單位球面上的映射分布,彩色部分(從藍(lán)色到紅色)表示fODF的映射值(從0到1),白線表示纖維的真實(shí)方向。 (a)交叉角=15°,μ=10; (b) 交叉角=15°,μ=20; (c) 交叉角=15°,μ=50; (d) 交叉角=90°,μ=10; (e) 交叉角=90°,μ=20; (f) 交叉角=90°,μ=50Fig.5 fODF mapped on unit sphere,where the color part (from blue to red) indicates the mapping of fODF value (from 0 to 1),white line represents the real direction. (a)Cross=15°, μ=10; (b)Cross=15°, μ=20; (c)Cross=15°, μ=50.(d)Cross=90°, μ=10; (e)Cross=90°, μ=20; (f)Cross=90°, μ=50

準(zhǔn)確重建fODF是纖維跟蹤的前提,而估計(jì)的纖維方向通常因?yàn)檎`差和計(jì)算精度等因素難以達(dá)到理想的纖維方向分布,如圖5紅色中心的鄰域周圍,這樣的誤差在很小范圍內(nèi)是允許的。在圖6中，分別對(duì)不同角度的纖維交叉重構(gòu),并計(jì)算了角度誤差(見圖7)。其中，Super-CSD[24]是Tournier[14]在CSD方法基礎(chǔ)上不斷優(yōu)化而形成的一套比較完善的球面反卷積方法,是目前效果最好的方法之一,并廣泛用于纖維跟蹤技術(shù)中; QBI方法具有較高的計(jì)算效率和較低的重建誤差,具有廣泛的應(yīng)用。

實(shí)驗(yàn)顯示，DB-SD在角度分辨率上高于其他兩種方法,對(duì)于35°以下的小角度纖維也能穩(wěn)定地重建,更重要的是所有的角度誤差均能夠保持在6°以下。隨著μ的提高，角度誤差在很小的范圍內(nèi)產(chǎn)生波動(dòng),這是因?yàn)槔w維信息比例的不斷增加使得角度誤差也愈加明顯,但仍然能夠保證在很小的角度誤差范圍內(nèi)。圖6 (d)顯示，由于階數(shù)的提高、待求解的未知系數(shù)增加，使得Super-CSD方法變得不穩(wěn)定；對(duì)于DB-SD,由于字典基函數(shù)加權(quán)系數(shù)的稀疏性,只需要少量的未知系數(shù)就能夠完整地重構(gòu),因此降低了因病態(tài)問題而引起的重構(gòu)不穩(wěn)定性。

圖6 模擬數(shù)據(jù)fODF的重建結(jié)果,其中SNR=40,從左到右依次表示模擬纖維交叉角度從15°～85°(間隔10°)。(a)DB-SD μ=20; (b)DB-SD μ=50; (c)Super-CSD 階數(shù)=12; (d) Super-CSD階數(shù)=16; (e)QBI階數(shù)=12; (f) QBI 階數(shù)=16Fig.6 The reconstruction results of fODF of simulated data, where SNR=40, from left to right were represented the fiber cross angle from 15 degree to 85 degree(interval 10 degree). (a)DB-SD μ=20; (b)DB-SD μ=50; (c)Super-CSD order=12; (d) Super-CSD order=16; (e)QBI order=12; (f)QBI order=16

圖7 模擬數(shù)據(jù)纖維方向分布重建結(jié)果的角度誤差(其中SNR=40)Fig.7 The angular error of the fiber orientation distribution reconstruction from simulated data. SNR=40

圖8是一組不同信噪比下的纖維方向錯(cuò)誤重構(gòu)概率的對(duì)比實(shí)驗(yàn)，其中Super-CSD和Qball方法的階數(shù)設(shè)置為14階; 實(shí)驗(yàn)中的DB-SD方法設(shè)置μ=20。實(shí)驗(yàn)分別對(duì)同一角度測(cè)試100次，累計(jì)這100次的重構(gòu)結(jié)果取平均值測(cè)試。實(shí)驗(yàn)顯示,本方法具有穩(wěn)定高效的纖維方向重構(gòu)效果,這是因?yàn)镈B-SD方法在提高角度分辨率的同時(shí)并未增加系數(shù)計(jì)算的復(fù)雜度,可以大大提高重建的穩(wěn)定性,在不同SNR下均能保持較高的重建效率。圖9是一組不同敏感系數(shù)和不同噪聲共同影響下的各角度重建誤差的實(shí)驗(yàn)。實(shí)驗(yàn)顯示，當(dāng)敏感系數(shù)較低時(shí)(見圖9中的(a)和(b)),信號(hào)本身重構(gòu)的難度增加,使得Q-ball方法的角度誤差普遍較大,這種情況下纖維的重建效果均不佳,纖維重建的結(jié)果并不可信; Super-CSD方法在纖維交叉角度低于45°時(shí)角度誤差較大,隨著角度的增加反映出了誤差逐漸下降,隨著誤差的增加也反映出了Super-CSD具有一定的抗噪性,但在100次的同等條件重建下,穩(wěn)定性仍然較低; DB-SD方法雖然在圖9(a)的35°時(shí)反映出了角度誤差約為10°,但總體角度誤差均較低,隨著噪聲等級(jí)的變化,依然能夠保證穩(wěn)定性與較低的角度誤差。當(dāng)敏感系數(shù)較高時(shí)(見圖9中(c)和(d)),Q-ball方法在交叉角度較大時(shí)重建的穩(wěn)定性較好; Super-CSD方法雖然在角度誤差上與DB-SD方法相當(dāng),但穩(wěn)定性卻弱于本方法。事實(shí)上,敏感稀疏與噪聲是臨床數(shù)據(jù)應(yīng)用時(shí)重要的影響參數(shù)之一,既要保證不同參數(shù)下重建的相對(duì)穩(wěn)定，更要保證誤差較低。

圖8 模擬數(shù)據(jù)下不同信噪比的纖維方向的錯(cuò)誤重構(gòu)概率。(a) SNR=10; (b) SNR=40Fig.8 The probability of false fiber detection of simulated reconstruction with different SNR. (a) SNR=10; (b) SNR=40

圖9 模擬數(shù)據(jù)下不同信噪比與敏感系數(shù)影響的纖維方向分布重建結(jié)果的角度誤差的均值和標(biāo)準(zhǔn)差。(a) b=1 000 s/mm2, SNR=10; (b) b=1 000 s/mm2, SNR=40; (c) b=3 000 s/mm2, SNR=10; (d) b=3 000 s/mm2, SNR=40Fig.9 The mean and standard deviation of the angular error in orientation reconstruction with different SNR and sensitivity coefficient from simulated data. (a) b=1 000 s/mm2, SNR=10; (b) b=1 000 s/mm2, SNR=40; (c) b=3 000 s/mm2, SNR=10; (d) b=3 000 s/mm2, SNR=40

圖10 扣帶與扣帶回皮質(zhì)相交區(qū)域的纖維方向重建實(shí)驗(yàn)對(duì)比,其中顏色表示纖維主方向在RGB規(guī)則上的映射。 (a)Super-CSD方法,球諧階數(shù)=12; (b)DB-SD方法,μ=20; (c)DB-SD方法,μ=50; (d)人腦冠狀面切片的DW-MRI圖,灰度顏色映射為FA值的大小,其中白色表示各向異性較大的區(qū)域Fig.10 The reconstruction of the intersection redion of cigulum and cingulate cortex with different methods,where the color representation mapping fiber main direction in RGB rules. (a)Super-CSD method,where SH order=12; (b)DB-SD method, where μ=20; (c)DB-SD method, where μ=50; (d) the human brain coronal sections in the lower left corner, gray value indicates the size of FA values, and white part said the area where a large anisotropy

圖11 下胼胝體與下扣帶相交區(qū)域的纖維方向重建實(shí)驗(yàn)對(duì)比,其中顏色表示纖維主方向在RGB規(guī)則上的映射。 (a)Super-CSD方法,球諧階數(shù)=12; (b)DB-SD方法,μ=20; (c)人腦冠狀面切片的DW-MRI圖,灰度顏色映射為FA值的大小,其中白色表示各向異性較大的區(qū)域Fig.11 The comparison of fiber direction reconstruction in the intersection region of the corpus callosun and cingulun,where the color representation mapping fiber main direction in RGB rules, ROI schematic representation the human brain coronal sections, gray value indicates the size of FA values. (a)Super-CSD method, where SH order=12 (b)DB-SD method, where μ=20;(c) The human brain coronal sections in the lower left corner, gray value indicates the size of FA values, and white part said the area where a large anisotropy

從本質(zhì)上來看，Super-CSD使用復(fù)雜的球諧函數(shù)作為基函數(shù),利用Tikhonov正則化方法來誘導(dǎo)負(fù)方向的消除。從計(jì)算效率來看,筆者使用簡(jiǎn)單的球面正矢函數(shù)作為基函數(shù),并未使用復(fù)雜的數(shù)值逼近算法。Super-CSD的角度分辨率的提高依賴于階數(shù)的增加,但從實(shí)驗(yàn)中可以看出,Super-CSD方法在高階非常不穩(wěn)定; 而QBI方法的角度分辨率較低,無法很好地重構(gòu)60°以下的纖維交叉,隨著階數(shù)的增加,也并未改善其重建精度。因此，從角度分辨率和計(jì)算效率來說，本方法都優(yōu)于目前常用的Super-CSD方法和QBI方法,具有很好的臨床應(yīng)用前景。

3.2臨床數(shù)據(jù)實(shí)驗(yàn)結(jié)果

由于QBI方法的角度分辨率低等原因，實(shí)際數(shù)據(jù)僅與重建效果較好的Super-CSD方法對(duì)比。給出了DB-SD方法與Super-CSD方法在分?jǐn)?shù)各向異性(fractional anisotropy,FA)背景下的纖維束交叉區(qū)域矢狀切面的纖維方向重建實(shí)驗(yàn)結(jié)果。其中,圖10的ROI區(qū)域?yàn)樯峡蹘c扣帶回皮質(zhì)的相交區(qū)域附近,圖11的ROI區(qū)域?yàn)橄码蓦阵w與下扣帶的相交區(qū)域附近。

由圖10中扣帶與扣帶回皮質(zhì)相交區(qū)域的多纖維交叉可見,不同方向的多纖維經(jīng)過區(qū)域,DB-SD方法都能夠很好地識(shí)別，尤其在綠色部分,纖維的方向保持一致,更符合纖維的真實(shí)分布；在圖11中，橫向扣帶穿過胼胝體纖維,DB-SD方法很好地還原了這一事實(shí),并且在下胼胝體與下扣帶相交區(qū)域,DB-SD方法表現(xiàn)出了高角度分辨率的優(yōu)勢(shì)。另外,隨著基函數(shù)的冪μ提高,重構(gòu)的角度誤差降低，角度分辨率提高,然而增加的計(jì)算量是否會(huì)提高重構(gòu)的時(shí)間,筆者做了一組對(duì)比實(shí)驗(yàn)。其中，?是采樣重建向量的分形次數(shù)(見表1)。實(shí)驗(yàn)分別對(duì)其中一層的256體素×256體素進(jìn)行實(shí)驗(yàn),計(jì)算單個(gè)體素重構(gòu)的平均時(shí)間(ms)。從表格可以看出,μ的提高對(duì)算法本身的重構(gòu)時(shí)間并未造成太大的影響,這是由于每個(gè)體素使用同樣的字典基函數(shù)分布,不同的是基函數(shù)所對(duì)應(yīng)的加權(quán)系數(shù),因此在對(duì)體素地群體重構(gòu)時(shí),并不需要反復(fù)計(jì)算每個(gè)體素的基函數(shù)分布,而只需一次計(jì)算。另外,由于基函數(shù)的加權(quán)稀疏的稀疏性,使得DB-SD方法通過幾個(gè)簡(jiǎn)單的基函數(shù)加權(quán)疊加就可以獲得很好的纖維方向分布。更值得一提的是，傳統(tǒng)方式提高階數(shù)增加了未知系數(shù)的維度，使得重構(gòu)的復(fù)雜度大大增加,而DB-SD方法的系數(shù)維度僅與球面基函數(shù)分布的個(gè)數(shù)相關(guān)。這樣,DB-SD相比其他方法可以大大地降低計(jì)算的復(fù)雜度,因此減少了重構(gòu)時(shí)間。那么,應(yīng)用時(shí)可以在保證穩(wěn)定的前提下，適當(dāng)?shù)靥岣擀虂碇貥?gòu)高角度分辨率的纖維。

表1每個(gè)體素中纖維方向重建的平均時(shí)間

Tab.1Themeantimeneededforfiberdirectionreconstructioninsinglevoxelwithdifferentmethod

分型次數(shù)平均重建時(shí)間/sDB-SDSuper-CSDμ=10μ=50μ=100階數(shù)=12?=313.513.413.635.7?=413.714.014.036.8?=515.615.915.943.5

4 討論和結(jié)論

近年來,隨著腦研究計(jì)劃在全球范圍內(nèi)的開展,大腦成像技術(shù)取得了巨大進(jìn)步。對(duì)于大腦微結(jié)構(gòu)和纖維群體成像而言,纖維跟蹤技術(shù)將臨床醫(yī)學(xué)與信息技術(shù)結(jié)合,旨在揭示復(fù)雜腦網(wǎng)絡(luò)的微結(jié)構(gòu)信息。隨著信息技術(shù)的大數(shù)據(jù)時(shí)代到來,給體素建模和數(shù)據(jù)分析帶來了更大的挑戰(zhàn)。纖維方向分布模型重建的精度與效率的平衡，是整個(gè)腦纖維成像的難點(diǎn)之一。針對(duì)這一問題,本研究設(shè)計(jì)了一種字典基函數(shù)框架下的纖維方向分布模型重建算法,在單位球下構(gòu)建基函數(shù)分布框架,實(shí)現(xiàn)了對(duì)不同的纖維方向分布，僅需少量的加權(quán)系數(shù)就可以快速重建。實(shí)驗(yàn)結(jié)果顯示,DB-SD方法確實(shí)能夠在角度分辨率和計(jì)算效率中展示出優(yōu)勢(shì)。目前,纖維方向分布模型重建有3種思路: 一種是通過Q空間的信號(hào)變換，還原纖維方向分布模型,如Q-ball; 另一種是通過球面反卷積方法,從信號(hào)中還原fODF,如CSD; 還有一種則是多峰函數(shù)擬合的方式。本研究的基本思路介于后兩種思路之間。首先,本研究采用了一種簡(jiǎn)單的雙葉函數(shù)，建立單位球上的字典基函數(shù)分布; 其次,本研究使用球面反卷積模型，計(jì)算字典基函數(shù)的加權(quán)系數(shù),通過簡(jiǎn)單的非負(fù)最小二乘方式，就可以獲得稀疏的字典基函數(shù)表示的纖維方向分布模型。從圖1的三維字典基函數(shù)分布可以看出,均勻的球面采樣所構(gòu)造出的字典基函數(shù)分布是均等的。為了突出纖維的主方向信息,本研究通過加冪的方式集中纖維方向(見圖5)，隨著冪的增加,fODF在單位球面上的映射分布更加集中于纖維的真實(shí)方向。雖然冪的增加帶來了一定的計(jì)算量,但DB-SD方法的加權(quán)系數(shù)是稀疏的(見圖3),這在計(jì)算中大大降低了基函數(shù)擬合的復(fù)雜度。本研究的一個(gè)特點(diǎn)是：利用簡(jiǎn)單而高效的方式，獲得高角度分辨率、低重構(gòu)誤差的纖維方向分布模型。從圖6中可以看出,DB-SD方法具有很好的穩(wěn)定性,并且對(duì)于小角度的纖維分布能夠很好地重構(gòu)。衡量纖維重構(gòu)算法的重要指標(biāo)之一是角度誤差,從圖7～9中可以看出,由于100次重構(gòu)中纖維角度產(chǎn)生了不同程度的變化,因此平均角度誤差也反映出DB-SD方法相比其他兩種方法具有更好的穩(wěn)定性。相比模擬數(shù)據(jù)的實(shí)驗(yàn)環(huán)境,實(shí)際數(shù)據(jù)具有更多的不確定因素和復(fù)雜性,如圖10、11所示。由于每個(gè)體素是獨(dú)立建模,在不考慮先驗(yàn)條件下,實(shí)際數(shù)據(jù)的高角度分辨率很難保證。但在實(shí)驗(yàn)中,DB-SD方法相對(duì)于Super-CSD方法,大大地提高了角度分辨率,更加適應(yīng)整體的纖維束走向。

同時(shí),本實(shí)驗(yàn)結(jié)果也表明:一方面,多峰函數(shù)擬合往往過于復(fù)雜，并且很難保證效率和穩(wěn)定性; 另一方面,纖維束本身具有一定的走向性約束,無論使用任何重構(gòu)方法都應(yīng)該符合整體的纖維束走向。但在圖10中,DB-SD方法在單條纖維束中也可能產(chǎn)生多纖維,這是由于實(shí)際數(shù)據(jù)中纖維往往具有不可預(yù)測(cè)的分支。

本研究提出了一種基于字典基函數(shù)框架的fODF估計(jì)方法,建立了字典基函數(shù)的球面反卷積模型。用球面字典基的簡(jiǎn)單加權(quán)形式，避免了偽峰與復(fù)雜的算法,避免了更高階的大規(guī)模病態(tài)逆問題的計(jì)算。纖維方向分布的稀疏性和非負(fù)性可以通過少量的基函數(shù)的加權(quán)來得到,實(shí)驗(yàn)測(cè)試表明，筆者提出的DB-SD方法更加逼近纖維的真實(shí)方向,角度分辨率大幅提高,以小范圍的角度誤差穩(wěn)定識(shí)別出小角度的纖維分叉。DB-SD方法在計(jì)算時(shí)間、角度分辨率和角度誤差等方面都優(yōu)于現(xiàn)有效果較好的Super-CSD方法和QBI方法。將進(jìn)一步促進(jìn)纖維跟蹤技術(shù)在臨床中應(yīng)用。

[1] 梁夏, 王金輝, 賀永. 人腦連接組研究: 腦結(jié)構(gòu)網(wǎng)絡(luò)和腦功能網(wǎng)絡(luò)[J]. 科學(xué)通報(bào), 2010 (16): 1565-1583.

[2] 朱麗君, 朱元貴, 曹河圻, 等. 全球腦研究計(jì)劃與展望[J]. 中國(guó)科學(xué)基金, 2013,27(6): 359-362.

[3] Mori S, Crain BJ, Chacko VP,etal. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[J]. Anuals of Neurology, 1999,45(2): 265-269.

[4] Mori S, van Zijl P. Fiber tracking: principles and strategies-a technical review[J]. NMR in Biomedicine, 2002,15(7-8): 468-480.

[5] Assaf Y, Pasternak O. Diffusion tensor imaging (DTI)-based white matter mapping in brain research: a review[J]. Journal of Molecular Neuroscience, 2008,34(1): 51-61.

[6] Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging[J]. Biophysical journal, 1994,66(1): 259-267.

[7] Wedeen VJ, Wang RP, Schmahmann JD,etal. Diffusion spectrum magnetic resonance imaging (DSI) tractography of crossing fibers[J].Neuroimage, 2008,41(4): 1267-1277.

[8] Cheng J, Ghosh A, Deriche R,etal. Model-free, regularized, fast, and robust analytical orientation distribution function estimation[C]//Proceedings of the 13th international conference on Medical image computing and computer-assisted intervention: Part I. Berlin: Springer-Verlag, 2010: 648-656.

[9] 李蓉, 馮遠(yuǎn)靜, 邵開來, 等. 磁共振擴(kuò)散高階張量成像的腦白質(zhì)纖維微結(jié)構(gòu)模型及特征提取算法[J]. 中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào), 2012,31(3): 365-373.

[10] Hlawitschka M, Scheuermann G. Tracking lines in higher order tensor fields[J]. Scientific Visualization: Advanced Concepts, 2010, 1: 124-135.

[11] Kaden E, Kn?sche TR, Anwander A. Parametric spherical deconvo lution: Inferring anatomical connectivity using diffusion MR imaging[J]. NeuroImage, 2007,37(2): 474-488.

[12] Jeurissen B, Leemans A, Tournier JD,etal. Estimating the number of fiber orientations in diffusion MRI voxels: a constrained spherical deconvolution study[C]//International Society Magnetic Resonance in Medicine. Stockholm: Wiley Online Library, 2010: 573.

[13] Tuch DS. Q-ball imaging[J]. Magnetic Resonance in Medicine, 2004,52(6): 1358-1372.

[14] Tournier J, Calamante F, Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution[J]. NeuroImage, 2007,35(4): 1459-1472.

[15] Parker GD, Marshall D, Rosin PL,etal. A pitfall in the reconstruct tion of fibre ODFs using spherical deconvolution of diffusion MRI data[J]. Neuroimage, 2013,65: 433-448.

[16] 張艷, 宋志堅(jiān). 腦白質(zhì)纖維束跟蹤算法的研究進(jìn)展[J]. 復(fù)旦學(xué)報(bào)(醫(yī)學(xué)版), 2014,41(1): 1-7.

[17] 吳占雄, 朱善安. 基于移動(dòng)最小二乘法的白質(zhì)纖維束走向跟蹤[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版), 2011,45(3): 458-461.

[18] Tournier J, Calamante F, Gadian DG,etal. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution[J]. NeuroImage, 2004,23(3): 1176-1185.

[19] Frank LR. Characterization of anisotropy in high angular resolution diffusion-weighted MRI[J]. Magnetic Resonance in Medicine, 2002,47(6): 1083-1099.

[20] Merlet S, Caruyer E, Deriche R. Parametric dictionary learning for modeling EAP and ODF in diffusion MRI[C]//Medical Image Computing and Computer Assisted Intervention (MICCAI).Beilin: Springer, 2012, 7512: 8 p.

[21] S?derman O, J?nsson B. Restricted diffusion in cylindrical geome try[J]. Journal of Magnetic Resonance, Series A, 1995,117(1): 94-97.

[22] Landman BA, Bogovic JA, Wan H,etal. Resolution of crossing fibers with constrained compressed sensing using diffusion tensor MRI[J]. NeuroImage, 2012,59(3): 2175-2186.

[23] Michailovich O, Rathi Y, Dolui S. Spatially regularized compressed sensing for high angular resolution diffusion imaging[J]. IEEE Transactions on Medical Imaging, 2011,30(5): 1100-1115.

[24] Calamante F, Tournier JD, Jackson GD,etal. Track-density imaging (TDI): super-resolution white matter imaging using whole-brain track-density mapping[J]. Neuroimage, 2010,53(4): 1233-1243.

A Novel Fiber Orientation Distribution Reconstruction Method Based on Dictionary Basis Function Framework

Wu Ye Feng Yuanjing*Li Fei Gao Chengfeng

(College of Information Engineering, Zhejiang University of Technology, Hangzhou 310023, China)

Brain WM fiber orientation reconstruction is the throat of fiber imaging that set the premise of fiber tractgraphy. In conventional approaches, constraint conditions are dependent on prior information of fiber orientation, limiting the improvement of computation efficiency and precision. In this work, we presented a dictionary basic function based fiber orientation distribution function (fODF) estimation method. First, the uniform distribution of dictionary basis functions on unit sphere was established; second, spherical deconvolution (SD) model based on the dictionary basis functions was established; third, the coefficients of basic functions were directly obtained using non-negative least square method. The proposed method successfully avoided the spurious peak, leases computation burden and eliminates higher order ill-posed inverse problem, which makes the sparsity and non-negativity needed for the representation straightforwardly represented by weighted sums of basic functions. Experiments on simulated and clinical brain data both verified that computational complexity of the proposed method is much less than that of Super-CSD and QBI under the same conditions, such as, computational time of the method is a third as long as that of Super-CSD; the angular resolution of the method increases to a range of 20°, while Super-CSD is 40° and QBI is 60°. Because the method is simple and efficient, angular error does not change with such parameters as angle, and with the addition of the improvement of angular resolution and angular error performance, the probability of error fiber orientation reconstruction decreases to less than 8%, which will provide accurate local fiber orientation for fiber tractgraphy.

dictionary basis function; spherical deconvolution; fiber orientation distribution; orientation distribution model

10.3969/j.issn.0258-8021. 2015. 03.006

2014-07-02， 錄用日期:2015-04-07

國(guó)家自然科學(xué)基金(61379020); 浙江省自然科學(xué)基金(LY13F030007)

TP391

A

0258-8021(2015) 03-0297-011

# 中國(guó)生物醫(yī)學(xué)工程學(xué)會(huì)會(huì)員(Member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)， E-mail: fyjing@zjut.edu.cn