高強 ，鮑園，趙俊

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

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

0 引言

冠心病是心血管疾病中發(fā)病率、死亡率最高的，據(jù)統(tǒng)計在西方國家每年約有30%的人死于該疾病，并且我國也處于不斷上升的趨勢[1]。冠心病的早期診斷具有重要價值[1-2]。冠狀動脈CT血管造影（computed tomography angiography，CTA）具有成本低、無創(chuàng)、準確性和敏感性高等特點，是繼數(shù)字減影血管造影（digital subtraction angiography,DSA）之后監(jiān)測和診斷心臟疾病的一種安全、可靠的方法[1]。CTA系統(tǒng)能成為檢測心臟疾病的有效方法，得益于高速的球管旋轉(zhuǎn)、寬大的探測器以及雙源掃描等方式，減少了心臟的掃描時間，例如當今超高端的CT，可在0.28 s內(nèi)完成360o的旋轉(zhuǎn)[3]，從而提高了系統(tǒng)時間分辨率，提升了圖像質(zhì)量。但是每次掃描重建的圖像，還可能會因心臟跳動頻率比CT探測器采集頻率高而產(chǎn)生運動偽影。

非剛性圖像配準方法是解決上述問題的方法之一[4]，它主要包括心臟冠狀動脈的特征點提取，生成運動向量場以及圖像的灰度映射三部分[5]。其中，在生成運動向量場方面，樣條插值方法可解決心臟冠狀動脈運動自由度高、不確定性高等問題，是消除心臟冠狀動脈偽影的主流方法[6～8]。薄板樣條、（thin plate spline，TPS）插值方法是運動向量場插值方法的一種選擇[8]，通過此方法可以找到通過所有控制點的最短路徑，從而模擬心臟冠脈的運動軌跡，但是TPS屬于全局性的配準方法，任意一個控制點發(fā)生變化都會對全局產(chǎn)生影響，配準精度并不是太高[6]。在執(zhí)行效率方面，該方法運行時間過長，無法達到實時處理的要求。心臟冠狀動脈運動是由心臟跳動與呼吸共同作用的結果，因此運動自由度比較高，需要對運動的局部性加以考慮。傳統(tǒng)的B樣條插值方法雖然在改變控制點時只會影響局部鄰域，能夠解決局部性變形的問題，但是隨著物體模型變形程度的增加，其運動向量場配準的精度和執(zhí)行速度也會下降。

我們提出了一種基于均勻網(wǎng)格多層次B樣條的CTA圖像的運動向量場插值方法，此方法可以從多尺度角度來控制網(wǎng)格的變化，一方面它從較稀疏的網(wǎng)格著手，能夠模擬全局的彈性形變；另一方面它能夠在上一層次的基礎之上，逐漸步入較細的網(wǎng)格層次，由于控制點之間的間距變小，可以模擬局部的變形，從而提高配準精度。另外，利用分割矩陣的數(shù)學方法，可將插值時大量造成冗余計算的共享節(jié)點進行過濾，不僅兼顧了原有算法的配準精度，同時也提高了算法的執(zhí)行效率。

1 方法

1.1 特征點集插值函數(shù)的構建方法

本研究提出的多層次B樣條算法[9]是一種多分辨率多尺度的插值方法，它允許從離散的特征點集中擬合出三維自由形變的曲面。令三自變量目標控制點集：

P屬于三維運動向量場V={(x,y,z,d)|0≤x＜m,0≤y＜n,0≤z＜l}內(nèi)部。為了逼近P，構建一個均勻的三-三次B樣條函數(shù)作為P的插值函數(shù)。令A是覆蓋V區(qū)域大小為(m+3)×(n+3)×(l+3)的三維空間，Aijk網(wǎng)格是A空間網(wǎng)格中的某一控制點，那么P的插值函數(shù)為：

其中0≤t＜1。

曲面的逼近問題是確定運動向量場中的所有控制點的值A，根據(jù)式（2）可知函數(shù)值f(x,y,z)與數(shù)據(jù)點(x,y,z)鄰域內(nèi)的64個控制點相關，并且按照LEE等[9]提出的最小二乘法的思想來保證待求的V中每個控制點Aabc都在A范圍內(nèi)，最終通過以下方式計算出所有的控制點：

其中wabc=Ba(s)Bb(t)Bc(u)，wefg=Be(s)Bf(t)Bg(u)

1.2 運動向量場迭代插值過程

包含于空間網(wǎng)格A中的V可從多個尺度逐層插值，并可產(chǎn)生一系列插值函數(shù)fk，而它們的和就是逼近目標的插值函數(shù)，計算式為：

最稀疏的網(wǎng)格是由參考圖像特征點與目標圖像特征點的差計算得到的，并且可利用此網(wǎng)格的格點得到比較粗略的近似插值函數(shù)。后續(xù)細化后的控制格點可以進一步逼近和消除殘差。當經(jīng)過數(shù)次迭代得到目標尺寸的網(wǎng)格時，插值即可完成，參見圖1的插值過程。

參考網(wǎng)格V代表運動向量場的初始位置，形變網(wǎng)格相當于運動向量場的位移量。起初網(wǎng)格比較稀疏，格點間的距離較大。經(jīng)過每次迭代插值維度成倍增加，格點間距離成倍縮小，并且使用每次迭代生成的向量場（即參考網(wǎng)格與形變網(wǎng)格之和）通過式（2）得到的結果與目標匹配點的差作為下一次迭代的輸入，這樣反復N次迭代，當達到目標大小時結束。

圖1 MBS插值流程圖Fig.1 Flow chart of MBS interpolation

1.3 分割矩陣插值方法

本研究采用了一種分割矩陣的數(shù)學方法，它允許對基于均勻網(wǎng)格的B樣條生成冠狀動脈的運動向量場進行迭代插值。TPS的插值方法忽視了冠狀動脈的有效取值范圍，對整個三維空間中所有坐標點進行全局插值，而本研究的方法可以只對提取的感興趣局部區(qū)域進行插值。

其基本方法是將原有的B樣條曲線分成左右兩條曲線段，B樣條矩陣[10]如下：

其中k=0,1,…,n-3且0≤t＜1，M為B樣條基函數(shù)系數(shù)。假設一個B樣條曲線P(t)是由控制點{P0,P1,P2,P3}構成的，我們將曲線二元細分，產(chǎn)生兩條曲線PL(t)和PR(t)，并將原有t的[0,1)的取值范圍劃分為[0,0.5)和[0.5,1)。然后對這兩條曲線段進行一系列的變換，最終得到基于控制點{P0,P1,P2,P3}插值點的表達式：

通過（7），我們可以知道PL(t)和PR(t)中存在被重復插值的點，即。這些重復點可以從任意一個曲線段中選取，最終可得到插值點公式：

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

過濾掉分割矩陣插值方法中產(chǎn)生的重復的插值點是提高算法執(zhí)行速度的關鍵。根據(jù)算法中的插值方法可以將待插值點的位置通過左右相鄰的兩個或者三個控制點通過式（8）計算得到，這樣可以不用考慮控制點的空間維度只需將控制點集P按照一維的排列方式表示并額外復制3組相同的點集，統(tǒng)一定義為P1、P2、P3以及P4，其中點集中的每個元素以D來表示，而、'以及'是待插值的點集（見圖2）。

圖2 待插值控制網(wǎng)格數(shù)據(jù)布局Fig.2 Layout of control grid data to be interpolated

每組點集按照同一點集中點與點之間的相鄰關系進行縱向排列，這樣可以對數(shù)據(jù)進行很好的并行處理，對于4組控制點集只需考慮彼此間對應位置即可。但是存在一個問題，那就是會出現(xiàn)大量插值點被重復運算導致計算資源的浪費（如圖2中的第一列P2和P3中的D2和D3以及第二列P1和P2中的D2和D3），這主要是和分割矩陣算法特點有關。為了解決此問題，不妨先只對進行插值，當完成了當前維度的運算之后，再對以及做相同維度數(shù)據(jù)間隔長度的插值。根據(jù)三維空間的幾何關系可知X方向插值數(shù)據(jù)間隔為1，Y方向插值數(shù)據(jù)間隔為X方向元素個數(shù)，Z方向插值數(shù)據(jù)間隔為XY平面的元素個數(shù)。通過上述方法得到如下結果：

（1）優(yōu)化前總操作次數(shù)=31×插值數(shù)據(jù)長度，其中31次操作包括：12次讀取操作、7次加法運算、7次乘法操作以及5次寫入操作（以上操作次數(shù)可根據(jù)式（8）計算得到）；

（2）優(yōu)化后總操作次數(shù)=31×插值數(shù)據(jù)長度－18×插值數(shù)據(jù)長度/數(shù)據(jù)重復周期，其中數(shù)據(jù)重復周期=維度長度×數(shù)據(jù)間隔。

因為維度長度和數(shù)據(jù)間隔只和當前插值數(shù)據(jù)的維度有關且相對不變，由此可見插值數(shù)據(jù)長度越大，提升性能就會越明顯。

2 實驗與結果

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

利用上海聯(lián)影醫(yī)療科技有限公司320排CT拍攝的兩套患者心臟的CTA數(shù)據(jù)對本研究提出的心臟冠狀動脈運動向量場三維插值算法進行實驗測試。本研究所使用的算法均在同一臺工作電腦上實現(xiàn)，具體配置見表1。

表1 實驗環(huán)境清單Tab.1 Experiment environment list

每套數(shù)據(jù)都包含一份CT重建后的投影數(shù)據(jù)，數(shù)據(jù)層厚為0.5 mm以及相關的從投影數(shù)據(jù)中提取并處理后的一套心臟冠狀動脈特征點坐標數(shù)據(jù)。測試數(shù)據(jù)是通過心電圖探測到心臟舒張期的信號后觸發(fā)CT掃描并進行濾波反投影變換后得到的，其中表2、表3和表4中出現(xiàn)的Phase1、Phase2、Phase3和Phase4分別表示CT每旋轉(zhuǎn)一定角度，通過上述方式所得到的數(shù)據(jù)。

表2 原始數(shù)據(jù)信息Tab.2 Original data information

表3 TPS與MBS執(zhí)行時間對比Tab.3 Time comparison between TPS and MBS

表4 TPS和MBS配準誤差對比Tab.4 MSE Comparison between TPS and MBS

2.2 實驗結果與討論

本實驗主要是對insight toolkit（ITK）提供的TPS算法與本研究提出的MBS算法在執(zhí)行速度、均方誤差（mean square error，MSE）以及校準后圖像效果三方面進行對比。

（1）從執(zhí)行速度來看，通過計算機性能計數(shù)器對兩個算法的執(zhí)行時間進行記錄得到表3的結果，從中可以發(fā)現(xiàn)TPS隨著匹配點的增加其執(zhí)行時間呈上升趨勢。而MBS算法的執(zhí)行時間變化不大，不受匹配點個數(shù)的影響，并且在執(zhí)行速度上明顯高于TPS算法。主要原因是MBS的執(zhí)行時間只和目標維度的大小以及迭代次數(shù)有關，即無論匹配點數(shù)是多少，都只需6次迭代從大小為16×16×12稀疏網(wǎng)格插值到516×516×284高密度網(wǎng)格。

（2）從配準精度來看，隨著匹配點的個數(shù)增多，兩種算法的配準精度都有所提高。薄板樣條插值TPS配準精度一般，這很可能和它的全局配準性質(zhì)有很大的關系，它不能像B樣條一樣進行局部性的微調(diào)。而MBS插值配準效果已經(jīng)達到了小于1個像素的程度，針對匹配點集中Phase3匹配點的6次迭代的配準結果（見圖 3），MBS首次迭代就已經(jīng)達到了小于1個像素的精度并且每次迭代其精度都能夠進一步提升。而TPS算法是通過線性方程組求解，在目標運動向量場的維度下進行插值得到的結果[11]，無法像MBS那樣通過迭代在多個尺度上進行累積插值。因此，MBS算法在空間局部性調(diào)整控制點以及多尺度樣條細化方面均具有優(yōu)勢。

圖3 MBS方法的迭代次數(shù)與配準精度關系圖Fig.3 Relationship between iteration times and accuracy of MBS

（3）從處理后的圖像效果來看（見圖4），基于TPS插值的CTA圖像中冠狀動脈部分沒有特別明顯的處理效果，并且在圖像整體上都有比較明顯的拉伸，甚至在心室部位帶來一定的偽影，這主要是因為TPS插值算法是在516×516×284整個圖像空間坐標上進行插值的；而MBS插值算法不僅消除了校準圖像中冠狀動脈部位的拖尾，而且沒有像TPS算法那樣對整幅圖像造成扭曲。并且通過圖 5和圖 6的運動向量場剖面圖也可以看到TPS的運動向量場整體上都處于扭曲的狀態(tài)，而MBS只有局部發(fā)生扭曲，其它部分沒有明顯的變化。

圖4 圖像對比圖Fig.4 Image comparison result

圖5 TPS運動向量場剖面圖Fig.5 Profile of motion vector field for TPS

圖6 MBS運動向量場剖面圖Fig.6 Profile of motion vector field for MBS

3 結論

本研究提出了一種基于均勻網(wǎng)格多層次B樣條CTA圖像的運動向量場快速插值方法，利用B樣條函數(shù)本身的空間局部性以及多尺度迭代的插值方式使運動向量場能夠產(chǎn)生較理想的形變，在配準誤差方面達到了遠小于1個像素的水平。另外，本方法采用分割矩陣的方式，將三維數(shù)據(jù)按照一維的方式進行插值，并利用插值點數(shù)值周期性重復，有效地過濾了重復的數(shù)據(jù)，極大地提高了算法的計算速度。實驗證明，本方法對CTA圖像去偽影有較好的效果，有較高的臨床應用價值。

此外，為了快速得到高質(zhì)量去偽影的CTA圖像，需要將CT從多個角度采集到的數(shù)據(jù)進行處理，這個時候需要使用提出的方法以并行的方式對生成的匹配點集進行插值，這涉及到多線程開發(fā)的方法，未來可在GPU上開發(fā)此算法，進一步提高診斷效率和有效性。