李世森，王熹芳

（天津大學(xué)，天津 300072）

凹包內(nèi)散亂點(diǎn)集Delaunay四面體角度剖分算法

李世森，王熹芳

（天津大學(xué)，天津 300072）

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

散亂點(diǎn)；Delaunay規(guī)則；空間三角形內(nèi)交；四面體

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

Delaunay三角剖分算法起源于1934年俄國數(shù)學(xué)家Delaunay提出的Delaunay準(zhǔn)則［2］。在應(yīng)用有限元方法解決各類問題時(shí)，Delaunay網(wǎng)格有嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)證明，能夠生成形態(tài)優(yōu)化的網(wǎng)格［3］，是網(wǎng)格劃分的主要方法之一。目前，Delaunay三角剖分算法廣泛地應(yīng)用于計(jì)算機(jī)圖形學(xué)、航天、地質(zhì)、土木工程等領(lǐng)域，體現(xiàn)了較強(qiáng)的實(shí)用價(jià)值。

Delaunay三角剖分算法大致可以分為三類：分而治之算法、逐點(diǎn)插入算法和三角網(wǎng)增長法［4］。雖然二維的Delaunay三角剖分算法已經(jīng)取得了一定的成果，但是三維空間的Delaunay四面體剖分算法，由于三維空間外包面的復(fù)雜性，目前的研究還不夠成熟。如文獻(xiàn)［5-6］需要先從任意多面體中剖分出一個(gè)凸空間，得到一個(gè)新的多面體，再把剖分出來的凸空間分割成多個(gè)四面體，再重復(fù)對新的多面體進(jìn)行剖分，直至剖分完畢；文獻(xiàn)［7-9］在對空間散亂點(diǎn)集進(jìn)行Delaunay四面體剖分時(shí)，需要先將散亂點(diǎn)集的外邊界進(jìn)行Delaunay三角剖分，由邊界面生成初始四面體，再通過交換或?qū)ふ規(guī)缀侮P(guān)系將初始四面體剖分更新為新的四面體。

本文提出了對邊界面已知的，外包面為凹的空間散亂點(diǎn)進(jìn)行Delaunay四面體剖分的算法。通過引入對于空間2個(gè)三角形是否內(nèi)交的判斷，保證了生成的四面體在凹包邊界面的里側(cè)，實(shí)現(xiàn)了凹包內(nèi)散亂點(diǎn)集的Delaunay四面體剖分。

1 凹包內(nèi)散亂點(diǎn)Delaunay四面體網(wǎng)格生成方法

1.1 相關(guān)概念

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

空間兩三角形內(nèi)交定義為在三維空間內(nèi)，當(dāng)一個(gè)三角形三條邊所圍成的區(qū)域與另一個(gè)三角形三條邊所圍成的區(qū)域的交集為一條線段時(shí)，則稱作兩三角形內(nèi)交。在本算法中，當(dāng)該線段與三角形的某一條邊重合時(shí)，不認(rèn)為兩三角形內(nèi)交。

1.2 基本思路

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

1.3 判斷空間兩個(gè)三角形內(nèi)交的數(shù)學(xué)方法

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

如果直接采用式（3）這一表達(dá)式進(jìn)行計(jì)算機(jī)編程，當(dāng)t2=0時(shí)（或者t2非常的接近于0時(shí)），會(huì)造成計(jì)算機(jī)數(shù)值溢出（或者計(jì)算結(jié)果喪失精度）。即式（3）在計(jì)算機(jī)中并不是總能求出具體的數(shù)值（但是式（3）在任何情況下，在數(shù)學(xué)上都可以做為交點(diǎn)坐標(biāo)的表達(dá)式），因此式（3）在后續(xù)的算法中一直保持表達(dá)式的形式。

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

通過判斷式（4）、（5）、（6）的正負(fù)來確定交點(diǎn)Q4與線段Q1Q2的位置關(guān)系。又因?yàn)楹停▁5-x4）2不影響式（4）、（5）、（6）的正負(fù)，所以只需要判斷-t1（t1-t2）（令M=-t1（t1-t2））的正負(fù)即可。當(dāng)M值為正時(shí)，則Q4內(nèi)分線段Q1Q2（可以說平面K1內(nèi)分線段Q1Q2）；當(dāng)M值為負(fù)時(shí)，則Q4外分線段Q1Q2（可以說平面K1外分線段Q1Q2）；當(dāng)M值等于0時(shí)，則線段Q1Q2至少有一個(gè)端點(diǎn)在平面K1上。

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

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

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

圖1 線段Q1Q2在平面內(nèi)Fig.1 Line segmentQ1Q2on the plane

2 凹包內(nèi)散亂點(diǎn)四面體剖分的算法

2.1 判斷空間兩個(gè)三角形內(nèi)交的算法

（1）輸入兩個(gè)三角形的6個(gè)點(diǎn)的點(diǎn)號及坐標(biāo)（xi,yi,zi），（i=1,2…6）。

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

（3）計(jì)算A，B，C，D，E，F(xiàn)，G，H的值。

（4）依次計(jì)算各個(gè)線段對應(yīng)的t1和t2的值。

（5）分別判斷各個(gè)線段對應(yīng)的-t1（t1-t2）的正負(fù)。當(dāng)-t1（t1-t2）值為正時(shí)，標(biāo)記值為1（存儲(chǔ)在數(shù)組中）；當(dāng)-t1（t1-t2）值為負(fù)時(shí)，標(biāo)記值為-1；當(dāng)-t1（t1-t2）等于0時(shí)，標(biāo)記值為0。

（6）檢驗(yàn)三角形Q1Q2Q3三條邊的標(biāo)記值。當(dāng)三角形Q1Q2Q3的三條邊，有任意一條邊標(biāo)記值為1時(shí)，則進(jìn)行步驟（8）；當(dāng)三條邊的三個(gè)標(biāo)記值均為0時(shí)，則進(jìn)行步驟（7）。其他情況，直接返回空間兩三角形不內(nèi)交。

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

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

（8）檢驗(yàn)另一三角形P1P2P3的三條邊的標(biāo)記值。方法同步驟（5），（6），（7）。

（9）檢驗(yàn)線段P4P5與Q4Q5的位置關(guān)系。當(dāng)R值為正時(shí)，標(biāo)記值為1；當(dāng)R值為負(fù)時(shí)，標(biāo)記值為-1；當(dāng)R值等于0時(shí)，標(biāo)記值為0（標(biāo)記值存儲(chǔ)在另一數(shù)組中）。

（10）判定空間兩三角形是否內(nèi)交。當(dāng)4個(gè)交點(diǎn)（Q4，Q5，P4，P5）標(biāo)記值有兩個(gè)為1或者4個(gè)交點(diǎn)的標(biāo)記值均為0時(shí)，返回空間兩三角形內(nèi)交；否則返回空間兩三角形不內(nèi)交。

2.2 算法運(yùn)行框圖及算法的輸入

本算法所需的輸入文件為邊界面輸入文件和散亂點(diǎn)輸入文件。邊界面的剖分格式為三角形，邊界面數(shù)據(jù)文件中每一行保存一個(gè)邊界三角形的三個(gè)點(diǎn)的點(diǎn)號，該數(shù)據(jù)文件的行數(shù)即為邊界面三角形的個(gè)數(shù)。散亂點(diǎn)數(shù)據(jù)文件中每一行保存一個(gè)散亂點(diǎn)的坐標(biāo)，順序?yàn)椤皒坐標(biāo)，y坐標(biāo)，z坐標(biāo)”，1號點(diǎn)保存在第一行，2號點(diǎn)保存在第二行，以此類推，行數(shù)即為空間散亂點(diǎn)的個(gè)數(shù)。

算法運(yùn)行框圖如圖3所示。

2.3 算法的具體步驟

算法具體分為以下幾個(gè)步驟：

（1）讀取邊界面三角形及空間散亂點(diǎn)。

（2）用一個(gè)變量q記錄已生成的四面體個(gè)數(shù)，用另一個(gè)變量p記錄已擴(kuò)展的四面體個(gè)數(shù)。p和q的初始值為0。

（3）對每一個(gè)邊界面如下運(yùn)算：

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

b）采用文獻(xiàn)［1］的算法，在散亂點(diǎn)中找到一個(gè)點(diǎn)與邊界面三角形生成一個(gè)四面體。

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

（4）p值增加1。

（5）比較p，q值。若p＞q，結(jié)束；否則，轉(zhuǎn)向步驟（6）。

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

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

（8）轉(zhuǎn)向步驟（4）。

圖3 程序運(yùn)行框圖Fig.3 Program operation diagram

圖4 1 000個(gè)空間散亂點(diǎn)的效果圖Fig.4 Sketch of 1 000 spatial scattered

3 凹包面效果圖

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

4 結(jié)論

（1）本文所述算法沿用了參考文獻(xiàn)［1］的算法，即根據(jù)四面體球缺角的特性在散亂點(diǎn)中尋找與之構(gòu)成Delaunay四面體的點(diǎn)的算法。

（2）在參考文獻(xiàn)［1］的算法的基礎(chǔ)上，本算法又提出了判斷空間兩三角形內(nèi)交的算法。

a）當(dāng)檢驗(yàn)空間兩三角形是否內(nèi)交時(shí)，若兩個(gè)需要檢驗(yàn)的三角形有公共邊（或重復(fù)），則不需要進(jìn)行內(nèi)交檢驗(yàn)。

b）在判斷空間兩三角形是否內(nèi)交的一般數(shù)學(xué)表達(dá)式中必然會(huì)出現(xiàn)除法運(yùn)算，而空間兩三角形的位置可以千變?nèi)f化，也就無法保證表達(dá)式中的分母不為0。為了得到通用的算法，本文推導(dǎo)出了最終的數(shù)學(xué)表達(dá)式，得到了不含有分母的通用判別式（M，R），可準(zhǔn)確判斷出任意兩個(gè)空間三角形是否相交。當(dāng)最終的數(shù)學(xué)表達(dá)式中對應(yīng)的分母為0時(shí)也對應(yīng)了一種空間兩三角形的位置關(guān)系。

c）通過對空間兩三角形內(nèi)交的算法的引入，本算法能夠?qū)Π及鼉?nèi)散亂點(diǎn)集進(jìn)行Delaunay四面體剖分。

（3）給出了整個(gè)算法實(shí)現(xiàn)的具體步驟。

（4）應(yīng)用本算法對空間散亂點(diǎn)個(gè)數(shù)為413的空間體進(jìn)行了Delaunay四面體剖分，并得到了效果圖。

參考文獻(xiàn)：

［1］邵鐵政,李世森.凸包內(nèi)空間散亂點(diǎn)集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］王建華,徐強(qiáng)尋,張銳.任意形狀三維物體的Delaunay網(wǎng)格生成算法［J］.巖石力學(xué)與工程學(xué)報(bào),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三角剖分算法的研究及應(yīng)用［D］.大連：大連海事大學(xué),2010.

［5］陳一民,李超,熊玉梅.任意多面體的四面體剖分算法［J］.計(jì)算機(jī)工程與應(yīng)用,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］李昌領(lǐng),張虹,朱海峰.一種任意多面體的四面體剖分的改進(jìn)算法［J］.計(jì)算機(jī)工程與應(yīng)用,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］關(guān)文革,武強(qiáng),賈麗萍,等.約束數(shù)據(jù)域Delaunay四面體網(wǎng)格生成算法［J］.華中科技大學(xué)學(xué)報(bào)：自然科學(xué)版,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程序設(shè)計(jì)［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-），男，河北冀縣人，副教授，主要從事港口航道與近海水流、泥沙研究。