張春雷，戴 麗，劉 宇，李 鶴

(東北大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽(yáng) 110819)

隨著計(jì)算機(jī)技術(shù)、圖像處理技術(shù)與臨床醫(yī)學(xué)的深度融合，手術(shù)導(dǎo)航系統(tǒng)應(yīng)運(yùn)而生，其能夠?qū)崿F(xiàn)手術(shù)工具對(duì)患者靶點(diǎn)的精確定位，并為醫(yī)生提供穩(wěn)定、可靠的可視化導(dǎo)航系統(tǒng)，目前已在諸多外科手術(shù)中展現(xiàn)出巨大應(yīng)用價(jià)值和市場(chǎng)潛力.患者配準(zhǔn)是手術(shù)導(dǎo)航系統(tǒng)的關(guān)鍵技術(shù)，其目標(biāo)是獲得患者實(shí)體空間與其三維虛擬模型所在圖像空間的剛性變換關(guān)系，以實(shí)現(xiàn)術(shù)前數(shù)據(jù)與術(shù)中數(shù)據(jù)的位置統(tǒng)一，配準(zhǔn)精度和效率在很大程度上決定了導(dǎo)航系統(tǒng)優(yōu)劣和手術(shù)成敗.實(shí)現(xiàn)患者配準(zhǔn)需借助分別在圖像空間和患者空間中描述的至少三組相同的特征基準(zhǔn)點(diǎn)[1]，再利用配準(zhǔn)算法計(jì)算完成.根據(jù)基準(zhǔn)點(diǎn)的配對(duì)情況，患者配準(zhǔn)方法可分為成對(duì)基準(zhǔn)法和非成對(duì)基準(zhǔn)法.

成對(duì)基準(zhǔn)法通常需將易在醫(yī)學(xué)影像中辨識(shí)和采集的醫(yī)學(xué)標(biāo)記物植入患者骨骼內(nèi)或粘于體表，其配準(zhǔn)精度通?？蛇_(dá)0.5～1.5 mm[2].Chen等[3]在口腔種植系統(tǒng)中將鈦釘作為醫(yī)學(xué)標(biāo)記物，得到的基準(zhǔn)配準(zhǔn)誤差和目標(biāo)配準(zhǔn)誤差分為1.12 mm和1.35 mm，但這種侵入式標(biāo)記物會(huì)使患者的痛苦和心理負(fù)擔(dān)倍增.Lin等[4]設(shè)計(jì)了一種自定義基準(zhǔn)標(biāo)記并將其置于患者頭部表面，提出了一種基于光學(xué)跟蹤系統(tǒng)的實(shí)時(shí)自動(dòng)患者配準(zhǔn)方法，平均配準(zhǔn)誤差約為0.7 mm，然而，采用這種方法時(shí)，一旦標(biāo)記發(fā)生移位或其重建模型出現(xiàn)缺損，就會(huì)造成較大的配準(zhǔn)誤差.

非成對(duì)基準(zhǔn)法通常通過(guò)點(diǎn)云數(shù)據(jù)描述患者空間和圖像空間中的共同解剖結(jié)構(gòu)特征，兩組點(diǎn)云非一一對(duì)應(yīng)關(guān)系.Chen等[5]利用顱頜面骨模型表面點(diǎn)云進(jìn)行配準(zhǔn)時(shí)的目標(biāo)誤差可達(dá)1mm左右.Kong等[6]以人工拾點(diǎn)方式在圖像空間中采集骨模型表面點(diǎn)云來(lái)大致描述骨模型的局部外廓，但配準(zhǔn)精度低、穩(wěn)定性差，且操作極為耗時(shí).最常用的三維點(diǎn)云配準(zhǔn)算法是由Besl等[7]提出的迭代最近點(diǎn)(iterative closet point, ICP)算法，其核心思想是基于最小二乘原理，迭代尋找兩組點(diǎn)云的最優(yōu)剛性變換關(guān)系.但I(xiàn)CP算法還存在諸多缺陷，若兩組點(diǎn)云初始位姿相差較大，且沒(méi)有良好的初始變換值，其配準(zhǔn)結(jié)果將極易陷入局部最優(yōu)解；最近點(diǎn)搜索過(guò)程嚴(yán)重影響計(jì)算效率，由此提出的基于多維二叉樹(shù)的最近點(diǎn)搜索策略[8]和加速ICP算法[9]雖然大幅提升了最近點(diǎn)搜索速度，但這會(huì)額外增加創(chuàng)建樹(shù)結(jié)構(gòu)的時(shí)間.此外，聚類ICP算法[10]、魯棒尺度ICP算法[11]和主成分分析法[12]等基于復(fù)雜點(diǎn)云幾何特征的配準(zhǔn)算法更適于解決規(guī)模大、輪廓特征明顯的點(diǎn)云匹配問(wèn)題.因受實(shí)際手術(shù)條件限制，很難在患者實(shí)體上采集和構(gòu)建可靠且?guī)缀翁卣髫S富的大規(guī)模點(diǎn)云，因此，傳統(tǒng)ICP算法及其變體算法均難以直接應(yīng)用于手術(shù)導(dǎo)航系統(tǒng)的患者配準(zhǔn)中，且需要找到一種兼顧精度高、效率高和可操作性好的患者配準(zhǔn)方法.

本文結(jié)合外科手術(shù)導(dǎo)航系統(tǒng)的應(yīng)用要求和ICP算法的特點(diǎn)，以提升配準(zhǔn)精度和可操作性為目標(biāo)，提出基于三點(diǎn)法初始配準(zhǔn)和ICP算法精確配準(zhǔn)的患者配準(zhǔn)方法.實(shí)驗(yàn)結(jié)果表明，該方法簡(jiǎn)便易行、結(jié)果穩(wěn)定，具備良好的配準(zhǔn)精度和配準(zhǔn)效果.

1 患者配準(zhǔn)問(wèn)題描述及數(shù)據(jù)獲取

手術(shù)導(dǎo)航系統(tǒng)主要分為3個(gè)子系統(tǒng)，即光學(xué)定位子系統(tǒng)、圖像子系統(tǒng)和患者子系統(tǒng).光學(xué)定位子系統(tǒng)(O)選用加拿大NDI公司Polaris Vega被動(dòng)式光學(xué)三維定位系統(tǒng)，其工作組件包括位置傳感器和被動(dòng)式定位工具(探針工具和定位剛體).其中，定位工具上固定4個(gè)位于同一平面但不共線的反射標(biāo)記球，位置傳感器能追蹤并返回每個(gè)標(biāo)記球的形心坐標(biāo)值或組合工具中心點(diǎn)的位姿信息.患者子系統(tǒng)(P)為治療對(duì)象，將定位剛體與患者保持固定來(lái)配合手術(shù)執(zhí)行.圖像子系統(tǒng)(I)用于患者的CT/MRI數(shù)據(jù)和三維模型的處理及可視化，基于3D Slicer醫(yī)學(xué)圖像處理平臺(tái)實(shí)現(xiàn).各子系統(tǒng)的坐標(biāo)系轉(zhuǎn)換關(guān)系如圖1所示.

各子系統(tǒng)坐標(biāo)系轉(zhuǎn)換矩陣的求解過(guò)程如下:

(1)

(2)

1) 將定位剛體固定于患者骨骼并一同進(jìn)行CT掃描，然后在3D Slicer中利用CT數(shù)據(jù)重建出其三維模型，得到圖像空間I，并提取完整的骨模型表面點(diǎn)云MI，樣本數(shù)通常可達(dá)幾萬(wàn)至幾十萬(wàn).

2) 在患者實(shí)體空間P中，利用探針工具尖端觸碰若干次患者骨骼表面，采集得到點(diǎn)云VP，樣本數(shù)通常僅能達(dá)到幾十個(gè)，這是由于在實(shí)際手術(shù)中會(huì)受到骨骼解剖結(jié)構(gòu)的可暴露范圍限制.

由此，即可將患者配準(zhǔn)問(wèn)題轉(zhuǎn)化為密集目標(biāo)點(diǎn)云MI(圖像空間)與稀疏源點(diǎn)云MP(患者空間)之間的配準(zhǔn)問(wèn)題，這兩組點(diǎn)云的樣本數(shù)相差較大，且輪廓和形狀等幾何特性差異懸殊.

2 配準(zhǔn)方法

2.1 基于三點(diǎn)法的初始配準(zhǔn)

將定位剛體(固定于患者骨骼)上的反光標(biāo)記球作為圖像空間I和患者空間P中的共同特征基準(zhǔn).

(3)

(4)

在光學(xué)定位子系統(tǒng)中，定位剛體上每個(gè)標(biāo)記球的球心坐標(biāo)均是在位置傳感器坐標(biāo)系O中實(shí)時(shí)更新和描述的.與上述方法同理，可獲得局部坐

(5)

設(shè)點(diǎn)F4在圖像空間I和患者空間P的齊次坐標(biāo)分別為fI和fP，則定義初始配準(zhǔn)誤差eInit為

eInit=|fI-TInit·fP|.

(6)

由于CT圖像存在偽影且三維重建時(shí)易導(dǎo)致標(biāo)記球模型發(fā)生畸變或缺損，而依賴于成對(duì)基準(zhǔn)點(diǎn)的配準(zhǔn)方法又對(duì)基準(zhǔn)點(diǎn)拾取精度要求較高，因此僅將上述方法用于初始配準(zhǔn).

2.2 基于ICP算法的精確配準(zhǔn)

設(shè)MP和MI的齊次坐標(biāo)矩陣分別為V′和U:

(7)

(8)

利用初始配準(zhǔn)矩陣TInit將源點(diǎn)云V′映射至圖像空間中，得到新的源點(diǎn)云V，即

V=TInit·V′.

(9)

執(zhí)行精確配準(zhǔn)的目標(biāo)是獲得源點(diǎn)云V相對(duì)于目標(biāo)點(diǎn)云U的轉(zhuǎn)換矩陣TICP，其核心思想是根據(jù)設(shè)定的目標(biāo)函數(shù)，基于最小二乘原理迭代估計(jì)點(diǎn)云之間的最優(yōu)匹配關(guān)系，方法如下:

1) 輸入源點(diǎn)云V={vi,i=1,2,…,m}和目標(biāo)點(diǎn)云U={uj,j=1,2,…,n}，其中m?n.

2) 目標(biāo)點(diǎn)云U的抽樣:點(diǎn)云U中樣本數(shù)龐大且分布密集，但其絕大部分與源點(diǎn)云V是不匹配的，從而成為干擾噪聲，不僅會(huì)降低配準(zhǔn)精度，還會(huì)嚴(yán)重增加最近點(diǎn)對(duì)搜索過(guò)程的計(jì)算耗時(shí).因此，對(duì)目標(biāo)點(diǎn)云U進(jìn)行合理抽樣簡(jiǎn)化是十分必要的.

以源點(diǎn)云V中各點(diǎn)vi為中心，建立以r為半徑的球區(qū)域Si；保留目標(biāo)點(diǎn)云U中所有位于區(qū)域Si中的樣本(nr個(gè))，得到更新后的目標(biāo)點(diǎn)云U(rh)(h為迭代計(jì)算次數(shù))，如圖3所示.如果r值設(shè)置過(guò)大，則目標(biāo)點(diǎn)云的樣本仍會(huì)較多；如果r值設(shè)置過(guò)小，則可能會(huì)造成目標(biāo)點(diǎn)云U中沒(méi)有任何樣本位于區(qū)域Si內(nèi)，從而造成關(guān)鍵樣本數(shù)據(jù)丟失，最終影響配準(zhǔn)精度.實(shí)際上，可將配準(zhǔn)誤差e作為球形抽樣區(qū)域半徑r的取值依據(jù)，即

r=k·e.

(10)

式中，k為可設(shè)定的比例系數(shù)(k> 0).

3) 源點(diǎn)云V的異常樣本剔除:將無(wú)目標(biāo)點(diǎn)云U樣本的球形區(qū)域所對(duì)應(yīng)的源點(diǎn)云樣本視為異常點(diǎn)，并予以剔除，得到更新后的源點(diǎn)云為V(r)，其樣本數(shù)為mr.若出現(xiàn)該問(wèn)題時(shí)，除r值設(shè)置過(guò)小的因素外，則很可能是在患者空間中采樣的操作失誤.

4) 最近點(diǎn)搜索配對(duì)與錯(cuò)誤剔除:對(duì)于源點(diǎn)云V(r)中的每個(gè)點(diǎn)，均在目標(biāo)點(diǎn)云U(r)中遍歷搜索與之歐式距離最小的點(diǎn)，以組成最近對(duì)應(yīng)點(diǎn)對(duì).為提升配準(zhǔn)精度，需刪除距離大于設(shè)定閾值d的異常對(duì)應(yīng)點(diǎn)對(duì)，并得到更新后的目標(biāo)點(diǎn)云P={pl}和源點(diǎn)云Q={ql}(l=1,2,…,N)，此時(shí)二者樣本數(shù)一致.

5) 建立目標(biāo)誤差函數(shù)E(TICP,h):

(11)

式中:h為迭代計(jì)算次數(shù)；N為最近對(duì)應(yīng)點(diǎn)對(duì)的數(shù)量；Rh為旋轉(zhuǎn)變換矩陣；Ph為平移變換向量.

6) 利用奇異值分析法[13]求解Rh和Ph，以使式(11)中E取得最小值.

7) 判斷迭代終止條件:如果e小于預(yù)設(shè)閾值ε或迭代次數(shù)h到達(dá)預(yù)設(shè)最大值hmax，則停止迭代，配準(zhǔn)完成；否則，重復(fù)步驟2)～6)，執(zhí)行第(h+1)次迭代，直至e<ε或h=hmax.

3 配準(zhǔn)實(shí)驗(yàn)案例

3.1 實(shí)驗(yàn)設(shè)置

以豬股骨和豬髂骨為實(shí)驗(yàn)對(duì)象進(jìn)行配準(zhǔn)實(shí)驗(yàn).本文涉及的配準(zhǔn)算法、數(shù)據(jù)操作及圖形可視化均基于Python與3D Slicer交互編程實(shí)現(xiàn)，計(jì)算機(jī)操作環(huán)境為Windows7 x64位系統(tǒng)、Intel(R) Core(TM) i3-6100(3.70 GHz)處理器和8 GB內(nèi)存(RAM).具體方法如下:

首先，將固定有定位剛體的股骨和髂骨經(jīng)CT掃描，并利用3D Slicer中重建三維模型；分別提取標(biāo)記球的三維虛擬模型和骨模型表面點(diǎn)云；獲得各標(biāo)記球心在光學(xué)位置傳感器和圖像空間中的位置；執(zhí)行初始配準(zhǔn)；利用探針工具分別在股骨和髂骨表面上的6個(gè)特定區(qū)域內(nèi)各采集4個(gè)點(diǎn)(如圖4所示)，獲得24個(gè)源點(diǎn)云樣本，圖像空間中的目標(biāo)點(diǎn)云(模型點(diǎn)云)樣本數(shù)分別為46 270和47 712個(gè)；分別利用本文提出的方法和傳統(tǒng)ICP算法執(zhí)行精確配準(zhǔn)；在患者空間中的每個(gè)采樣區(qū)內(nèi)，保證采樣數(shù)相同，僅改變采樣位置，共重復(fù)執(zhí)行30組實(shí)驗(yàn).其中，設(shè)置關(guān)鍵參數(shù):hmax=150，r=6 mm，d=3 mm.

3.2 點(diǎn)云配準(zhǔn)效果驗(yàn)證

以股骨配準(zhǔn)為例，通過(guò)預(yù)設(shè)實(shí)驗(yàn)數(shù)據(jù)的方式比較本文方法和傳統(tǒng)ICP算法的點(diǎn)云配準(zhǔn)效果.如圖5所示，預(yù)先給定一個(gè)確定的剛性變換矩陣Tdef，并將初始目標(biāo)點(diǎn)云A0變換為點(diǎn)云A1；在點(diǎn)云A1中提取子集源點(diǎn)云B0，以模擬通過(guò)探針工具在患者空間采集的點(diǎn)云數(shù)據(jù)；給定初始配準(zhǔn)矩陣T0，分別利用本文方法和傳統(tǒng)ICP算法對(duì)目標(biāo)點(diǎn)云A0和源點(diǎn)云B0進(jìn)行匹配，在本文方法中，經(jīng)初始配準(zhǔn)后(第一次迭代時(shí))的源點(diǎn)云為B1，此時(shí)點(diǎn)云A0經(jīng)抽樣后得到點(diǎn)云A0S，其樣本數(shù)減少至2 473個(gè)；利用最終的配準(zhǔn)矩陣將點(diǎn)云A1映射至初始目標(biāo)點(diǎn)云A0，并得到驗(yàn)證點(diǎn)云A10；觀察點(diǎn)云A0和點(diǎn)云A10的全局匹配效果，即表示目標(biāo)點(diǎn)云A0與源點(diǎn)云B0的實(shí)際配準(zhǔn)效果.

圖6為本文方法與傳統(tǒng)ICP算法的配準(zhǔn)效果.由圖6可知，當(dāng)采用本文方法時(shí)，配準(zhǔn)后的目標(biāo)點(diǎn)云與源點(diǎn)云的位置和姿態(tài)趨于一致，重合度較高；但利用傳統(tǒng)ICP算法時(shí)，兩組點(diǎn)云的整體匹配效果明顯較差，配準(zhǔn)結(jié)果陷入局部最優(yōu).

3.3 結(jié)果分析

利用配準(zhǔn)后兩組點(diǎn)云中各對(duì)應(yīng)點(diǎn)對(duì)的距離平均值作為配準(zhǔn)誤差度量，各組實(shí)驗(yàn)結(jié)果如圖7所示.當(dāng)采用本文方法對(duì)豬股骨和豬髂骨執(zhí)行配準(zhǔn)時(shí)，其最大配準(zhǔn)誤差分別為1.08mm和1.05mm，平均配準(zhǔn)誤差分別為(0.83±0.10) mm和(0.86±0.09) mm，由此可見(jiàn)，配準(zhǔn)誤差平均值及標(biāo)準(zhǔn)差均較小，各組結(jié)果分布較平穩(wěn)；此外，配準(zhǔn)運(yùn)算的耗時(shí)較短，其平均值分別僅為(0.034±0.003) s和(0.035±0.004) s.然而，當(dāng)利用傳統(tǒng)ICP算法進(jìn)行配準(zhǔn)時(shí)，其最大誤差分別達(dá)4.31 mm和4.18 mm，平均誤差分別為(3.35±0.42) mm和(3.28±0.44) mm，且運(yùn)算耗時(shí)均在0.8s左右，顯然，傳統(tǒng)ICP算法對(duì)應(yīng)的配準(zhǔn)誤差及各組結(jié)果波動(dòng)范圍均較大，且運(yùn)算耗時(shí)明顯更長(zhǎng).由此可見(jiàn)，相較于傳統(tǒng)ICP算法，本文所提方法在實(shí)現(xiàn)準(zhǔn)確、穩(wěn)定和高效的患者配準(zhǔn)應(yīng)用中更具適用性和可行性.

為進(jìn)一步說(shuō)明采用本文方法執(zhí)行患者配準(zhǔn)的結(jié)果穩(wěn)定且具有可重復(fù)性，分別在股骨表面前側(cè)、旁側(cè)和后側(cè)設(shè)置患者空間采樣區(qū)，保持其他條件設(shè)置相同，對(duì)每種采樣配置各執(zhí)行10組配準(zhǔn)實(shí)驗(yàn)，結(jié)果如圖8所示.股骨前側(cè)、旁側(cè)和后側(cè)采樣配置對(duì)應(yīng)的平均配準(zhǔn)誤差分別為(0.85±0.09) mm,(0.89±0.09) mm和(0.81±0.10) mm，可見(jiàn)，誤差平均值均小于0.90 mm，標(biāo)準(zhǔn)差相近均未超過(guò)0.10 mm.由實(shí)驗(yàn)數(shù)據(jù)可知，在患者空間中不同的采樣配置對(duì)使用本文方法所得的配準(zhǔn)結(jié)果影響較小，從而能在一定程度上表明本文方法具有較好的穩(wěn)定性和可重復(fù)性.

實(shí)際上，配準(zhǔn)精度會(huì)受到多方面因素的影響，例如球形抽樣半徑r、最近點(diǎn)對(duì)的距離閾值d、迭代次數(shù)h和目標(biāo)點(diǎn)云樣本密度等參數(shù)，以及源點(diǎn)云采樣誤差、三維模型重建誤差和計(jì)算誤差等因素.

本文通過(guò)控制變量試驗(yàn)，僅考察r值對(duì)配準(zhǔn)精度的影響.保持股骨配準(zhǔn)實(shí)驗(yàn)的其余設(shè)置不變，僅改變r(jià)值，得到其與配準(zhǔn)誤差e的關(guān)系，如圖9所示.可見(jiàn)，在不同初始配準(zhǔn)誤差eInit下的e隨r的變化趨勢(shì)基本相同.以eInit=1.5 mm為例進(jìn)行詳細(xì)說(shuō)明，當(dāng)r<2 mm時(shí)，e大致呈遞減趨勢(shì)，此時(shí)目標(biāo)點(diǎn)云中所保留的樣本數(shù)較少，e更易受到eInit的影響；當(dāng)2 mmeIint.同時(shí)，r值上限還受到配準(zhǔn)運(yùn)算耗時(shí)和預(yù)期配準(zhǔn)精度的約束.若要求縮短配算耗時(shí)，則r值應(yīng)盡量減??；在算法實(shí)現(xiàn)中，由于運(yùn)算耗時(shí)很短，則可令r值按照較小增量Δr變化并遍歷計(jì)算出各自所對(duì)應(yīng)的實(shí)際配準(zhǔn)誤差，最后再?gòu)闹泻Y選出滿足預(yù)設(shè)配準(zhǔn)精度的r值范圍.

4 結(jié) 論

1) 采用本文所提方法進(jìn)行患者配準(zhǔn)具有比傳統(tǒng)ICP算法更佳的配準(zhǔn)精度和穩(wěn)定性.

2) 合理抽樣簡(jiǎn)化圖像空間點(diǎn)云利于降低點(diǎn)云噪聲對(duì)配準(zhǔn)精度產(chǎn)生的不利影響和縮短運(yùn)算耗時(shí)，因此應(yīng)根據(jù)初始配準(zhǔn)誤差和運(yùn)算耗時(shí)等約束條件合理設(shè)定球形抽樣區(qū)域的半徑r.

3) 通過(guò)本文方法解決患者配準(zhǔn)問(wèn)題的可操作性強(qiáng)，臨床應(yīng)用潛力大，既無(wú)需在圖像空間和患者空間執(zhí)行規(guī)模大且繁瑣的采樣，也利于在不侵入患者體內(nèi)植入標(biāo)記物的前提下保證較高的配準(zhǔn)精度.