1、 PAGE 9基于曲率球的參數曲面自適應三角化方法研究楊偉民1,2，樓天良1，林志偉2，沈洪垚21.義烏工商職業技術學院機電信息分院，義烏，322000 2.浙江大學現代制造工程研究所，杭州，310027摘要：參數曲面的三角化算法是有限元分析和數控機床路徑規劃的基礎。文中提出了一種新的基于分層遞推思想的三角化算法，直接在曲面上生成三角網格。根據曲面微分幾何定義相對保守的曲率球概念，將局部曲面近似成曲率球面，在球面上建立三角形邊和面片的誤差模型和遞推關系，計算三角形各個頂點在歐式空間以及參數空間的坐標。然后采用逐層遞推的方法，從曲線邊界開始，逐層生成三角形網格，直至滿足收斂準則。文章最后將此算法

2、用C+程序語言編程并通過NURBS曲面進行實例驗證，結果表明該算法能夠較好的對參數曲面進行三角網格劃分。關鍵詞：參數曲面；三角化；曲率球；誤差控制中圖分類號：TP391.73Parametric surface triangulation based on curvature sphereWeimin Yang1,2，Tianliang Lou1，Zhiwei Lin2，Hongyao Shen21.Yiwu Industrial & Commercial College，Yiwu，322000 2.Zhejiang University，Hangzhou，310027Abstract：Par

3、ametric surface triangulation algorithm is the basis of both the finite element method and the tool path generation in CNC machining. In this paper, a new method is proposed to directly generate triangular meshes on the surface through a layered recursive method. The concept of curvature sphere is d

4、efined. The local surface is approximated to be the spherical surface of the curvature sphere. On the curvature surface, the error models of the triangular mesh and the original surface are established and the coordinates of the triangular vertexes, both in the parametric domain and the Euclidean sp

5、ace, are calculated. Then the layered recursive method is applied to generate the triangular meshes layer by layer from the boundary of the surface until the convergent condition is met. The proposed method is implemented in C+ programing language and is tested by a NURBS surface example. Results sh

6、ow that the proposed method can deal well with the parametric surface triangulation problem.Keywords：Parametric surface; triangulation; curvature sphere; error control0 引言三角化是指對已知曲面方程的參數曲面進行三角形離散劃分。它是有限元分析和數控機床路徑規劃的基礎。三角化算法的優劣直接關系到這些計算的計算量和準確度。首先要求保證三角化離散的精度，即保證自由曲面離散后的三角面片能夠逼近原始曲面；其次要求離散具有自適應性，即生成的

7、三角面片的疏密程度能夠隨著曲面曲率變化而變化，在曲面曲率大的地方生成的三角面片致密一些，而在曲面曲率小的地方生成的三角面片相對稀疏一些。由于路徑規劃中干涉檢測算法需要遍歷曲面中所有三角面片，所以在保證精度的前提下，曲面三角面片數量越少越好。對參數曲面而言，最簡單、最直觀的三角化方法是映射法。這種方法先在參數空間實現平面三角化，如Delaunay平面三角剖分 ADDIN EN.CITE ADDIN EN.CITE.DATA HYPERLINK l _ENREF_1 o Cignoni, 1998 #252 1, HYPERLINK l _ENREF_2 o 余杰, 2010 #9 2，然后將得到

8、的三角形映射到歐式空間。這種方法由于參數空間和歐式空間的映射沒有空間形態的必然聯系，得到的三角面片沒有自適應性，難以評價三角化精度 ADDIN EN.CITE 王青2003131117王青程艷旗林海鮑虎軍浙江大學國家重點實驗室,浙江大學國家重點實驗室,浙江大學國家重點實驗室,浙江大學國家重點實驗室 浙江杭州310027 ,浙江杭州310027 ,浙江杭州310027 ,浙江杭州310027隱函數曲面的實時采樣與三角化技術浙江大學學報(理學版)浙江大學學報(理學版)621-62706隱函數曲面多邊形化曲面離散20031008-949733-1246/NCnki HYPERLINK l _ENRE

9、F_3 o 王青, 2003 #1 3。徐松等在傳統映射法的基礎上，在平面三角化算法中環邊統一處理，采用自適應三角網格加密法，有效地處理帶有特征約束條件的任意曲面三角剖分問題 ADDIN EN.CITE 徐松2000242217徐松王劍英中國科學技術大學計算機系!合肥230027,中國科學技術大學計算機系!合肥230027曲面的自適應三角網格剖分計算機輔助設計與圖形學學報計算機輔助設計與圖形學學報267-27104三角剖分特征約束自適應網格加密曲面200011-2925/TP11-2925/TPCnki HYPERLINK l _ENREF_4 o 徐松, 2000 #2 4。關振群等提出一種

10、基于映射法的復雜三維組合曲面的有限元網格全自動生成方法，通過引入虛邊界解決了閉合曲面在參數域中邊界不完整的問題，并通過調節虛邊界提高了復雜組合曲面網格生成的質量 ADDIN EN.CITE 關振群2003353317關振群隋曉峰顧元憲李云鵬大連理工大學工業裝備結構分析國家重點實驗室,大連理工大學工業裝備結構分析國家重點實驗室,大連理工大學工業裝備結構分析國家重點實驗室,大連理工大學工業裝備結構分析國家重點實驗室 大連116024 ,大連116024 ,大連116024 ,大連116024復雜三維組合曲面的有限元網格生成方法計算力學學報計算力學學報409-41604網格剖分閉合曲面虛邊界組合曲面

11、20031007-470821-1373/O3Cnki HYPERLINK l _ENREF_5 o 關振群, 2003 #3 5。吳寶海等采用改進波前法對曲面進行三角化細分，從曲面兩條邊界向曲面中心同時推進，避免了常規波前法中由于曲面角點的不良形態導致網格規劃的失敗以及生成低質量的網格 ADDIN EN.CITE 吳寶海2005464417吳寶海王尚錦西安交通大學能源與動力工程學院,西安交通大學能源與動力工程學院 西安710049 ,西安710049參數曲面網格生成的改進波前法計算機輔助設計與圖形學學報計算機輔助設計與圖形學學報1686-169008參數曲面三角剖分波前法20051003-9

12、77511-2925/TPCnki HYPERLINK l _ENREF_6 o 吳寶海, 2005 #4 6。梁義等結合映射法和前沿推進技術，通過曲率計算自動識別曲率特征并創建密度場，對組合參數曲面生成幾何自適應的網格 ADDIN EN.CITE 梁義2010676617梁義,陳建軍,陳立崗,鄭耀,幾何自適應參數曲面網格生成計算機輔助設計與圖形學學報計算機輔助設計與圖形學學報327-33520102曲面網格生成 自適應 Delaunay三角化 黎曼度量 中軸20101003-977597390XVIP HYPERLINK l _ENREF_7 o 梁義, 2010 #6 7。李李等利用均勻采

13、樣點來估算曲面曲率，針對裁剪曲面進行三角剖分，提高剖分效率和精度 ADDIN EN.CITE 李李2006585517李李 王亞平 裁剪曲面自適應三角化剖分計算機應用計算機應用2006年S1期TP391.41針對CAD/CAM領域中裁剪曲面的三角剖分問題,從提高效率和滿足剖分精度的角度出發,提出了一種較為實用的自適應離散方法,與傳統的自適應方法相比,三角片數量有所減少,速度有明顯的提高。2006 HYPERLINK l _ENREF_8 o 李李, 2006 #5 8。王飛等基于B樣條曲面的適應性分割，將背齊爾曲面片的三角化處理擴展到整張曲面，消除相鄰面間邊界處可能出現的裂縫 ADDIN EN

14、.CITE 王飛1997797717王飛,董國輝,基于b樣條曲面適應性分割的三角化算法北京郵電大學學報北京郵電大學學報95-98202計算機圖形學 B樣條曲面 三角化算法19971007-532191520AVIP HYPERLINK l _ENREF_9 o 王飛, 1997 #7 9。Shimada等則提出用BMG（Bubbul Mesh Generation）對曲面進行三角化 ADDIN EN.CITE Shimada19988108817Shimada, KenjiGossard, David C.Automatic triangular mesh generation of trim

15、med parametric surfaces for finite element analysisComputer Aided Geometric DesignComputer Aided Geometric Design199-222153TriangulationDelaunay triangulationVoronoi polygonsSphere packingPhysically based modelBubble meshingDynamic simulationFinite element analysis19980167-8396/science/article/pii/S

16、016783969700037X/10.1016/S0167-8396(97)00037-X HYPERLINK l _ENREF_10 o Shimada, 1998 #8 10。在本文中，我們結合曲面的局部微分幾何，提出曲率球概念，在局部范圍內，用曲率球近似曲面，在曲率圓上直接近似生成三角面片，并利用分層遞推思想，將生成的的三角片面覆蓋至整張曲面。1 參數曲面的微分幾何對于給定的參數曲面，只要知道參數域內的任意一點坐標值，就可以根據相應的公式計算出曲面在歐式空間內對應的點，以及曲面在該點的一階、二階偏導、和單位法向量。曲面在點處的法曲率可以表示為： （1）其中：，為曲面第一基本型；，為曲面

17、第二基本型。根據式(1)可以得到曲面在該點的主曲率： （2）其中：；。圖1 曲面的微觀幾何特性不妨設，根據、在數軸上取值正負不同，可以把曲面的微觀幾何特性分成5種情況，如圖1所示。圖1b所示為橢圓拋物面，對應圖1a中的第1和第5兩種情況，即、同時為正或同時為負。當、同時為正時，橢圓拋物面開口朝上；當、同時為負時，橢圓拋物面開口朝下。圖1c所示為拋物面，對應圖1a中的第2和第4兩種情況，即、兩者中其中有一個為零。當為零時，拋物面開口朝下；當為零時，拋物面開口朝下。圖1d所示為猴鞍面，對應圖1a中的第3中情況，即為正，為負。為了簡化討論和計算，引入曲率球的概念。曲率球的定義如下：設、為曲面上處的兩

18、個主曲率，取，把球心在過的法線上，半徑為的球稱為曲面在點的曲率球。由曲率球的定義可知，曲率球是以曲面上某點的極大曲率為半徑建立的球面，因此在將曲面附近的點轉移到曲率球上考慮，從誤差角度上講，是偏安全的。2 單個三角面片生成及誤差控制上節定義了曲率球的概念，本文涉及到的誤差都是將曲面局部近似到曲率球面上分析的。本節詳述在給定的精度條件下曲面上單個三角面片的生成。2.1 曲線上弦生成在考慮在曲面上獲得一條符合精度要求的三角形邊之前，先考慮在曲線上獲得一條符合精度要求的弦。對空間參數曲線及其曲線上的一個已知點，在曲線上點附近求點，使得弦和原始曲線段之間的弓高誤差為不大于給定的弦容差。近似地，本文在點

19、的曲率圓上考慮問題。如圖2所示，設曲線上弦長度為，近似用圓上的弦代替。在圓中，為中點，連接交圓于點，則線段就是弦的弓高，也即原始曲線上弦的弓高。設曲線上點曲率圓半徑為，在中，易得：圖2 在曲線上得到一條符合精度要求的弦 （3）對式（3）忽略2次小量，得：。由上可知，只要保證，就能保證弦的精度。臨界情況下，弦的長度可以表示為： （4）下面根據式（4）中得到的值反向近似求解點參數值。設一質點沿著曲線運動，在某點速率為，方向為曲線切線方向。曲線的參數是時間的變量，寫成，隨著時間遞增，因此。由于： （5）所以： （6）令，它們分別代表，點參數值。對在處一階泰勒展開并將（6）式代入得： （7）為質點走過

20、的路程，近似認為該路程等于弦長度，即：。將式（4）代入（7）式，得： （8）根據式（8）即可求出點的參數值，在曲線上的位置可完全確定。運用上述方法可以求得任意曲線上一條符合精度要求的弦。2.2 曲面上弦生成下面考慮在曲面上獲得一條符合精度要求的弦，即三角形的一條邊。設是曲面上一點，其參數坐標為；球為點處的曲率球；是曲面上領域內的任意一點。當趨近，即時，近似認為點落在球的球面上；那么曲面上弦的弓高就是弦的兩個端點都落在球上時的最大弓高，即在球最大圓上的弓高。圖3 在曲面上得到一條符合精度要求的弦如圖3所示，考慮過點的、向曲線，由式（8）可在、向曲線上分別找到符合弓高誤差的弦和弦： （9） （10

21、）式（9）、（10）中為曲率球的半徑，為弦、容許誤差，、分別為曲面對、的一階偏導。式（9）、（10）中用曲率球半徑求出的、值是偏安全的。如圖3所示，分析點的誤差，點坐標可以表示為。令為，的夾角，記弦長度為，、長度分別為、，則由圖3可知： （11）式（11）中， 長度隨著變化而變化。顯然，在曲率球面上，當最長時，弓高誤差最大。因此在最壞的情況下，即時，此時。而由式（4）可知，因此可得曲面上弦的最大弓高誤差為： （12）要保證弦的弓高誤差小于給定弦容差，即，只需保證： （13）即在計算式（9）、（10）中的、時的容許誤差要小于給定弦容差的四分之一。2.3 曲面上三角形生成本文中提出的生成三角形面片

22、的算法是遞推的，因此認為的一條邊的兩個端點在三維空間的坐標以及在參數域內值已經在前一步中計算出來。這里要考慮的問題是如何根據這條已知邊在曲面上找到第三個點，使得點分別和，連起來之后得到的盡量能夠符合精度要求。 如圖4所示，用2.1中的方法，對點進行一次運算，在參數域內得到點坐標；同理對點進行運算，在參數域內得到點坐標，這里用是讓盡量靠近，從而提高精度。取中點，其在參數域內坐標為：由此得到的就是所要求的三角形，如圖4。下面對進行誤差分析。在微觀情況下，由于三邊是在同一點鄰域內產生的弦，在這點鄰域內曲面的曲率變化很小，近似認為這三邊所在的曲率球的半徑相同。由于三邊弓高一樣，所以的三條邊長相同，認為

23、是等邊三角形。圖4 在曲面上得到一個符合精度要求的三角形的三個頂點都落在曲面上，近似情況下以三角形其中一個頂點的曲率球代替附近曲面。如圖5所示，垂直于所在平面交球面于，則所在平面到曲面的弓高就近似為長度。圖5三角面片在曲率球上的近似弓高圖6 三角形面片弓高和三角形邊弓高的關系如圖6所示，圓弧為球大圓上的圓弧，球半徑為，即為在2.2中弦最大弓高誤差。則在直角中： （14）在直角中，： （15）在直角中： （16）令，則： （17）忽略平方項，得： （18）所以要想把的弓高誤差限定在，只要在2.2中保證： （19）綜合式（13）、（19），若要保證三角面片生成精度，只需在計算式（9）、（10）中的

24、、時的容許誤差要小于給定弦容差的四分之一。3 三角網格生成第2節中介紹了曲面上單個三角面片的生成以及誤差控制方法，下面闡述曲面三角網格生成。三角網格生成時遵循“點對點，邊對邊”原則，即每個三角形的頂點必須對應另外三角形的頂點，每個三角形的邊對應另外三角形的邊；不允許出現某個三角形的頂點對應另外三角形的邊，否則生成的三角網格會出現裂縫。本文提出一種分層遞推生成三角面片的方法，直接在參數曲面上生成三角網格，生成的三角網格無裂縫而且符合精度要求。對一張參數曲面，其中，如圖7所示，其邊界線分別是線，線，線，線。按照本文提出的分層遞推方法對曲面逐層進行三角網格劃分，假設第層劃分結束后形成的的圖形如圖7所

25、示，其中左邊有三角網格的是已經進行三角劃分過的曲面，右邊的是沒有進行三角劃分的曲面，它們之間的分界折線是由當前一層三角形的邊首尾依次連接得到的，記這一層折線為第層折線。下面對第層折線進行輪廓更行，同時在曲面上生成三角形，得到第層三角形面片。圖7 第k層三角面片生成圖8 折線段輪廓更新分析第層折線段、，首先考慮，夾角。當時，則直接連接得到，同時刪除、，成為新的，成為，下面的線段序號依次減1。下一個分析對象為。當時，則對線段和夾角考慮。當時，則直接連接得到，同時刪除、，成為新的，下面的線段序號都不變。下一個分析對象為。當，那么對根據兩端點，按照2.3節中方法計算出曲面上一點，得到的，同時刪除，成為

26、新的，成為新的，成為新的。下一個分析對象是。圖8所示的情況由于得到的新三角形面片都在折線右側，即在曲面未三角化區域，此時不存在問題。但當得到的新三角形面片在折線左側，也就是在曲面已經三角化的區域，此時新三角形會和原來的三角形重疊起來，這種情況應該避免。事實上，圖8中的夾角是根據向量，得到的： （20）式（8）無法計算當線段和的夾角大于180時的情況。而新三角形出現在左邊正是因為沒有檢測到和夾角大于180的情況。為此引入凹角以及凸角，直觀的，相對于右邊還沒有三角化了的區域，當角是圖8中的情況時，認為是凹角，它是凹向未三角化區域的；當是圖9的情況時，認為是凸角，它是凸向未三角化區域的。圖9 凸角定

27、義以上根據觀察給出了凹凸角的簡單定義，下面給出凹凸角精確的數學定義。點位于曲面上，如圖8所示，點法向量可由得到。令向量，理論上當生成的三角面片無窮小的時候，都是位于的切平面上的，此時有和方向完全一致。而實際上因為生成的三角形不可能無窮小，和方向不可能真正完全一致。人為規定一個，向量方向一致的判據。設向量和向量的夾角，人為規定當時，和同向；當時，和反向。而當和同向時，是凸角；當和反向時，是凹角。所以在前述的算法中，加上對的判斷，如果是凸角，則跳出直接去判斷和的夾角是凹角還是凸角。如果是凹角，則進入判斷的情況。3.1 邊界處理以上討論中沒有涉及到曲面邊界，下面針對邊界進行討論。對第層折線的第一條線

28、段，它的起點所在的曲面邊界可能有3種情況。如圖10所示，分別是在線上，在線上，在線上。這三種情況可以簡單地從點的參數坐標判斷出來。分析在線上這種情況，其他也是類似的。在圖10中，為線和線在的切線的夾角。當時，按照2.1節中方法在線上得到一段符合精度要求的弦,得到的三角形作為一個新的三角形，同時刪除，將變為，其他都保持不變。而當時，則在考慮和的所成角的凹凸情況。如果是凸角則直接根據兩端點、按2.3節中方法形成新三角形，同時處理；如果和的所成角是凹角，計算和的角度：如果該角度大于90則按照2.3方法根據、形成新三角形，同時處理，否則連接、以及端點得到新的三角形，同時處理。對層折線段的最后一條線處理

29、和第一條線是類似的。也分為三種情況，的終點分別和線，線，線相交；處理也和類似，都是遵循先和前一條線段求角，然后再考慮分別和曲面界面的切線求角，這里不再贅述。圖10 網格生成邊界處理 圖11 網格生成的收斂條件3.2 網格的起始和收斂遞推法對曲面進行三角化處理依賴于第一層折線段的生成。本文中第一層折線段取自曲面邊界，用2.1節中方法在邊界曲線上連續取出一系列折線段。下面分析網格生成的收斂條件。可以想象當網格生成快要結束的時候，對應的那一層折線段會碰到曲面的邊界，如圖11所示。每一層折線段的起點和終點都落在曲面的邊界上，在每一層折線段對每條線段生成三角形之前，先對這條的終點進行判斷，如果的終點在邊

30、界上，那么的起點也必然在曲面邊界上，跳過和，下一個要分析的是；如果的終點不在邊界上，則照常處理。對這樣一層處理結束之后對新得到逐個掃描，判斷終點是否落在邊界上，并記錄終點落在邊界上的個數，記為。當滿足： （21）則可以判斷曲面三角網格生成結束。最后再針對生成最后一層三角形，結束網格生成。4 實例演示本文提出的算法用C+語言編程，以下用一張NURBS曲面驗證本文算法。該曲面的控制頂點信息如下：節點矢量分別為：權值系數為： 。三角化結果如圖12所示。由圖12可知，本文算法能夠較好地對參數曲面進行三角化，在曲面曲率小的地方，能夠生成較疏的三角網格；而在曲率大的地方，生成相對密的三角網格。圖12 NURBS曲面三角化實例5 總結三角化算法是有限元分析和數控機床路徑規劃的基礎。本文提出了一種新的針對自由曲面三角化的算法。定義了曲率球的概念，將復雜的曲面誤差模型轉化到曲率球上考慮，大大簡化了誤差計算。先后在曲面上獲得符合精度要求的弦和三角形面片。在此基礎上，利用遞推算法，從曲面邊界開始，逐層遞推生成三角網格，并確定了三角網格生成的收斂條件。本文算法經過實例驗證，表明算法確實有效。參考文獻 ADDIN EN.REFLIST 1.Cignoni P, Montani C, Scopigno R (1998) DeWall: A fast divide and c