凹包內散亂點集Delaunay四面體角度剖分算法

李世森，王熹芳

（天津大學，天津 300072）

凹包內散亂點集Delaunay四面體角度剖分算法

李世森，王熹芳

（天津大學，天津 300072）

在邵鐵政［1］三維空間散亂點集Delaunay四面體剖分算法的基礎上，提出了一種不含有除法運算（不存在被0除或喪失計算精度的情形）的通用的判定空間兩三角形內交的算法，可以實現凹包內散亂點集的Delaunay四面體剖分。該算法已經通過Fortran語言編程實現并且給出了算例。

散亂點；Delaunay規則；空間三角形內交；四面體

Biography：LI Shi?sen（1969-）,male,associate professor.

Delaunay三角剖分算法起源于1934年俄國數學家Delaunay提出的Delaunay準則［2］。在應用有限元方法解決各類問題時，Delaunay網格有嚴謹的數學證明，能夠生成形態優化的網格［3］，是網格劃分的主要方法之一。目前，Delaunay三角剖分算法廣泛地應用于計算機圖形學、航天、地質、土木工程等領域，體現了較強的實用價值。

Delaunay三角剖分算法大致可以分為三類：分而治之算法、逐點插入算法和三角網增長法［4］。雖然二維的Delaunay三角剖分算法已經取得了一定的成果，但是三維空間的Delaunay四面體剖分算法，由于三維空間外包面的復雜性，目前的研究還不夠成熟。如文獻［5-6］需要先從任意多面體中剖分出一個凸空間，得到一個新的多面體，再把剖分出來的凸空間分割成多個四面體，再重復對新的多面體進行剖分，直至剖分完畢；文獻［7-9］在對空間散亂點集進行Delaunay四面體剖分時，需要先將散亂點集的外邊界進行Delaunay三角剖分，由邊界面生成初始四面體，再通過交換或尋找幾何關系將初始四面體剖分更新為新的四面體。

本文提出了對邊界面已知的，外包面為凹的空間散亂點進行Delaunay四面體剖分的算法。通過引入對于空間2個三角形是否內交的判斷，保證了生成的四面體在凹包邊界面的里側，實現了凹包內散亂點集的Delaunay四面體剖分。

1 凹包內散亂點Delaunay四面體網格生成方法

1.1 相關概念

為了有效地解決凹包內散亂點Delaunay四面體剖分，本算法給出了空間兩三角形內交的定義。

空間兩三角形內交定義為在三維空間內，當一個三角形三條邊所圍成的區域與另一個三角形三條邊所圍成的區域的交集為一條線段時，則稱作兩三角形內交。在本算法中，當該線段與三角形的某一條邊重合時，不認為兩三角形內交。

1.2 基本思路

［1］在Delaunay準則定義［3］、球缺角定義［1］、四面體球缺角定義及特性［1］的基礎上，給出了凸包內散亂點集Delaunay四面體剖分算法，本算法在此基礎上又加入了判斷空間兩三角形內交的算法。給定一個空間散亂點集以及該散亂點集的外包面（以空間三角形面組成），首先用參考文獻［1］中的方法對散亂點集進行Delaunay四面體劃分，然后判斷生成四面體的各個面是否與已知的邊界面內交。若內交繼續尋找下一點構成四面體，若不內交，則生成一個Delaunay四面體，這樣就使得參考文獻［1］中的方法能夠處理凹包面內的散亂點。

1.3 判斷空間兩個三角形內交的數學方法

為了判斷這兩三角形是否內交，需要首先分別判斷一個三角形的三條邊所在的直線與另一三角形所在平面的位置關系：相交、平行（交點在無窮遠處相交）、在面內（有無窮多個交點）。以判斷三角形Q1Q2Q3的邊Q1Q2所在直線與平面K1的位置關系為例，有聯立方程組

如果直接采用式（3）這一表達式進行計算機編程，當t2=0時（或者t2非常的接近于0時），會造成計算機數值溢出（或者計算結果喪失精度）。即式（3）在計算機中并不是總能求出具體的數值（但是式（3）在任何情況下，在數學上都可以做為交點坐標的表達式），因此式（3）在后續的算法中一直保持表達式的形式。

求解完直線Q1Q2與平面K1的交點表達式后，需要進一步判斷交點與線段Q1Q2的位置關系。設直線Q1Q2與平面K1的交點為Q4（x7,y7,z7），將x7,y7,z7的表達式依次帶入 （x4-x7）×（x7-x5），（y4-y7）×（y7-y5），（z4-z7）×（z7-z5）中，則可以得到

通過判斷式（4）、（5）、（6）的正負來確定交點Q4與線段Q1Q2的位置關系。又因為和（x5-x4）2不影響式（4）、（5）、（6）的正負，所以只需要判斷-t1（t1-t2）（令M=-t1（t1-t2））的正負即可。當M值為正時，則Q4內分線段Q1Q2（可以說平面K1內分線段Q1Q2）；當M值為負時，則Q4外分線段Q1Q2（可以說平面K1外分線段Q1Q2）；當M值等于0時，則線段Q1Q2至少有一個端點在平面K1上。

至此，不需要求出式（3）具體數值，而采用無除法表達式M的正負，即可判斷線段與平面的位置關系。

當三角形Q1Q2Q3的某一條邊被平面K1內分時，則三角形Q1Q2Q3與平面K1的交集為一條線段（設為Q4Q5）。當三角形Q1Q2Q3的三條邊對應的-t1（t1-t2）值均為0時（則必有一條邊對應的t2值為0），則可判斷為t2值為0的線段在平面K1內，即三角形Q1Q2Q3與平面K1的交集為線段Q4Q5（即Q1Q2），如圖1所示。除了上述兩種情況之外，空間兩三角形不可能發生內交。

同理可得，三角形P1P2P3與平面K2的交集（設為線段P4P5）。兩條線段的位置關系（兩條線段必然在同一直線上）可能存在如圖2所示的5種情況。判斷兩條交線段位置關系的方法，與上文判斷交點與線段位置關系的方法類似（即通過判斷交點Q4，Q5與線段P4P5的位置關系和交點P4，P5與線段Q4Q5的位置關系來判定兩條交線段的位置關系）。

圖1 線段Q1Q2在平面內Fig.1 Line segmentQ1Q2on the plane

2 凹包內散亂點四面體剖分的算法

2.1 判斷空間兩個三角形內交的算法

（1）輸入兩個三角形的6個點的點號及坐標（xi,yi,zi），（i=1,2…6）。

（2）判斷兩個三角形是否有公共邊（或者完全相同）。若存在公共邊，則直接返回空間兩三角形不內交；否則，進行步驟（3）。

（3）計算A，B，C，D，E，F，G，H的值。

（4）依次計算各個線段對應的t1和t2的值。

（5）分別判斷各個線段對應的-t1（t1-t2）的正負。當-t1（t1-t2）值為正時，標記值為1（存儲在數組中）；當-t1（t1-t2）值為負時，標記值為-1；當-t1（t1-t2）等于0時，標記值為0。

（6）檢驗三角形Q1Q2Q3三條邊的標記值。當三角形Q1Q2Q3的三條邊，有任意一條邊標記值為1時，則進行步驟（8）；當三條邊的三個標記值均為0時，則進行步驟（7）。其他情況，直接返回空間兩三角形不內交。

圖2 兩交線段的位置關系圖Fig.2 Positional figure of two intersecting line segments

（7）檢驗三條邊對應t2的值。找出三條邊中t2值為0（或者t2非常的接近于0）的邊（該邊即為線段Q4Q5）。

（8）檢驗另一三角形P1P2P3的三條邊的標記值。方法同步驟（5），（6），（7）。

（9）檢驗線段P4P5與Q4Q5的位置關系。當R值為正時，標記值為1；當R值為負時，標記值為-1；當R值等于0時，標記值為0（標記值存儲在另一數組中）。

（10）判定空間兩三角形是否內交。當4個交點（Q4，Q5，P4，P5）標記值有兩個為1或者4個交點的標記值均為0時，返回空間兩三角形內交；否則返回空間兩三角形不內交。

2.2 算法運行框圖及算法的輸入

本算法所需的輸入文件為邊界面輸入文件和散亂點輸入文件。邊界面的剖分格式為三角形，邊界面數據文件中每一行保存一個邊界三角形的三個點的點號，該數據文件的行數即為邊界面三角形的個數。散亂點數據文件中每一行保存一個散亂點的坐標，順序為“x坐標，y坐標，z坐標”，1號點保存在第一行，2號點保存在第二行，以此類推，行數即為空間散亂點的個數。

算法運行框圖如圖3所示。

2.3 算法的具體步驟

算法具體分為以下幾個步驟：

（1）讀取邊界面三角形及空間散亂點。

（2）用一個變量q記錄已生成的四面體個數，用另一個變量p記錄已擴展的四面體個數。p和q的初始值為0。

（3）對每一個邊界面如下運算：

a）判斷該邊界面是否出現在已生成的q個四面體中。若未出現，則轉到步驟（b）；否則，處理下一個邊界面。

b）采用文獻［1］的算法，在散亂點中找到一個點與邊界面三角形生成一個四面體。

c）檢驗（b）中生成的四面體的每一個面與所有已知邊界面及已生成的四面體的各個面均不內交時，將四面體4個點的點號按順序保存在數據文件中，q值增加1；否則，從全部點集中臨時刪除此點，回到步驟（b）。

（4）p值增加1。

（5）比較p，q值。若p＞q，結束；否則，轉向步驟（6）。

（6）對第p個四面體的4個三角形進行是否可擴展的判斷：若該三角形是邊界面，或者是兩個四面體的公共面，則不可擴展；否則進行步驟（7）。

（7）對第p個四面體的每一個可擴展的三角形采用文獻［1］的算法生成四面體，并且保證該四面體的各個面不與已生成的四面體的各個面內交。將每一個生成的四面體4個點的點號按順序保存在數據文件中，q值增加1。

（8）轉向步驟（4）。

圖3 程序運行框圖Fig.3 Program operation diagram

圖4 1 000個空間散亂點的效果圖Fig.4 Sketch of 1 000 spatial scattered

3 凹包面效果圖

本文運用上述算法，對空間散亂點個數為413的空間體進行了Delaunay四面體剖分，得到如圖4所示的整體效果圖及x，y，z方向的截面圖，由于整體效果圖中不能很好的顯示內部的網格，故而又給出了x，y，z方向的截面圖。此外，本文還給出了空間散亂點個數為1 000的L型空間體效果圖，如圖4所示。可以看到，本算法可以對凹包面內散亂點進行Delaunay四面體剖分，這對空間散亂點集Delaunay四面體剖分算法的研究有一定的意義。

4 結論

（1）本文所述算法沿用了參考文獻［1］的算法，即根據四面體球缺角的特性在散亂點中尋找與之構成Delaunay四面體的點的算法。

（2）在參考文獻［1］的算法的基礎上，本算法又提出了判斷空間兩三角形內交的算法。

a）當檢驗空間兩三角形是否內交時，若兩個需要檢驗的三角形有公共邊（或重復），則不需要進行內交檢驗。

b）在判斷空間兩三角形是否內交的一般數學表達式中必然會出現除法運算，而空間兩三角形的位置可以千變萬化，也就無法保證表達式中的分母不為0。為了得到通用的算法，本文推導出了最終的數學表達式，得到了不含有分母的通用判別式（M，R），可準確判斷出任意兩個空間三角形是否相交。當最終的數學表達式中對應的分母為0時也對應了一種空間兩三角形的位置關系。

c）通過對空間兩三角形內交的算法的引入，本算法能夠對凹包內散亂點集進行Delaunay四面體剖分。

（3）給出了整個算法實現的具體步驟。

（4）應用本算法對空間散亂點個數為413的空間體進行了Delaunay四面體剖分，并得到了效果圖。

參考文獻：

［1］邵鐵政,李世森.凸包內空間散亂點集Delaunay四面體角度剖分算法［J］.水道港口,2013,146（1）：95-98.

SHAO T Z，LI S S.Delaunay Angle Algorithm of Spatial Scattered Point Set Delaunay Triangulation for Convex Hull［J］.Journal of Waterway and Harbor,2013,146（1）：95-98.

［2］Lawson C L.Software for C1 Surface Interpolation,Mathematical Software III［M］.New York：Academic Press,1977：161-194.

［3］王建華,徐強尋,張銳.任意形狀三維物體的Delaunay網格生成算法［J］.巖石力學與工程學報,2003,22（5）：717-720.

WANG J H,XU Q X,ZHANG R.Delaunay algorithm and related procedure to generate the tetrahedron mesh for an object with ar?bitrary boundary［J］.Chinese Journal of Rock Mechanics and Engineering，2003,22（5）：717-720.

［4］李麗.三維空間Delaunay三角剖分算法的研究及應用［D］.大連：大連海事大學,2010.

［5］陳一民,李超,熊玉梅.任意多面體的四面體剖分算法［J］.計算機工程與應用,2003（30）：69-93.

CHEN Y M,LI C,XIONG Y M.Algorithm of Dividing an Arbitrary Polyhedron to Tetrahedrons［J］.Computer Engineering and Applications,2003（30）：69-93.

［6］李昌領,張虹,朱海峰.一種任意多面體的四面體剖分的改進算法［J］.計算機工程與應用,2012（25）：20-38.

LI C L,ZHANG H,ZHU H F.Improved Algorithm of Dividing an Arbitrary Polyhedron into Tetrahedrons［J］.Computer Engineer?ing and Applications,2012（25）：20-38.

［7］關文革,武強,賈麗萍,等.約束數據域Delaunay四面體網格生成算法［J］.華中科技大學學報：自然科學版,2005（5）：67-69.

GUAN W G，WU Q，JIA L P，et al.Algorithm of mesh generation of Delaunay tetrahedral in constrained domain［J］.Journal of Huazhong University of Science and Technology,2005（5）：67-69.

［8］XU H,WU Q.Design&implementation of visualization for 3d sandwich geological bodies［J］.Computer Applications,2001,21（12）：56-60.

［9］WU Q,XU H.An approach to computer modeling and visualization of geological faults in 3d［J］.Computers and Geoscience,2003,29（4）：503-509.

［10］彭國倫.Fortran95程序設計［M］.北京：中國電力出版社，2010.

Delaunay angle algorithm of scattered point set delaunay triangulation for concave hull

LI Shi?sen,WANG Xi?fang
（Tianjin University，Tianjin300072，China）

Based onDelaunay Angle Algorithm of Spatial Scattered Point Set Delaunay Triangulation for Con?vex Hulldefined by SHAO Tie?zheng,a new common algorithm of judging two triangles intersection which did not contain the division（without the case of divided by zero or loss of accuracy）was proposed in this paper,and the algo?rithm could solve spatial scattered point set Delaunay triangulation for concave hull.The actual programming opera?tion of the new algorithm was also carried out by Fortran,and an example was given.

scattered points;Delaunay rules;spatial triangles intersection;tetrahedron

O 182.2

A

1005-8443（2014）02-0180-05

2013-04-02；

2013-06-20

李世森（1969-），男，河北冀縣人，副教授，主要從事港口航道與近海水流、泥沙研究。