高強(qiáng) ,鮑園,趙俊
1 上海交通大學(xué) 生物醫(yī)學(xué)工程學(xué)院,上海市,200240
2 上海聯(lián)影醫(yī)療科技有限公司,上海市,201807
冠心病是心血管疾病中發(fā)病率、死亡率最高的,據(jù)統(tǒng)計(jì)在西方國(guó)家每年約有30%的人死于該疾病,并且我國(guó)也處于不斷上升的趨勢(shì)[1]。冠心病的早期診斷具有重要價(jià)值[1-2]。冠狀動(dòng)脈CT血管造影(computed tomography angiography,CTA)具有成本低、無(wú)創(chuàng)、準(zhǔn)確性和敏感性高等特點(diǎn),是繼數(shù)字減影血管造影(digital subtraction angiography,DSA)之后監(jiān)測(cè)和診斷心臟疾病的一種安全、可靠的方法[1]。CTA系統(tǒng)能成為檢測(cè)心臟疾病的有效方法,得益于高速的球管旋轉(zhuǎn)、寬大的探測(cè)器以及雙源掃描等方式,減少了心臟的掃描時(shí)間,例如當(dāng)今超高端的CT,可在0.28 s內(nèi)完成360o的旋轉(zhuǎn)[3],從而提高了系統(tǒng)時(shí)間分辨率,提升了圖像質(zhì)量。但是每次掃描重建的圖像,還可能會(huì)因心臟跳動(dòng)頻率比CT探測(cè)器采集頻率高而產(chǎn)生運(yùn)動(dòng)偽影。
非剛性圖像配準(zhǔn)方法是解決上述問(wèn)題的方法之一[4],它主要包括心臟冠狀動(dòng)脈的特征點(diǎn)提取,生成運(yùn)動(dòng)向量場(chǎng)以及圖像的灰度映射三部分[5]。其中,在生成運(yùn)動(dòng)向量場(chǎng)方面,樣條插值方法可解決心臟冠狀動(dòng)脈運(yùn)動(dòng)自由度高、不確定性高等問(wèn)題,是消除心臟冠狀動(dòng)脈偽影的主流方法[6~8]。薄板樣條、(thin plate spline,TPS)插值方法是運(yùn)動(dòng)向量場(chǎng)插值方法的一種選擇[8],通過(guò)此方法可以找到通過(guò)所有控制點(diǎn)的最短路徑,從而模擬心臟冠脈的運(yùn)動(dòng)軌跡,但是TPS屬于全局性的配準(zhǔn)方法,任意一個(gè)控制點(diǎn)發(fā)生變化都會(huì)對(duì)全局產(chǎn)生影響,配準(zhǔn)精度并不是太高[6]。在執(zhí)行效率方面,該方法運(yùn)行時(shí)間過(guò)長(zhǎng),無(wú)法達(dá)到實(shí)時(shí)處理的要求。心臟冠狀動(dòng)脈運(yùn)動(dòng)是由心臟跳動(dòng)與呼吸共同作用的結(jié)果,因此運(yùn)動(dòng)自由度比較高,需要對(duì)運(yùn)動(dòng)的局部性加以考慮。傳統(tǒng)的B樣條插值方法雖然在改變控制點(diǎn)時(shí)只會(huì)影響局部鄰域,能夠解決局部性變形的問(wèn)題,但是隨著物體模型變形程度的增加,其運(yùn)動(dòng)向量場(chǎng)配準(zhǔn)的精度和執(zhí)行速度也會(huì)下降。
我們提出了一種基于均勻網(wǎng)格多層次B樣條的CTA圖像的運(yùn)動(dòng)向量場(chǎng)插值方法,此方法可以從多尺度角度來(lái)控制網(wǎng)格的變化,一方面它從較稀疏的網(wǎng)格著手,能夠模擬全局的彈性形變;另一方面它能夠在上一層次的基礎(chǔ)之上,逐漸步入較細(xì)的網(wǎng)格層次,由于控制點(diǎn)之間的間距變小,可以模擬局部的變形,從而提高配準(zhǔn)精度。另外,利用分割矩陣的數(shù)學(xué)方法,可將插值時(shí)大量造成冗余計(jì)算的共享節(jié)點(diǎn)進(jìn)行過(guò)濾,不僅兼顧了原有算法的配準(zhǔn)精度,同時(shí)也提高了算法的執(zhí)行效率。
本研究提出的多層次B樣條算法[9]是一種多分辨率多尺度的插值方法,它允許從離散的特征點(diǎn)集中擬合出三維自由形變的曲面。令三自變量目標(biāo)控制點(diǎn)集:
P屬于三維運(yùn)動(dòng)向量場(chǎng)V={(x,y,z,d)|0≤x<m,0≤y<n,0≤z<l}內(nèi)部。為了逼近P,構(gòu)建一個(gè)均勻的三-三次B樣條函數(shù)作為P的插值函數(shù)。令A(yù)是覆蓋V區(qū)域大小為(m+3)×(n+3)×(l+3)的三維空間,Aijk網(wǎng)格是A空間網(wǎng)格中的某一控制點(diǎn),那么P的插值函數(shù)為:
其中0≤t<1。
曲面的逼近問(wèn)題是確定運(yùn)動(dòng)向量場(chǎng)中的所有控制點(diǎn)的值A(chǔ),根據(jù)式(2)可知函數(shù)值f(x,y,z)與數(shù)據(jù)點(diǎn)(x,y,z)鄰域內(nèi)的64個(gè)控制點(diǎn)相關(guān),并且按照LEE等[9]提出的最小二乘法的思想來(lái)保證待求的V中每個(gè)控制點(diǎn)Aabc都在A范圍內(nèi),最終通過(guò)以下方式計(jì)算出所有的控制點(diǎn):
其中wabc=Ba(s)Bb(t)Bc(u),wefg=Be(s)Bf(t)Bg(u)
包含于空間網(wǎng)格A中的V可從多個(gè)尺度逐層插值,并可產(chǎn)生一系列插值函數(shù)fk,而它們的和就是逼近目標(biāo)的插值函數(shù),計(jì)算式為:
最稀疏的網(wǎng)格是由參考圖像特征點(diǎn)與目標(biāo)圖像特征點(diǎn)的差計(jì)算得到的,并且可利用此網(wǎng)格的格點(diǎn)得到比較粗略的近似插值函數(shù)。后續(xù)細(xì)化后的控制格點(diǎn)可以進(jìn)一步逼近和消除殘差。當(dāng)經(jīng)過(guò)數(shù)次迭代得到目標(biāo)尺寸的網(wǎng)格時(shí),插值即可完成,參見(jiàn)圖1的插值過(guò)程。
參考網(wǎng)格V代表運(yùn)動(dòng)向量場(chǎng)的初始位置,形變網(wǎng)格相當(dāng)于運(yùn)動(dòng)向量場(chǎng)的位移量。起初網(wǎng)格比較稀疏,格點(diǎn)間的距離較大。經(jīng)過(guò)每次迭代插值維度成倍增加,格點(diǎn)間距離成倍縮小,并且使用每次迭代生成的向量場(chǎng)(即參考網(wǎng)格與形變網(wǎng)格之和)通過(guò)式(2)得到的結(jié)果與目標(biāo)匹配點(diǎn)的差作為下一次迭代的輸入,這樣反復(fù)N次迭代,當(dāng)達(dá)到目標(biāo)大小時(shí)結(jié)束。
圖1 MBS插值流程圖Fig.1 Flow chart of MBS interpolation
本研究采用了一種分割矩陣的數(shù)學(xué)方法,它允許對(duì)基于均勻網(wǎng)格的B樣條生成冠狀動(dòng)脈的運(yùn)動(dòng)向量場(chǎng)進(jìn)行迭代插值。TPS的插值方法忽視了冠狀動(dòng)脈的有效取值范圍,對(duì)整個(gè)三維空間中所有坐標(biāo)點(diǎn)進(jìn)行全局插值,而本研究的方法可以只對(duì)提取的感興趣局部區(qū)域進(jìn)行插值。
其基本方法是將原有的B樣條曲線分成左右兩條曲線段,B樣條矩陣[10]如下:
其中k=0,1,…,n-3且0≤t<1,M為B樣條基函數(shù)系數(shù)。假設(shè)一個(gè)B樣條曲線P(t)是由控制點(diǎn){P0,P1,P2,P3}構(gòu)成的,我們將曲線二元細(xì)分,產(chǎn)生兩條曲線PL(t)和PR(t),并將原有t的[0,1)的取值范圍劃分為[0,0.5)和[0.5,1)。然后對(duì)這兩條曲線段進(jìn)行一系列的變換,最終得到基于控制點(diǎn){P0,P1,P2,P3}插值點(diǎn)的表達(dá)式:
通過(guò)(7),我們可以知道PL(t)和PR(t)中存在被重復(fù)插值的點(diǎn),即。這些重復(fù)點(diǎn)可以從任意一個(gè)曲線段中選取,最終可得到插值點(diǎn)公式:
過(guò)濾掉分割矩陣插值方法中產(chǎn)生的重復(fù)的插值點(diǎn)是提高算法執(zhí)行速度的關(guān)鍵。根據(jù)算法中的插值方法可以將待插值點(diǎn)的位置通過(guò)左右相鄰的兩個(gè)或者三個(gè)控制點(diǎn)通過(guò)式(8)計(jì)算得到,這樣可以不用考慮控制點(diǎn)的空間維度只需將控制點(diǎn)集P按照一維的排列方式表示并額外復(fù)制3組相同的點(diǎn)集,統(tǒng)一定義為P1、P2、P3以及P4,其中點(diǎn)集中的每個(gè)元素以D來(lái)表示,而、'以及'是待插值的點(diǎn)集(見(jiàn)圖2)。
圖2 待插值控制網(wǎng)格數(shù)據(jù)布局Fig.2 Layout of control grid data to be interpolated
每組點(diǎn)集按照同一點(diǎn)集中點(diǎn)與點(diǎn)之間的相鄰關(guān)系進(jìn)行縱向排列,這樣可以對(duì)數(shù)據(jù)進(jìn)行很好的并行處理,對(duì)于4組控制點(diǎn)集只需考慮彼此間對(duì)應(yīng)位置即可。但是存在一個(gè)問(wèn)題,那就是會(huì)出現(xiàn)大量插值點(diǎn)被重復(fù)運(yùn)算導(dǎo)致計(jì)算資源的浪費(fèi)(如圖2中的第一列P2和P3中的D2和D3以及第二列P1和P2中的D2和D3),這主要是和分割矩陣算法特點(diǎn)有關(guān)。為了解決此問(wèn)題,不妨先只對(duì)進(jìn)行插值,當(dāng)完成了當(dāng)前維度的運(yùn)算之后,再對(duì)以及做相同維度數(shù)據(jù)間隔長(zhǎng)度的插值。根據(jù)三維空間的幾何關(guān)系可知X方向插值數(shù)據(jù)間隔為1,Y方向插值數(shù)據(jù)間隔為X方向元素個(gè)數(shù),Z方向插值數(shù)據(jù)間隔為XY平面的元素個(gè)數(shù)。通過(guò)上述方法得到如下結(jié)果:
(1)優(yōu)化前總操作次數(shù)=31×插值數(shù)據(jù)長(zhǎng)度,其中31次操作包括:12次讀取操作、7次加法運(yùn)算、7次乘法操作以及5次寫(xiě)入操作(以上操作次數(shù)可根據(jù)式(8)計(jì)算得到);
(2)優(yōu)化后總操作次數(shù)=31×插值數(shù)據(jù)長(zhǎng)度-18×插值數(shù)據(jù)長(zhǎng)度/數(shù)據(jù)重復(fù)周期,其中數(shù)據(jù)重復(fù)周期=維度長(zhǎng)度×數(shù)據(jù)間隔。
因?yàn)榫S度長(zhǎng)度和數(shù)據(jù)間隔只和當(dāng)前插值數(shù)據(jù)的維度有關(guān)且相對(duì)不變,由此可見(jiàn)插值數(shù)據(jù)長(zhǎng)度越大,提升性能就會(huì)越明顯。
利用上海聯(lián)影醫(yī)療科技有限公司320排CT拍攝的兩套患者心臟的CTA數(shù)據(jù)對(duì)本研究提出的心臟冠狀動(dòng)脈運(yùn)動(dòng)向量場(chǎng)三維插值算法進(jìn)行實(shí)驗(yàn)測(cè)試。本研究所使用的算法均在同一臺(tái)工作電腦上實(shí)現(xiàn),具體配置見(jiàn)表1。
表1 實(shí)驗(yàn)環(huán)境清單Tab.1 Experiment environment list
每套數(shù)據(jù)都包含一份CT重建后的投影數(shù)據(jù),數(shù)據(jù)層厚為0.5 mm以及相關(guān)的從投影數(shù)據(jù)中提取并處理后的一套心臟冠狀動(dòng)脈特征點(diǎn)坐標(biāo)數(shù)據(jù)。測(cè)試數(shù)據(jù)是通過(guò)心電圖探測(cè)到心臟舒張期的信號(hào)后觸發(fā)CT掃描并進(jìn)行濾波反投影變換后得到的,其中表2、表3和表4中出現(xiàn)的Phase1、Phase2、Phase3和Phase4分別表示CT每旋轉(zhuǎn)一定角度,通過(guò)上述方式所得到的數(shù)據(jù)。
表2 原始數(shù)據(jù)信息Tab.2 Original data information
表3 TPS與MBS執(zhí)行時(shí)間對(duì)比Tab.3 Time comparison between TPS and MBS
表4 TPS和MBS配準(zhǔn)誤差對(duì)比Tab.4 MSE Comparison between TPS and MBS
本實(shí)驗(yàn)主要是對(duì)insight toolkit(ITK)提供的TPS算法與本研究提出的MBS算法在執(zhí)行速度、均方誤差(mean square error,MSE)以及校準(zhǔn)后圖像效果三方面進(jìn)行對(duì)比。
(1)從執(zhí)行速度來(lái)看,通過(guò)計(jì)算機(jī)性能計(jì)數(shù)器對(duì)兩個(gè)算法的執(zhí)行時(shí)間進(jìn)行記錄得到表3的結(jié)果,從中可以發(fā)現(xiàn)TPS隨著匹配點(diǎn)的增加其執(zhí)行時(shí)間呈上升趨勢(shì)。而MBS算法的執(zhí)行時(shí)間變化不大,不受匹配點(diǎn)個(gè)數(shù)的影響,并且在執(zhí)行速度上明顯高于TPS算法。主要原因是MBS的執(zhí)行時(shí)間只和目標(biāo)維度的大小以及迭代次數(shù)有關(guān),即無(wú)論匹配點(diǎn)數(shù)是多少,都只需6次迭代從大小為16×16×12稀疏網(wǎng)格插值到516×516×284高密度網(wǎng)格。
(2)從配準(zhǔn)精度來(lái)看,隨著匹配點(diǎn)的個(gè)數(shù)增多,兩種算法的配準(zhǔn)精度都有所提高。薄板樣條插值TPS配準(zhǔn)精度一般,這很可能和它的全局配準(zhǔn)性質(zhì)有很大的關(guān)系,它不能像B樣條一樣進(jìn)行局部性的微調(diào)。而MBS插值配準(zhǔn)效果已經(jīng)達(dá)到了小于1個(gè)像素的程度,針對(duì)匹配點(diǎn)集中Phase3匹配點(diǎn)的6次迭代的配準(zhǔn)結(jié)果(見(jiàn)圖 3),MBS首次迭代就已經(jīng)達(dá)到了小于1個(gè)像素的精度并且每次迭代其精度都能夠進(jìn)一步提升。而TPS算法是通過(guò)線性方程組求解,在目標(biāo)運(yùn)動(dòng)向量場(chǎng)的維度下進(jìn)行插值得到的結(jié)果[11],無(wú)法像MBS那樣通過(guò)迭代在多個(gè)尺度上進(jìn)行累積插值。因此,MBS算法在空間局部性調(diào)整控制點(diǎn)以及多尺度樣條細(xì)化方面均具有優(yōu)勢(shì)。
圖3 MBS方法的迭代次數(shù)與配準(zhǔn)精度關(guān)系圖Fig.3 Relationship between iteration times and accuracy of MBS
(3)從處理后的圖像效果來(lái)看(見(jiàn)圖4),基于TPS插值的CTA圖像中冠狀動(dòng)脈部分沒(méi)有特別明顯的處理效果,并且在圖像整體上都有比較明顯的拉伸,甚至在心室部位帶來(lái)一定的偽影,這主要是因?yàn)門(mén)PS插值算法是在516×516×284整個(gè)圖像空間坐標(biāo)上進(jìn)行插值的;而MBS插值算法不僅消除了校準(zhǔn)圖像中冠狀動(dòng)脈部位的拖尾,而且沒(méi)有像TPS算法那樣對(duì)整幅圖像造成扭曲。并且通過(guò)圖 5和圖 6的運(yùn)動(dòng)向量場(chǎng)剖面圖也可以看到TPS的運(yùn)動(dòng)向量場(chǎng)整體上都處于扭曲的狀態(tài),而MBS只有局部發(fā)生扭曲,其它部分沒(méi)有明顯的變化。
圖4 圖像對(duì)比圖Fig.4 Image comparison result
圖5 TPS運(yùn)動(dòng)向量場(chǎng)剖面圖Fig.5 Profile of motion vector field for TPS
圖6 MBS運(yùn)動(dòng)向量場(chǎng)剖面圖Fig.6 Profile of motion vector field for MBS
本研究提出了一種基于均勻網(wǎng)格多層次B樣條CTA圖像的運(yùn)動(dòng)向量場(chǎng)快速插值方法,利用B樣條函數(shù)本身的空間局部性以及多尺度迭代的插值方式使運(yùn)動(dòng)向量場(chǎng)能夠產(chǎn)生較理想的形變,在配準(zhǔn)誤差方面達(dá)到了遠(yuǎn)小于1個(gè)像素的水平。另外,本方法采用分割矩陣的方式,將三維數(shù)據(jù)按照一維的方式進(jìn)行插值,并利用插值點(diǎn)數(shù)值周期性重復(fù),有效地過(guò)濾了重復(fù)的數(shù)據(jù),極大地提高了算法的計(jì)算速度。實(shí)驗(yàn)證明,本方法對(duì)CTA圖像去偽影有較好的效果,有較高的臨床應(yīng)用價(jià)值。
此外,為了快速得到高質(zhì)量去偽影的CTA圖像,需要將CT從多個(gè)角度采集到的數(shù)據(jù)進(jìn)行處理,這個(gè)時(shí)候需要使用提出的方法以并行的方式對(duì)生成的匹配點(diǎn)集進(jìn)行插值,這涉及到多線程開(kāi)發(fā)的方法,未來(lái)可在GPU上開(kāi)發(fā)此算法,進(jìn)一步提高診斷效率和有效性。