基于移動(dòng)頂點(diǎn)的肝臟模型漸進(jìn)切割算法

張峰峰， 黃軻， 于凌濤， 詹蔚

(1.蘇州大學(xué) 機(jī)電工程學(xué)院，江蘇 蘇州 215006； 2.蘇州大學(xué) 蘇州納米科技協(xié)同創(chuàng)新中心，江蘇 蘇州 215123； 3.哈爾濱工程大學(xué) 機(jī)電工程學(xué)院，黑龍江 哈爾濱 150001； 4.蘇州大學(xué)附屬第一醫(yī)院，江蘇 蘇州 215000)

原發(fā)性肝癌是常見的消化系統(tǒng)惡性腫瘤之一，其中約85%為肝細(xì)胞癌。目前，肝癌已成為全球第二大癌癥相關(guān)死亡原因，每年新增肝癌患者30萬～100萬例[1-2]。肝癌的治療復(fù)雜，部分肝臟切除手術(shù)是避免患者死亡和肝癌轉(zhuǎn)移最好的方法，但手術(shù)過程中醫(yī)生無法直觀地觀察到肝臟內(nèi)部的組織結(jié)構(gòu)。應(yīng)用增強(qiáng)現(xiàn)實(shí)(augment reality，AR)技術(shù)能夠?qū)⑻摂M肝臟模型疊加到真實(shí)圖像或視頻中，為醫(yī)生提供更多手術(shù)信息[3]。通過圖像配準(zhǔn)可以將肝臟虛擬模型與真實(shí)肝臟的位置重疊，疊加的虛擬模型能夠清晰地顯示肝臟的內(nèi)部結(jié)構(gòu)，為醫(yī)生提供直觀的指導(dǎo)。為了保證增強(qiáng)現(xiàn)實(shí)指導(dǎo)的準(zhǔn)確性，疊加的虛擬模型必須與真實(shí)場(chǎng)景的肝臟保持一致，虛擬模型需要實(shí)時(shí)進(jìn)行相應(yīng)的切割和變形。

肝臟模型的切割既要滿足一定的準(zhǔn)確性,又要滿足一定的實(shí)時(shí)性。實(shí)時(shí)性、穩(wěn)定性和真實(shí)性是衡量切割效果的關(guān)鍵因素[4]，直接決定了切割的效果。實(shí)時(shí)性與切割的刷新率有關(guān)，為了達(dá)到流暢的切割效果，通常要求刷新率不低于25 Hz。這就需要算法滿足一定的時(shí)間復(fù)雜度和空間復(fù)雜度要求[5-7]。

用于增強(qiáng)現(xiàn)實(shí)的虛擬模型需要由患者CT圖像或MRI圖像重建得到，目前三維重建的方法主要有2大類：體繪制、面繪制。體繪制效果更好，能更好地反映患者體內(nèi)真實(shí)情況，但體繪制模型計(jì)算量巨大，很難滿足實(shí)際使用的需求。面繪制模型計(jì)算量小，便于操作，是目前主流的三維重建方法[8]。

針對(duì)面繪制的模型，切割方法又主要有面模型切割法和體模型切割法。面模型切割法由頂點(diǎn)連接成大量三角形構(gòu)成模型的表面網(wǎng)格。切割過程中主要通過對(duì)三角形單元進(jìn)行操作，改變模型網(wǎng)格的拓?fù)浣Y(jié)構(gòu)。體模型切割法需要將網(wǎng)格模型進(jìn)行多面體劃分，使模型內(nèi)部以多面體單元填充，切割過程中主要對(duì)多面體單元進(jìn)行操作，常用的多面體有四面體、六面體和超六面體[9]。目前大多數(shù)學(xué)者集中于四面體切割法的研究，根據(jù)不同的剖分策略可將四面體切割的方法分為去除法和體元剖分法[10]。去除法是由Bronielsen等[11]提出的，去除法的基本思想是通過碰撞檢測(cè)找到與手術(shù)器械發(fā)生碰撞的四面體單元，將其從模型中刪除。但其真實(shí)感差，切口邊緣會(huì)產(chǎn)生嚴(yán)重的鋸齒。體元剖分法真實(shí)感強(qiáng)，切口細(xì)膩，但其計(jì)算量巨大，隨著切割的進(jìn)行，單元體數(shù)量急劇增加，很難滿足實(shí)時(shí)性要求[12]。另外Mor[13]提出了一種累進(jìn)切割法，其基本思想是在手術(shù)刀完全掠過某四面體單元之前提前將其分裂。存在計(jì)算量大的問題。Nienhuys[14]和Serby等[15]提出了一種移動(dòng)頂點(diǎn)的方法，在切割過程中搜索切割路徑附近的頂點(diǎn)，將其移動(dòng)到切割線上。但這種移動(dòng)模型頂點(diǎn)容易產(chǎn)生退化單元。

目前針對(duì)面模型切割的研究較少，各種方法分別存在計(jì)算量大和網(wǎng)格質(zhì)量差的問題，且大多數(shù)方法只對(duì)表面切割進(jìn)行了研究，沒有對(duì)切割面進(jìn)行研究。本文針對(duì)上述問題對(duì)面模型切割的頂點(diǎn)移動(dòng)法進(jìn)行了研究和優(yōu)化，并按照Delaunay三角剖分準(zhǔn)則對(duì)局部三角面片進(jìn)行重構(gòu)，最后對(duì)切割面的構(gòu)造進(jìn)行了研究。

1 切割算法

如圖1所示，本文算法的核心思路為：

1)進(jìn)行碰撞檢測(cè)，獲得切割點(diǎn)具體位置；

2)搜索最近頂點(diǎn)，將搜索到的最近點(diǎn)移動(dòng)到切割點(diǎn)的位置；

3)優(yōu)化局部面片，檢測(cè)最近點(diǎn)附近的退化三角形并進(jìn)行優(yōu)化；

4)分裂網(wǎng)格，將切割線兩側(cè)的頂點(diǎn)分別向不同的方向移動(dòng)一定的距離以形成切口；

5)構(gòu)造切割面，連接相應(yīng)頂點(diǎn)形成切割面。

圖1 移動(dòng)頂點(diǎn)算法思路Fig.1 Block diagram of moving vertex algorithm

以上所有計(jì)算在一幀中完成，稱為一個(gè)計(jì)算幀，即每一個(gè)計(jì)算幀從檢測(cè)到切割操作開始，切割面構(gòu)造完結(jié)束。

1.1 碰撞檢測(cè)

碰撞檢測(cè)用于決定何時(shí)發(fā)生切割并得到切割位置，是后續(xù)一切操作的基礎(chǔ)。如圖2所示，常用的碰撞檢測(cè)方法為包圍盒法，其中包圍盒又分為包圍球、軸對(duì)齊包圍盒(AABB包圍盒)、方向包圍盒(OBB包圍盒)和離散方向凸包圍盒(k-Dops包圍盒)[16-17]。包圍盒算法簡單，計(jì)算量小，但無論哪種包圍盒都無法完全貼合模型。為了能夠精確地檢測(cè)切割位置，本文針對(duì)面模型切割提出了一種新的碰撞檢測(cè)方法。

圖2 包圍盒算法對(duì)比示意Fig.2 Comparison of bounding box algorithms

在模型切割仿真中，切割工具的簡化形式有3種：1)將切割工具視為與網(wǎng)格表面的一系列碰撞點(diǎn)；2)使用簡單的平面或直線表示切割工具;3)渲染出完整的切割工具，但在切割計(jì)算時(shí)使用其中軸線進(jìn)行計(jì)算[18]。本文使用第3種方法(如圖3)，在手術(shù)刀刀刃上固定2個(gè)點(diǎn)Se和St，使用2點(diǎn)之間的連線進(jìn)行切割計(jì)算。

圖3 本文所采用手術(shù)刀具示意Fig.3 Schematic diagram of surgical knife used in this article

在切割過程中從點(diǎn)Se向點(diǎn)St投射一條光線，如果光線被肝臟模型遮擋無法到達(dá)St即發(fā)生切割。檢測(cè)到切割后查找遮擋光線的三角面片，記為切割三角形ΔT1T2T3。該三角面片所在平面可以由三點(diǎn)式平面方程得到：

(1)

過點(diǎn)Se和點(diǎn)St的直線由兩點(diǎn)式直線方程等到：

(2)

由式(1)、(2)即可得到切割點(diǎn)的精確位置Ci(xi,yi,zi)。

1.2 搜索最近點(diǎn)

在本文方法中準(zhǔn)確地獲取到最近點(diǎn)是算法的關(guān)鍵，其結(jié)果決定了切割能否成功以及切割的最終效果。本文使用與傳統(tǒng)模型切割不同的計(jì)算頻率，本文方法中手術(shù)刀每移動(dòng)一定距離，進(jìn)行一次碰撞檢測(cè)和切割計(jì)算，該方法不僅具有更好的延遲反饋，還能通過設(shè)置不同的步長來調(diào)整計(jì)算頻率，在精度和效率之間取得平衡。

通過計(jì)算ΔT1T2T3的頂點(diǎn)與切割點(diǎn)Ci距離的方式搜索最近點(diǎn)，但最近點(diǎn)并不總是取計(jì)算得出距離最近的點(diǎn)，在切割過程中應(yīng)該避免圖4所示的特殊情況。

圖4 特殊的最近點(diǎn)Fig.4 Special closest point

在圖4中，當(dāng)相臨2次計(jì)算得到的最近點(diǎn)Mi-1和Mi為一個(gè)四邊形的2個(gè)相對(duì)頂點(diǎn)時(shí)，Mi-1和Mi沒有直接連接，無法形成切割線，三角面片不能正常分裂。故每次計(jì)算得到臨時(shí)最近點(diǎn)Mi后，需要對(duì)其和上一個(gè)最近點(diǎn)Mi-1進(jìn)行判斷，判斷條件為：點(diǎn)Mi-1屬于當(dāng)前切割三角形ΔT1T2T3或點(diǎn)Mi屬于上一個(gè)切割三角形ΔT1T2T3。滿足該條件說明點(diǎn)Mi和點(diǎn)Mi-1屬于同一個(gè)三角形，2個(gè)點(diǎn)由三角形的一條邊相連，能夠正常形成切割線。若不滿足判斷條件，取距離次之的點(diǎn)做為臨時(shí)最近點(diǎn)再進(jìn)行判斷，如果如圖5所示點(diǎn)Mi-1和Mi仍然不能直接連接,則取第3個(gè)點(diǎn)做為臨時(shí)最近點(diǎn)。如果仍然不滿足條件，表示步長設(shè)置過大，2次切割檢測(cè)跨過了若干三角形，導(dǎo)致點(diǎn)Mi和點(diǎn)Mi-1所屬的2個(gè)三角形沒有公共邊或公共點(diǎn)，此時(shí)應(yīng)適當(dāng)減小檢測(cè)步長。

圖5 特殊的次最近點(diǎn)Fig.5 Special second closest point

1.3 優(yōu)化面片

狹長或小面積的退化三角形會(huì)影響后續(xù)的切割計(jì)算，還會(huì)導(dǎo)致模型顯示質(zhì)量下降。切割使用的原始模型遵循Delaunay三角剖分原則，不包含退化三角形。本文的移動(dòng)頂點(diǎn)法每次只移動(dòng)一個(gè)頂點(diǎn)，故不會(huì)出現(xiàn)小面積三角形，但并不能避免出現(xiàn)狹長三角形。因?yàn)楸疚牡姆椒看沃灰苿?dòng)一個(gè)最近點(diǎn)Mi，狹長三角形也只會(huì)出現(xiàn)在包含點(diǎn)Mi的三角形中。為了消除狹長三角形，只需要對(duì)包含點(diǎn)Mi的三角形進(jìn)行檢測(cè)優(yōu)化即可。

Delaunay三角剖分是一種由給定點(diǎn)集V形成優(yōu)質(zhì)表面網(wǎng)格的算法。Delaunay 三角剖分包含2個(gè)重要準(zhǔn)則，一個(gè)是空?qǐng)A特性：在Delaunay三角網(wǎng)中任一三角形的外接圓中不存在V中的其他頂點(diǎn)；另一個(gè)是最大化最小角特性：在給定點(diǎn)集V可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大[19]。

圖6 Delaunay三角剖分Fig.6 Delaunay triangulation

在本文方法中，模型頂點(diǎn)點(diǎn)集V并不是固定不變的，隨著切割的進(jìn)行，總會(huì)有部分頂點(diǎn)坐標(biāo)發(fā)生變化，不能直接使用Delaunay三角剖分。但在4個(gè)頂點(diǎn)構(gòu)成的2個(gè)相臨三角形中，如果2個(gè)三角形的外接圓范圍內(nèi)均不包含第4個(gè)頂點(diǎn)，則這4個(gè)頂點(diǎn)滿足空?qǐng)A特性，否則總是可以連接四邊形的另一個(gè)對(duì)角線使其滿足空?qǐng)A特性。且滿足空?qǐng)A特性連接方式的最小角大于不滿足空?qǐng)A特性連接方式的最小角。即滿足空?qǐng)A特性的連接方式也滿足最大化最小角特性。

圖7 交換對(duì)角線Fig.7 Exchange diagonal

本文方法對(duì)每個(gè)包含最近點(diǎn)Mi的三角形和其相臨三角形組成的點(diǎn)集進(jìn)行空?qǐng)A特性判斷，將不滿足空?qǐng)A特性的點(diǎn)集按照上述方式重新連接。保證每次移動(dòng)最近點(diǎn)后三角面片中沒有退化三角形。

1.4 構(gòu)造切割面

在網(wǎng)格模型的切割研究中更多的研究集中于模型表面的切割。不同于先生成完整的切割線再計(jì)算切口的方法，漸進(jìn)式切割的計(jì)算是隨著切割動(dòng)作逐漸進(jìn)行的，切割面的構(gòu)造必須也是逐漸進(jìn)行的。對(duì)于本文的方法即在每次頂點(diǎn)移動(dòng)后對(duì)切割面進(jìn)行計(jì)算。另外由于在漸進(jìn)式切割中不能確定當(dāng)前計(jì)算幀切割是否結(jié)束，所以在構(gòu)造切割面時(shí)每一個(gè)計(jì)算幀都要按照切割結(jié)束來生成閉合的切割面。圖8以切割線左側(cè)部分為例說明切割面的構(gòu)造方法。

圖8 切割面的構(gòu)造過程Fig.8 Cutting surface construction process

為了使切割真實(shí)性更好，通過在每一個(gè)計(jì)算幀中記錄手術(shù)刀尖端Sti的坐標(biāo)的方式得到切割深度。然后依次連接點(diǎn)Mi、Sti-1、Mi-1和點(diǎn)Mi-1、Sti-1、Sti-2構(gòu)成ΔMiSti-1Mi-1和ΔMi-1Sti-1Sti-2即可形成完整的切割面。

由于每次碰撞檢測(cè)為一個(gè)計(jì)算幀，因此手術(shù)刀接觸模型的瞬間即為第一個(gè)計(jì)算幀，此時(shí)表面最近M1與手術(shù)刀尖端點(diǎn)St1距離很近，可認(rèn)為點(diǎn)M1與點(diǎn)St1完全重合，ΔM2St1M1為狹長的退化三角形，故只保留點(diǎn)M1，刪除點(diǎn)St1和ΔM2St1M1。另外，在每一個(gè)計(jì)算幀中，切割線左右部分共用點(diǎn)Mi、點(diǎn)Sti和點(diǎn)Sti-1，即共用ΔMiStiSti-1，所以不需要構(gòu)造ΔMiStiSti-1。但在下一個(gè)切割循環(huán)中點(diǎn)Mi變?yōu)镸i-1并在三角面片的分裂過程中被復(fù)制一份連接切割線右側(cè)三角形，左右部分不再共用ΔMiStiSti-1，所以該三角形需要在下一個(gè)計(jì)算幀中構(gòu)造(如圖9)。

圖9 切割面上的特殊三角形Fig.9 Special triangle on the cutting surface

2 實(shí)驗(yàn)結(jié)果與分析

使用上述針對(duì)面模型切割的移動(dòng)頂點(diǎn)法對(duì)肝臟模型進(jìn)行切割實(shí)驗(yàn)。實(shí)驗(yàn)硬件環(huán)境為：Intel Core i3-8100(主頻3.6 GHz)，8 G內(nèi)存，NVIDIA Geforce GTX1050Ti顯卡，在Unity 2018.2和Visual Studio 2015中使用C#語言進(jìn)行編程實(shí)現(xiàn)。使用的肝臟模型是經(jīng)患者CT掃描分割后三維重建得到，由8 706個(gè)頂點(diǎn)和15 872個(gè)三角形單元組成。

圖10展示了針對(duì)面模型切割的移動(dòng)頂點(diǎn)法的切割效果，圖10(a)為直接將移動(dòng)頂點(diǎn)法應(yīng)用于面模型切割的效果。在點(diǎn)M1移動(dòng)后，出現(xiàn)退化三角形概率較大，雖然隨著切割的進(jìn)行，隨后的頂點(diǎn)移動(dòng)有可能使退化三角形正?；匦伦兂烧Ｈ切危荒鼙ＷC每一個(gè)退化三角形都能正?；?。圖10(b)為使用本文提出的方法優(yōu)化的結(jié)果。在原移動(dòng)頂點(diǎn)法出現(xiàn)退化三角形的位置經(jīng)過交換對(duì)角線的操作后成功消除了退化三角形，每次頂點(diǎn)移動(dòng)后都進(jìn)行檢測(cè)優(yōu)化可以確保每一個(gè)退化三角形都能正?；?。圖10(c)為本文方法構(gòu)造出的切割面。本文致力于肝臟面模型切割的拓?fù)浣Y(jié)構(gòu)變化研究，不涉及生物力學(xué)和變形，肝臟切口大小通過一個(gè)參數(shù)控制。為了更好的顯示切割面，圖10(c)將切口大小參數(shù)進(jìn)行了適當(dāng)?shù)姆糯蟆．?dāng)切口較小時(shí)切口邊緣頂點(diǎn)移動(dòng)距離小不會(huì)出現(xiàn)退化三角形，如圖10(c)切口增大時(shí)切口邊緣將出現(xiàn)退化三角形。在實(shí)際操作中不能預(yù)測(cè)切口大小，故需要對(duì)切口邊緣的三角形進(jìn)行優(yōu)化，使用與最近點(diǎn)相關(guān)三角形優(yōu)化相同的方法對(duì)切口邊緣三角形優(yōu)化結(jié)果如圖10(d)。

圖10 模型切割效果Fig.10 Model cutting effect

圖11 多路徑、交叉切割示意Fig.11 Schematic diagram of multi-path, cross cutting

圖12 交叉切割實(shí)驗(yàn)刷新率Fig.12 Cross cutting experiment refresh rate

另外本文對(duì)切割過程中出現(xiàn)切割路徑交叉和多次切割的情況進(jìn)行了實(shí)驗(yàn)，如圖11，并記錄了實(shí)驗(yàn)過程的刷新率。實(shí)驗(yàn)顯示本文的方法具有良好的穩(wěn)定性。切割過程中，圖像的平均刷新率達(dá)到了119 Hz(如圖12)，為模擬組織變形的計(jì)算留出了足夠的空間。

3 結(jié)論

1)本文提出光線投射算法進(jìn)行精確的切割檢根據(jù)Delaunay三角剖分準(zhǔn)則對(duì)頂點(diǎn)移動(dòng)后的局部面片進(jìn)行了優(yōu)化和重構(gòu)以消除原移動(dòng)頂點(diǎn)法可能出現(xiàn)的退化的畸形三角形。

2)本文還提出了一種針對(duì)面模型切割的切割面構(gòu)造方法，能夠在漸進(jìn)切割過程中逐步生成完整閉合的切割面，切面真實(shí)度大大提高。

3)本文提出的方法很好地解決表面切割模型的實(shí)時(shí)性問題并具有良好的穩(wěn)定性，且實(shí)時(shí)刷新率能夠穩(wěn)定在119±2 Hz，具有很好的實(shí)時(shí)性。