(1 北京航空航天大學(xué)計算機學(xué)院，北京 100036； 2 青島大學(xué)附屬醫(yī)院，山東省數(shù)學(xué)醫(yī)學(xué)與計算機手術(shù)重點實驗室，山東省高等學(xué)校數(shù)學(xué)醫(yī)學(xué)臨床診療與營養(yǎng)健康協(xié)同創(chuàng)新中心)

手術(shù)導(dǎo)航(IGS)是指醫(yī)生在術(shù)前利用醫(yī)學(xué)影像設(shè)備和計算機圖形學(xué)的方法，對患者多模態(tài)的圖像數(shù)據(jù)進行3維重建和可視化處理，獲得3維模型，制定合理、定量的手術(shù)計劃，開展術(shù)前模擬；在術(shù)中通過配準操作，把3維模型與患者的實際體位、空間中手術(shù)器械的實時位置統(tǒng)一在一個坐標(biāo)系下，并利用3維定位系統(tǒng)對手術(shù)器械在空間中的位置實時采集并顯示，醫(yī)生通過觀察3維模型中手術(shù)器械與病變部位的相對位置關(guān)系，對病人進行導(dǎo)航手術(shù)治療。IGS通?？梢詭椭t(yī)生更安全地執(zhí)行手術(shù)，并且對病人造成的損傷更小，使得病人可以更快速恢復(fù)。在現(xiàn)今醫(yī)療條件下，IGS已經(jīng)被廣泛應(yīng)用于多種手術(shù)當(dāng)中，例如在經(jīng)皮冠狀動脈介入治療[1-3](PCI)手術(shù)中利用X線造影圖像幫助醫(yī)生執(zhí)行手術(shù)，定位病變位置。然而在PCI等微創(chuàng)介入手術(shù)過程中，通常導(dǎo)管是可見的，血管并不可見。外科醫(yī)生需要知道導(dǎo)管相對于具體目標(biāo)的位置。為了達到該目的，傳統(tǒng)的做法是使用造影劑，使得血管可見。然而造影劑具有毒性，不可以持續(xù)使用[4-5]。同時使用造影劑還會存在外滲等問題[6-7]。因此使用從術(shù)中的X線圖像中提取出來的導(dǎo)管與術(shù)前從病人的CT血管造影(CTA)數(shù)據(jù)中提取的血管樹配準來定位導(dǎo)管成為了解決上述問題的一種方法[8]。目前有兩種新型的導(dǎo)管配準方法，基于隱式馬爾科夫模型(HMM)的導(dǎo)管配準方法[9]和基于形狀相似度的導(dǎo)管配準方法[10]。為了探討兩種配準算法配準精度以及時間代價，本研究對荷蘭Rotterdam大學(xué)醫(yī)學(xué)中心、法國Henri Mondor大學(xué)及意大利Circolo醫(yī)院聯(lián)合提供的開源數(shù)據(jù)集分別進行了實驗。

1 資料與方法

1.1 資料來源

選取荷蘭Rotterdam大學(xué)醫(yī)學(xué)中心、法國Henri Mondor大學(xué)及意大利Circolo醫(yī)院聯(lián)合提供的開源數(shù)據(jù)集。該數(shù)據(jù)集開源了從一次真實微創(chuàng)介入手術(shù)中提取出來的血管樹數(shù)據(jù)以及在該次介入手術(shù)中導(dǎo)管在每幀X線圖像上的位置信息，共74幀。該數(shù)據(jù)集開源的血管樹是一棵具有728個采樣點、102個分支段的復(fù)雜血管樹。這棵血管樹具有采樣點數(shù)量眾多、分支情況復(fù)雜等特點。因此具備足夠的復(fù)雜性來測試基于HMM的導(dǎo)管配準算法和基于形狀相似度的導(dǎo)管配準算法的精度與時間消耗。同時通過該數(shù)據(jù)集提供的導(dǎo)管在每幀X線圖像上的位置信息，來模擬導(dǎo)管在血管中移動情況。實驗在一臺CPU為Intel Core i5-3210M，內(nèi)存為8 GB的筆記本電腦上進行。對于兩種配準算法，我們均使用Powell優(yōu)化算法求解。

1.2 基于HMM的導(dǎo)管配準算法

HMM是一個具有N個狀態(tài)的系統(tǒng)S={s1，s2…sN}(圖1)。HMM會根據(jù)狀態(tài)遷移概率以及當(dāng)前的觀測值在離散的時間點t上發(fā)生的狀態(tài)改變。aij表示HMM從狀態(tài)si遷移到sj的概率，且aij≥0，∑jaij=1。

圖1 HMM

根據(jù)RABINER[11]提出的Viterbi算法，該算法在時間t=t1時選取從時間t=0到t1時間段內(nèi)的一條最佳路徑，又稱為Viterbi路徑。Viterbi路徑是指從任意狀態(tài)開始到達狀態(tài)si可能性最高的路徑。在時間t內(nèi)到達狀態(tài)si的Viterbi路徑的概率為δt(i)。對于HMM，每一個狀態(tài)都有一個初始概率，所有的初始概率構(gòu)成一個集合Π={π1，π2…πN}，且∑πi=1。因此有

δ0(i)=πi

而時間t內(nèi)的Viterbi路徑的概率δt(i)可以遞歸計算

其中Ot(i)為HMM在時間t以及狀態(tài)為Si的概率。

在基于HMM的導(dǎo)管配準算法中，狀態(tài)Si表示導(dǎo)管在血管樹中的位置。狀態(tài)遷移概率aij表示導(dǎo)管從血管樹的一個位置pi移動到pj的概率。在一個時間段內(nèi)，導(dǎo)管移動到pi附近的位置的概率應(yīng)該大于移動到較遠位置的概率。因此遷移概率是一個與pipj之間距離D(pi，pj)有關(guān)的量，可以通過以下方式計算，

在基于HMM的導(dǎo)管配準算法中Ot(i)通過將2D的導(dǎo)管投影位置數(shù)據(jù)與3D血管樹位置數(shù)據(jù)進行一次3D-2D的配準來計算。

其中E(Ct，V)是時刻t導(dǎo)管Ct與血管樹的配準結(jié)果。

我們選取所有可能的血管段Vi與導(dǎo)管Ct配準，并將最好的匹配作為導(dǎo)管與血管的配準結(jié)果。

M(Ct，Vi)是導(dǎo)管Ct與血管段Vi的配準結(jié)果。

導(dǎo)管Ct與血管段Vi配準的數(shù)學(xué)本質(zhì)是一個二次優(yōu)化問題，

求一個變換T，使得配準結(jié)果M(Ct，Vi)最小。其中F(c，Vi，T)是配準度量項，

其中，v是血管段Vi中距離導(dǎo)管點c最近的血管點。變換T將v做一個剛體變換并投影到成像平面上。

1.3 基于形狀相似度的導(dǎo)管配準算法

基于形狀相似度的導(dǎo)管配準算法分為兩步。首先，使用一個形狀相似度度量去尋找與導(dǎo)管形狀最相似的血管段。第二步執(zhí)行導(dǎo)管與該血管段的剛性配準。血管樹被表示成一個點集以及一個邊集合，G=(P，ε)，其中P是血管點的集合，ε是血管樹中邊的集合。血管段V(p)={p，p1，…pn}表示從點p沿血管路徑到血管根結(jié)點的所有血管點的集合。對于兩條曲線，我們認為兩條曲線在相同位置點上的切向量越接近，則兩條曲線越相似。因此，對于血管中的一點p∈P，相似性度量定義為，

其中Cl是2D導(dǎo)管的長度，Tgt(C(u))表示導(dǎo)管在位置u處的切向量，Tgt(Vproj(p，u))表示血管段p的投影在位置u處的切向量。從上式可以看出S(p)∈[0，Cl]。因此，Cl為相似度的最大值。我們將計算所有從血管根結(jié)點到血管葉子結(jié)點l的血管段與導(dǎo)管的相似度。這種血管段稱為葉血管段，葉血管段與導(dǎo)管的相似度定義為葉血管段中任意一段血管與導(dǎo)管相似度的最大值，

我們選擇相似度最大的k段葉血管段進行后面的導(dǎo)管與血管段的配準。

配準方程度量導(dǎo)管中的點與葉血管段中的最近點距離之和。導(dǎo)管尖端C1與葉血管段l投影中的最近點匹配，其中這些葉血管段是之前選取的相似度最高的那幾段葉血管段。

其中p是葉血管段投影中距離導(dǎo)管尖端的最近點。后續(xù)的每一個導(dǎo)管點Ci在與前一個導(dǎo)管點Ci-1匹配的血管點pi-1一定距離h內(nèi)選取距離Ci最近的血管點pi作為匹配點。

最終的配準結(jié)果可以通過以下方式計算，

D(Ci,l,T,pi-1)

其中W(x)是一個權(quán)重項。L(ci，c1)是導(dǎo)管點ci沿導(dǎo)管路徑到導(dǎo)管尖端c1的距離。因為在配準過程中，越靠近導(dǎo)管尖端的點，配準精度要求越高，也越重要。因此我們根據(jù)導(dǎo)管點沿導(dǎo)管路徑到導(dǎo)管尖端的距離給出一個權(quán)重項，

這個配準度量M足夠的快速，因為它僅僅尋找特定領(lǐng)域中的最近點。最后，最佳變換T是那些血管段與導(dǎo)管配準的結(jié)果中相似度最高的配準結(jié)果中的變換T。

1.4 觀察指標(biāo)

從荷蘭Rotterdam大學(xué)醫(yī)學(xué)中心、法國Henri Mondor大學(xué)及意大利Circolo醫(yī)院聯(lián)合提供的74幀開源導(dǎo)管數(shù)據(jù)集中選取第一幀導(dǎo)管數(shù)據(jù)，修改不同的配準參數(shù)σs、σa進行配準實驗，并記錄配準的計算時間及成對點間的平均距離。

2 結(jié) 果

2.1 基于HMM的導(dǎo)管配準算法實驗結(jié)果

調(diào)整配準參數(shù)σa、σs得到如下實驗結(jié)果，見表1。根據(jù)實驗1～5，調(diào)整配準參數(shù)σa、σs并不會影響成對點平均距離，即不會影響最后的配準精度。根據(jù)實驗1～3，σa會影響配準時間，在本數(shù)據(jù)集的實驗中σa的最佳值是8.0。根據(jù)實驗1、4、5，參數(shù)σs對配準時間的影響非常小，最終時間在60 ms左右。

2.2 基于形狀相似度的導(dǎo)管配準算法實驗結(jié)果

調(diào)整葉血管段數(shù)、領(lǐng)域范圍得到以下實驗結(jié)果，見表2。減少葉血管段數(shù)會明顯降低配準時間。降低領(lǐng)域范圍也會明顯降低配準時間，但是會對配準精度產(chǎn)生影響。

表1 基于HMM的導(dǎo)管配準算法實驗結(jié)果

表2 基于形狀相似度的導(dǎo)管配準算法實驗結(jié)果

2.3 連續(xù)幀情況下兩種配準算法的實驗結(jié)果

我們在連續(xù)的20幀情況下，分別針對兩種配準算法進行了導(dǎo)管配準實驗，并統(tǒng)計了算法的平均運行時間和成對對應(yīng)點的平均距離。根據(jù)我們的實驗結(jié)果表明，基于HMM的導(dǎo)管配準算法平均運行時間為73 ms，而基于形狀相似度的導(dǎo)管配準算法平均運行時間為832 ms。基于HMM的導(dǎo)管配準中成對點平均距離為3.042 31 mm，基于形狀相似度的導(dǎo)管配準算法中成對點的平均距離為2.506 33 mm。在配準精度方面，基于HMM的導(dǎo)管配準算法的配準精度要低于基于形狀相似度的導(dǎo)管配準算法。

3 討 論

IGS系統(tǒng)[12-17]是當(dāng)前醫(yī)學(xué)圖像處理的熱門研究領(lǐng)域，它融合了醫(yī)學(xué)[18-20]、圖像處理[21-22]及計算機圖形學(xué)[23-24]等各個方面的知識，吸引了中外許多的研究者。近年來IGS系統(tǒng)的研究者越來越多地將研究重心放在IGS系統(tǒng)中的配準技術(shù)上[25-26]。配準是IGS系統(tǒng)中的關(guān)鍵技術(shù)，負責(zé)手術(shù)器械以及病變器官的定位[27-29]。一般情況下，在執(zhí)行導(dǎo)管配準之前，必須先執(zhí)行分割算法，將導(dǎo)管從2D X線圖像中分割出來。這也是一個非常具有挑戰(zhàn)性的任務(wù)。HEIBEL等[30]提出了一種從肝臟手術(shù)中的X線圖像中提取導(dǎo)管的算法。WAGNER等[31]提出了一種從連續(xù)X線圖像幀中提取導(dǎo)管的算法。不同的導(dǎo)管提取算法適用于不同的手術(shù)場景。導(dǎo)管提取算法精度上的差異也會對后續(xù)導(dǎo)管配準算法的精度產(chǎn)生一定的影響。

本研究首先對基于HMM的導(dǎo)管配準方法和基于形狀相似度的導(dǎo)管配準算法進行實驗，測試出能使兩種算法到達最佳結(jié)果的參數(shù)，再分別以最佳參數(shù)進行連續(xù)幀實驗，并記錄在連續(xù)幀情況下兩種算法的平均運行時間和成對點間的平均距離。本研究結(jié)果顯示，基于HMM的導(dǎo)管配準方法具有配準結(jié)果穩(wěn)定、運算速度快等特點?；贖MM的導(dǎo)管配準方法會充分考慮之前幀的配準情況，考慮到導(dǎo)管在幀間的運動穩(wěn)定性和連貫性，因此在連續(xù)幀情況下同樣可以保持較快的運算速度、極少的運行時間、較高的配準精度和配準穩(wěn)定性。本研究結(jié)果顯示，基于形狀相似度的導(dǎo)管配準算法比基于HMM的導(dǎo)管配準算法具有更高的配準精度，但是以此為代價是消耗更多的計算時間。本研究結(jié)果同時顯示，在連續(xù)幀的情況下，基于形狀相似度的導(dǎo)管配準算法所需要的計算時間約為基于HMM的導(dǎo)管配準算法的12倍?？梢钥闯鲈诰C合情況下，基于形狀相似度的導(dǎo)管配準算法比基于HMM的導(dǎo)管配準算法具有更高的配準精度，而基于HMM的導(dǎo)管配準算法比基于形狀相似度的導(dǎo)管配準算法擁有更快的計算速度。

綜上所述，雖然基于形狀相似度的導(dǎo)管配準算法具有更高的配準精度，然而其過長的配準時間難以達到在臨床中實時進行導(dǎo)管配準的要求；而基于HMM的導(dǎo)管配準算法雖然在精度上比基于形狀相似度的導(dǎo)管配準算法要弱，但是精度上基本滿足了臨床要求，而其配準速度也可以滿足在臨床中導(dǎo)管實時配準這個需求。因此基于HMM的導(dǎo)管配準算法具有更高的臨床價值，值得臨床推廣應(yīng)用。