高強(qiáng) ，鮑園，趙俊

1 上海交通大學(xué) 生物醫(yī)學(xué)工程學(xué)院，上海市，200240

2 上海聯(lián)影醫(yī)療科技有限公司，上海市，201807

0 引言

冠心病是心血管疾病中發(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í)行效率。

1 方法

1.1 特征點(diǎn)集插值函數(shù)的構(gòu)建方法

本研究提出的多層次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)

1.2 運(yùn)動(dòng)向量場(chǎng)迭代插值過(guò)程

包含于空間網(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

1.3 分割矩陣插值方法

本研究采用了一種分割矩陣的數(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)公式：

1.4 插值算法優(yōu)化方法

過(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ì)越明顯。

2 實(shí)驗(yàn)與結(jié)果

2.1 測(cè)試數(shù)據(jù)與環(huán)境

利用上海聯(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

2.2 實(shí)驗(yàn)結(jié)果與討論

本實(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

3 結(jié)論

本研究提出了一種基于均勻網(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)一步提高診斷效率和有效性。