腦組織牽拉變形矯正的XFEM 建模及其驗(yàn)證

王巍偉 李 萍 章琛曦 宋志堅(jiān)*1(復(fù)旦大學(xué)數(shù)字醫(yī)學(xué)研究中心，上海 200032)

2(上海市醫(yī)學(xué)圖像處理與計(jì)算機(jī)輔助手術(shù)重點(diǎn)實(shí)驗(yàn)室，上海 200032)

引言

以CT、MRI 等醫(yī)學(xué)圖像為基礎(chǔ)的神經(jīng)外科手術(shù)導(dǎo)航系統(tǒng)(image guided navigation system，IGNS)，對(duì)于輔助醫(yī)生準(zhǔn)確識(shí)別與切除腦部病變組織、避免損傷重要神經(jīng)或功能區(qū)等具有重要的臨床意義［1］。但是，在神經(jīng)外科手術(shù)中，由于開顱后顱內(nèi)壓的改變、腦脊液的流出、牽拉、腫瘤切除等因素的影響，腦組織會(huì)產(chǎn)生較為嚴(yán)重的變形。這導(dǎo)致由術(shù)前數(shù)據(jù)建立的圖像引導(dǎo)系統(tǒng)與實(shí)際術(shù)中腦組織情況產(chǎn)生較大偏差，大幅降低了IGNS 的定位精度。腦組織變形的產(chǎn)生可分為3 種類型:第1 種是硬腦膜打開后產(chǎn)生，第2 種是手術(shù)牽拉產(chǎn)生，第3 種是病灶切除產(chǎn)生。牽拉變形是外科醫(yī)生為尋找到達(dá)大腦皮層下目標(biāo)(如腫瘤)的路徑，用腦壓板撥開腦組織導(dǎo)致的腦組織變形［2］。目前，國(guó)內(nèi)外的研究主要集中在第1 種腦組織變形上，較少考慮牽拉和腫瘤切除引起的第2、3 種腦組織變形，而這兩種腦組織變形正是降低IGNS 術(shù)中導(dǎo)航精度的主要因素。矯正手術(shù)牽拉和切除引起腦組織變形，主要在于都會(huì)導(dǎo)致腦組織拓?fù)浣Y(jié)構(gòu)改變，且矯正牽拉變形是矯正切除變形的基礎(chǔ)，因此矯正腦組織牽拉變形顯得尤為重要。

目前主要有兩類矯正腦組織變形的方法:術(shù)中影像和生物力學(xué)模型矯正。雖然術(shù)中影像如術(shù)中CT(intraoperative CT，iCT)和術(shù)中MR(intraoperative MR，iMR)能夠較準(zhǔn)確地矯正腦組織變形，但其缺點(diǎn)是成本高、操作不方便、容易引起感染且延長(zhǎng)了手術(shù)時(shí)間，因此限制了術(shù)中影像矯正的臨床應(yīng)用。生物力學(xué)模型的矯正方法通過(guò)建立生物力學(xué)模型，然后利用有限元方法(finite element method，F(xiàn)EM)的求解來(lái)矯正腦組織變形，從而克服了術(shù)中影像方法的缺點(diǎn)。Ferrant 使用線彈性生物力學(xué)模型建模，通過(guò)刪除牽拉路徑上所有網(wǎng)格，實(shí)現(xiàn)模型的更新［3］。Miga 應(yīng)用固結(jié)理論，建立線性多孔彈性生物力學(xué)模型，沿牽拉路徑切割網(wǎng)格，產(chǎn)生兩倍的網(wǎng)格節(jié)點(diǎn)，實(shí)現(xiàn)模型的更新［4］。動(dòng)物實(shí)驗(yàn)表明，該模型可矯正平均66%的腦組織位移誤差，降低80%的定位誤差［5］。Sun 同樣使用線性多孔彈性生物力學(xué)模型，通過(guò)FEM 計(jì)算，模擬牽拉腦組織引起的腦組織變形［6］。該方法利用立體攝像機(jī)獲取腦壓板的運(yùn)動(dòng)位移場(chǎng)，并以此為邊界條件對(duì)模型進(jìn)行求解。此方法獲取大約75%的腦皮層運(yùn)動(dòng)。

以上研究者都利用FEM 來(lái)處理腦組織拓?fù)浣Y(jié)構(gòu)的改變，概括起來(lái)有兩種處理方式:完全重新網(wǎng)格化或者將拓?fù)浣Y(jié)構(gòu)改變區(qū)域局部網(wǎng)格化［7］。不管FEM 采用哪種處理方法，都需要改變網(wǎng)格原有的拓?fù)浣Y(jié)構(gòu)，會(huì)破壞原有網(wǎng)格的連續(xù)性，而擴(kuò)展有限元(extended finite element method，XFEM)［8］則克服了這個(gè)缺點(diǎn)。它處理腦組織變形首先不考慮拓?fù)浣Y(jié)構(gòu)變化產(chǎn)生的不連續(xù)性，而是直接將腦組織網(wǎng)格化，再根據(jù)裂縫的具體形狀，沿著裂縫經(jīng)過(guò)的網(wǎng)格節(jié)點(diǎn)加上新的虛自由度。這樣處理腦組織牽拉產(chǎn)生的裂縫，能在不破壞原有網(wǎng)格連續(xù)性的基礎(chǔ)上準(zhǔn)確地表達(dá)裂縫產(chǎn)生的拓?fù)洳贿B續(xù)性。Vigneron 首次提出將XFEM 應(yīng)用于線彈性生物力學(xué)模型中，并利用基于iMR 圖像的非剛性配準(zhǔn)技術(shù)，獲取牽拉前后腦組織的位移，以此作為邊界條件求解模型［2，9-10］。此方法使用Canny 算子對(duì)配準(zhǔn)前后的圖像進(jìn)行邊緣檢測(cè)，結(jié)果表明配準(zhǔn)后圖像的Hausdorff 距離減小。但該方法僅得到定性結(jié)果，沒(méi)有給出XFEM 在矯正牽拉時(shí)引起腦組織變形有效性的定量性結(jié)果。

本研究提出，將XFEM 應(yīng)用于線彈性生物力學(xué)模型，用來(lái)矯正牽拉變形的系統(tǒng)框架。選取腦組織體模作為試驗(yàn)對(duì)象，用iCT 圖像作為評(píng)估的金標(biāo)準(zhǔn);利用植入的不銹鋼標(biāo)記物，定量分析了XFEM 在矯正腦組織牽拉變形的有效性。同時(shí)，歸納比較不同尺寸網(wǎng)格對(duì)于XFEM 矯正牽拉變形的影響。

1 牽拉變形矯正系統(tǒng)框架

基于XFEM 矯正腦組織牽拉變形的系統(tǒng)框架如圖1 所示。步驟1 和步驟2 為打開硬腦膜后的腦變形狀態(tài)(以下稱為術(shù)前狀態(tài))，得到硬膜打開后已矯正腦組織變形的腦組織網(wǎng)格;步驟3 通過(guò)iCT/iMR 掃描得到手術(shù)中牽拉產(chǎn)生的腦組織變形圖像;步驟4 將牽拉變形圖像通過(guò)點(diǎn)配準(zhǔn)與術(shù)前圖像配準(zhǔn)，以確定在術(shù)前圖像坐標(biāo)系中腦壓板初始位置和腦壓板牽拉位移;步驟5 將獲得的腦壓板初始位置和腦壓板牽拉位移引入到基于術(shù)前圖像的XFEM線彈性模型中，作為邊界條件;步驟6 根據(jù)步驟5 求得的邊界條件和已知的零位移邊界條件求解XFEM線彈性模型方程，據(jù)此，得到牽拉變形網(wǎng)格;步驟7更新腦組織術(shù)前圖像，得到矯正腦組織牽拉變形后的圖像;步驟8 將更新后的圖像顯示給用戶。

1.1 XFEM 原理

FEM 無(wú)法很好地處理拓?fù)浣Y(jié)構(gòu)改變的不連續(xù)情況，這是由于決定網(wǎng)格節(jié)點(diǎn)生物力學(xué)屬性的形函數(shù)(nodal shape function，NSF)是連續(xù)函數(shù)。XFEM擴(kuò)展了FEM 的功能，通過(guò)增加裂縫處網(wǎng)格節(jié)點(diǎn)的虛自由度來(lái)更好地模擬裂縫等不連續(xù)情況，克服了FEM 處理裂縫時(shí)需要重新進(jìn)行網(wǎng)格化的缺點(diǎn)。XFEM 可表達(dá)為

式中，uXFEM為XFEM 位移場(chǎng)。

式(1)等號(hào)右側(cè)第1 部分是FEM 網(wǎng)格節(jié)點(diǎn)位移場(chǎng)，I 表示FEM 網(wǎng)格節(jié)點(diǎn)集合，φi是FEM 的NSF，ui是網(wǎng)格節(jié)點(diǎn)FEM 自由度。與FEM 不同，XFEM 在網(wǎng)格點(diǎn)集J 和K 兩個(gè)集合上新增加了虛自由度，用這些虛自由度來(lái)表述裂縫的不連續(xù)性。點(diǎn)集J 和K的關(guān)系可表達(dá)為

圖1 應(yīng)用XFEM 建模矯正牽拉變形的系統(tǒng)框架Fig．1 XFEM framework for retraction correction

如圖2 空心圓點(diǎn)所示，當(dāng)裂縫將網(wǎng)格單元的所有連線都完全分割時(shí)，此類受影響網(wǎng)格節(jié)點(diǎn)組成點(diǎn)集J 。式(1)等號(hào)右側(cè)第2 部分新增加Heaviside 函數(shù)來(lái)表達(dá)裂縫對(duì)點(diǎn)集J 產(chǎn)生的不連續(xù)性。單位階躍函數(shù)H 用公式表示為

式中，x 為網(wǎng)格節(jié)點(diǎn)，x*為裂縫上與距離x 點(diǎn)最近點(diǎn)，en為x*處裂縫法向量，aj表示點(diǎn)集J 的虛自由度。

如圖2 實(shí)心方點(diǎn)所示，當(dāng)裂縫末端進(jìn)入某網(wǎng)格單元但沒(méi)有將該單元節(jié)點(diǎn)所有連線完全分割時(shí)，此類受影響網(wǎng)格節(jié)點(diǎn)組成受裂尖影響的點(diǎn)集K 。

式(1)等號(hào)右側(cè)第3 部分表達(dá)了點(diǎn)集K 的不連續(xù)性，包括半徑和角度的計(jì)算，Clk 為點(diǎn)集K 的新增虛自由度，F(xiàn)l為點(diǎn)集K 的新增不連續(xù)函數(shù)，可表示為

式中，r 和θ 分別為點(diǎn)信K 中某點(diǎn)x 與裂尖形成的法向量和切向量的夾角。

1.2 腦組織分割與網(wǎng)格化

使用圖像分割算法將術(shù)前采集圖像中的顱骨和皮膚等組織分割出來(lái)，僅保留包含腦組織部分的圖像，并完成對(duì)腦組織的三維重建。在此基礎(chǔ)上，利用類八叉樹算法將腦組織的三維幾何模型劃分為僅由均勻的六面體單元組成的網(wǎng)格。

圖2 XFEM 原理示意Fig．2 The schematic diagram of XFEM

1.3 XFEM 應(yīng)用于線性彈性模型

線彈性模型是應(yīng)用在腦組織變形矯正中最簡(jiǎn)單有效的生物力學(xué)模型［11］，它將腦組織看作是一個(gè)各向同性的線性彈性物體，用一系列準(zhǔn)靜態(tài)偏微分方程［12］表示。牽拉變形矯正的核心問(wèn)題就是使用XFEM 方法求解線性方程，即

式中，K 是網(wǎng)格節(jié)點(diǎn)XFEM 剛度矩陣，P 是網(wǎng)格節(jié)點(diǎn)受力向量，a 是網(wǎng)格節(jié)點(diǎn)位移向量。

K 由兩個(gè)參數(shù)決定:楊氏模量E 和泊松比v，選取楊氏模量E = 3 000 Pa 和泊松比v = 0．45［13］。

1.4 邊界條件獲取與應(yīng)用

為了求解式(5)中a 的唯一解，需要邊界條件約束，即已知部分網(wǎng)格節(jié)點(diǎn)的位移，包括受裂縫影響網(wǎng)格節(jié)點(diǎn)位移和零位移網(wǎng)格節(jié)點(diǎn)。受裂縫影響網(wǎng)格節(jié)點(diǎn)包括腦組織牽拉的初始位置，即切開腦組織位置，可以從術(shù)前圖像獲取;還有這些網(wǎng)格節(jié)點(diǎn)受裂縫影響產(chǎn)生的位移，可以通過(guò)對(duì)比牽拉前后掃描影像上的腦壓板坐標(biāo)變化。零位移網(wǎng)格節(jié)點(diǎn)為腦干部位網(wǎng)格節(jié)點(diǎn)，可以將兩種邊界條件引入式(5)進(jìn)行迭代運(yùn)算，求得牽拉變形后網(wǎng)格。

1.5 牽拉變形結(jié)果更新

用XFEM 計(jì)算得到的牽拉變形網(wǎng)格更新術(shù)前圖像獲得變形后的圖像，以此改進(jìn)了反插值算法［12］。原有算法只能處理網(wǎng)格節(jié)點(diǎn)的位移信息，能正確更新部分腦組織(非牽拉區(qū)域)圖像，無(wú)法正確更新?tīng)坷瓍^(qū)域圖像。改進(jìn)的算法能夠正確處理牽拉變形網(wǎng)格中的裂縫信息，所以能更新包括牽拉區(qū)域在內(nèi)的所有牽拉變形后的腦組織圖像。

從變形后的六面體單元出發(fā)，尋找出單元內(nèi)所有像素點(diǎn);根據(jù)單元內(nèi)8 個(gè)節(jié)點(diǎn)位移場(chǎng)計(jì)算出這些像素點(diǎn)的變形位移場(chǎng);再根據(jù)位移場(chǎng)，找到像素點(diǎn)在未變形前圖像中的位置。利用八鄰域內(nèi)三線性插值獲得像素點(diǎn)的灰度值，將此作為變形后該點(diǎn)的灰度值。經(jīng)過(guò)對(duì)所有變形后的單元進(jìn)行上述處理，從而獲得變形后的三維數(shù)據(jù)場(chǎng)，并對(duì)其進(jìn)行三維顯示。

2 試驗(yàn)驗(yàn)證

2.1 腦組織體模概述

為驗(yàn)證所提出的系統(tǒng)框架，選用腦組織體模，并采用iCT 圖像作為驗(yàn)證金標(biāo)準(zhǔn)。腦組織體模如圖3 所示，重1 ～2 kg，由6%水凝膠PVA-C 材料制備，外形接近大腦形態(tài)，楊氏模量E ?3 000 Pa，變形在20%內(nèi)，應(yīng)力應(yīng)變曲線接近線性。

圖3 腦組織體模Fig．3 The brain phantom

將腦組織體模固定于有機(jī)玻璃容器中，在容器外側(cè)不同高度和底部表面粘貼共10 個(gè)2 mm 的塑料標(biāo)記物，用于配準(zhǔn)獲取邊界條件。在雙側(cè)體模的額葉、頂葉和枕葉部位，分別植入4、13 和7 枚直徑從1 ～2 mm 不等的不銹鋼標(biāo)記物。其中，頂葉為牽拉區(qū)域，13 枚標(biāo)志物都位于較淺的裂縫周圍組織中，額葉和枕葉除了標(biāo)記物23 和24，其他均為淺部標(biāo)記物。這些標(biāo)記物隨體模的運(yùn)動(dòng)而運(yùn)動(dòng)，因此可以對(duì)比這些標(biāo)記物在體模內(nèi)牽拉變形前后的位移變化，對(duì)比XFEM 框架矯正前后圖像中這些標(biāo)記物的位移變化，評(píng)估模型的精度。

2.2 圖像分割與網(wǎng)格化

在術(shù)前圖像的感興趣區(qū)域設(shè)定種子點(diǎn)，通過(guò)區(qū)域增長(zhǎng)算法尋求其最大連通域，并手動(dòng)填充未選中的感興趣區(qū)域;再使用腐蝕算法處理，將腦組織與有機(jī)玻璃容器之間殘存的細(xì)微連接去除;最后，重建獲得僅包含腦組織的三維幾何模型。為便于XFEM 建模，因此對(duì)三維模型表面進(jìn)行平滑，以此來(lái)簡(jiǎn)化模型表面。然后進(jìn)行六面網(wǎng)格化。

2.3 邊界條件獲取與應(yīng)用

選擇螺旋模式，0.5 mm 層厚，無(wú)間隔掃描的CT(Siemens SOMATOM Definition AS)序列。首先，人工在體模上切開一道裂縫，將此設(shè)定為腦壓板初始位置，進(jìn)行掃描作為術(shù)前圖像;然后，插入腦壓板分別向右側(cè)牽拉5.7 mm，向左側(cè)牽拉4.8 mm，總計(jì)10.5 mm，進(jìn)行術(shù)中掃描。兩次掃描都重建為(512×512 × 421)像素的斷層圖像。最后，利用容器上的標(biāo)記物，將術(shù)中的圖像配準(zhǔn)到術(shù)前圖像坐標(biāo)系中，獲得從術(shù)中圖像變換到術(shù)前圖像的變換關(guān)系，找到術(shù)中圖像中的腦壓板在術(shù)前空間中的坐標(biāo)，比較變形前后圖像中的腦壓板坐標(biāo)變化，獲得腦壓板的牽拉位移。此外，零位移節(jié)點(diǎn)從腦干部位網(wǎng)格表面手工選取。

在體模的幾何模型進(jìn)行六面體網(wǎng)格化的過(guò)程中，為研究不同尺寸網(wǎng)格對(duì)XFEM 計(jì)算矯正牽拉變形的影響，采用5 種尺寸網(wǎng)格來(lái)比較矯正牽拉變形的結(jié)果，分別為2.5 mm (70 370個(gè)節(jié)點(diǎn)，65 044個(gè)單元)、3.5 mm (28 899 個(gè)節(jié)點(diǎn)，26 098個(gè)單元)、4.75 mm (13 221 個(gè) 節(jié) 點(diǎn)，11 648 個(gè) 單 元)，5.0 mm(11 614個(gè)節(jié)點(diǎn)，10 183個(gè)單元)和7.5 mm (4 373個(gè)節(jié)點(diǎn)，3 689個(gè)單元)。

建立XFEM 框架的線性方程組，利用工作站(Intel Xeon E5540 處理器，12GB 內(nèi)存)進(jìn)行迭代運(yùn)算，求得網(wǎng)格節(jié)點(diǎn)的位移和與裂縫相關(guān)信息。

2.4 圖像更新與驗(yàn)證結(jié)果評(píng)估

XFEM 計(jì)算得到變形網(wǎng)格后，用改進(jìn)的反插值算法更新術(shù)前圖像，得到牽拉變形后圖像。

模型的預(yù)測(cè)誤差由XFEM 框架矯正的圖像與實(shí)際變形圖像中不銹鋼標(biāo)記物坐標(biāo)差計(jì)算得出，模型的矯正精度由矯正后圖像中不銹鋼標(biāo)記物的位移與真實(shí)位移的百分比得出，有

預(yù)測(cè)誤差 =‖B - C‖2

矯正精度 =(‖B - A‖2- ‖C - A‖2)× 100%式中，A 為術(shù)前圖像，B 為矯正變形后圖像，C 為實(shí)際術(shù)中圖像。

3 試驗(yàn)結(jié)果

使用XFEM 矯正10.5 mm 的牽拉變形，植入體模的24 個(gè)不銹鋼標(biāo)記物，實(shí)際變形大小為0 ～6 mm，平均(2.3 ±1.4)mm。

在利用基于XFEM 計(jì)算得到的變形后網(wǎng)格對(duì)術(shù)前圖像進(jìn)行更新時(shí)，改進(jìn)的反插值算法能夠識(shí)別節(jié)點(diǎn)位移信息和裂縫信息。如圖4 所示，分別將改進(jìn)的反插值方法和傳統(tǒng)的反插值法更新得到的圖像與實(shí)際牽拉圖像進(jìn)行對(duì)比;(a)是利用傳統(tǒng)的反插值算法更新得到的圖像，明顯看出牽拉裂縫處的形態(tài)無(wú)法準(zhǔn)確更新，(b)是改進(jìn)了反插值算法更新得到的圖像，可看出已準(zhǔn)確地更新了牽拉裂縫圖像;(c)是實(shí)際牽拉圖像?？梢钥闯觯?a)與(c)在圖像裂縫處的形態(tài)相差較大，而(b)與(c)在裂縫處的形態(tài)吻合較好。

圖4 術(shù)中牽拉圖像的更新結(jié)果。(a)傳統(tǒng)反插值法更新圖像的結(jié)果;(b)改進(jìn)的反插值法更新圖像結(jié)果;(c)實(shí)際牽拉圖像Fig．4 The images updated by two algorithms． (a)The warped images updated by traditional backinterpolation algorithm; (b) The warped images updated by modified back- interpolation algorithm;(c)Real retraction images

表1 列出采用XFEM 框架矯正相同牽拉變形時(shí)網(wǎng)格尺寸對(duì)模型矯正精度的影響，網(wǎng)格尺寸分別為2.5、3.5、4.75、5.0、7.5 mm。尺寸最小的2.5 mm網(wǎng)格無(wú)法計(jì)算出牽拉變形的結(jié)果;3.5 mm 網(wǎng)格的預(yù)測(cè)誤差為0.4 mm，矯正精度為(87.1 ± 16.4%);4.75 mm 網(wǎng)格的預(yù)測(cè)誤差為0.4 mm，矯正精度為(85.8 ± 16.2)%;5.0 mm 網(wǎng)格的預(yù)測(cè)誤差為0.3 mm，矯正精度為(87.5 ±16.0)%;尺寸最大的7.5 mm 網(wǎng)格的預(yù)測(cè)誤差為0.5 mm，矯正精度為(83.1±20.3)%。從預(yù)測(cè)上看，5.0 mm 網(wǎng)格的預(yù)測(cè)誤差最小，7.5 mm 網(wǎng)格的預(yù)測(cè)誤差最大，平均預(yù)測(cè)誤差為0.4 mm。從矯正精度上看，5.0 mm 網(wǎng)格的矯正精度最高，7.5 mm 網(wǎng)格的矯正精度最低，平均矯正精度為85.9%。從計(jì)算時(shí)間上看，3.5 mm 網(wǎng)格的計(jì)算時(shí)間為4 440 s，遠(yuǎn)遠(yuǎn)長(zhǎng)于其他幾種網(wǎng)格的計(jì)算時(shí)間;7.5 mm 網(wǎng)格的計(jì)算時(shí)間最短，僅僅60 s。

表1 不同尺寸網(wǎng)格對(duì)XFEM 矯正牽拉變形結(jié)果的影響Tab． 1 The correction accuracy and calculating time of different size meshes

圖5 是尺寸為5.0 mm 的網(wǎng)格矯正牽拉變形后圖像中不銹鋼標(biāo)記物的位移和實(shí)際術(shù)中圖像中不銹鋼標(biāo)記物的位移對(duì)比統(tǒng)計(jì)結(jié)果。在圖5(a)中，菱形表示尺寸為5.0 mm 的網(wǎng)格矯正牽拉后圖像中不銹鋼標(biāo)記物位移，正方形表示實(shí)際術(shù)中圖像中不銹鋼標(biāo)記物位移，三角形表示對(duì)應(yīng)每一個(gè)不銹鋼標(biāo)記物的預(yù)測(cè)誤差。圖5(b)為尺寸為5.0 mm 網(wǎng)格中每一個(gè)不銹鋼標(biāo)記物的矯正精度，所有的均大于50%;標(biāo)記物5 ～18 正好位于頂葉牽拉區(qū)域，其矯正精度處于波動(dòng)中;最大位移出現(xiàn)在標(biāo)記物13，零位移出現(xiàn)在枕葉標(biāo)記物24 和25。

圖6 為5.0 mm 網(wǎng)格建立的XFEM 框架矯正牽拉變形的三維結(jié)果。圖6(a)是術(shù)前掃描得到的網(wǎng)格，包含了裂縫位置，即腦壓板牽拉的初始位置。圖6(b)是XFEM 框架矯正牽拉變形網(wǎng)格的云紋圖，紅色表示變形最大位置，藍(lán)色代表無(wú)變形位置;可以看出，受到牽拉影響，網(wǎng)格裂縫周圍出現(xiàn)的大小不同的位移，最大位移出現(xiàn)在裂縫周圍。圖6(c)是XFEM 矯正牽拉變形更新得到的圖像三維重建結(jié)果。

圖5 尺寸5.0 mm 的網(wǎng)格矯正牽拉變形后圖像中不銹鋼標(biāo)記物位移與實(shí)際術(shù)中圖像中標(biāo)記物位移的統(tǒng)計(jì)結(jié)果。(a)5.0 mm 網(wǎng)格矯正變形后圖像中不銹鋼標(biāo)記物位移，實(shí)際術(shù)中圖像中不銹鋼標(biāo)記物位移和預(yù)測(cè)誤差;(b)網(wǎng)格中每一個(gè)不銹鋼標(biāo)記物的矯正精度Fig．5 The statistical results of metal markers in the updated images with mesh size 5.0 mm comparing with markers in intraoperative images．(a)The displacement of metal markers in the updated images and intraoperative images，and the correction errors of each marker;(b)The correction accuracy of each metal marker

從腦壓板末端形態(tài)與受裂縫影響網(wǎng)格節(jié)點(diǎn)分布關(guān)系，探討不同網(wǎng)格尺寸對(duì)XFEM 矯正牽拉變形的影響。如圖7 所示，黑色半圓形曲線表示腦壓板末端的形態(tài)，黑色網(wǎng)格為局部網(wǎng)格，深黑色點(diǎn)為插入腦壓板后受其產(chǎn)生裂縫影響的網(wǎng)格節(jié)點(diǎn)。圖7(a)是3.5 mm 網(wǎng)格，節(jié)點(diǎn)的分布完全能夠反映腦壓板末端形態(tài);圖7(b)是4.75 mm 網(wǎng)格，節(jié)點(diǎn)的分布大致能夠反映出腦壓板末端形態(tài);圖7(c)是5.0 mm 網(wǎng)格，節(jié)點(diǎn)的分布大致能夠反映腦壓板末端形態(tài);圖7(d)是7.5 mm 網(wǎng)格，節(jié)點(diǎn)的分布無(wú)法正確反映腦壓板末端形態(tài)。

4 討論

筆者提出了一種將擴(kuò)展有限元(XFEM)應(yīng)用于線彈性生物力學(xué)模型來(lái)矯正牽拉變形的系統(tǒng)框架，并利用腦組織體模，實(shí)際驗(yàn)證了該牽拉變形矯正框架的有效性。驗(yàn)證結(jié)果表明，該系統(tǒng)框架能夠有效地矯正術(shù)中的腦組織牽拉變形，提高IGNS 術(shù)中導(dǎo)航的精確性及可靠性。

圖6 XFEM 矯正框架的工作流程。(a)術(shù)前網(wǎng)格;(b)矯正牽拉后變形網(wǎng)格;(c)矯正牽拉變形后圖像三維重建結(jié)果Fig．6 Illustration of the workflow of our XFEM framework． (a)Pre-deformed mesh;(b)The deformed mesh updated by our framework;(c)The three-dimensional brain phantom visualized by ray casting algorithm

在利用XFEM 模擬計(jì)算出的變形網(wǎng)格更新圖像時(shí)，改進(jìn)的反插值算法與原算法相比，能準(zhǔn)確地更新受牽拉部位的腦組織圖像，盡管圖像的外輪廓較為平滑，但并不影響牽拉變形矯正結(jié)果評(píng)估的客觀性。這是因?yàn)槟X壓板是對(duì)腦組織內(nèi)部施加力，而腦組織表面的位移是皮層下運(yùn)動(dòng)傳導(dǎo)的結(jié)果，外輪廓平滑不會(huì)影響腦組織內(nèi)部的位移。此外，模型矯正精度的評(píng)估主要依靠植入腦組織內(nèi)部的不銹鋼標(biāo)記物，腦組織表面平滑不會(huì)影響皮層下標(biāo)記物的位置。

驗(yàn)證試驗(yàn)中植入腦組織體模的不銹鋼標(biāo)記物主要位于較淺腦組織，尤其是頂葉部位。由于受牽拉的腦組織一般較淺，因此可以最大程度地捕獲到腦組織牽拉對(duì)周圍組織的影響。Vigneron 并沒(méi)有采用植入標(biāo)記物方法，而是通過(guò)矯正前后圖像的Hausdorff 距離來(lái)評(píng)估矯正結(jié)果，這不能定量評(píng)價(jià)矯正結(jié)果的好壞，同時(shí)評(píng)價(jià)結(jié)果受到非剛性配準(zhǔn)精度的影響。

圖7 選取不同尺寸網(wǎng)格，利用XFEM 框架矯正牽拉變形，受裂縫影響網(wǎng)格節(jié)點(diǎn)分布與腦壓板末端形態(tài)之間的關(guān)系。(a)3.5mm 網(wǎng)格;(b)4.75mm 網(wǎng)格;(c)5.0mm 網(wǎng)格;(d)7.5 mm 網(wǎng)格Fig． 7 The relationship between nodes fully intersected by the crack and the end of the retractors in different sizes of meshes． (a)3.5mm mesh; (b ) 4.75mm mesh; (c)5.0mm mesh;(d)7.5mm mesh

在網(wǎng)格化的過(guò)程中，網(wǎng)格尺寸的選取應(yīng)該考慮到牽拉變形的大小。在試驗(yàn)中，單側(cè)腦壓板牽拉位移平均為5.2 mm，通過(guò)5 種不同尺寸網(wǎng)格模型矯正結(jié)果的對(duì)比，發(fā)現(xiàn)與牽拉位移最接近的5.0 mm 網(wǎng)格矯正精度最高。小于牽拉位移的3.5 mm 網(wǎng)格的矯正精度與5.0 mm 網(wǎng)格的矯正精度相近，而計(jì)算時(shí)間遠(yuǎn)長(zhǎng)于5.0 mm 網(wǎng)格的結(jié)果;選用遠(yuǎn)小于牽拉位移的2.5 mm 網(wǎng)格，無(wú)法計(jì)算出牽拉變形結(jié)果;大于牽拉位移的7.5 mm 網(wǎng)格的矯正精度最低。因此，網(wǎng)格尺寸相對(duì)于牽拉位移過(guò)小時(shí)，較多的節(jié)點(diǎn)和單元成倍地延長(zhǎng)了計(jì)算時(shí)間，同時(shí)線性方程組可能無(wú)法迭代達(dá)到收斂。在網(wǎng)格相對(duì)于牽拉位移過(guò)大時(shí)，雖然系統(tǒng)計(jì)算負(fù)擔(dān)減小，但會(huì)導(dǎo)致網(wǎng)格節(jié)點(diǎn)剛度矩陣K 相對(duì)于網(wǎng)格節(jié)點(diǎn)受力矩陣P 過(guò)大;在求解線性方程組中網(wǎng)格節(jié)點(diǎn)位移向量a 時(shí)，分母遠(yuǎn)大于分子，導(dǎo)致計(jì)算得到的網(wǎng)格節(jié)點(diǎn)位移肯定小于真實(shí)位移。

網(wǎng)格尺寸的選擇，不僅與牽拉變形的大小有關(guān)，而且與該尺寸網(wǎng)格中受腦壓板影響網(wǎng)格節(jié)點(diǎn)形態(tài)和腦壓板末端的形態(tài)一致性相關(guān)。在選取不同尺寸網(wǎng)格矯正牽拉變形時(shí)，雖然網(wǎng)格尺寸對(duì)矯正精度影響較小，但是隨著網(wǎng)格尺寸增大，模型的矯正精度有下降趨勢(shì)，初步認(rèn)為這可能與插入腦壓板后受裂縫影響的網(wǎng)格節(jié)點(diǎn)分布有關(guān)。網(wǎng)格中受裂縫影響網(wǎng)格節(jié)點(diǎn)的分布越接近腦壓板末端形態(tài)，對(duì)這些節(jié)點(diǎn)施加力就越接近真實(shí)腦壓板末端作用于腦組織的力。但是，當(dāng)網(wǎng)格尺寸變大時(shí)，網(wǎng)格中受裂縫影響節(jié)點(diǎn)的分布越來(lái)越不能反映腦壓板末端形態(tài)，對(duì)節(jié)點(diǎn)施加的力就越不能真實(shí)地反映腦壓板對(duì)腦組織的作用力。

在試驗(yàn)中，使用腦組織體模來(lái)驗(yàn)證線彈性XFEM 框架矯正腦組織牽拉變形的有效性，但這與開顱手術(shù)中腦組織的牽拉變形情況還是存在差距。后續(xù)工作準(zhǔn)備將矯正系統(tǒng)框架應(yīng)用于動(dòng)物試驗(yàn)和臨床試驗(yàn)，以進(jìn)一步驗(yàn)證該矯正系統(tǒng)的可靠性和實(shí)用性。此外，XFEM 框架使用iCT 來(lái)獲得術(shù)中腦組織牽拉變形圖像，雖然得到的邊界條件非常準(zhǔn)確，但是iCT 存在輻射且費(fèi)用高，因此后續(xù)將利用簡(jiǎn)單可靠的術(shù)中設(shè)備(如三維激光掃描儀或術(shù)中超聲)來(lái)代替iCT/iMR，在不降低本系統(tǒng)框架矯正精度的前提下，增加臨床的實(shí)用性。

［1］ Zhang C，Wang M，Song Z. A brain-deformation framework based on a linear elastic model and evaluation using clinical data［J］. IEEE Trans Biomed Eng，2011，58(1):191 -199.

［2］ Vigneron LM，Duflot MP，Robe PA，et al. 2D XFEM-based modeling of retraction and successive resections for preoperative image update［J］. Comput Aided Surg，2009，14(1 -3):1 -20.

［3］ Ferrant M，Nabavi A，Macq B，et al. Registration of 3D intraoperative MR images of the brain using a finite-element biomechanical model［J］. IEEE Trans Med Imaging，2001，20(12):1384 -1397.

［4］ Miga MI，Roberts DW，Kennedy FE，et al. Modeling of retraction and resection for intraoperative updating of images［J］.Neurosurgery，2001，49(1):75 -84，84 -85.

［5］ Lamprich BK，Miga MI. Analysis of model-updated MR images to correct for brain deformation due to tissue retraction［C］ //Medical Imaging 2003:Visualization，Image-Guided Procedures，and Display. San Diego:SPIE，2003:552 -560.

［6］ Sun H，Kennedy F，Carlson E，et al. Modeling of Brain Tissue Retraction Using Intraoperative Data［C］ //Barillot C，Haynor D，Hellier P，eds. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2004. Berlin Heidelberg:Springer，2004:3217，225 -233.

［7］ Duflot M，Nguyen-Dang H. A meshless method with enriched weight functions for fatigue crack growth［J］. International Journal for Numerical Methods in Engineering，2004，59(14):1945 -1961.

［8］ Mo?s N. Dolbow J，Belytschko T. A finite element method for crack growth without remeshing［J］. International Journal for Numerical Methods in Engineering，1999，46(1):131 -150.

［9］ Vigneron LM，Noels L，Warfield SK，et al. Serial FEM/XFEMbased update of preoperative brain Images using intraoperative MRI［J］. Int J Biomed Imaging，2012，2012:872783.

［10］ Vigneron LM，Warfield SK，Robe PA，et al. 3D XFEM-based modeling of retraction for preoperative image update ［J］.Comput Aided Surg，2011，16(3):121 -134.

［11］ Wittek A，Hawkins T，Miller K. On the unimportance of constitutive models in computing brain deformation for imageguided surgery［J］. Biomech Model Mechanobiol，2009，8(1):77 -84.

［12］ Liu Y，Song Z. A robust brain deformation framework based on a finite element model in IGNS［J］. Int J Med Robot，2008，4(2):146 -157.

［13］ Platenik LA，Miga MI，Roberts DW，et al. In vivo quantification of retraction deformation modeling for updated image-guidance during neurosurgery［J］. IEEE Trans Biomed Eng，2002，49(8):823 -835.