鄭偉,周楊,李文華,劉帥奇,張曉丹,馬澤鵬
(1.河北大學(xué) 電子信息工程學(xué)院,河北 保定 071002;2.河北省數(shù)字醫(yī)療工程重點實驗室,河北 保定 071002;3.河北省機器視覺中心,河北 保定 071002;4.河北大學(xué)附屬醫(yī)院 CT-MRI診斷科,河北 保定 071000)
彌散張量成像(diffision tensor imaging,DTI)是20世紀90年代新興的一種大腦結(jié)構(gòu)成像方法,它利用分子在大腦中的運動情況來繪成圖像,可以有效揭示阿爾茨海默病、精神分裂[1-2]、閱讀障礙等引起的大腦組織結(jié)構(gòu)反常變化,引導(dǎo)醫(yī)療人員進行大腦手術(shù).在臨床上,對同一個病人不同時期DTI圖像研究分析和診斷具有十分重要的臨床價值.在做這些類型的比較之前,必須要進行DTI圖像配準.DTI圖像配準是指將2組或者多組DTI圖像進行空間幾何變換,使各個像素或體素在相同解剖結(jié)構(gòu)上能夠相對應(yīng),而且DTI圖像配準的結(jié)果將直接影響纖維束追蹤的效果.基于光流場理論[3]的demons算法由于其算法的高效性和理論基礎(chǔ)的完整性而被廣泛應(yīng)用于非剛性圖像配準中[4],原始的demons算法由Thifion等在1998年提出[4],其原理是利用參考圖像的梯度信息以及浮動圖像與參考圖像的差作為形變向量驅(qū)動浮動圖像朝著參考圖像的方向發(fā)生形變.原始的demons算法收斂速度慢,對于大形變的配準取得的效果也不理想.針對以上問題,一些學(xué)者往驅(qū)動力公式中加入新的信息項或參數(shù),例如幾何形狀約束項[5],李果霖等[6]將系數(shù)f(θ)引入驅(qū)動力公式解決了參考圖像和浮動圖像梯度差異過大引起的過配準問題.蔣晨嬌等[7]為了解決傳統(tǒng)DTI圖像配準方法容易忽略空間信息及高維特性的問題,提出基于改進雞群優(yōu)化算法的AD患者DTI圖像配準,提高了配準精度,取得較好的配準結(jié)果.為了保持配準前后解剖結(jié)構(gòu)的一致性,Guo等[8]提出DTI配準的高維空間標準化算法,先對待配準圖像進行空間標準化,然后對DTI圖像進行高維標準化.通過對比實驗驗證了該算法的有效性.張桂梅等[9]提出自適應(yīng)分數(shù)階demons算法解決了圖像容易陷入局部最優(yōu)的問題.許磊等[10]針對原始demons算法不適合配準大形變圖像的缺點對demons驅(qū)動力公式進行了重新定義.針對多模態(tài)醫(yī)學(xué)圖像配準效果差精度低的問題,Simona等[11]提出結(jié)合稀疏和稠密特征改進腦DTI圖像多模態(tài)配準方法,使用直方圖梯度的互信息作為配準標準,提高了多模態(tài)圖像配準的精度.
Wang等[12]提出active demons算法,該算法對于形變較大的圖像有良好的配準效果,收斂速度和精度相對demons算法有很大提高,但圖像的拓撲結(jié)構(gòu)容易發(fā)生改變.對于均化系數(shù)α的取值問題,α不能兼顧圖像的大小形變,而在配準過程中大小形變是同時存在的,針對這個問題,薛鵬等[13]提出基于平衡系數(shù)的active demons配準算法,引入新的參數(shù)平衡系數(shù)k,固定均化系數(shù),手動選取最優(yōu)值,得到最優(yōu)平衡系數(shù)k.王昌等[14]提出基于局部聯(lián)合熵梯度的active demons算法,在驅(qū)動力公式中加入新的信息項局部聯(lián)合熵梯度,并且引入了多分辨率的策略.薛文靜等[15]也將局部熵引入active demons算法,將改進了active demons算法的驅(qū)動力公式應(yīng)用于多模態(tài)的醫(yī)學(xué)圖像配準.徐曉瑩等[16]將平衡系數(shù)k當做一個變量引入多通道DTI圖像配準,提出變參數(shù)active demons算法,令k取不同值,可以動態(tài)改變形變向量u的大小,其中手動調(diào)整均化系數(shù)α,系數(shù)的最優(yōu)值通過實驗測得,但是收斂速度依然不是很高.
本文在此基礎(chǔ)上提出了基于變慣性系數(shù)active demons算法的DTI多通道配準,將慣性系數(shù)b引入變參數(shù)active demons算法中,并且改變b的取值結(jié)合平衡系數(shù)一起循環(huán)迭代,實驗結(jié)果顯示與變參數(shù)active demons算法相比進一步提高了算法的收斂速度并且得到了最好的配準精度,同時圖像的拓撲結(jié)構(gòu)得到了很好的保持.
水分子在純凈液體中向各個方向彌散的程度相同,把這種性質(zhì)稱作各向同性.但是大腦中有各種各樣的組織會阻礙水分子的運動,水分子向各方向的彌散存在差異,這種差異稱作各向異性.由于各向異性的存在,水分子在大腦中的移動可以看作橢球形,這種性質(zhì)可以用數(shù)學(xué)上的對稱二階張量來表示,如式(1)所示:
(1)
對稱矩陣有9個元素,獨立的元素只有6個,因此至少需要在6個不同的方向上測量彌散張量.
連續(xù)運動的2幀圖像存在亮度的變化即光流現(xiàn)象,Thirion等[17]將光流的概念引入到圖像配準中,把參考圖像和浮動圖像看作連續(xù)運動的圖像,由于光流現(xiàn)象的存在,這2幅圖像存在灰度梯度的差異,利用這種差異構(gòu)造一個驅(qū)動力來驅(qū)動浮動圖像向參考圖像發(fā)生形變來進行圖像的配準.假定2幅圖像的亮度不變,則二維圖像的光流場方程如式(2)所示:
(2)
式中,u表示點(x,y)處的形變向量;M(x,y)和S(x,y)分別表示浮動圖像和參考圖像的灰度值;S(x,y)表示參考圖像在(x,y)處的梯度.為防止S(x,y)趨于0或不存在時方程失去意義,將式(2)改為式(3):
(3)
式中,α為均化系數(shù).
由demons算法的原理可知連續(xù)的2幀圖像形變向量較小,所以該算法不適用于大形變圖像配準,且收斂速度慢.
針對以上demons算法出現(xiàn)的問題,Rogelj等[18]提出active demons算法把浮動圖像的梯度信息也引入形變方程見式4,參考圖像的梯度信息和浮動圖像梯度信息的合力作為驅(qū)動浮動圖像發(fā)生形變的demons力,同時改變均化系數(shù)α的取值也可以改變驅(qū)動力的大小,該算法有效解決了單一依靠參考圖像梯度信息的缺陷,收斂速度和精度得到提高.
(4)
均化系數(shù)α越大,形變向量u越小,收斂速度變慢,配準時間較長,但精度更高,適合小形變區(qū)域的圖像配準;α越小,形變向量u越大,收斂速度更快,配準時間較短,但精度低,適合大形變區(qū)域的圖像配準.
對均化系數(shù)α分析可知,α不能同時兼顧大形變和小形變的圖像配準.薛鵬等[13]受參數(shù)α的啟發(fā),在分母的梯度模值平方項上加入一個新的參數(shù)平衡系數(shù)k,平衡系數(shù)k和均化系數(shù)α聯(lián)合調(diào)整驅(qū)動力強度,如式(5)所示:
(5)
如果固定α,k值減小,u增大,形變向量增大,收斂速度加快;k值增大,u減小,形變向量減小,收斂速度變慢.
將6個獨立的分量Dxx、Dyy、Dzz、Dxy、Dxz、Dyz分別作為算法的輸入疊加求和再取平均值作為驅(qū)動參考圖像發(fā)生形變的demons力.對于DTI圖像中的任意體素其位移的迭代如式(6)所示:
(6)
針對active demons算法收斂速度慢的問題,在un前加入慣性系數(shù)b,根據(jù)物理學(xué)中的慣性概念[19],慣性系數(shù)b為上次迭代對本次迭代的影響.當上次結(jié)果對本次迭代影響較小時,令b取較小的值,反之則增大b的取值.通過動態(tài)改變慣性系數(shù)的取值來影響驅(qū)動力強度是本文研究的核心.變慣性系數(shù)active demons算法的驅(qū)動力迭代公式如式(7)所示.
(7)
利用式(7)計算浮動圖像的形變向量,位移場的平滑使用高斯濾波器,插值采用三線性插值,利用主方向保留法(presreving principle diffusion,PPD)對張量的方向進行重定向,具體流程如圖1所示:
圖1 變慣性系數(shù)active demons算法配準流程Fig.1 Active demons algorithm with variable inertia coefficient registration process
本文使用的DTI數(shù)據(jù)來自標準數(shù)據(jù)庫ADNI(https://ida.loni.usc.edu/login.jsp?Project=AD-NI),數(shù)據(jù)來自同一AD患者在2個不同時期采用相同設(shè)備和參數(shù)采集的DTI數(shù)據(jù),設(shè)備為3T的GE醫(yī)療系統(tǒng),重復(fù)時間(TR)為9 050 ms,回波時間(TE)為61.8 ms,翻轉(zhuǎn)角(FA)為90°,圖像矩陣大小為256×256×60,厚度為2.7 mm,反映附加梯度場快慢的參數(shù)b為1 000 s/mm2.
采用FSL(http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Fsllnstallation)對DTI原始實驗數(shù)據(jù)做預(yù)處理,原始數(shù)據(jù)格式轉(zhuǎn)換后先對DTI數(shù)據(jù)做渦流矯正輸出.nii.gz格式的DTI數(shù)據(jù),再去除非腦組織得到全腦mask文件,最后計算FA、MD等全張量數(shù)據(jù).原始的軸向圖見圖2.
引入慣性系數(shù)b,單一的b值在取值大于1的情況下可以暫時增加算法的收斂速度,但是隨著迭代次數(shù)的增加,配準精度降低,不能得到好的實驗結(jié)果,本文的改進是動態(tài)改變慣性系數(shù)的取值,使配準前期的b取一個較大的值隨著迭代的增加b逐漸減小,配合平衡系數(shù)k一起調(diào)節(jié)驅(qū)動力的大小.下面先討論單一的b值對實驗結(jié)果影響.令均化系數(shù)α=0.007,k[16]的取值范圍為0.5 a.軸向參考FA圖;b.軸向浮動FA圖.圖2 原始軸向圖像Fig.2 Original axial image 圖3 單一b值的MSE圖Fig.3 MSE diagram for a single b value 縱坐標MSE(mean square error,均方誤差)為均方誤差,橫坐標為迭代次數(shù).當b的取值小于1時,圖像的收斂速度和精度明顯降低,這是因為配準前期k的取值較小,這樣可以配準大形變量,b的取值如果也較小會減小前期的形變量.因此b的取值小于1會降低收斂速度和精度.當b等于1時就等價于變參數(shù)active demons 的迭代方程.當b大于1時,慣性系數(shù)b與平衡系數(shù)k一同加快收斂速度,隨著迭代的進行,大形變得到配準,可是對于小形變出現(xiàn)過配準,圖像拓撲結(jié)構(gòu)發(fā)生改變,主要原因是b值偏大,不能兼顧小形變.從圖3可以看出無論b值大于1還是小于1都不能達到理想的效果甚至出現(xiàn)過配準現(xiàn)象. 本文提出的基于變慣性系數(shù)的active demons算法可以解決以上出現(xiàn)的問題,前期通過讓b取較大的值加快收斂速度,隨著迭代的進行,大形變區(qū)域得到配準,再降低b值,配準小形變區(qū)域.本實驗中根據(jù)文獻[16],α=0.007,高斯平滑的標準差σ采用Park等[20]在demons算法下設(shè)定的11、9、5、3.標準差σ隨迭代次數(shù)的增加在迭代次數(shù)每1/4處減小,平衡系數(shù)k隨著σ按0.5、0.7、0.9、1的順序依次增大,b隨著標準差按1.15、1.05、1.01、1的順序依次減小.實驗結(jié)果如圖4所示. 圖4 3種算法的對比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算法取得了較好的實驗結(jié)果,見表1. 表1 3種算法在相同迭代次數(shù)下的MSE值比較Tab.1 Comparison of MSE values of the three algorithms with the same number of iterations 配準后的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算法的配準效果要遠遠好于demons算法,并且相對于其他2種算法也取得了最小的一個MSE值. 在變參數(shù)active demons算法基礎(chǔ)上做改進,引入慣性系數(shù)b,討論b的取值對實驗結(jié)果的影響,得出結(jié)論單一b值無法取得理想的結(jié)果,b在一定范圍內(nèi)取合適的變值可以提高算法的收斂速度和精度.比較配準后圖像FA圖,發(fā)現(xiàn)配準后圖像的拓撲結(jié)構(gòu)得到很好的保持.本文算法也存在如下不足,由于改進算法引入了新參數(shù),在進行配準實驗時,為了達到高配準精度,針對不同的圖像類型,可能需要進行一定的參數(shù)取值實驗工作,以選擇最優(yōu)的配準結(jié)果.3 總結(jié)