王玉慧， 張玉茹

（北京航空航天大學(xué)機(jī)械工程及自動(dòng)化學(xué)院，北京 100191）

用Catmull-Clark細(xì)分及網(wǎng)格調(diào)整方法重構(gòu)牙齒曲面

王玉慧， 張玉茹

（北京航空航天大學(xué)機(jī)械工程及自動(dòng)化學(xué)院，北京 100191）

首先根據(jù)牙齒表面測(cè)量數(shù)據(jù)點(diǎn)，計(jì)算出其長方體包圍盒；并據(jù)此構(gòu)造細(xì)分曲面的初始網(wǎng)格；采用矩陣對(duì)角化方法，推導(dǎo)Catmull-Clark細(xì)分極限點(diǎn)的表達(dá)式，計(jì)算初始網(wǎng)格的頂點(diǎn)經(jīng)過細(xì)分后的極限點(diǎn)；按照極限點(diǎn)逼近數(shù)據(jù)點(diǎn)的原則移動(dòng)控制網(wǎng)格頂點(diǎn)，經(jīng)過逐次再細(xì)分、再調(diào)整網(wǎng)格，使各級(jí)網(wǎng)格在數(shù)據(jù)點(diǎn)的“引導(dǎo)”下逐步變形，使網(wǎng)格逐步逼近牙齒表面的測(cè)量數(shù)據(jù)點(diǎn)集合，實(shí)現(xiàn)牙齒表面模型的三維重建。

計(jì)算機(jī)輔助幾何設(shè)計(jì)；曲面重構(gòu)；細(xì)分造型；牙齒模型

在虛擬現(xiàn)實(shí)環(huán)境中，虛擬場(chǎng)景建模是一項(xiàng)重要工作，直接影響系統(tǒng)的性能。本文研究力覺交互虛擬現(xiàn)實(shí)牙科手術(shù)培訓(xùn)中由測(cè)量數(shù)據(jù)點(diǎn)重構(gòu)牙齒表面模型。

關(guān)于牙齒表面的重構(gòu)，LI Zhong[1]采用雙三次貝塞爾曲面進(jìn)行了牙齒表面模型的重構(gòu)，針對(duì)每一條輪廓線，用若干三次貝塞爾曲線段拼接得到；然后再對(duì)不同輪廓上的對(duì)應(yīng)點(diǎn)進(jìn)行貝塞爾曲線插值，得到由G2連續(xù)的雙三次貝塞爾曲面片拼接而成的表面模型。文獻(xiàn)[2]采用B-樣條曲面進(jìn)行牙齒表面曲面模型的重構(gòu)，由于B-樣條反算要解較大的線性方程組，計(jì)算量較大。Mikrogeorgis G[3]利用人類第六顆牙齒的斷層圖像得到斷層輪廓數(shù)據(jù)點(diǎn)，通過人機(jī)交互的方式進(jìn)行基于三角片的牙齒表面模型重構(gòu)。Isaac Newton Lima da Silva[4]由牙齒的斷層掃描圖像得到牙齒的斷層數(shù)據(jù)點(diǎn)，然后構(gòu)建牙齒的多面體模型，再利用商用軟件Pro/E得到牙齒的實(shí)體模型。

用曲面擬合進(jìn)行物體表面模型的重構(gòu)計(jì)算復(fù)雜、計(jì)算量較大；采用直接給物體表面的數(shù)據(jù)點(diǎn)建立拓?fù)潢P(guān)系的方法，所建立的物體表面模型的質(zhì)量依賴于測(cè)量數(shù)據(jù)點(diǎn)的測(cè)量精度及數(shù)據(jù)點(diǎn)的分布情況。本文嘗試采用細(xì)分造型構(gòu)建牙齒表面模型。

細(xì)分方法是通過將多邊形網(wǎng)格中的每個(gè)多邊形按照一定的規(guī)則分成幾個(gè)子多邊形，從而得到更光滑的網(wǎng)格。其算法簡單、直觀，適用于構(gòu)造復(fù)雜曲面，Chaikin[5]于 1974年提出了通過重復(fù)割去多邊形的角，最終得到光滑的極限曲線的方法，后來被證明該極限曲線是以多邊形為控制多邊形的二次 B-樣條曲線。Catmull和 Clark[6]提出基于四邊形的細(xì)分方法，其細(xì)分極限曲面是雙三次 B-樣條曲面，文中指出曲面在規(guī)則點(diǎn)處（關(guān)聯(lián)邊數(shù)為4）能達(dá)到曲率連續(xù)，而在異常點(diǎn)處曲面是切線連續(xù)。Doo D, Sabin[7]采用離散付立葉變換的方法證明了Catmull和Clark的結(jié)論。A A Ball[8-9]分別采用矩陣的特征值性質(zhì)和離散付立葉變換的方法證明了以上結(jié)論。Loop[10]提出了一種基于三角形網(wǎng)格的細(xì)分方法，較之基于四邊形網(wǎng)格的細(xì)分方法算法更簡單。Suzuki[11]提出了一種基于 loop細(xì)分和網(wǎng)格調(diào)整的曲面重構(gòu)算法。本文采用 Catmull-Clark細(xì)分方法，并針對(duì)數(shù)據(jù)點(diǎn)有噪聲，利用改進(jìn)的網(wǎng)格調(diào)整方法進(jìn)行牙齒曲面的三維重構(gòu)，力圖使重構(gòu)的網(wǎng)格在逼近牙齒數(shù)據(jù)點(diǎn)的同時(shí)具有較好的光順性。

1 算法原理

1.1 Catmull-Clark細(xì)分極限點(diǎn)表達(dá)式

由文獻(xiàn)[9]，經(jīng)過 Catmull-Clark計(jì)算細(xì)分后新的點(diǎn)點(diǎn)、邊點(diǎn)和面點(diǎn)的公式為

圖1 Catmull-Clark細(xì)分算法示意圖

求得矩陣的特征值和特征向量如表1所示。

表1 矩陣的特征值和特征向量

采用與n=4時(shí)同樣的求細(xì)分極限點(diǎn)的方法求得，當(dāng)n=3時(shí)

1.2 細(xì)分網(wǎng)格的構(gòu)造與調(diào)整

首先，計(jì)算牙齒表面數(shù)據(jù)點(diǎn)的包圍盒，將包圍盒略作放大，然后建立如圖2左所示的初始網(wǎng)格；對(duì)初始網(wǎng)格進(jìn)行一次 Catmull-Clark細(xì)分得到圖2右所示的網(wǎng)格；然后用式(3)、式(4)求出網(wǎng)格上每個(gè)頂點(diǎn)的細(xì)分極限點(diǎn)。

設(shè)牙齒數(shù)據(jù)點(diǎn)集合為 Pi(i = 0 ,1,… ,s )，網(wǎng)格頂點(diǎn)集合為 Vi(i = 0 ,1,… ,m)，其細(xì)分極限點(diǎn)集為( i = 0 ,1,… ,m)，其中 m為數(shù)據(jù)點(diǎn)個(gè)數(shù)。針對(duì) n =4和 n =3兩種情況，用式(6)和式(7)計(jì)算每一個(gè)極限點(diǎn)。再對(duì)于每一個(gè)極限點(diǎn)，在數(shù)據(jù)點(diǎn)集合中求出與它距離最近的點(diǎn) Pj。

δi則表示 Vi的對(duì)應(yīng)細(xì)分極限點(diǎn)與其最近數(shù)據(jù)點(diǎn) Pj之間的向量差。

則網(wǎng)格中所有頂點(diǎn)距離其最近數(shù)據(jù)點(diǎn)的平均距離為

接下來進(jìn)行網(wǎng)格調(diào)整，調(diào)整的原則是移動(dòng)控制網(wǎng)格頂點(diǎn)V使得其所對(duì)應(yīng)的新的細(xì)分極限點(diǎn)V∞'重合于其在牙齒表面數(shù)據(jù)點(diǎn)中的最近距離點(diǎn)。

如果要根據(jù)牙齒表面數(shù)據(jù)點(diǎn)按照網(wǎng)格頂點(diǎn)的調(diào)整原則計(jì)算出所有新的控制頂點(diǎn)V，則需要解龐大的線性方程組，本文參考[11]的方法，用每次僅移動(dòng)V而不改變iE，iF進(jìn)行迭代的方法，即，在考慮使得 'V 對(duì)應(yīng)的細(xì)分極限點(diǎn)逼近數(shù)據(jù)點(diǎn)時(shí)，只移動(dòng)V，而不考慮其周圍的關(guān)聯(lián)點(diǎn)。等到一個(gè)點(diǎn)作為V移動(dòng)完畢后，再更換一個(gè)點(diǎn)作為V，這時(shí)就僅考慮當(dāng)前V點(diǎn)的移動(dòng)。通過逐次迭代，使得調(diào)整后的網(wǎng)格點(diǎn)對(duì)應(yīng)的細(xì)分極限點(diǎn)逐步接近其對(duì)應(yīng)的最近點(diǎn)，直到滿足距離誤差ε＜E為止。

根據(jù)網(wǎng)格頂點(diǎn)的調(diào)整原則，原則上應(yīng)使得細(xì)分的極限點(diǎn)重合于牙齒表面數(shù)據(jù)點(diǎn)中距其最近的數(shù)據(jù)點(diǎn)，即使得

其中 α為細(xì)分極限點(diǎn)逼近牙齒表面數(shù)據(jù)點(diǎn)逼近項(xiàng)的系數(shù)；β為反映網(wǎng)格調(diào)整過程中，當(dāng)一個(gè)控制網(wǎng)格上的頂點(diǎn)被調(diào)整過程中，其周圍的點(diǎn)對(duì)該點(diǎn)的制約因素。這兩個(gè)參數(shù)的選取視具體應(yīng)用問題的要求，通過實(shí)驗(yàn)選定。

調(diào)整網(wǎng)格控制頂點(diǎn)算法步驟如下：

（1） 根據(jù)牙齒表面測(cè)量數(shù)據(jù)點(diǎn)建立初始網(wǎng)格；

（2） 分別用式(6)和式(7)計(jì)算網(wǎng)格中 4=n和 3=n 的每個(gè)控制頂點(diǎn)的 Catmull-Clark細(xì)分極限點(diǎn)；

（3） 求出牙齒表面測(cè)量數(shù)據(jù)點(diǎn)中距離每個(gè)極限點(diǎn)最近的點(diǎn)；

（4） 由式(9)計(jì)算網(wǎng)格上某一頂點(diǎn)的細(xì)分極限點(diǎn)與其最近數(shù)據(jù)點(diǎn)的距離；

（5） 由式(10)計(jì)算網(wǎng)格上所有點(diǎn)的細(xì)分極限點(diǎn)到其最近數(shù)據(jù)點(diǎn)的平均距離；

（6） 設(shè)定誤差閾值ε，如果ε＜E，則不再調(diào)整網(wǎng)格；

（7） 如果ε＞E，則對(duì)于網(wǎng)格中所有頂點(diǎn)用式(13)計(jì)算網(wǎng)格頂點(diǎn)V經(jīng)過一次移動(dòng)后的新位置 'V，將V移動(dòng)到 'V；

（8） 將調(diào)整后的網(wǎng)格進(jìn)行一次 Catmull-Clark細(xì)分；

（9） 用細(xì)分得到的新網(wǎng)格代替上一級(jí)網(wǎng)格，然后重復(fù)(2)～(9)的過程。

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

圖2 網(wǎng)格建立及調(diào)整過程

表2 初始網(wǎng)格參數(shù)

將初始網(wǎng)格進(jìn)行一次調(diào)整后的網(wǎng)格如圖2(b)所示，調(diào)整后再細(xì)分一次的網(wǎng)格如圖2(c)所示，在此基礎(chǔ)上再次調(diào)整后的網(wǎng)格如圖2(d)所示，圖2(e)、圖 2(f)分別為第二次調(diào)整網(wǎng)格后再進(jìn)行一次、二次細(xì)分后的網(wǎng)格。圖2(g)為牙齒模型的渲染圖。

圖3 平均誤差與迭代次數(shù)的關(guān)系

3 結(jié) 論

本文把細(xì)分算法應(yīng)用于力覺交互的虛擬現(xiàn)實(shí)牙科醫(yī)生手術(shù)培訓(xùn)系統(tǒng)的虛擬場(chǎng)景建模中，構(gòu)造出牙齒曲面的不同精確程度的網(wǎng)格，并且針對(duì)原始數(shù)據(jù)點(diǎn)云有噪聲的問題，在網(wǎng)格調(diào)整中加入了光順項(xiàng)使得構(gòu)造的牙齒表面網(wǎng)格在滿足精度要求的同時(shí)具備一定的光順性。且細(xì)分算法簡單，易于實(shí)現(xiàn)。

[1]LI Zhong, MA Li-zhuang, TAN Wu-zheng, et al.Reconstruction from contour lines based on bi-cubic Bézier spline surface [J]. Journal of Zhejiang University SCIENCE A, 2006, 7(7)：1241-1246.

[2]吳 雯. 人工牙的三維重建及其交互實(shí)現(xiàn)[D]. 北京：中國科學(xué)院, 2000.

[3]Mikrogeorgis G, Lyroudia K L, Nikopoulos N, et al.3D computer-aided reconstruction of six teeth with morphological abnormalities [J]. International Endodontic Journal, 1999, 32(2)：88-93.

[4]Isaac Newton Lima da Silva, Gustavo Frainer Barbosa,Rodrigo Borowski Grecco Soares, et al. Creating three-dimensional tooth models from tomographic images [J]. Stomatologija, Baltic Dental and Maxillofacial Journal, 2008, (10)：67-71.

[5]Chaikin G. An algorithm for high speed curve generation [J]. Computer Graphics and Image Processing, 1974, (3)：346-349.

[6]Catmull E, Clark J. Recursively generated B-spline surfaces on arbitrary topological surfaces [J].Computer-Aided Design, 1978, 10(6)：350-355.

[7]Doo D, Sabin M. Behaviour of recursive division surfaces near extraordinary points [J].Computer-Aided Design, 1978, 10(6) ：356-360.

[8]Ball A A, Storry D J T. Conditions for tangent plane continuity over recursively generated B-spline surfaces [J].ACM Trans. Graph., 1988, 7(2)：83-102.

[9]Ball A A, Storry D J T. A matrix approach to the analysis of recursively generated B-spline surfaces [J].Computer-Aided Design, 1986, 18(8)：437-442.

[10]Charles Teorell Loop. Smooth subdivision surfaces based on triangles [D]. Master's thesis, University of Utah, Department of Mathematics, 1987.

[11]Hiromasa Suzuki, Shingo Takeuchi, Fumihiko Kimura,et al. Subdivision surface fitting to a range of points[C]//Proceedings of the 7th Pacific Conference on Computer Graphics and Applications. Seoul Korea, 1999：158-167.

Tooth Surface Reconstruction by Catmull-Clark Subdivision and Mesh Modification

WANG Yu-hui, ZHANG Yu-ru
( School of Mechanical Engineering and Automation, Beijing University of Aeronautics and Astronautics, Beijing 100191, China )

Bounding box of data points are obtained according to data points of tooth surface, based on which initial control meshes are constructed. The limit points of Catmull-Clark subdivision are calculated by the expression derived by means of matrix diagonalization. Vertexes of control mesh are moved in order that each limit point approximates to its data point. Then subdivision and modification of the meshes are applied repeatedly so that the mesh at each level becomes deformed gradually under the “guide” of points to approximate to data points of tooth surface. Thus the three-dimensional reconstruction of tooth model is achieved.

computer aided geometric design; surface reconstruction; subdivision; tooth model

TP 391

A

1003-0158(2010)06-0056-06

2009-04-15

王玉慧（1964-），女，河南鄭州人，副教授，博士，主要研究方向?yàn)橛?jì)算機(jī)圖形圖像處理。