胡 巖 王田苗 王君臣 胡 磊

(北京航空航天大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院,北京 100191)

基于 Kriging算法的手術(shù)導(dǎo)航三維形變技術(shù)

胡 巖 王田苗 王君臣 胡 磊

(北京航空航天大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院,北京 100191)

針對(duì)計(jì)算機(jī)輔助手術(shù)三維導(dǎo)航技術(shù)中無(wú)法利用術(shù)前 CT重建模型進(jìn)行病人個(gè)體化手術(shù)的問(wèn)題,提出一種基于克里金(Kriging)算法的三維骨骼形變(bone morphing)技術(shù).利用術(shù)中采集和醫(yī)學(xué)圖像模型的兩片點(diǎn)云進(jìn)行配準(zhǔn),在 CT重建模型上選取 4個(gè)合適的控制點(diǎn),通過(guò)求解 Kriging算法建立的方程組,獲得 CT重建模型上所有點(diǎn)相對(duì)于實(shí)際病例變形后的位置,實(shí)現(xiàn)醫(yī)學(xué)圖像模型向個(gè)體化的變形.最后,以人體股骨 CT模型和塑料模型為對(duì)象進(jìn)行了形變的精度分析和關(guān)節(jié)手術(shù)導(dǎo)航評(píng)估實(shí)驗(yàn).結(jié)果表明:這種形變方法簡(jiǎn)單可靠,可以用于前交叉韌帶(ACL,Anterior Cruciate Ligament)重建手術(shù)的導(dǎo)航評(píng)估.

膝關(guān)節(jié);手術(shù)導(dǎo)航;克里金;配準(zhǔn);骨骼形變

計(jì)算機(jī)輔助外科手術(shù)三維導(dǎo)航技術(shù)目的在于實(shí)現(xiàn)外科手術(shù)過(guò)程可視化,減少手術(shù)創(chuàng)傷和輻射,同時(shí)提高手術(shù)精度[1].利用機(jī)器人和導(dǎo)航技術(shù),將醫(yī)學(xué)圖像和手術(shù)過(guò)程中的實(shí)際解剖組織關(guān)聯(lián)起來(lái),提供給醫(yī)生一個(gè)可視化、交互式、精確的手術(shù)環(huán)境.臨床上通過(guò)常用的醫(yī)學(xué)圖像重建出手術(shù)對(duì)象(骨骼、軟組織、血管等)的三維形狀[2].醫(yī)學(xué)圖像形變技術(shù),就是將 CT重建模型和實(shí)際病例對(duì)象聯(lián)系起來(lái),實(shí)現(xiàn)針對(duì)個(gè)體病例的變形,這樣才能將實(shí)際病例的操作映射到 CT模型,引導(dǎo)醫(yī)生完成手術(shù).因此形變技術(shù)是手術(shù)導(dǎo)航技術(shù)的關(guān)鍵,形變的精度直接影響到手術(shù)的成敗.

醫(yī)學(xué)圖像的形變分為二維和三維圖像的形變,前者包括多邊形形狀漸變和圖像自然漸變,后者包括剛體和軟組織的形變.本文只針對(duì)三維剛體的形變.圖像變形技術(shù)是研究圖像中點(diǎn)到點(diǎn)的空間映射和變換問(wèn)題,該技術(shù)可以應(yīng)用于眾多手術(shù)適應(yīng)癥(全膝關(guān)節(jié)置換、前交叉韌帶重建)的無(wú)圖像導(dǎo)航.文獻(xiàn)[3-4]提出了基于模型的技術(shù),通過(guò)骨骼形變(bone morphing)技術(shù)獲得的骨骼表面點(diǎn)推理出病人個(gè)體化的三維表面,利用基于模型的統(tǒng)計(jì)主元分析法(PCA,Principal Component Analysis)實(shí)現(xiàn)通用模型到術(shù)中數(shù)據(jù)的形變和姿態(tài)的優(yōu)化.文獻(xiàn)[5]運(yùn)用了最近點(diǎn)迭代(ICP,Iterative Closed Point)法進(jìn)行形變和姿態(tài)的優(yōu)化.文獻(xiàn)[6-7]利用 PCA獲得的模型,根據(jù)病人的身高、體重和年齡,推理出了所有可能的模型.法國(guó)Praxim公司在文獻(xiàn)[5]研究的基礎(chǔ)上推出了相應(yīng)的產(chǎn)品 Surgetics系統(tǒng)[8-9].國(guó)外的方法一般都需要進(jìn)行統(tǒng)計(jì)學(xué)分析,算法的實(shí)現(xiàn)周期比較長(zhǎng).而本文提出的方法簡(jiǎn)單實(shí)用,無(wú)需對(duì)病體進(jìn)行統(tǒng)計(jì)分析.該技術(shù)可以應(yīng)用到骨科的多個(gè)領(lǐng)域,應(yīng)用時(shí)無(wú)需任何 X片和(或)CT等影像,避免 X射線輻射,減少費(fèi)用支出,并且可以針對(duì)個(gè)體化的病例進(jìn)行手術(shù),將會(huì)提高手術(shù)精度.

1 三維模型形變的數(shù)學(xué)描述

骨骼形變技術(shù)是將給定的原數(shù)字化的幾何模型 S變換到目標(biāo)數(shù)字幾何模型 T,T是由 CT重建的三維模型,S是術(shù)中在光學(xué)定位儀下采集的數(shù)字化病例.所謂 S到T的配準(zhǔn),就是找到一種變形函數(shù),使得 S中的點(diǎn)變換到 T的相應(yīng)元素上:

用一種通用的流程來(lái)描述所有的醫(yī)學(xué)模型變形問(wèn)題:

1)分別抽取 S和 T中的特征元素集合,記為FS和 FT,(ps)j和(pt)j分別為特征集合中對(duì)應(yīng)點(diǎn)的元素;

2)定義目標(biāo)函數(shù):

式中,D為廣義距離函數(shù);N對(duì)應(yīng)變換點(diǎn)的數(shù)量;a和 b為待求參數(shù).

3)加入邊界條件并使用優(yōu)化算法求解目標(biāo)函數(shù)的全局極小點(diǎn),得到參數(shù) a和 b.

在手術(shù)三維導(dǎo)航中,一般需要將術(shù)前 CT/MRI三維重建模型與手術(shù)空間中采集的數(shù)字化特征點(diǎn)進(jìn)行對(duì)應(yīng)變形.但是,變形前還必須進(jìn)行配準(zhǔn),因此這里的 S和T都是笛卡爾坐標(biāo)系,而廣義距離就是三維坐標(biāo)點(diǎn)之間的歐氏距離.

2 基于 Kriging算法的模型形變

Kriging算法在揭示空間數(shù)據(jù)非均勻、各向異性分布規(guī)律方面具有優(yōu)勢(shì),因而用于空間散亂點(diǎn)數(shù)據(jù)的插值.在膝關(guān)節(jié)的三維模型中,設(shè)第 j個(gè)控制點(diǎn)變形前后的坐標(biāo)為

(Ps)j即手術(shù)中測(cè)的坐標(biāo)值.建立方程組

式中

其中 c為常數(shù).

式(2)右邊可以看作仿射變換和彈性變形的組合,彈性變形部分需要邊界條件:

因此,參數(shù) a和 b可通過(guò)以下方程組解:

為保證式(2)有唯一解,控制點(diǎn)個(gè)數(shù)應(yīng)不少于 4,采用股骨的生理特征點(diǎn)作為控制點(diǎn).手術(shù)中,在股骨髁上取 4個(gè)股骨表面生理標(biāo)志點(diǎn),如圖1所示.

圖1 圖像坐標(biāo)系與股骨的生理標(biāo)志點(diǎn)

由于這 4個(gè)點(diǎn)的實(shí)際位置均由醫(yī)生的經(jīng)驗(yàn)確定,選取時(shí)要盡量對(duì)應(yīng),與模型點(diǎn)的對(duì)應(yīng)關(guān)系必然存在誤差.本文采用了文獻(xiàn)[10]中的改進(jìn)方法:

式中,I是單位矩陣;λ是歸一化參數(shù),調(diào)節(jié)了仿射變換和彈性變形的權(quán)重,λ越大,彈性變形所占的比重越小.在 z方向上,參數(shù)求解方程組可以寫(xiě)成如下的矩陣形式:

求出參數(shù) a和 b后,再將它們代回到式(9),設(shè) M為模型上點(diǎn)的總個(gè)數(shù),任意點(diǎn)變形后的 z坐標(biāo)都可以式(9)表示,x和 y坐標(biāo)可以用同樣的方法求解.

3 基于表面特征點(diǎn)的形變方法

手術(shù)中需要抽取圖像和實(shí)際病例的對(duì)應(yīng)點(diǎn)進(jìn)行配準(zhǔn),而后通過(guò)選取圖像上的特征點(diǎn),結(jié)合上節(jié)的算法就可以完成變形,變形過(guò)程如圖 2所示.

圖2 變形過(guò)程示意圖

3.1 數(shù)字化 CT模型重建

術(shù)前 CT圖像是具有幾何坐標(biāo)的醫(yī)學(xué)圖像格式,是對(duì)連續(xù)空間進(jìn)行采樣得到的三維離散灰度場(chǎng).利用 Marchingcubes算法[11]重建出關(guān)節(jié)的三維模型,它是面繪制模型,其幾何描述是基于掃描CT數(shù)據(jù)時(shí)定義的坐標(biāo)系.采用可視化工具包(VTK,Visualization Toolkit)對(duì)三維模型進(jìn)行交互式操作,例如模型分割、三角網(wǎng)格優(yōu)化、光線投射拾取單個(gè)表面點(diǎn)、大面積拾取表面點(diǎn)等.醫(yī)學(xué)圖像上的特征點(diǎn)云可以通過(guò)先將圖像數(shù)據(jù)重建成面繪制模型,然后利用光線投射法交互式地拾取重建模型上特征區(qū)域表面上的三維點(diǎn)集,稱(chēng)其為點(diǎn)云T,如圖 3a所示.這一過(guò)程稱(chēng)為數(shù)字化圖像空間特征曲面.

圖3 特征點(diǎn)云

3.2 數(shù)字化個(gè)體模型獲取

術(shù)中采用加拿大 NDI公司的光學(xué)定位儀 Polaris Spectrum,標(biāo)稱(chēng)定位精度小于 0.25mm.在手術(shù)對(duì)象上安裝可被光學(xué)定位儀識(shí)別的 Marker,它是術(shù)中實(shí)際空間坐標(biāo)系.利用探針在手術(shù)對(duì)象(股骨髁和脛骨平臺(tái))表面劃動(dòng),由定位儀實(shí)時(shí)記錄探針末端點(diǎn)在 Marker坐標(biāo)系下的坐標(biāo),稱(chēng)其為點(diǎn)云 S,如圖 3b所示.這一過(guò)程稱(chēng)為數(shù)字化實(shí)際空間特征曲面.

3.3 三維配準(zhǔn)及形變

運(yùn)用 ICP算法將 CT重建模型和術(shù)中采集病例信息進(jìn)行三維配準(zhǔn).需要注意的是特征區(qū)域不必要連續(xù),每片點(diǎn)云都可以是對(duì)空間中若干個(gè)不同區(qū)域數(shù)字化的結(jié)果.然后,運(yùn)用第 2節(jié)的算法對(duì)CT模型進(jìn)行變形.三維形變的流程如圖 4所示.

圖4 手術(shù)三維導(dǎo)航的變形流程

4 形變精度及導(dǎo)航應(yīng)用實(shí)驗(yàn)

4.1 骨骼形變的精度

本文采用自主研發(fā)的光電導(dǎo)航系統(tǒng)為平臺(tái)進(jìn)行仿真評(píng)估實(shí)驗(yàn),以塑料關(guān)節(jié)模型骨為實(shí)驗(yàn)對(duì)象,如圖 5所示.實(shí)驗(yàn)前對(duì)志愿者下肢關(guān)節(jié)進(jìn)行 CT掃描,層間距為 0.5 mm,x和 y方向間距為0.586mm.實(shí)驗(yàn)中使用灰度值 300對(duì)關(guān)節(jié) CT進(jìn)行三維重建.

圖5 實(shí)驗(yàn)平臺(tái)、關(guān)節(jié)模型骨及探針

前交叉韌帶(ACL,Anterior Cruciate Ligament)的植入點(diǎn)在股骨髁間窩和脛骨平臺(tái),因此這兩個(gè)部位的變形最重要,選取其作為導(dǎo)航測(cè)試區(qū)域.利用 OpenGL工具包建立三維測(cè)試環(huán)境,進(jìn)行ACL重建的仿真評(píng)估.采用 NDI定位儀標(biāo)定過(guò)的探針模擬手術(shù)器械,其尖點(diǎn)在 Marker坐標(biāo)系下的空間坐標(biāo)可實(shí)時(shí)得到,利用該探針在塑料骨的股骨髁和脛骨平臺(tái)表面滑動(dòng),采集關(guān)節(jié)骨骼表面特征.圖 6所示為變形前后的股骨模型.

圖6 變形前后的CT模型和術(shù)中點(diǎn)云

骨骼形變的精度取決于計(jì)算機(jī)的自動(dòng)控制和人為交互,求解變形后的 CT重建模型和醫(yī)生進(jìn)行數(shù)字化的骨面特征之間的均方根誤差.如果它大于 1mm,則此次形變是不合格的,需要重新采集點(diǎn)云.本文進(jìn)行了 300次實(shí)驗(yàn),隨機(jī)抽取其中的20次數(shù)據(jù)分析,如表 1所示,其中僅有 2次測(cè)得的誤差超過(guò) 1mm.剔除這 2次結(jié)果重新實(shí)驗(yàn),得到的結(jié)果是 0.803mm,證明超過(guò) 1mm的誤差是醫(yī)生的操作引起的,與骨骼形變算法無(wú)關(guān).由于骨骼形變精度主要取決于覆蓋骨骼表面的點(diǎn)云,點(diǎn)云數(shù)量需要至少 300個(gè)點(diǎn),且均勻分布于骨骼表面,否則形變誤差較大.將探針?lè)旁诠趋辣砻娌煌恢?計(jì)算機(jī)中顯示探針針尖,計(jì)算尖點(diǎn)和變形后點(diǎn)的歐幾里德距離,股骨最大為 0.838mm,脛骨最大為 0.916mm.產(chǎn)生誤差的原因包括定位器誤差、探針誤差、剛體的位移、配準(zhǔn)和重建的數(shù)學(xué)誤差、Kriging算法的插值誤差等.

表 1 骨骼形變的精度

4.2 導(dǎo)航應(yīng)用實(shí)驗(yàn)

驗(yàn)證韌帶重建手術(shù)是否成功的重要標(biāo)準(zhǔn)[12]主要是兩點(diǎn).利用該文提出的形變技術(shù)進(jìn)行了針對(duì) ACL重建的術(shù)后仿真和評(píng)估.圖 4所示為 CT模型的形變流程,將變形后的三維模型保存成STL文件,讀入自行開(kāi)發(fā)的“無(wú)圖像導(dǎo)航下的 ACL重建仿真評(píng)估系統(tǒng)”.利用探針在塑料模型骨上點(diǎn)取植入點(diǎn)位置,重建虛擬 ACL,脛骨繞股骨按人體運(yùn)動(dòng)方式旋轉(zhuǎn),觀察 ACL是否與骨骼面碰撞以及不等距值的變化,如圖 7所示.不等距值不發(fā)生變化,表示重建的 ACL符合不等距性;反之則不符合.從圖 7中可以發(fā)現(xiàn)在脛骨運(yùn)動(dòng)時(shí),雖然ACL有不同程度的延長(zhǎng)或收縮,但重建 ACL的不等距值保持不變,為 2.360 039mm,證明重建的ACL符合不等距性,變形后的 CT模型可以用于ACL重建的導(dǎo)航評(píng)估.

圖7 ACL術(shù)后導(dǎo)航評(píng)估

5 結(jié) 論

本文提出了一種基于 Kriging算法的三維骨骼形變技術(shù),將該技術(shù)用于 ACL重建的導(dǎo)航評(píng)估手術(shù),取得了良好的效果,術(shù)中形變經(jīng)過(guò)多次實(shí)驗(yàn)均在 1min以?xún)?nèi).在本文的實(shí)驗(yàn)條件下形變后精度在1mm以?xún)?nèi),證明該方法可以用于膝關(guān)節(jié)ACL重建術(shù)后導(dǎo)航評(píng)估,并且流程簡(jiǎn)單,不需要復(fù)雜的圖像分割,具有一定的通用性.該技術(shù)能夠提高手術(shù)質(zhì)量,減輕醫(yī)生手術(shù)的疲勞程度,降低醫(yī)患輻射傷害和手術(shù)成本.但是,該技術(shù)只是初步用于模型骨的導(dǎo)航評(píng)估功能實(shí)驗(yàn),其精度還有待于進(jìn)一步的提高,并應(yīng)用于 ACL的重建手術(shù)規(guī)劃中,例如可以幫助醫(yī)生進(jìn)行植入點(diǎn)的規(guī)劃、顯示植入隧道等.總之,骨骼技術(shù)對(duì)于實(shí)時(shí)三維顯示骨骼缺陷、手術(shù)鉆孔規(guī)劃、移植物的尺寸選擇、根據(jù)骨骼冠狀面的移植物位姿的選擇等將會(huì)有重要且深遠(yuǎn)的意義.

References)

[1]Nolte L P,Thomas B.Basic principles of CAOS[J].Care Injured,2004,35(1):6-16

[2]王君臣,王田苗,徐源,等.基于 ICP算法的手術(shù)導(dǎo)航三維配準(zhǔn)技術(shù)[J].北京航空航天大學(xué)學(xué)報(bào),2009,35(4):434-438 Wang Junchen,Wang Tianmiao,Xu Yuan,et al.Registration method based on ICP algorithm for 3D surgical navigation[J].Journal of Beijing University of Aeronautics and Astronautics,2009,35(4):434-438(in Chinese)

[3]Fleute M,Lavallee.Building a complete surface model from sparse data using statistical shape models[J].MICCAI,1998:879-887

[4]Stindel E,Briard J L,Merloz P,et al.Bone morphing:3D morphological data for tool knee arthroplasty[J].Computer Aided Surgery,2002:156-168

[5]Chan C S,Edwards P J,Hawkes D J.Integration of ultrasound based registration with statistical shape models for computer-assisted orthopaedic surgery[J].SPIE,Medical Imaging,2003:414-424

[6]Kumar T R,Johannes H,Nolte L P,et al.Bone morphing with statistical shape models for enhanced visualization[J].Medical Imaging,2004,5367:122-130

[7]Kumar T R,Martin A S,Haydar T,et al.Statistical deformable models for robust 3D surface extrapolation from sparse data[J].Medical Image Analysis,2007,11(2):99-109

[8]Plaweski S,Pearle A,Granchi C,et al.PRAX IM ACL navigation system using bonemorphing:Navigation and MIS in orthopedic surgery[M].Berlin:Springer,2007:315-323

[9]Pena E,Calvo B,Martinez M A,et al.Influence of the tunnel angle in ACL reconstructionson the biomechanics of the knee joint[J].Clin Biomech,2006,21(5):508-516

[10]Mǜller J,Mang A,Buzug T M.A template-deformation method for facial reproduction[C]//Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis.Croatia:Institute of Electrical＆Electronics Engine,2005:359-364

[11]Timothy SN,Hong Y.A survey of the marching cubes algorithm[J].Computer＆Graphics,2006,30:854-879

[12]Markus F,Stéphane L,Rém i J.Incorporating a statistically based shape model into a system for computer-assisted anterior cruciate ligament surgery[J].Medical Image Analysis,1999,3(3):209-222

(編 輯:文麗芳)

3Dmorphing method based on Kriging algorithm for surgical navigation

Hu Yan Wang Tianmiao Wang Junchen Hu Lei

(School of Mechanical Engineering and Automation,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)

A setof preoperative CT reconstructed model can not beused in operating on individual patient in 3D surgical navigation systems.According to the problem above,a 3D bone morphing method based on Kriging algorithm was presented.Two point sets,one of which was acquired from physical space by the optical localizer and the other was from 3Dmodel reconstructed by the marching cube method,were used to calculate the transformation matrix between the preoperative CT image space and the intra operative marker space.Then 4 appropriate reference points were chosen to solve the equations based on Kriging algorithm,and the positions of points on CTmodel after deformation respect to the physical case were acquired.At last,precision analysis of morphing and evaluation experiment of knee jointnavigation was carried out with CT and plastic model.The result shows that under this condition the final morphing accuracy is less than 1mm,and this technology can be used in operation evaluation of the anterior cruciate ligament navigation reconstruction.

knee joint;surgical navigation;Kriging;registration;bone morphing

TP 391.9

A

1001-5965(2010)05-0555-05

2009-04-07

國(guó)家杰出青年科學(xué)基金資助項(xiàng)目(60525314);國(guó)家科技支撐計(jì)劃資助項(xiàng)目(2006BAI03A 16);北京市科技計(jì)劃資助項(xiàng)目(D 020602500001093)

胡 巖(1979-),男,山東濟(jì)南人,博士生,buaarobot@126.com.