劉秀玲，陳棟，董亞龍，董斌，王洪瑞

（1.河北大學(xué) 電子信息工程學(xué)院 醫(yī)工交叉研究中心，河北 保定 071002；2.河北大學(xué)附屬醫(yī)院 計算機(jī)中心，河北 保定 071000）

隨著計算機(jī)圖形學(xué)技術(shù)、數(shù)字圖形處理技術(shù)、虛擬現(xiàn)實技術(shù)等技術(shù)的快速發(fā)展，虛擬手術(shù)作為一個多學(xué)科交叉的新領(lǐng)域，越來越廣泛地應(yīng)用到醫(yī)學(xué)手術(shù)的訓(xùn)練和相關(guān)醫(yī)學(xué)的研究中.虛擬手術(shù)訓(xùn)練要達(dá)到與傳統(tǒng)訓(xùn)練同樣的效果，就必須保證虛擬仿真過程的實時性和模型的真實感，因而建立一個快速、高效、精確的軟組織模型是虛擬形變仿真中最為關(guān)鍵的一部分［1－2］.

三角網(wǎng)格模型因其具有簡單的模型結(jié)構(gòu)、較好的通用性及網(wǎng)格穩(wěn)定性，廣泛應(yīng)用在真實感圖形模型的構(gòu)建中.三維網(wǎng)格細(xì)分曲面模型是在原始三角網(wǎng)格曲面基礎(chǔ)上經(jīng)過網(wǎng)格細(xì)分算法獲取的曲面模型，使得經(jīng)過處理之后的三角網(wǎng)格模型更準(zhǔn)確地描述三維物體的幾何特性和物理特征.國內(nèi)外的學(xué)者對模型細(xì)分的方法進(jìn)行了大量研究，并提出了大量的著名方法，如Loop細(xì)分法［3－4］，Buterfly細(xì)分算法［5］.Sederberg等［6］在B 樣條構(gòu)建方法的基礎(chǔ)上引入了非均勻節(jié)點區(qū)間的概念，提出了非均勻細(xì)分曲面的方法.Zorin在經(jīng)典的Butterfly算法的基礎(chǔ)上進(jìn)行改進(jìn)使得網(wǎng)格模型具有更為光滑的特性［7］.張文明等［8］提出了局域Delaunay剖分的動態(tài)點單元一體化三維網(wǎng)格的自適應(yīng)生成算法，使得模型可以準(zhǔn)確地反映出三維物體的幾何特征［8］.然而，常見的三角網(wǎng)格加密技術(shù)是對整體模型處理，雖然可以達(dá)到較為真實的形變效果，但是大幅增加了計算形變所需的時間.

本文采用了Loop三角網(wǎng)格細(xì)分算法對軟組織模型形變區(qū)域進(jìn)行局部處理，使得形變區(qū)域三角網(wǎng)的密度增加，在保證形變快速性的基礎(chǔ)上，提高了軟組織發(fā)生形變時的形變精確度和視覺沉浸感；同時對該區(qū)域的三角網(wǎng)格進(jìn)行拉普拉斯平滑優(yōu)化處理，從而進(jìn)一步提高了形變區(qū)域形變效果的真實程度.實驗結(jié)果表明該算法是準(zhǔn)確有效的，能夠較好地兼顧軟組織形變對準(zhǔn)確性和快速性的要求.

1 Loop網(wǎng)格細(xì)分算法

三維模型的表面網(wǎng)格大多數(shù)采用的是三角形網(wǎng)格的形式，網(wǎng)格中三角面片的大小、形狀直接影響模型的逼真程度.經(jīng)過網(wǎng)格細(xì)分之后的模型三角網(wǎng)格能夠精確地表達(dá)復(fù)雜模型表面的形體特征.細(xì)分方法具有一些主要的特點，可以概括為5個方面：拓?fù)涞倪m應(yīng)性；可伸縮性；網(wǎng)格唯一性；數(shù)值穩(wěn)定性；實現(xiàn)簡易性.將網(wǎng)格細(xì)分后的軟組織模型應(yīng)用到虛擬手術(shù)訓(xùn)練系統(tǒng)中，不僅可以反映軟組織的幾何特性，還可以在軟組織受力形變的仿真過程中更精確地體現(xiàn)出組織表面形變的效果.

本文對軟組織進(jìn)行網(wǎng)格自適應(yīng)細(xì)分的方法就是基于經(jīng)典的Loop網(wǎng)格細(xì)分算法.Loop細(xì)分算法采用的拓?fù)湟?guī)則是在每一個三角形網(wǎng)格的每一條邊上插入1個新頂點，再將各個頂點連接起來，將1個三角面片細(xì)分為4個，新生成的網(wǎng)格數(shù)量是原始網(wǎng)格數(shù)量的4倍，這種方式只生成了2種頂點，即邊生成的新頂點（E－頂點）和原來的頂點（V－頂點）［9］.

Loop算法的頂點獲得的規(guī)則如下：

1）內(nèi)部邊生成的點E－頂點：設(shè)內(nèi)部邊的2 個頂點分別為ν0，ν1，由此邊組成的三角形為Δν0ν1ν2和Δν0ν1ν3.則E－頂點定義為

2）邊界邊生成的點E－頂點，設(shè)該邊是由頂點ν0，ν1組成，則

3）對于每個內(nèi)部頂點V－頂點，設(shè)其相鄰頂點為ν0，ν1，…，νn，那么得到的頂點V－頂點為νv＝（1－nβn）ν＋其中βn 為鄰點權(quán)值，定義為

Loop細(xì)分算法生成的4類頂點如圖1所示.

圖1 Loop細(xì)分算法Fig.1 Loop subdivision algorithm

圖2 基于Loop的自適應(yīng)網(wǎng)格細(xì)分及優(yōu)化Fig.2 Subdivision and optimization based on the Loop subdivision algorithm

2 自適應(yīng)網(wǎng)格細(xì)分算法的設(shè)計與實現(xiàn)

軟組織形變模擬時用到的模型網(wǎng)格與一般實體模型網(wǎng)格所要達(dá)到的目的有所不同.通常模型網(wǎng)格劃分的目的是為了更準(zhǔn)確地描述出三維物體的幾何特性和物理特性，而形變模型不僅要反應(yīng)出軟組織的外形特征，還要針對外力的作用對形變區(qū)域進(jìn)行自適應(yīng)的網(wǎng)格加密，才能更準(zhǔn)確地模擬形變效果.

2.1 自適應(yīng)網(wǎng)格細(xì)分算法概述

通常自適應(yīng)網(wǎng)格細(xì)分技術(shù)的應(yīng)用是為了描述幾何體的幾何特征，對幾何體局部網(wǎng)格進(jìn)行細(xì)分.而本文在軟組織形變仿真中，自適應(yīng)網(wǎng)格細(xì)分技術(shù)應(yīng)用在針對形變區(qū)域的精細(xì)化操作中，目的是為了保證形變的真實性和良好的視覺效果.

曲面的細(xì)分是從原始的三角網(wǎng)格數(shù)據(jù)出發(fā)，根據(jù)預(yù)先設(shè)定好的細(xì)分迭代規(guī)則，在給定區(qū)域三角網(wǎng)格內(nèi)插入新的頂點，然后將這些行新頂點連接形成新的三角網(wǎng)格數(shù)據(jù)，實現(xiàn)對原始三角網(wǎng)格的細(xì)分效果.

首先要確定網(wǎng)格細(xì)分區(qū)域，然后為避免在新細(xì)分出的三角網(wǎng)格內(nèi)出現(xiàn)“較大”的三角面片，對需要進(jìn)行進(jìn)一步細(xì)分的三角面片制定了一個細(xì)分規(guī)則，最后將滿足要求的三角網(wǎng)格重新用Delaunay三角剖分方法進(jìn)行重構(gòu)，使得經(jīng)Delaunay三角剖分之后的三角網(wǎng)格質(zhì)量飽滿并且具有唯一性，最終達(dá)到理想的自適應(yīng)網(wǎng)格細(xì)分效果.算法實現(xiàn)的步驟如下：1）識別細(xì)分區(qū)域邊界，利用Loop細(xì)分算法添加新的頂點.2）區(qū)域內(nèi)的所有頂點進(jìn)行Delaunay三角剖分，構(gòu)建新的三角網(wǎng)格.3）判斷三角網(wǎng)格中的三角面片是否滿足要求，對不滿足要求的三角面片繼續(xù)進(jìn)一步細(xì)分操作.

自適應(yīng)網(wǎng)格算法的設(shè)計流程如圖2所示.

2.2 網(wǎng)格細(xì)分區(qū)域范圍的劃定

與以往自適應(yīng)網(wǎng)格劃分區(qū)域的目的不同，本文細(xì)分區(qū)域的劃定是為了保證形變區(qū)域在形變計算后形變效果的逼真，所以，網(wǎng)格加密區(qū)域范圍的劃定就是根據(jù)外力的值對模型表面形變區(qū)域的確定過程.當(dāng)外力作用在模型表面的某一位置時，根據(jù)質(zhì)點彈簧模型的力學(xué)方程式中，F(xiàn)ij表示質(zhì)點i，j之間的彈簧力；K 為彈簧的彈性系數(shù)；Lij為2個質(zhì)點之間靜止時原有的距離；pi，pj分別為質(zhì)點i，j的三維坐標(biāo)向量.求得質(zhì)點間的相互作用力，以及力在模型表面的傳導(dǎo)范圍，當(dāng)某2個質(zhì)點的力小于0.01N 時，視為外力到達(dá)作用邊界，不會再向外圍繼續(xù)傳導(dǎo).如圖3所示，假設(shè)外力作用位置為圖中點1的位置.

圖3 模型表面三角網(wǎng)格Fig.3 Triangle mesh of surface model

1）假設(shè)該質(zhì)點周圍的質(zhì)點是固定的，根據(jù)公式（1）計算出該點的虛擬位移量；

2）目標(biāo)質(zhì)點的移動會對周圍6個質(zhì)點產(chǎn)生彈簧拉力作用，根據(jù)這個拉力逐個計算出與目標(biāo)質(zhì)點相連的其他6個質(zhì)點的虛擬位移量.返回1），重復(fù)計算，直到包括目標(biāo)質(zhì)點在內(nèi)的7個質(zhì)點相對平衡為止；

3）分別以2，3，4，5，6，7作為目標(biāo)質(zhì)點，逐個計算與其相連質(zhì)點的虛擬位移量，進(jìn)行力的傳導(dǎo).返回1），直到這些質(zhì)點都到達(dá)力的平衡位置.

在這種方法中，力的計算只是為了獲取作用力的影響區(qū)域，質(zhì)點的位移并沒有顯示在界面上，因而不會對模型的渲染造成影響.最終獲得的這個影響區(qū)域就是需要進(jìn)行網(wǎng)格細(xì)分的區(qū)域.如圖4所示，中心深色區(qū)域即外力為10N 時的需要進(jìn)行細(xì)分區(qū)域.

2.3 基于幾何特性的細(xì)分規(guī)則制定

由于原始三角網(wǎng)格中的每個三角面片大小未必一致，為達(dá)到更好的三角網(wǎng)質(zhì)量，要使經(jīng)過細(xì)分之后生成網(wǎng)格的三角面片具有較為一致的幾何形狀，為此，需要對細(xì)化區(qū)域內(nèi)的三角面片是否需要進(jìn)一步細(xì)化進(jìn)行判斷.

設(shè)細(xì)分區(qū)域內(nèi)的某一個三角形邊長為l1i，l2i，l3i，當(dāng)max（l1i，l2i，l3i）＞α?xí)r，則該三角面片需要進(jìn)一步細(xì)分處理.α為設(shè)定的一個固定值，用來控制細(xì)分后三角網(wǎng)格中每個三角形的最大邊長.

3 實驗結(jié)果及分析

在CPU：Pentium E2200 2.20GHz，內(nèi)存2G，顯卡NVIDIA Geforce 9700GT，顯存512 MB 的PC 機(jī)上，以VTK 可視化工具包為基礎(chǔ)，用VS2005作為開發(fā)環(huán)境，實現(xiàn)本文提出的自適應(yīng)網(wǎng)格生成及優(yōu)化算法.

如圖4所示，在確定細(xì)分區(qū)域之后，將本文算法應(yīng)用到該區(qū)域的軟組織模型表面上進(jìn)行網(wǎng)格細(xì)分，得到的三角網(wǎng)格如圖5所示.

圖4 網(wǎng)格細(xì)分區(qū)域的劃定Fig.4 Subdivision area of the model

圖5 原始的三角網(wǎng)格與細(xì)化后的三角網(wǎng)格Fig.5 Original mesh and subdivision mesh

進(jìn)一步驗證本文算法生成的模型網(wǎng)格在形變仿真時的實時性，將原始的軟組織表面模型整體進(jìn)行Loop網(wǎng)格細(xì)分，與本文方法進(jìn)行對比.將外力作用在同一位置上，隨著作用力逐漸增大，受力質(zhì)點達(dá)到穩(wěn)定所消耗的時間也逐漸增大，記錄數(shù)據(jù)得到表1.繪制出作用力與消耗時間的對比圖，如圖6所示.由實驗數(shù)據(jù)可以得出，本文提出的自適應(yīng)網(wǎng)格細(xì)分算法，具有良好的仿真實時性，大幅縮短了采用傳統(tǒng)加密方法時消耗的時間.

表1 形變仿真耗時Tab.1 Time of deformation simulation

圖6 變形仿真耗時柱狀圖Fig.6 Time of deformation simulation

軟組織表面在形變過程中，經(jīng)過三角網(wǎng)格細(xì)化后的模型對于形變效果的表達(dá)能力有了一定的提高，但模型表面依然會出現(xiàn)“棱角”現(xiàn)象，如圖7a所示.為消除這種現(xiàn)象，對軟組織表面模型進(jìn)行拉普拉斯平滑優(yōu)化處理，經(jīng)過優(yōu)化后的模型效果如圖7b所示.對模型采取平滑優(yōu)化處理之后這種棱角得到了很好的解決，增強(qiáng)了變形的真實性，保證了虛擬手術(shù)訓(xùn)練對沉浸感的要求.

圖7 軟組織表面模型形變效果Fig.7 Deformation effect of soft tissue surface model

4 結(jié)論

本文提出的三角網(wǎng)格自適應(yīng)加密算法，較好地保證了形變結(jié)果的精確性，同時大大縮短了形變計算的時間，并且應(yīng)用平滑優(yōu)化技術(shù)使得模型表面更加光滑、平順，提升了視覺的沉浸感.為了進(jìn)一步提高虛擬手術(shù)訓(xùn)練系統(tǒng)的整體性能，將對形變模型的結(jié)構(gòu)特點和受力特點進(jìn)行更深一步的研究；完善網(wǎng)格細(xì)分的規(guī)則，使細(xì)分出的三角網(wǎng)格質(zhì)量更均勻；采用CUDA 并行技術(shù)提高形變計算的速度；對碰撞檢測程序進(jìn)行修改，保證碰撞檢測的快速和精確.

［1］ SEMERARO F，BERGAMASCO M，F(xiàn)RISOLI A，et al.Virtual reality prototype in healthcare simulation training［J］.Resuscitation，2008，77（1）：S60－S61.

［2］ ROMERO G，MAROTO J，F(xiàn)ELEZ J，et al.Virtual reality applied to a full simulator of electrical sub－stations［J］.Electric Power Systems Research，2008，78（3）：409－417.

［3］ LOOP C.Smooth subdivision surfaces based on triangles［D］.Salt lake City：University of Utah，1987.

［4］ CHENG F，F(xiàn)AN F，LAI S.Loop subdivision surface based progressive interpolation［J］.Journal of Computer Science and Technology，2009，24（1）：39－46.

［5］ DYN N，LEVIN D，GREGORY J A.A butterfly subdivision scheme or surface interpolatory with tension control［J］.ACM Transactions on Graphics，1990，9（2）：160－169.

［6］ SEDERBERG T，ZHENG J，SWELL D，et al.Non－uniform recursive subdivision surfaces［Z］.In Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques，New York，1998.

［7］ 賈太峰.基于任意拓?fù)渚W(wǎng)的細(xì)分改進(jìn)算法研究［D］.無錫：江南大學(xué)，2008.JIA Taifeng.The subdivision research on improved algorithm over arbitrary topological meshes［D］.Wuxi：Jiangnan University，2008.

［8］ 張文明，劉彬，徐剛.三維實體網(wǎng)格自適應(yīng)劃分算法［J］.機(jī)械工程學(xué)報，2009，45（11）：266－270.ZHANG Wenming，LIU Bin，XU Gang.Three dimensional entity mesh generation algorithm［J］.Journal of Mechanical Engineering，2009，45（11）：266－270.

［9］ 潘青，徐國良.常平面曲率細(xì)分曲面的構(gòu)造［J］.計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，2011，23（6）：964－970.PAN Qing，XU Guoliang.Construction of constant mean curvature subdivision surfaces［J］.Journal of Computer－aided Design ＆Computer Graphics，2011，23（6）：964－970.