基于變慣性系數(shù)active demons算法的DTI多通道配準(zhǔn)研究

鄭偉，周楊，李文華，劉帥奇，張曉丹，馬澤鵬

(1.河北大學(xué) 電子信息工程學(xué)院，河北 保定 071002；2.河北省數(shù)字醫(yī)療工程重點(diǎn)實(shí)驗(yàn)室，河北 保定 071002；3.河北省機(jī)器視覺(jué)中心，河北 保定 071002；4.河北大學(xué)附屬醫(yī)院 CT-MRI診斷科，河北 保定 071000)

彌散張量成像(diffision tensor imaging,DTI)是20世紀(jì)90年代新興的一種大腦結(jié)構(gòu)成像方法，它利用分子在大腦中的運(yùn)動(dòng)情況來(lái)繪成圖像,可以有效揭示阿爾茨海默病、精神分裂[1-2]、閱讀障礙等引起的大腦組織結(jié)構(gòu)反常變化，引導(dǎo)醫(yī)療人員進(jìn)行大腦手術(shù).在臨床上，對(duì)同一個(gè)病人不同時(shí)期DTI圖像研究分析和診斷具有十分重要的臨床價(jià)值.在做這些類型的比較之前，必須要進(jìn)行DTI圖像配準(zhǔn).DTI圖像配準(zhǔn)是指將2組或者多組DTI圖像進(jìn)行空間幾何變換，使各個(gè)像素或體素在相同解剖結(jié)構(gòu)上能夠相對(duì)應(yīng)，而且DTI圖像配準(zhǔn)的結(jié)果將直接影響纖維束追蹤的效果.基于光流場(chǎng)理論[3]的demons算法由于其算法的高效性和理論基礎(chǔ)的完整性而被廣泛應(yīng)用于非剛性圖像配準(zhǔn)中[4]，原始的demons算法由Thifion等在1998年提出[4]，其原理是利用參考圖像的梯度信息以及浮動(dòng)圖像與參考圖像的差作為形變向量驅(qū)動(dòng)浮動(dòng)圖像朝著參考圖像的方向發(fā)生形變.原始的demons算法收斂速度慢，對(duì)于大形變的配準(zhǔn)取得的效果也不理想.針對(duì)以上問(wèn)題，一些學(xué)者往驅(qū)動(dòng)力公式中加入新的信息項(xiàng)或參數(shù)，例如幾何形狀約束項(xiàng)[5]，李果霖等[6]將系數(shù)f(θ)引入驅(qū)動(dòng)力公式解決了參考圖像和浮動(dòng)圖像梯度差異過(guò)大引起的過(guò)配準(zhǔn)問(wèn)題.蔣晨嬌等[7]為了解決傳統(tǒng)DTI圖像配準(zhǔn)方法容易忽略空間信息及高維特性的問(wèn)題，提出基于改進(jìn)雞群優(yōu)化算法的AD患者DTI圖像配準(zhǔn)，提高了配準(zhǔn)精度，取得較好的配準(zhǔn)結(jié)果.為了保持配準(zhǔn)前后解剖結(jié)構(gòu)的一致性，Guo等[8]提出DTI配準(zhǔn)的高維空間標(biāo)準(zhǔn)化算法,先對(duì)待配準(zhǔn)圖像進(jìn)行空間標(biāo)準(zhǔn)化，然后對(duì)DTI圖像進(jìn)行高維標(biāo)準(zhǔn)化.通過(guò)對(duì)比實(shí)驗(yàn)驗(yàn)證了該算法的有效性.張桂梅等[9]提出自適應(yīng)分?jǐn)?shù)階demons算法解決了圖像容易陷入局部最優(yōu)的問(wèn)題.許磊等[10]針對(duì)原始demons算法不適合配準(zhǔn)大形變圖像的缺點(diǎn)對(duì)demons驅(qū)動(dòng)力公式進(jìn)行了重新定義.針對(duì)多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)效果差精度低的問(wèn)題，Simona等[11]提出結(jié)合稀疏和稠密特征改進(jìn)腦DTI圖像多模態(tài)配準(zhǔn)方法，使用直方圖梯度的互信息作為配準(zhǔn)標(biāo)準(zhǔn)，提高了多模態(tài)圖像配準(zhǔn)的精度.

Wang等[12]提出active demons算法，該算法對(duì)于形變較大的圖像有良好的配準(zhǔn)效果，收斂速度和精度相對(duì)demons算法有很大提高，但圖像的拓?fù)浣Y(jié)構(gòu)容易發(fā)生改變.對(duì)于均化系數(shù)α的取值問(wèn)題，α不能兼顧圖像的大小形變，而在配準(zhǔn)過(guò)程中大小形變是同時(shí)存在的，針對(duì)這個(gè)問(wèn)題，薛鵬等[13]提出基于平衡系數(shù)的active demons配準(zhǔn)算法，引入新的參數(shù)平衡系數(shù)k，固定均化系數(shù)，手動(dòng)選取最優(yōu)值，得到最優(yōu)平衡系數(shù)k.王昌等[14]提出基于局部聯(lián)合熵梯度的active demons算法，在驅(qū)動(dòng)力公式中加入新的信息項(xiàng)局部聯(lián)合熵梯度，并且引入了多分辨率的策略.薛文靜等[15]也將局部熵引入active demons算法，將改進(jìn)了active demons算法的驅(qū)動(dòng)力公式應(yīng)用于多模態(tài)的醫(yī)學(xué)圖像配準(zhǔn).徐曉瑩等[16]將平衡系數(shù)k當(dāng)做一個(gè)變量引入多通道DTI圖像配準(zhǔn)，提出變參數(shù)active demons算法，令k取不同值，可以動(dòng)態(tài)改變形變向量u的大小，其中手動(dòng)調(diào)整均化系數(shù)α，系數(shù)的最優(yōu)值通過(guò)實(shí)驗(yàn)測(cè)得，但是收斂速度依然不是很高.

本文在此基礎(chǔ)上提出了基于變慣性系數(shù)active demons算法的DTI多通道配準(zhǔn)，將慣性系數(shù)b引入變參數(shù)active demons算法中，并且改變b的取值結(jié)合平衡系數(shù)一起循環(huán)迭代，實(shí)驗(yàn)結(jié)果顯示與變參數(shù)active demons算法相比進(jìn)一步提高了算法的收斂速度并且得到了最好的配準(zhǔn)精度，同時(shí)圖像的拓?fù)浣Y(jié)構(gòu)得到了很好的保持.

1 變慣性系數(shù)active demons算法及在DTI多通道配準(zhǔn)中的應(yīng)用

水分子在純凈液體中向各個(gè)方向彌散的程度相同,把這種性質(zhì)稱作各向同性.但是大腦中有各種各樣的組織會(huì)阻礙水分子的運(yùn)動(dòng)，水分子向各方向的彌散存在差異,這種差異稱作各向異性.由于各向異性的存在,水分子在大腦中的移動(dòng)可以看作橢球形，這種性質(zhì)可以用數(shù)學(xué)上的對(duì)稱二階張量來(lái)表示，如式(1)所示:

(1)

對(duì)稱矩陣有9個(gè)元素，獨(dú)立的元素只有6個(gè)，因此至少需要在6個(gè)不同的方向上測(cè)量彌散張量.

連續(xù)運(yùn)動(dòng)的2幀圖像存在亮度的變化即光流現(xiàn)象，Thirion等[17]將光流的概念引入到圖像配準(zhǔn)中，把參考圖像和浮動(dòng)圖像看作連續(xù)運(yùn)動(dòng)的圖像，由于光流現(xiàn)象的存在,這2幅圖像存在灰度梯度的差異，利用這種差異構(gòu)造一個(gè)驅(qū)動(dòng)力來(lái)驅(qū)動(dòng)浮動(dòng)圖像向參考圖像發(fā)生形變來(lái)進(jìn)行圖像的配準(zhǔn).假定2幅圖像的亮度不變，則二維圖像的光流場(chǎng)方程如式(2)所示:

(2)

式中，u表示點(diǎn)(x,y)處的形變向量；M(x,y)和S(x,y)分別表示浮動(dòng)圖像和參考圖像的灰度值；S(x,y)表示參考圖像在(x,y)處的梯度.為防止S(x,y)趨于0或不存在時(shí)方程失去意義，將式(2)改為式(3):

(3)

式中，α為均化系數(shù).

由demons算法的原理可知連續(xù)的2幀圖像形變向量較小，所以該算法不適用于大形變圖像配準(zhǔn)，且收斂速度慢.

針對(duì)以上demons算法出現(xiàn)的問(wèn)題，Rogelj等[18]提出active demons算法把浮動(dòng)圖像的梯度信息也引入形變方程見(jiàn)式4，參考圖像的梯度信息和浮動(dòng)圖像梯度信息的合力作為驅(qū)動(dòng)浮動(dòng)圖像發(fā)生形變的demons力，同時(shí)改變均化系數(shù)α的取值也可以改變驅(qū)動(dòng)力的大小，該算法有效解決了單一依靠參考圖像梯度信息的缺陷,收斂速度和精度得到提高.

(4)

均化系數(shù)α越大，形變向量u越小，收斂速度變慢，配準(zhǔn)時(shí)間較長(zhǎng)，但精度更高，適合小形變區(qū)域的圖像配準(zhǔn)；α越小，形變向量u越大，收斂速度更快，配準(zhǔn)時(shí)間較短，但精度低，適合大形變區(qū)域的圖像配準(zhǔn).

對(duì)均化系數(shù)α分析可知，α不能同時(shí)兼顧大形變和小形變的圖像配準(zhǔn).薛鵬等[13]受參數(shù)α的啟發(fā)，在分母的梯度模值平方項(xiàng)上加入一個(gè)新的參數(shù)平衡系數(shù)k，平衡系數(shù)k和均化系數(shù)α聯(lián)合調(diào)整驅(qū)動(dòng)力強(qiáng)度，如式(5)所示:

(5)

如果固定α，k值減小，u增大，形變向量增大，收斂速度加快；k值增大，u減小，形變向量減小，收斂速度變慢.

將6個(gè)獨(dú)立的分量Dxx、Dyy、Dzz、Dxy、Dxz、Dyz分別作為算法的輸入疊加求和再取平均值作為驅(qū)動(dòng)參考圖像發(fā)生形變的demons力.對(duì)于DTI圖像中的任意體素其位移的迭代如式(6)所示:

(6)

針對(duì)active demons算法收斂速度慢的問(wèn)題，在un前加入慣性系數(shù)b，根據(jù)物理學(xué)中的慣性概念[19]，慣性系數(shù)b為上次迭代對(duì)本次迭代的影響.當(dāng)上次結(jié)果對(duì)本次迭代影響較小時(shí)，令b取較小的值，反之則增大b的取值.通過(guò)動(dòng)態(tài)改變慣性系數(shù)的取值來(lái)影響驅(qū)動(dòng)力強(qiáng)度是本文研究的核心.變慣性系數(shù)active demons算法的驅(qū)動(dòng)力迭代公式如式(7)所示.

(7)

利用式(7)計(jì)算浮動(dòng)圖像的形變向量，位移場(chǎng)的平滑使用高斯濾波器，插值采用三線性插值，利用主方向保留法(presreving principle diffusion，PPD)對(duì)張量的方向進(jìn)行重定向,具體流程如圖1所示:

圖1 變慣性系數(shù)active demons算法配準(zhǔn)流程Fig.1 Active demons algorithm with variable inertia coefficient registration process

2 實(shí)驗(yàn)及結(jié)果分析

2.1 實(shí)驗(yàn)數(shù)據(jù)及預(yù)處理

本文使用的DTI數(shù)據(jù)來(lái)自標(biāo)準(zhǔn)數(shù)據(jù)庫(kù)ADNI(https://ida.loni.usc.edu/login.jsp?Project=AD-NI)，數(shù)據(jù)來(lái)自同一AD患者在2個(gè)不同時(shí)期采用相同設(shè)備和參數(shù)采集的DTI數(shù)據(jù)，設(shè)備為3T的GE醫(yī)療系統(tǒng)，重復(fù)時(shí)間(TR)為9 050 ms,回波時(shí)間(TE)為61.8 ms,翻轉(zhuǎn)角(FA)為90°,圖像矩陣大小為256×256×60，厚度為2.7 mm，反映附加梯度場(chǎng)快慢的參數(shù)b為1 000 s/mm2.

采用FSL(http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Fsllnstallation)對(duì)DTI原始實(shí)驗(yàn)數(shù)據(jù)做預(yù)處理，原始數(shù)據(jù)格式轉(zhuǎn)換后先對(duì)DTI數(shù)據(jù)做渦流矯正輸出.nii.gz格式的DTI數(shù)據(jù)，再去除非腦組織得到全腦mask文件,最后計(jì)算FA、MD等全張量數(shù)據(jù).原始的軸向圖見(jiàn)圖2.

2.2 變慣性系數(shù) active demons算法下的多通道DTI配準(zhǔn)結(jié)果與分析

引入慣性系數(shù)b，單一的b值在取值大于1的情況下可以暫時(shí)增加算法的收斂速度，但是隨著迭代次數(shù)的增加，配準(zhǔn)精度降低，不能得到好的實(shí)驗(yàn)結(jié)果，本文的改進(jìn)是動(dòng)態(tài)改變慣性系數(shù)的取值，使配準(zhǔn)前期的b取一個(gè)較大的值隨著迭代的增加b逐漸減小，配合平衡系數(shù)k一起調(diào)節(jié)驅(qū)動(dòng)力的大小.下面先討論單一的b值對(duì)實(shí)驗(yàn)結(jié)果影響.令均化系數(shù)α=0.007，k[16]的取值范圍為0.5

a．軸向參考FA圖；b.軸向浮動(dòng)FA圖.圖2 原始軸向圖像Fig.2 Original axial image

圖3 單一b值的MSE圖Fig.3 MSE diagram for a single b value

縱坐標(biāo)MSE(mean square error,均方誤差)為均方誤差，橫坐標(biāo)為迭代次數(shù).當(dāng)b的取值小于1時(shí)，圖像的收斂速度和精度明顯降低，這是因?yàn)榕錅?zhǔn)前期k的取值較小，這樣可以配準(zhǔn)大形變量，b的取值如果也較小會(huì)減小前期的形變量.因此b的取值小于1會(huì)降低收斂速度和精度.當(dāng)b等于1時(shí)就等價(jià)于變參數(shù)active demons 的迭代方程.當(dāng)b大于1時(shí)，慣性系數(shù)b與平衡系數(shù)k一同加快收斂速度，隨著迭代的進(jìn)行，大形變得到配準(zhǔn)，可是對(duì)于小形變出現(xiàn)過(guò)配準(zhǔn)，圖像拓?fù)浣Y(jié)構(gòu)發(fā)生改變，主要原因是b值偏大，不能兼顧小形變.從圖3可以看出無(wú)論b值大于1還是小于1都不能達(dá)到理想的效果甚至出現(xiàn)過(guò)配準(zhǔn)現(xiàn)象.

本文提出的基于變慣性系數(shù)的active demons算法可以解決以上出現(xiàn)的問(wèn)題，前期通過(guò)讓b取較大的值加快收斂速度，隨著迭代的進(jìn)行，大形變區(qū)域得到配準(zhǔn)，再降低b值，配準(zhǔn)小形變區(qū)域.本實(shí)驗(yàn)中根據(jù)文獻(xiàn)[16]，α=0.007，高斯平滑的標(biāo)準(zhǔn)差σ采用Park等[20]在demons算法下設(shè)定的11、9、5、3.標(biāo)準(zhǔn)差σ隨迭代次數(shù)的增加在迭代次數(shù)每1/4處減小，平衡系數(shù)k隨著σ按0.5、0.7、0.9、1的順序依次增大，b隨著標(biāo)準(zhǔn)差按1.15、1.05、1.01、1的順序依次減小.實(shí)驗(yàn)結(jié)果如圖4所示.

圖4 3種算法的對(duì)比MSE曲線Fig.4 Comparison of MSE curves of three algorithms

AD(active demons),IAD(improve active demons),IIAD(improve inertia coefficient active demons)分別為active demons算法、變參數(shù)active demons算法、變慣性系數(shù)active demons算法的MSE曲線，由圖4可知，本文算法迭代27次的結(jié)果與變參數(shù)active demons算法迭代50次的結(jié)果相同，本文算法的MSE值比變參數(shù)active demons算法提高13%，比active demons算法的MSE值提高19%,變慣性系數(shù)active demons算法相較于active demons算法和變參數(shù)active demons算法取得了較好的實(shí)驗(yàn)結(jié)果，見(jiàn)表1.

表1 3種算法在相同迭代次數(shù)下的MSE值比較Tab.1 Comparison of MSE values of the three algorithms with the same number of iterations

配準(zhǔn)后的FA圖及FA差值圖如圖5所示.

a.原始demons算法；b.active demons算法；c.變參數(shù)active demons算法；d.變慣性系數(shù)active demons算法.圖5 4種算法的FA圖和FA差值圖Fig.5 FA images and FA difference graphs of the four algorithms

可以看出變慣性系數(shù)active demons算法的配準(zhǔn)效果要遠(yuǎn)遠(yuǎn)好于demons算法，并且相對(duì)于其他2種算法也取得了最小的一個(gè)MSE值.

3 總結(jié)

在變參數(shù)active demons算法基礎(chǔ)上做改進(jìn)，引入慣性系數(shù)b，討論b的取值對(duì)實(shí)驗(yàn)結(jié)果的影響，得出結(jié)論單一b值無(wú)法取得理想的結(jié)果，b在一定范圍內(nèi)取合適的變值可以提高算法的收斂速度和精度.比較配準(zhǔn)后圖像FA圖，發(fā)現(xiàn)配準(zhǔn)后圖像的拓?fù)浣Y(jié)構(gòu)得到很好的保持.本文算法也存在如下不足，由于改進(jìn)算法引入了新參數(shù),在進(jìn)行配準(zhǔn)實(shí)驗(yàn)時(shí),為了達(dá)到高配準(zhǔn)精度,針對(duì)不同的圖像類型,可能需要進(jìn)行一定的參數(shù)取值實(shí)驗(yàn)工作,以選擇最優(yōu)的配準(zhǔn)結(jié)果.