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

手術(shù)導(dǎo)航(IGS)是指醫(yī)生在術(shù)前利用醫(yī)學(xué)影像設(shè)備和計(jì)算機(jī)圖形學(xué)的方法，對患者多模態(tài)的圖像數(shù)據(jù)進(jìn)行3維重建和可視化處理，獲得3維模型，制定合理、定量的手術(shù)計(jì)劃，開展術(shù)前模擬；在術(shù)中通過配準(zhǔn)操作，把3維模型與患者的實(shí)際體位、空間中手術(shù)器械的實(shí)時(shí)位置統(tǒng)一在一個(gè)坐標(biāo)系下，并利用3維定位系統(tǒng)對手術(shù)器械在空間中的位置實(shí)時(shí)采集并顯示，醫(yī)生通過觀察3維模型中手術(shù)器械與病變部位的相對位置關(guān)系，對病人進(jì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)的位置。為了達(dá)到該目的，傳統(tǒng)的做法是使用造影劑，使得血管可見。然而造影劑具有毒性，不可以持續(xù)使用[4-5]。同時(shí)使用造影劑還會存在外滲等問題[6-7]。因此使用從術(shù)中的X線圖像中提取出來的導(dǎo)管與術(shù)前從病人的CT血管造影(CTA)數(shù)據(jù)中提取的血管樹配準(zhǔn)來定位導(dǎo)管成為了解決上述問題的一種方法[8]。目前有兩種新型的導(dǎo)管配準(zhǔn)方法，基于隱式馬爾科夫模型(HMM)的導(dǎo)管配準(zhǔn)方法[9]和基于形狀相似度的導(dǎo)管配準(zhǔn)方法[10]。為了探討兩種配準(zhǔn)算法配準(zhǔn)精度以及時(shí)間代價(jià)，本研究對荷蘭Rotterdam大學(xué)醫(yī)學(xué)中心、法國Henri Mondor大學(xué)及意大利Circolo醫(yī)院聯(lián)合提供的開源數(shù)據(jù)集分別進(jìn)行了實(shí)驗(yàn)。

1 資料與方法

1.1 資料來源

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

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

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

圖1 HMM

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

δ0(i)=πi

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

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

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

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

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

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

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

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

求一個(gè)變換T，使得配準(zhǔn)結(jié)果M(Ct，Vi)最小。其中F(c，Vi，T)是配準(zhǔn)度量項(xiàng)，

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

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

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

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

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

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

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

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

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

其中W(x)是一個(gè)權(quán)重項(xiàng)。L(ci，c1)是導(dǎo)管點(diǎn)ci沿導(dǎo)管路徑到導(dǎo)管尖端c1的距離。因?yàn)樵谂錅?zhǔn)過程中，越靠近導(dǎo)管尖端的點(diǎn)，配準(zhǔn)精度要求越高，也越重要。因此我們根據(jù)導(dǎo)管點(diǎn)沿導(dǎo)管路徑到導(dǎo)管尖端的距離給出一個(gè)權(quán)重項(xiàng)，

這個(gè)配準(zhǔn)度量M足夠的快速，因?yàn)樗鼉H僅尋找特定領(lǐng)域中的最近點(diǎn)。最后，最佳變換T是那些血管段與導(dǎo)管配準(zhǔn)的結(jié)果中相似度最高的配準(zhǔn)結(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ù)，修改不同的配準(zhǔn)參數(shù)σs、σa進(jìn)行配準(zhǔn)實(shí)驗(yàn)，并記錄配準(zhǔn)的計(jì)算時(shí)間及成對點(diǎn)間的平均距離。

2 結(jié) 果

2.1 基于HMM的導(dǎo)管配準(zhǔn)算法實(shí)驗(yàn)結(jié)果

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

2.2 基于形狀相似度的導(dǎo)管配準(zhǔn)算法實(shí)驗(yàn)結(jié)果

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

表1 基于HMM的導(dǎo)管配準(zhǔn)算法實(shí)驗(yàn)結(jié)果

表2 基于形狀相似度的導(dǎo)管配準(zhǔn)算法實(shí)驗(yàn)結(jié)果

2.3 連續(xù)幀情況下兩種配準(zhǔn)算法的實(shí)驗(yàn)結(jié)果

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

3 討 論

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

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

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