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

(東北大學 機械工程與自動化學院，遼寧 沈陽 110819)

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

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

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

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

1 患者配準問題描述及數(shù)據(jù)獲取

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

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

(1)

(2)

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

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

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

2 配準方法

2.1 基于三點法的初始配準

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

(3)

(4)

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

(5)

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

eInit=|fI-TInit·fP|.

(6)

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

2.2 基于ICP算法的精確配準

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

(7)

(8)

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

V=TInit·V′.

(9)

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

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

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

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

r=k·e.

(10)

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

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

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

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

(11)

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

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

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

3 配準實驗案例

3.1 實驗設(shè)置

以豬股骨和豬髂骨為實驗對象進行配準實驗.本文涉及的配準算法、數(shù)據(jù)操作及圖形可視化均基于Python與3D Slicer交互編程實現(xiàn)，計算機操作環(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中重建三維模型；分別提取標記球的三維虛擬模型和骨模型表面點云；獲得各標記球心在光學位置傳感器和圖像空間中的位置；執(zhí)行初始配準；利用探針工具分別在股骨和髂骨表面上的6個特定區(qū)域內(nèi)各采集4個點(如圖4所示)，獲得24個源點云樣本，圖像空間中的目標點云(模型點云)樣本數(shù)分別為46 270和47 712個；分別利用本文提出的方法和傳統(tǒng)ICP算法執(zhí)行精確配準；在患者空間中的每個采樣區(qū)內(nèi)，保證采樣數(shù)相同，僅改變采樣位置，共重復執(zhí)行30組實驗.其中，設(shè)置關(guān)鍵參數(shù):hmax=150，r=6 mm，d=3 mm.

3.2 點云配準效果驗證

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

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

3.3 結(jié)果分析

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

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

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

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

4 結(jié) 論

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

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

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