黃香，程筱勝，戴寧，閆國(guó)棟，董光雷

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

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

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

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

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

1 骨架驅(qū)動(dòng)模型變形原理

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

圖1 軸線驅(qū)動(dòng)變形Fig 1 Axial curve drive deformation

整個(gè)變形的流程如圖2所示:

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

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

2 “骨架線”的構(gòu)建

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

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

3 網(wǎng)格變形

3． 1 定義局部標(biāo)架

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

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

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

圖3 起始標(biāo)架計(jì)算Fig 3 The creation of the first local frame

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

3． 2 劃分區(qū)域

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

圖4 模型點(diǎn)、標(biāo)架原點(diǎn)與單位切矢位置關(guān)系Fig 4 The relationship among the model vertex，the origin of frame and the unit tangent

依據(jù)上述判別公式，程序每循環(huán)1次，需要遍歷的網(wǎng)格頂點(diǎn)數(shù)就減少一部分，這就有效地提高了計(jì)算速度。我們?yōu)榫植繕?biāo)架及其對(duì)應(yīng)的網(wǎng)格點(diǎn)索引值建立一個(gè)結(jié)構(gòu)體，用C++語(yǔ)言定義如下:

struct Pt_Mesh

{

std::vector＜HPoint＞base_pt;/*表示離散后的曲線上的點(diǎn)的坐標(biāo)即 Pi，HPoint 是一個(gè)三維坐標(biāo)點(diǎn)類(lèi)*/

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

std::vector＜float＞ind_rm;/*網(wǎng)格上的每個(gè)點(diǎn)所對(duì)應(yīng)的一個(gè)軸向參數(shù)，在變形過(guò)程中保持不變，會(huì)在接下來(lái)的章節(jié)中進(jìn)行描述*/

};

3． 3 模型變形

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

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

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

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

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

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

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

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

曲線上標(biāo)架的數(shù)量決定了模型變形結(jié)果的平滑程度。如果標(biāo)架的數(shù)量較少，則變形效果有可能令人不十分滿(mǎn)意。但如果標(biāo)架數(shù)量過(guò)多的話(huà)，又會(huì)延長(zhǎng)計(jì)算時(shí)間。因此，我們需要在兩者之間找到一個(gè)平衡點(diǎn)。

4 實(shí)驗(yàn)與分析

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

4． 1 切牙整體變形

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

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

4． 2 磨牙整體變形

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

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

4． 3 雙尖牙局部變形

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

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

選取切牙整體繞X軸順時(shí)針旋轉(zhuǎn)10°作為考察對(duì)象，如表1所示。

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

由表1我們可以看出:當(dāng)標(biāo)架數(shù)量增多時(shí)，劃分區(qū)域的時(shí)間、局部標(biāo)架時(shí)間和變形時(shí)間都會(huì)相應(yīng)變長(zhǎng)，但變形效果越來(lái)越好。當(dāng)標(biāo)架數(shù)量為200和500時(shí)，牙齒表面形態(tài)的效果相差無(wú)幾?？紤]到程序的計(jì)算速度以及變形效果的好壞程度，選取每條曲線標(biāo)架數(shù)量為200。

5 結(jié)束語(yǔ)

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

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

航空航天大學(xué)，2006．

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

［3］金樹(shù)人，姚月玲，高勃，等．應(yīng)用快速成型法制作磨牙樹(shù)脂全冠［J］．第四軍醫(yī)大學(xué)學(xué)報(bào)，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］徐崗，汪國(guó)昭，陳小雕．自由變形技術(shù)及其應(yīng)用［J］．計(jì)算機(jī)研究與發(fā)展，2010，47(2):344-352．

［9］魏斌，袁修干．基于NURBS曲面的軸變形方法［J］．北京航空航天大學(xué)學(xué)報(bào)，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］施法中．計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理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．