高云凱，張婷婷，楊肇通，羅 鍵

（同濟大學 汽車學院，上海201804）

用一定的測量手段對實物或模型進行數(shù)據(jù)測量，然后重構實物的CAD（computer aided design）模型，實現(xiàn)產(chǎn)品設計或制造過程，即為逆向工程中實物的逆向重構過程［1］.曲面重建是采用逆向工程方法生產(chǎn)CAD模型的關鍵步驟［2］.近年來，隨著計算機輔助幾何設計CGAD（computer graphic aided design）及測量技術的不斷發(fā)展，曲面重建技術在機械產(chǎn)品設計制造數(shù)字化過程中發(fā)揮著舉足輕重的作用.因此，如何高效獲得較高精度的重建曲面已經(jīng)成為越來越多生產(chǎn)設計者關注的重點.

實現(xiàn)曲面重建較為傳統(tǒng)的方法是基于特征面的區(qū)域分割法，具體流程如圖1a所示.該方法首先分析了曲面的微分信息，提取相應的特征邊或特征面并進行數(shù)據(jù)分塊，最后根據(jù)不同曲面的特性利用相關多項式進行擬合，完成曲面的逆向重構.國內(nèi)外許多學者在這方面研究較為深入，如Huang等［3］提出先基于網(wǎng)格邊的方向曲率來識別邊界，然后利用區(qū)域增長對網(wǎng)格面片進行分組的區(qū)域分割法；董明曉等［4］提出基于數(shù)據(jù)點曲率變化的區(qū)域分割法；劉勝蘭［5］提出基于三角網(wǎng)格模型的四邊域重構法等.這些算法結(jié)構簡單，對特征較為明顯或曲面較為規(guī)則的模型能夠較好實現(xiàn)曲面重構的質(zhì)量，但一旦涉及復雜程度較高的車身曲面或自由曲面，其適用性就會不足，原因是這些復雜曲面的特征相比一般規(guī)則區(qū)域并不明顯，且較為分散.

在傳統(tǒng)曲面重構算法的基礎上，為了適應復雜曲面的重構，一些學者提出了基于曲率線追蹤的曲面重構法.如圖1b所示，該算法的核心是依據(jù)獲取的主方向通過數(shù)值積分的方式［6］對曲率線進行追蹤.如Alliez等［7］提出在曲面參數(shù)空間求解微分方程來確定追蹤步的曲率線追蹤法，但該算法在全局參數(shù)化的基礎上進行追蹤映射過程中精度難以保證，而且計算量大；龐旭芳等［8］提出利用鄰域點信息建立曲率場并通過點對應的曲面微分信息來實現(xiàn)追蹤，但曲率線在曲面特征處變化不明顯；Marinov等［9］提出基于三角網(wǎng)格面片的局部參數(shù)化法來追蹤曲率線，Kalogerakis等［10］提出在離散點云數(shù)據(jù)中直接迭代步長的空間免映射追蹤法，這2種方法數(shù)值積分的過程較為繁瑣，且需要復雜的數(shù)據(jù)結(jié)構進行網(wǎng)格密度的控制.

本文基于曲率線追蹤的思想提出一種基于網(wǎng)格遍歷的曲面網(wǎng)格重劃算法，該算法簡化了曲率線追蹤過程，提出的鄰域點剔除法能夠靈活控制重劃網(wǎng)格密度，保證網(wǎng)格分布均勻可靠，實現(xiàn)較好的曲面網(wǎng)格重劃.算法的實現(xiàn)基于逆向過程中的三角網(wǎng)格模型，先采用局部參數(shù)化的一般二次曲面法來估算網(wǎng)格曲面的微分信息，并在此基礎上建立模型的半邊數(shù)據(jù)結(jié)構及提取模型的自由邊界特征，最終采用基于網(wǎng)格遍歷的曲率線算法采樣.

圖1 曲面重建算法流程Fig.1 Flowchart of surface reconstruction algorithm

1 三角網(wǎng)格頂點數(shù)據(jù)預處理

1.1 頂點主曲率場的建立

通常來講，曲面的一階偏導量為切向量和法向量，二階偏導量為曲率等相關信息，曲率信息主要指主曲率及主方向.三角網(wǎng)格頂點的主曲率為該點所有法曲率中的最大值與最小值，主方向則為主曲率對應法平面上的切向量，它們都是曲面重要的微分幾何特性.因此，確定三角網(wǎng)格頂點的離散微分量是網(wǎng)格重劃的基礎.

本文采用局部一般二次曲面法［5］來確定網(wǎng)格頂點的主曲率場.空間任意復雜曲面的局部形狀都可由一般二次曲面片S（u，v）=（u，v，h（u，v））來近似擬合，其中S（u，v）是一般二次曲面片在局部坐標系下的參數(shù)形式，且二次多項式h（u，v）的數(shù)學表達式為

式中：u，v，h分別為局部坐標系puvh中u，v，h軸對應的坐標值；A，B，C，D，E，F(xiàn)為二次多項式系數(shù).

對于離散網(wǎng)格的任意頂點P=｛pi｝，pi∈R3，i為網(wǎng)格頂點編號，i∈｛1，…，N｝，N為網(wǎng)格頂點數(shù).本文采用質(zhì)心面積夾角法［11］來估算頂點的法矢Npi.以Npi為局部坐標系的h軸建立一組單位正交基｛u，v｝，uv平面垂直于h軸，形成局部坐標系puvh.利用三角網(wǎng)格的拓撲關系，分別選取頂點對應的二階鄰域點集，投影至建立好的局部坐標系puvh中，作為曲面擬合的樣本點.

在獲取鄰域點的局部信息后，二次曲面片就可通過加權最小二乘法來擬合.

式中：i為當前網(wǎng)格頂點的編號；j為對應的鄰域點編號；uj，vj，hj為以點i為原點的局部坐標系puvh下j的坐標值；ωij為權值.

權值的分配按照鄰域點pj與中心點pi的距離來確定，即

式中：dpipj為鄰域點pj距頂點pi的距離，max（dpipj）為所有鄰域點距頂點pi的最大距離.二次多項式的系數(shù)可由法方程來求得.

在獲取離散頂點pi的局部最佳擬合多項式后，可用原點p（u=0，v=0）的微分特性來求解離散頂點的微分量，即通過擬合多項式的系數(shù)和微分幾何的相關定義來確定主曲率及主方向.

1.2 自由邊界的識別

自由邊界是半封閉模型的重要特征之一，它既是曲率線追蹤的終止條件，又是網(wǎng)格重劃的起點或終點.自由邊界不但可以保證網(wǎng)格重劃的完整性，而且在曲面性質(zhì)、受力分析等方面也起到了重要作用.因此，該特征需要從模型中提取出來，進而為后期的算法做準備.

根據(jù)半邊結(jié)構的特性，利用初始化的數(shù)據(jù)結(jié)構可提取模型的自由邊界.首先，遍歷模型的半邊數(shù)組，判斷該半邊的對邊是否存在；若存在，則該邊就不是自由邊界；若不存在，則該邊為自由邊界，將其標記并存儲在自由邊數(shù)組中.

在查找完所有網(wǎng)格的自由邊界后，同樣利用半邊結(jié)構的特性對自由邊數(shù)組進行排序.如圖2所示，設a為自由邊數(shù)組中的第1個元素，a的下一個半邊為b，b的對邊為d，f為d的下一個半邊，由此找到自由邊數(shù)組中a的下一個元素f，以此類推，直至回到a，最后完成當前數(shù)組的排序.每組自由邊界均為閉環(huán).

圖2 自由邊界的識別Fig.2 The detection of free boundary

2 網(wǎng)格遍歷的曲率線追蹤

曲率線分為最大曲率線及最小曲率線.由定義可知，最大曲率線上每個點的切向量均與當前點的最大曲率的曲面主方向共線；同理，最小曲率線上每個點的切向量均與當前點最小曲率的曲面主方向共線［7］.最大曲率線和最小曲率線相交成直角，構成曲面上的一張正交網(wǎng)［12］，若對曲率線兩兩求交，則連接交點就會形成曲面的重劃網(wǎng)格.基于曲率線追蹤的曲面重構算法效率高、效果好，能夠?qū)崿F(xiàn)較少的人機交互，適用范圍廣.

本文提出基于網(wǎng)格遍歷的曲率線追蹤算法，算法的核心是：沿追蹤點的主方向?qū)⒃擖c投影至二維參數(shù)平面來尋找下一個追蹤點，且所有的追蹤點均位于三角網(wǎng)格面片的頂點或邊.

曲率線追蹤的過程分為沿最大、最小主方向的追蹤.在追蹤最大主方向時，首先設定所有的網(wǎng)格頂點為曲率線段的種子點集，然后在追蹤每一組最大的曲率線段時將標記過的網(wǎng)格頂點從種子點集中剔除，最后直至種子點集為空才完成最大曲率線段的采樣.同理，重新設定所有的網(wǎng)格頂點為曲率線段的種子點集，按上述步驟完成最小曲率線段的采樣.

2.1 曲率線段追蹤起始點的采樣

曲率比值δi=Ki，max/Ki，min記錄了曲面在這一點沿不同方向的彎曲程度，其中Ki，max與Ki，min分別為點i的法曲率的最大值與最小值.顯然，各向異性區(qū)域越廣，主曲率方向越明顯，反之則亦然.在臍點處，曲面沿各個方向的彎曲程度相同，即δi為零，此時的主曲率方向為任意方向［13］.因此，應尋找種子點集中曲率比值最大的點為每組曲率線段的追蹤起始點.這里需要注意的是，追蹤起始點不能為自由邊界上的點，原因就在于追蹤過程涉及一階鄰域三角形，而自由邊界上點的一階鄰域三角形不完整，會導致無法追蹤.

2.2 曲率線段的追蹤

在確定曲率線段的追蹤起始點后，可根據(jù)種子點的曲率信息展開各向異性區(qū)域的追蹤.由于追蹤點可能為網(wǎng)格頂點或網(wǎng)格邊上的點，所以將追蹤的起始位置分為2種情況，即基于點的追蹤和基于邊的追蹤.

基于點的追蹤采用重心坐標法來進行局部參數(shù)化［14］.首先，標記追蹤點為P，Qj為其一階鄰域點，分別計算點P到Qj之間的距離dj，找到距離最大的點Q0并將其投影至過點p與法矢NP垂直的切平面上，以Q0的投影點Q′0為基準，沿著鄰域三角網(wǎng)格面片的排列順序（逆時針）將點P的一階鄰域點及一階鄰域邊按順序展開，形成如圖3a所示的分塊區(qū)域.其切平面的示意如圖3b所示，θi為原三角網(wǎng)格面片在P處的內(nèi)角，參數(shù)α滿足公式（4）：

然后，投影點P處的主方向T至切平面并標記為T′，根據(jù)T′與向量PQ′0的夾角θ′來判定T′所屬的分塊區(qū)域，這樣就可以確定下一個追蹤點所屬的三角網(wǎng)格面片.

最后，將上一步確認的三角網(wǎng)格面片投影至切平面形成△PQ′iQ′i+1，i∈（0，n），如圖4所示，求T′與該投影三角形的交點K′及K′到其所在邊的2個頂點的距離比dQ′2K′/dK′Q′1映射關系來獲取下一個追蹤點K所在三角網(wǎng)格邊上的位置.

圖3 重心坐標法的局部參數(shù)化示意Fig.3 Local parameterization using barycenter coordinates method

圖4 當前追蹤點P的鄰域三角網(wǎng)格投影示意Fig.4 The projection diagram of current tracking point on neighborhood triangular meshes

基于邊的追蹤采用基平面展開法來進行局部參數(shù)化［15］.基平面是指當前追蹤步的前一個追蹤線段所在的三角網(wǎng)格面片，即第n步追蹤的曲率線段所在的網(wǎng)格面片為基平面，則第n+1步追蹤的網(wǎng)格面片就為目標平面，具體示意如圖5所示.假設點K為當前某一邊上的追蹤點，點P為點K的前一個追蹤點，那么線段KP所在的三角網(wǎng)格就是基平面，這里記為△P0P1P2；與此同時，與點K所在邊的相鄰三角形就為目標平面，即△P1P2P3.點K的下一個追蹤點的求解方式也在切平面上展開，首先，將點P1，P2，P3及K的主方向TK依次投影至過點K的切平面上，從而獲得目標映射平面△P′1P′2P′3及投影后的主方向T′K；然后，在二維平面△P′1P′2P′3內(nèi)過點K沿T′K的方向求下一條邊上的映射追蹤點L′；最后，根據(jù)點L′在邊P′3P′2的位置來逆映射出追蹤點L在邊P3P2上的位置.

圖5 基平面展開法的局部參數(shù)化示意Fig.5 The base plane parameterization diagram using barycenter coordinates method

追蹤點微分信息的求解可分為2種情況，即若追蹤點為網(wǎng)格頂點，則其主曲率及主方向可直接獲取；若追蹤點為網(wǎng)格邊上的點，則可用線性插值法來確定.

2.3 曲率線段的分布

曲率線段的追蹤流程如圖6所示.首先，對種子點集進行點采樣，用以確定曲率線段的追蹤起始點并對其進行標記；然后，展開各向異性區(qū)域的曲率線分布并分析追蹤點位置，從而有針對性地展開基于點、邊的追蹤算法，以便確定下一個追蹤點的位置.這里需要注意，若下一個追蹤點與鄰近網(wǎng)格頂點的距離小于用戶設定的閾值（設L為當前追蹤點所在三角網(wǎng)格的邊長，本文取0.05L），則將當前追蹤點替換為鄰近的網(wǎng)格頂點.

圖6 曲率線段的追蹤流程Fig.6 Flowchart of curvature line tracking

當每組曲率線段的追蹤遇到下述2種情況時，追蹤終止.

（1）追蹤到敏感區(qū)域.所謂敏感區(qū)域，是指已被追蹤過的區(qū)域，其微觀表現(xiàn)為數(shù)據(jù)結(jié)構中的敏感性.一個三角網(wǎng)格面片有2個方向的敏感性，即最大主曲率方向敏感與最小主曲率方向敏感.如果被定義在最大主曲率方向上敏感，則當最大主曲率線段追蹤到該三角網(wǎng)格面片時子追蹤步應立即停止；同理，當最小主曲率線段追蹤到最小主曲率敏感的三角網(wǎng)格面片時，該子追蹤步也應當停止.簡單地說，每追蹤到一個三角網(wǎng)格面片時，都要對該網(wǎng)格面片進行標記.

（2）追蹤線段到達模型的自由邊界.在曲率線段追蹤完畢后，需要設定曲率線段的密度來控制網(wǎng)格分布的大小.曲率線段的密度是指2條相鄰的曲率線段間的空間距離大小，其分布將會影響網(wǎng)格重構的質(zhì)量.因此，在進行曲率線段的追蹤時應考慮它在給定位置處的理想密度［4］.本文提出一種動態(tài)控制曲率線段密度的方法，即實時鄰域點剔除法.該算法的核心是：在每條追蹤線段完畢后，將追蹤線段所在的三角網(wǎng)格m階鄰域點進行標記（m為用戶自定義，本文選取m=5），標記過的點不會成為追蹤起始點.

如圖7所示，假設頂點A為某三角網(wǎng)格模型中曲率比值最大的點，B，C，D緊隨其后；3組線段為追蹤的曲率線段.那么，當完成以點A為追蹤起始點的子追蹤步后，接下來應選取除A外的曲率比值最大點，則B應作為第2組曲率線段的追蹤起始點，但由于點B，C均處在第1組曲率追蹤線段所在網(wǎng)格面片的m階鄰域范圍內(nèi)，所以下一組曲率線段的追蹤起始點為不在范圍內(nèi)的點D.

圖7 曲率線段的密度控制示意Fig.7 The diagram of curvature line density control

按上述追蹤方法沿網(wǎng)格頂點的最大、最小主方向循環(huán)進行曲率線段的追蹤，直至形成正交的主曲率網(wǎng)格.對網(wǎng)格上所有的節(jié)點進行求交，這些交點便是重劃網(wǎng)格的采樣點，而重劃的網(wǎng)格邊則為采樣點的主曲率線段，由此來獲取逆向過程中曲面模型的重劃網(wǎng)格.

3 算法應用實例分析與比較

通過一些三角網(wǎng)格模型的實例，如管道、車門、車身翼子板等來測試曲面網(wǎng)格的重劃算法，從而實現(xiàn)自曲率估算→自由邊界特征的提取→曲率線段的追蹤→網(wǎng)格重劃等一系列過程.圖8～10為算法應用實例，圖11及12為相關算法與本文算法實例比較.

本文采用絕對值分組法來顯示曲率的估算結(jié)果，根據(jù)人眼對顏色的分辨能力，按照色帶由紅到綠再到藍，亮度由暗到亮再到暗的變化方式選出10種最具代表性的顏色來繪制曲率，相應的顏色及RGB（red，green，blue）值為：紫（0.5，0，0）、紅（1.0，0，0）、橙（1.0，0.5，0）、淺紅（1.0，0.5，1.0）、淺黃（1.0，1.0，0.5）、綠（0，1.0，0）、淺藍（0，1.0，1.0）、藍綠（0.5，0.5，1.0）、藍（0，0，1.0）、青（0，0.5，0.5）.曲率值由大到小等分為10組，對應上述10種顏色由亮到暗變化.

圖8 管道的網(wǎng)格重劃過程Fig.8 The surface remeshing of pipe

圖9 車門的網(wǎng)格重劃過程Fig.9 The surface remeshing of car door

如圖8b，9b，10b所顯示的曲率估算結(jié)果來看，在圓角、孔洞等過渡曲面處色彩鮮明，曲率值明顯較大，而在平面或曲面較為平緩區(qū)域曲率渲染色彩位于曲率值較小的范圍內(nèi)，顏色較暗，模型各區(qū)域的曲率渲染色彩分布充分反映曲面曲率值大小，因此曲率估算結(jié)果較接近真實值.如圖8c，9c，10c所示，管道的上下端的自由邊界，車門及翼子板的外圍自由邊界都能夠完全被識別，因此利用半邊數(shù)據(jù)結(jié)構查找自由邊的方法為后期曲率追蹤終止條件的設定奠定堅實的基礎.最后，如圖8a，8d，8e，8f，9a，9d，9e，9f，10a，10d，10e，10f所示，最大、最小曲率線沿曲率線的主方向展開形成規(guī)則的曲率線網(wǎng)，網(wǎng)格的大小分布合理；曲率較大位置（圓角或倒角處）的網(wǎng)格分布較密；曲面平坦處分布較疏.所以，對特征不明顯或較為分散的復雜曲面，該算法能夠提供理想的曲面網(wǎng)格重劃結(jié)果.

圖10 車身翼子板的網(wǎng)格重劃過程Fig.10 The surface remeshing of fender

如圖11所示，維納斯模型在嘴角右下方有凹陷，頭部發(fā)型復雜，脖頸曲面較為規(guī)則.相比文獻［8］算法，本文算法不僅在凹陷處實現(xiàn)較為規(guī)則的網(wǎng)格重劃，且對不同復雜程度的曲面都具有良好的適應性.

圖11 文獻［8］算法與本文算法的比較Fig.11 Comparison with reference［8］

如圖12所示，本文算法已經(jīng)達到了通過追蹤來分割離散曲面的目標，兔子網(wǎng)格密度能按照曲面的變化合理分布，兔子后背處網(wǎng)格較頭部更密集，體現(xiàn)了其肌肉的紋理變化.但是要獲得高質(zhì)量的結(jié)果，仍需要增加網(wǎng)格簡化（mesh simplification）、凸分解（convex decom-position）等后處理步驟.

圖12 文獻［9］算法與本文算法的比較Fig.12 Comparison with reference［9］

4 結(jié)語

主要研究了逆向工程中曲面網(wǎng)格的重劃算法.首先，采用局部一般二次曲面法建立了三角網(wǎng)格頂點的主曲率場，并提出了基于半邊結(jié)構的自由邊界識別法，該方法能夠完全識別出模型的自由邊界；然后，提出了一種基于網(wǎng)格遍歷的曲率線段追蹤法，該方法采用了將目標三角網(wǎng)格投影至切平面的二維參數(shù)化法，準確、高效地追蹤到了曲率線段；最后，在追蹤過程中提出了基于鄰域點剔除思想的曲率線段動態(tài)密度控制法.

本文提出的方法能夠在特征不明顯或結(jié)構復雜的網(wǎng)格曲面上進行曲率線段的追蹤，并實現(xiàn)密度適中的曲率線分布，從而最大限度減小了后期重構曲面與原始網(wǎng)格模型間的差距，為三角網(wǎng)格模型的曲面重建提供了保障.

［1］ 邢健，付大鵬，郝德成.基于逆向工程的汽輪機葉片型面CAD建模方法的研究［J］.機械設計與制造，2011（5）：223.XING Jian，F(xiàn)U Dapeng，HAO Decheng.Study of CAD modeling for turbine blade profile based on reverseengineering［J］.Machinery Design &Manufacture，2011（5）：223.

［2］ 盛忠起，蔡光起.逆向工程及曲面重建技術［J］.機械設計與制造，2001，12（6）：102.SHENG Zhongqi，CAI Guangqi.Reverse engineering and surface reconstruction［J］.Machinery Design &Manufacture，2001，12（6）：102.

［3］ Huang J B，Menq C H.Automatic data segmentation for geometric feature extraction from unorganized 3-D coordinate points［J］.IEEE Transactions on Robotics and Automation，2001，17（3）：268.

［4］ 董明曉，鄭康平，姚斌.曲面重構中點云數(shù)據(jù)的區(qū)域分割研究［J］.中國圖象圖形學報，2005，10（5）：575.DONG Mingxiao，ZHENG Kangping，YAO Bin.Research on point cloud data segmentation in surface reconstruction，2005，10（5）：575.

［5］ 劉勝蘭.逆向工程中自由曲面與規(guī)則曲面重建關鍵技術研究［D］南京：南京航空航天大學，2004.LIU Shenglan.Research on key technology inreconstruction of free-form ®ular surfaces in reverseengineering［D］.Nanjing：Nanjing University of Aeronautics and Astronautics，2004.

［6］ Allize P，Ucelli G，Gotsman C，et al.Recent advances in remeshing of surface［M］.Mathematics and Visualization.Berlin：Springer-Verlag，2008.

［7］ Alliez P，Cohen S D，Devillers O，et al.Anisotropic polygon remeshing［J］.ACM Transaction on Graphics（TOG），2003，22（3）：485.

［8］ 龐旭芳，宋展.散亂點云曲率場的計算［J］.中國科學院先進技術研究通報，2010，4（7）：39.PANG Xufang，SONG Zhan.High-throughput multi-antigen immunoassays on a microfluidicchip［J］.Bulletin of Adavanced Technology Research，2010，4（7）：39.

［9］ Marinov M，Kobbelt L.Direct anisotropic quad-dominant remeshing［C］//Proceedings of Pacific Graphics.Seoul：IEEE，2004：207-216.

［10］ Kalogerakis E，Nowrouzezahrai D，Simari P，et al.Extracting lines of curvature from noisy point clouds［J］.Computer-Aided Design，2009，41（4）：282.

［11］ 齊寶明.三角網(wǎng)格離散曲率估計和Taubin方法改進［D］.大連：大連理工大學，2008.QI Baoming.Curvatures estimation and the improvementof Taubin’s method on triangular mesh［D］.Dalian：Dalian University of Technology，2008.

［12］ 朱心雄.自由曲線曲面造型技術［M］.北京：科學技術出版社，2000.ZHU Xinxiong.Curve and curved surface model technology［M］.Beijing：Scientific and Technical Publishers，2000.

［13］ 朱為鵬，高成英，羅笑南.全局各向異性四邊形主導網(wǎng)格重建方法［J］.軟件學報，2012，23（5）：1305.ZHU Weipeng， GAO Chengying，LUO Xiaonan.Global anisotropic quad-dominant remeshing［J］.Journal of Software，2012，23（5）：1305.

［14］ Welch W，Witkin A.Free-form shape design using triangulated surfaces［C］//SIGGRAPH94 Conference Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Techniques.Orlando：SIGGRAPH，1994：247-256.

［15］ Sorkine O，Cohen D，Goldenthal R，et al.Bounded-distortion piecewise mesh parameterization ［C］//Proceedings of the Conference on Visualization’02.Washington D C：IEEE，2002：355-362.