一種非負(fù)高階張量擬牛頓搜索的纖維方向分布估計(jì)方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及醫(yī)學(xué)成像、圖像處理、數(shù)值分析、三維重建、計(jì)算機(jī)科學(xué)、神經(jīng)解剖學(xué)等 領(lǐng)域，尤其是一種用于腦白質(zhì)纖維成像的非負(fù)纖維方向分布的高階笛卡爾擴(kuò)散張量成像方 法。
【背景技術(shù)】
[0002] 腦白質(zhì)纖維成像是一種非侵入式的獲取人腦白質(zhì)區(qū)域的體素神經(jīng)纖維微結(jié)構(gòu)信 息并進(jìn)行三維重構(gòu)和展示的信息技術(shù)。其主要方法是對原始的擴(kuò)散加權(quán)磁共振數(shù)據(jù)進(jìn)行體 素建模，獲得每個(gè)體素內(nèi)的纖維方向分布情況，形成具有解剖學(xué)意義的纖維空間微結(jié)構(gòu)。通 過大量研究，該領(lǐng)域的學(xué)者在腦纖維微結(jié)構(gòu)重構(gòu)算法、不確定信息處理等方面已取得了一 系列成果。然而，當(dāng)前流行的纖維重構(gòu)算法大都無法保證所求得的纖維方向分布函數(shù)的非 負(fù)性。這顯然不符合擴(kuò)散信號(hào)的實(shí)際物理意義，同時(shí)，它還會(huì)在一定程度上促成偽峰的形 成，從而影響纖維成像的結(jié)果，極大地阻礙了該技術(shù)在臨床醫(yī)學(xué)中的應(yīng)用。

【發(fā)明內(nèi)容】

[0003] 為了克服現(xiàn)有纖維成像方法無法保證擴(kuò)散函數(shù)非負(fù)性的問題，本發(fā)明提出一種以 高階張量為導(dǎo)向的基于擬牛頓搜索的非負(fù)纖維方向分布估計(jì)方法。
[0004] 本發(fā)明所采用的技術(shù)方案如下所示：
[0005] -種非負(fù)高階張量擬牛頓搜索的纖維方向分布估計(jì)方法，所述的纖維方向分布估 計(jì)方法包括以下步驟：
[0006] (1)數(shù)據(jù)預(yù)處理:讀取腦部擴(kuò)散加權(quán)磁共振數(shù)據(jù)，獲取施加梯度方向g時(shí)的磁共振 信號(hào)S(g)和未施加梯度方向時(shí)的磁共振信號(hào)So,以及相應(yīng)的梯度方向數(shù)據(jù)，選取感興趣區(qū) 域，并計(jì)算該區(qū)域的擴(kuò)散衰減信號(hào)S(g)/S〇;
[0007] (2)將感興趣區(qū)域內(nèi)每個(gè)體素中的擴(kuò)散衰減信號(hào)逐個(gè)建模為具有擴(kuò)散形態(tài)的橢球 分布模型，建模過程如下：
[0008] 2.1)體素微結(jié)構(gòu)建模:將擴(kuò)散衰減信號(hào)假設(shè)為沿重建向量v的單條纖維信號(hào)響應(yīng) 函數(shù)R (v，g)與擴(kuò)散函數(shù)D (v)在球面上的卷積：
[0009]
[0010] 其中，糾,.4) = 1?-近似看作一個(gè)高斯分布函數(shù)，g= {giERlx31 i = l，· · ·，n}為梯 度方向，v={vPeRlx3|p = l，. . .，K}為在單位球面上采樣的重建向量，R1X3表示維度為1X3 的實(shí)數(shù)域矩陣空間，η和K分別表示梯度方向和重建向量的個(gè)數(shù)，μ = ε b是表征擴(kuò)散效率ε與 擴(kuò)散敏感系數(shù)b共同影響的一個(gè)參數(shù)，擴(kuò)散函數(shù)D(v)的表達(dá)式如下：
[0011]
[0012 ] drs表示單項(xiàng)式v;V:. V丨ts的系數(shù)，1為高階張量的階數(shù)，r，s分別表示重建向量v = (Vx， Vy，VZ)的基方向VX和Vy的指數(shù)，表示第j個(gè)張量的張量系數(shù)，j = l，. . .，m，m表示張量的個(gè)數(shù)
，./? = ν:? -表示第j個(gè)張量單項(xiàng)式，并滿足p + + + 1，c是由m個(gè)張 量系數(shù)組成的系數(shù)向量，表達(dá)式為〇=|>1，\2，.人]1^(>) = [；^(>)彳2(>)，...；^(>)]由1]1個(gè) 張量單項(xiàng)式構(gòu)成；
[0013] 2.2)數(shù)學(xué)模型：
[0014] 擴(kuò)散加權(quán)磁共振信號(hào)有η個(gè)擴(kuò)散梯度方向gi，i = l，...，n，并且沿重建向量v進(jìn)行 重建，那么系數(shù)向量c通過最小化下面的代價(jià)函數(shù)J(c)求得：
[0015]
[0016] 其中，Ei = S(gi)/So是第i個(gè)擴(kuò)散梯度方向gi上的衰減信號(hào)； 是一個(gè)mXm維的矩陣，其值只與擴(kuò)散梯度方向gl、重建向量v以及參數(shù)μ有關(guān)，對于每一個(gè)擴(kuò) 散梯度方向gi，都有一個(gè)Qi矩陣與之對應(yīng)；
[0017] (3)計(jì)算張量系數(shù)向量c，得到擴(kuò)散函數(shù)D(v)，再計(jì)算每個(gè)采樣點(diǎn)處的擴(kuò)散函數(shù)值， 最后將擴(kuò)散函數(shù)值擬合成擴(kuò)散模型，搜索極值并計(jì)算纖維方向。
[0018] 進(jìn)一步，所述步驟(3)中，所述張量系數(shù)向量c的計(jì)算包括以下步驟：
[0019] 3.1)在單位半球面上均勻采樣321個(gè)離散的點(diǎn)，以球心為原點(diǎn)獲取這321個(gè)重建向 量v，計(jì)算單條纖維響應(yīng)函數(shù)R(v，g)的值，設(shè)定高階張量模型的階數(shù)1，計(jì)算單項(xiàng)式矩陣F (v)，進(jìn)而計(jì)算出步驟2.2)中的矩陣&;
[0020] 3.2)使用BFGS擬牛頓搜索算法迭代求解2.2)中的最小化問題，步驟如下：
[0021]步驟3.2.1已知代價(jià)函數(shù)J(c)，選取一個(gè)初始點(diǎn)(^作為第一次迭代的搜索起始點(diǎn)， 迭代次數(shù)計(jì)為k=l，并設(shè)置最大迭代次數(shù)k_max，計(jì)算J(c)的梯度向量表達(dá)式：
[0022]
[0023] 步驟3.2.2對于第k次迭代，令表示第k次迭代的梯度向量，計(jì)算擬牛頓方
向dk = _Hktk作為第k次迭代的搜索方向，以ck為起點(diǎn)，沿方向dk進(jìn)行一維搜索，求得本次搜索 的可接受步長ak，矩陣Hk表示的是第k次迭代時(shí)⑶點(diǎn)處Hesse矩陣的逆的近似矩陣，其更新方 法如下-
[0024]
[0025] 其中，I是單位陣，Sk = Ck+1-Ck是相鄰兩次迭代的解的差，0k = tk+1_tk是相鄰兩次迭 代的梯度向量之差；
[0026] 步驟3.2.3更新張量系數(shù)向量Ck+1 = Ck+akdk，計(jì)算對應(yīng)的代價(jià)函數(shù)值J(ck+1)，當(dāng)達(dá) 到最大迭代次數(shù)或相鄰兩次迭代的代價(jià)函數(shù)值滿足時(shí)，則終止迭代，否則， 返回步驟3.2.2，進(jìn)入下一次迭代;參數(shù)σ為一個(gè)較小的常數(shù)，用于判斷算法是否收斂到一個(gè) 局部極小值。
[0027]再進(jìn)一步，所述步驟(3)還包括3.3)將得到的張量系數(shù)用于擬合擴(kuò)散函數(shù)，獲取纖 維方向分布函數(shù)模型，搜索極值并計(jì)算纖維方向，步驟如下：
[0028]步驟3.3.1對十二面體進(jìn)行5次細(xì)分，得到10242個(gè)球面上的相鄰等距點(diǎn)，以球心為 原點(diǎn)得到相應(yīng)個(gè)數(shù)的重建向量V，通過3.2)得到的張量系數(shù)向量可以求得擴(kuò)散函數(shù)
[0029]
[0030] 步驟3.3.2由前面得到的擴(kuò)散函數(shù)求得在10242個(gè)重建向量上的擴(kuò)散函數(shù)值，即纖 維方向分布函數(shù)的值，通過搜索纖維方向分布函數(shù)值中的極值點(diǎn)來獲取纖維的主方向，極 值點(diǎn)的搜索方法如下：
[0031 ]對每一個(gè)重建向量vq，q=l，…，10242,在10242個(gè)重建向量中搜索出與vq的夾角小 于Θ的所有向量，比較Vq與這些向量所對應(yīng)的纖維方向分布函數(shù)值的大小，若Vq所對應(yīng)的值 最大，則判斷v q為該體素的一個(gè)極值方向；依次遍歷所有重建向量，最后得到N個(gè)極值方向， 這些方向即為纖維的主方向。
[0032] 本發(fā)明的有益效果體現(xiàn)在:本發(fā)明采用了用高階笛卡爾張量擬合纖維方向分布函 數(shù)平方根的方式，保證了所得纖維方向分布函數(shù)的非負(fù)性，且角度分辨率較高，實(shí)驗(yàn)效果 好。
【具體實(shí)施方式】
[0033] 下面對本發(fā)明做進(jìn)一步說明。
[0034] -種非負(fù)高階張量擬牛頓搜索的纖維方向分布估計(jì)方法，所述的纖維方向分布估 計(jì)方法包括以下步驟：
[0035] (1)數(shù)據(jù)預(yù)處理:讀取腦部擴(kuò)散加權(quán)磁共振數(shù)據(jù)，獲取施加梯度方向g時(shí)的磁共振 信號(hào)S(g)和未施加梯度方向時(shí)的磁共振信號(hào)So,以及相應(yīng)的梯度方向數(shù)據(jù)，選取感興趣區(qū) 域，并計(jì)算該區(qū)域的擴(kuò)散衰減信號(hào)S(g)/So;
[0036] (2)將感興趣區(qū)域內(nèi)每個(gè)體素中的擴(kuò)散衰減信號(hào)逐個(gè)建模為具有擴(kuò)散形態(tài)的橢球 分布模型，建模過程如下：
[0037] 2.1)體素微結(jié)構(gòu)建模:將擴(kuò)散衰減信號(hào)假設(shè)為沿重建向量v的單條纖維信號(hào)響應(yīng) 函數(shù)R (v，g)與擴(kuò)散函數(shù)D (v)在球面上的卷積：
[0038]
[0039] 其中，j?:卜，@ '「近似看作一個(gè)高斯分布函數(shù)，g= {gi￡Rlx31 i = l，. . .，n}為梯 度方向，v={vPeRlx3|p = l，. . .，K}為在單位球面上采樣的重建向量，R1X3表示維度為1X3 的實(shí)數(shù)域矩陣空間，η和K分別表示梯度方向和重建向量的個(gè)數(shù)，μ = ε b是表征擴(kuò)散效率ε與 擴(kuò)散敏感系數(shù)b共同影響的一個(gè)參數(shù)，擴(kuò)散函數(shù)D(v)的表達(dá)式如下：
[0040]
[0041 ] drs表示單項(xiàng)式v; V, vf Μ的系數(shù)，1為高階張量的階數(shù)，r，s分別表示重建向量v = (vx，Vy，vz)的基方向vx和Vy的指數(shù)，λ」表示第j個(gè)張量的張量系數(shù)，j = l，. . .，m，m表示張量的 個(gè)數(shù)^(/ +彳+ 2)- ./>) = v:v:v廣表示第j個(gè)張量單項(xiàng)式，并滿足/ = f-I^ + s + 1，: c是由m 個(gè)張量系數(shù)組成的系數(shù)向量，表達(dá)式為= 由m個(gè)張量單項(xiàng)式構(gòu)成；
[0042] 2.2)數(shù)學(xué)模型：
[0043]擴(kuò)散加權(quán)磁共振信號(hào)有η個(gè)擴(kuò)散梯度方向gi，i = l，...，n，并且沿重建向量v進(jìn)行 重建，那么系數(shù)向量c通過最小化下面的代價(jià)函數(shù)J(c)求得：
[0044]
[0045] 其中，Ei = S(gi)/S〇是第i個(gè)擴(kuò)散梯度方向gi上的衰減信號(hào)；約>)> 是一個(gè)mXm維的矩陣，其值只與擴(kuò)散梯度方向gl、重建向量v以及參數(shù)μ有關(guān)，對于每一個(gè)擴(kuò) 散梯度方向gi，都有一個(gè)Qi矩陣與之對應(yīng)；
[0046] (3)計(jì)算張量系數(shù)向量c，得到擴(kuò)散函數(shù)D(v)，再計(jì)算每個(gè)采樣點(diǎn)處的擴(kuò)散函數(shù)值， 最后將擴(kuò)散函數(shù)值擬合成擴(kuò)散模型，搜索極值并計(jì)算纖維方向。<