基于骨架驅動的牙齒形態(tài)設計

黃香，程筱勝，戴寧，閆國棟，董光雷

(南京航空航天大學機電學院，江蘇南京 210016)

隨著科技的飛速發(fā)展，計算機輔助設計(CAD)與計算機輔助制造(CAM)技術進入了口腔修復學領域，這不僅提高了制作口腔修復體的效率，而且也減輕了患者的痛苦。目前已有以下著名的口腔CAD/CAM修復系統(tǒng)［1］:德國 Sirona公司的 CEREC 3D系統(tǒng)、德國Kavo公司的Everest系統(tǒng)和丹麥3Shape A/S公司的3SHAPE系統(tǒng)等。

將空間幾何變形應用于牙齒形態(tài)的設計是這些系統(tǒng)中的一個關鍵技術。國內的口腔CAD相關研究一般集中在北京大學口腔醫(yī)學院［2］、第四軍醫(yī)大學［3］等口腔醫(yī)學院，研究的平臺還停留在利用逆向工程軟件作一些面向應用的二次開發(fā)，缺少核心技術。

1984年Barr［4］首次提出了整體和局部變形的思想。該變形手段只適用于特定的幾類幾何形狀，要產生任意的、比較復雜的形狀變形則有很大的難度，不適合牙齒表面形態(tài)的設計。1986年Sederberg等［5］提出了自由變形(FFD)的概念。該方法適用于柔性物體動畫，變形范圍廣，但工作量較為繁重。Celniker等［6］提出了基于物理模型的曲面造型方法，Lazarus等［7］則提出了軸變形(AxDf)的方法。依據(jù)所使用變形工具的不同，可將上述變形技術分為4類［8］:基于體的變形、基于曲面的變形、基于曲線的變形、基于點的變形。其中，曲線驅動變形技術具有多方面優(yōu)勢，既操作簡單(點驅動變形優(yōu)勢)、易于實現(xiàn)，同時也可以實現(xiàn)比較復雜的變形(體、面驅動變形優(yōu)勢)。本研究中的骨架線由參數(shù)曲線來表示。

Lazarus等［7］在構建標架時，若曲線某一點處的曲率不存在，則該點處標架構建失敗。建立模型與曲線的映射關系時，可以選擇模型上的點與線上的最近點建立關系。這些最近點可以用曲線的凸包性和Oslo算法［9］確定，也可以使用離散樣條曲線來確定。基于上述研究，作者提出了一種快速實現(xiàn)牙齒修復體表面形態(tài)設計的方法。

1 骨架驅動模型變形原理

骨架線驅動模型變形如圖1所示，變形原理為:先為模型M構造近似骨架C，建立近似骨架C與模型M之間的映射關系;移動曲線上的控制點使曲線從C變形到C';依據(jù)曲線在離散點處的新局部坐標，更新模型上的點坐標。這一變形過程可以用一個簡單的關系式:F:Ω3→C→Ω3來表示［10］。

圖1 軸線驅動變形Fig 1 Axial curve drive deformation

整個變形的流程如圖2所示:

圖2 變形流程圖Fig 2 Flowchart of deformatio n

在建立映射關系時，模型需要被劃分成一個個子空間，這些子空間的連續(xù)性由曲線的連續(xù)性來保證。由于曲線是連續(xù)且光滑的，并且模型的子空間被劃分的很密，因此，可以近似地認為這些子空間也是連續(xù)且光滑的。在變形時，可以保證模型的曲面具有較好的光順性。

2 “骨架線”的構建

參數(shù)曲線的表現(xiàn)形式有:貝塞爾(Bezier)曲線、均勻B樣條曲線和非均勻有理B樣條(NURBS)曲線等。其中Bezier曲線缺乏局部特性;均勻B樣條曲線雖然局部性較好，但是沒有保留Bezier曲線端點處的幾何性質;NURBS曲線同時可以克服上述的缺點，但由于引入了權因子，一旦權因子選取的不合適，可能會導致很壞的參數(shù)化，從而使問題變得復雜化［11］。綜合多方面的因素，我們采用3次準均勻B樣條構造參數(shù)曲線，它具備良好的端點性質和局部性質，支持曲線局部編輯，便于對模型的局部區(qū)域實施變形。

勾畫曲線的交互手段類似于文獻［12］中所使用的方法。給定一個原始牙齒模型，用戶依據(jù)牙齒的表面形態(tài)，對該模型進行勾畫，選定興趣區(qū)域，這個興趣區(qū)域被參數(shù)曲線首末兩控制點構成的切平面所限定。勾畫時，用戶可以自由地調節(jié)屏幕視點，使操作更加便利。勾畫完成后所得到的曲線位于模型的表面。用戶可以通過平移、旋轉等操作調整曲線的位置，使其控制點落于模型之內。把控制點位于模型內部的三維(3D)參數(shù)曲線定義為“骨架線”。

3 網(wǎng)格變形

3． 1 定義局部標架

為明確映射關系，首先需要在參數(shù)曲線(“骨架線”)上建立局部坐標系。在建立坐標系之前，需要將曲線均勻離散化，并把離散點數(shù)據(jù)存儲到容器中。構造曲線上的局部活動標架有多種方法。其中弗朗內特(Frenet)標架是一種最常見的活動標架，但該標架在曲線曲率消失的地方無法進行定位，而且，在曲線的拐點處，該標架的主法矢和副法矢有可能會反向或產生不必要的旋轉，進而可能會導致模型的暴力扭曲［13］。為避免上述的問題發(fā)生，采用如下方法來建立坐標標架。

將活動標架F表示為一個四元組(t，T，N，B)，其中t為標架原點對應曲線C的差值參數(shù)，T=C'(t)/‖C'(t)‖，為C上參數(shù)t點處的單位切向量，‖·‖為歐式范數(shù)，N、B也是標架原點處的單位向量，T、N、B兩兩正交構成F的三個坐標軸。

首先構造一個t=0處的起始標架。令T0為該處的單位切矢T0=C'(0)/‖C'(0)‖，然后將T0投影到XY平面得到T0p，將T0p旋轉90°后再正規(guī)化即為N0，如圖3所示。B0則由T0和N0的叉積計算得到，即B0=T0×N0。如果T0的方向與世界坐標系下的向量(0，0，1)的方向一致，則該計算過程將會失效。為了解決這個問題，我們檢驗T0及其投影向量的合理性，選擇世界坐標系下的任意一個平面(XY、YZ或者XZ)，保證T0投影到該平面之后的向量不是一個零向量。

圖3 起始標架計算Fig 3 The creation of the first local frame

沿曲線在t∈［0，1］的定義域內構建其余標架。已知當前t=ti處曲線上點的坐標為Pi、單位切矢Ti，位于該點處的平面 πi。將向量 Ni－1投影到平面 πi上，并將其正規(guī)化即為Ni，則第三個向量Bi為Bi=Ti×Ni。在計算過程中要保證 α ＜90°，α 為 Bi與 Bi－1之間的夾角大小即這種算法簡單高效，避免了曲線暴力扭曲的發(fā)生。

3． 2 劃分區(qū)域

在變形之前，需要構建曲線與原始網(wǎng)格上的所有點之間的映射關系。我們采用文獻［14］的方法來確定網(wǎng)格與曲線之間的關系。如圖4所示，〈·|·〉表示兩個向量之間的點積，vm表示網(wǎng)格上的點，Pi表示曲線上點的坐標即局部坐標原點，Ti即單位切矢，Pi和Ti確定平面πi。劃分區(qū)域時，所用的判定式定義如下:

圖4 模型點、標架原點與單位切矢位置關系Fig 4 The relationship among the model vertex，the origin of frame and the unit tangent

依據(jù)上述判別公式，程序每循環(huán)1次，需要遍歷的網(wǎng)格頂點數(shù)就減少一部分，這就有效地提高了計算速度。我們?yōu)榫植繕思芗捌鋵木W(wǎng)格點索引值建立一個結構體，用C++語言定義如下:

struct Pt_Mesh

{

std::vector＜HPoint＞base_pt;/*表示離散后的曲線上的點的坐標即 Pi，HPoint 是一個三維坐標點類*/

std::vector＜int＞ index; /*Pi所對應的那一塊網(wǎng)格數(shù)據(jù)點的索引值*/

std::vector＜float＞ind_rm;/*網(wǎng)格上的每個點所對應的一個軸向參數(shù)，在變形過程中保持不變，會在接下來的章節(jié)中進行描述*/

};

3． 3 模型變形

定義 πk－1和 πk分別是由向量(Nk－1，Bk－1)和(Nk，Bk)所張成的空間平面。位于這兩個平面之間的點集vj可以由 πk－1平面上的標架來確定它的局部坐標。

設OXYZ為全局坐標系，位于平面πk－1和πk之間的網(wǎng)格上的點為 vm，Pk－1為標架 Fk－1的原點，而 Tk－1、Nk－1、Bk－1構成了該標架的 3個坐標軸，則 vm在局部標架Fk－1下的坐標為(ˉx，ˉy，ˉz)，它們之間的變換關系如下:

這樣就建立了全局坐標系與局部坐標系之間變換公式。令vm沿著Tk－1的方向投影到平面πk－1和πk上的點分別為vmp和vmn，vm和軸向參數(shù)rm通過等式

關聯(lián)起來，如圖5所示。dist是兩個點之間的歐式距離，這個參數(shù)在變形過程中保持不變。

圖5 網(wǎng)格點vm分別投影到面πk－1和面πk，得到vmp和vmnFig 5 vmand its projections vmpand vmnon πk －1and πk，respectively

當曲線C變形后，曲線上的離散點處的標架將進行重新計算。設變化后的曲線為C'。由于在計算曲線C的起始標架時，T0投影到世界坐標平面時具有隨機性，若采用相同的方法來計算C'的起始標架，則有可能會出現(xiàn)兩次標架差別很大的情況。因此，我們采用如下的方法來計算變化后的起始標架［9］。曲線C在 t=0 處的標架為(T0，N0，B0)，T'0為變化后的曲線在t=0的單位切矢。令ω為T0和T'0之間的夾角，向量R=T0×T'0，N0和B0沿著R的方向旋轉ω角度得到N'0和B'0。起始標架計算完成后，其余標架的計算方法同原始曲線在t=tk處計算標架的方法。從而vm新的坐標值為

但是這一公式并沒有考慮到控制點在移動曲線上時，曲線產生軸向拉伸的情況。為了把這種拉伸效果傳遞到網(wǎng)格，將v'm的公式變形為

其中，v'mp為變化后的網(wǎng)格點在 π'k－1上的投影，v'mp與向量T'k－1所確定的直線與平面π'k的交點為v'mn，如圖5所示。

曲線上標架的數(shù)量決定了模型變形結果的平滑程度。如果標架的數(shù)量較少，則變形效果有可能令人不十分滿意。但如果標架數(shù)量過多的話，又會延長計算時間。因此，我們需要在兩者之間找到一個平衡點。

4 實驗與分析

開發(fā)實現(xiàn)平臺:Windows XP系統(tǒng)，Visual C++．net，Hoops圖形顯示包。作者通過切牙、磨牙的整體變形和雙尖牙的牙尖變形來驗證算法的有效性與可行性。由于曲線采用了3次準均勻B樣條，因此，曲線上的某個點發(fā)生變化時，其附近區(qū)域也會有相應的改變。圖6～8分別表示3類牙的變形效果。

4． 1 切牙整體變形

圖6 切牙整體變形Fig 6 The global deformation of anterior tooth

圖6顯示了切牙整體旋轉的過程:a圖是切牙原始模型;b圖是調整后的參數(shù)曲線，位于原始模型內部;c圖是變形后的切牙牙齒形態(tài)，其中虛線表示原始“骨架線”，實線表示變化后的“骨架線”，箭頭方向表示變形方向。

4． 2 磨牙整體變形

圖7 磨牙整體變形Fig 7 The global deformation of molar tooth

圖7顯示了磨牙整體拉伸的過程:a圖是磨牙原始模型;b圖是磨牙調整后的“骨架線”，位于原始模型內部;c圖是拉伸以后的磨牙牙齒形態(tài)，箭頭方向表示變形方向。

4． 3 雙尖牙局部變形

圖8 雙尖牙局部變形Fig 8 The local deformation of bicuspid tooth

圖8顯示了雙尖牙某一牙尖局部拉伸的過程:a圖是雙尖牙原始牙尖形態(tài);b圖是對變形區(qū)域進行曲線勾畫，以及調整參數(shù)曲線后的形態(tài)，曲線位于模型內部;c圖是變形后的雙尖牙牙尖形態(tài)，其中虛線表示原始“骨架線”，實線表示變化后的“骨架線”，箭頭方向表示變形方向。

選取切牙整體繞X軸順時針旋轉10°作為考察對象，如表1所示。

表1 切牙變形結果的比對(繞X軸旋轉10°)Tab 1 The comparison of anterior tooth's deformation results

由表1我們可以看出:當標架數(shù)量增多時，劃分區(qū)域的時間、局部標架時間和變形時間都會相應變長，但變形效果越來越好。當標架數(shù)量為200和500時，牙齒表面形態(tài)的效果相差無幾。考慮到程序的計算速度以及變形效果的好壞程度，選取每條曲線標架數(shù)量為200。

5 結束語

現(xiàn)今的口腔CAD/CAM修復，一般是調用數(shù)據(jù)庫中的標準牙冠作為修復體，該修復體是一個統(tǒng)計模型，盡管具備了牙齒形態(tài)的大部分特征，但每一個患者的牙齒形態(tài)特征還是存在一定差異性的。我們依據(jù)鄰牙、頜牙牙面的約束關系，可以對修復體進行整體變形，也可以對修復體的局部特征(例如牙尖)進行修改。我們提出的基于骨架線驅動的牙齒形態(tài)設計其優(yōu)點如下:(A)在牙齒表面勾畫曲線時，可以隨意地刪除、增加曲線上的控制點，若對整體都不滿意則可以進行重繪。這提高了操作的交互性與簡便性。(B)在建立骨架上的局部標架時，我們所使用的方法克服了Frenet標架的缺點，骨架線既可以是直線也可以是曲線。也就是說，即使線上的曲率不存在也可以在該點處定義標架;同時也克服了因曲線出現(xiàn)暴力扭曲而導致的不滿意的牙齒形態(tài)變形。(C)牙齒模型區(qū)域劃分時所用的算法計算簡單，實用性較好，可以有效地提高計算效率。(D)牙齒模型上的點從局部坐標映射回全局坐標時，我們考慮到曲線變形會出現(xiàn)軸向拉伸的情況，因此，為了使牙齒形態(tài)的變形具有更好的光順性，對等式(5)做了一些改進，進而將這種拉伸傳遞到網(wǎng)格。

該方法還有一些問題亟待解決，如多條參數(shù)曲線如何實現(xiàn)控制變形，以及在變形過程中如何通過相應算法消除局部自交。這些將是我們今后研究的重點。

航空航天大學，2006．

［2］呂培軍，李彥生，王勇，等．國產口腔修復CAD-CAM系統(tǒng)的研究與開發(fā)［J］．中華口腔醫(yī)學雜志，2002，37(5):367-370．

［3］金樹人，姚月玲，高勃，等．應用快速成型法制作磨牙樹脂全冠［J］．第四軍醫(yī)大學學報，2003，24(8):700-702．

［4］BARR A H．Global and local deformations of solid primitives［J］．Computers Graphics，1984，18(3):21-34．

［5］SEDERBERG T W，PARRY R．Free-form deformation of solid geometric model［J］．Computer Graphics，1986，20(4):151-160．

［6］CELNIKER G，GOSSARD D．Deformable curve and surface finite-elements for free-form shape design［J］．Computer Graphics，1991，25(4):257-266．

［7］LAZARUS F，COQUILLART S，JANCENE P．Axial deformation:an intuitive technique［J］．Computer Aided Design，1994，26(8):607-613．

［8］徐崗，汪國昭，陳小雕．自由變形技術及其應用［J］．計算機研究與發(fā)展，2010，47(2):344-352．

［9］魏斌，袁修干．基于NURBS曲面的軸變形方法［J］．北京航空航天大學學報，1997，23(5):546-550．

［10］JAMES G，DOMINIQUE B．A survey of spatial deformation from a user-centered perspective［J］．ACM Transactions on Graphics，2008，27(4):1-21．

［11］施法中．計算機輔助幾何設計與非均勻有理B樣條［M］．北京:高等教育出版社，2001:309．

［12］KHO Y，GARLAND M．Sketching mesh deformations［C］//Proceedings of the 2005 Symposium on Interactive 3D Graphics．Washington D C:ACM，2005:147-154．

［13］ BLOOMENTHAL J．Graphics Gems［M］．San Diego:Academic Press Professional，1990:567-571．

［14］FORSTMANN S，OHYA J，KROHN-GRIMBERGHE A，et al．Deformation styles for spline-based skeletal animation［C］//Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on computer animation，San Diego:Eurographics Association，2007:141-150．