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

邵鐵政，李世森

（天津大學，天津300072）

隨著計算機能力的不斷強大，三維空間內的有限元計算被越來越廣泛的使用，因此三維空間內的網格劃分也隨之不斷地被研究和改進。三維問題所采用的網格劃分單元一般是四面體、六面體。而四面體網格具有靈活性較高及能夠更好的逼近邊界的特點。四面體網格生成的算法已經有許多重要的進展，如局部變換法、Delaunay 算法、四叉樹/八叉樹算法、單元生長法、網格前沿法等，其中以Delaunay 算法應用最廣，因為其具有嚴謹的數學理論證明作為基礎，能夠由初始散亂點生成形態優化的網格，故稱為Delaunay 網格［1］。

相對于二維領域而言，目前Delaunay 算法在三維領域內的研究還不夠成熟，因為三維空間內，點、邊、面的管理及單元之間的關系更加復雜，單元重疊的可能性較大，而且不易檢驗。經研究大部分Delaunay 算法［1-3］，散亂點集的外邊界進行Delaunay 三角剖分,而且需要先在散亂點內部形成一個初始的四面體剖分，然后通過交換或者尋找幾何關系將內部四面體的剖分更新為新的Delaunay 四面體剖分，這樣使得網格劃分過程的前期工作量較大。本文適用于外包面為凸的散亂點四面體剖分，不需要內部初始四面體劃分，提高了四面體網格劃分的效率。

1 三維散亂點Delaunay 四面體網格生成方法

1.1 相關概念

Delaunay 算法源于1934 年B.Delaunay 提出的Delaunay 準則［4］。

三維Delaunay 準則的定義如下：

（1）定義1:對于三維空間內由一堆散點組成的四面體網格，若每一個四面體的外接球均不包含除此四面體頂點以外的網格點，則此四面體稱為Delaunay 四面體，由這些四面體形成的網格稱為Delaunay 網格。在已知大量初始散點的情況下，若不存在多點共球的情況，則由這些散點形成的Delaunay 網格是唯一的［5-6］。

為了方便尋找Delaunay 四面體，本文引入了一個新的角度概念——球缺角。由二維Delaunay 三角形圓心角判定方法推出三維空間下的部分相關定義如下。

為了闡述方便本文引入了平面域內、域外定義：

（2）定義2:如圖1-a 以平面ABC 為例，右手展開拇指與其余四指垂直，四指沿A-B-C 方向握拳，拇指所指方向的空間，即平面ABC 上部O 點所在空間為域內，反方向的空間及平面ABC 則為域外。

四面體球缺角定義如下：

（3）定義3：如圖1-a，ΔABC 為在三維空間內一個三角形，D 為在空間內與之構成四面體的點，O 為A、B、C、D 四點的外接球球心，半徑為R。OA 為O 與ΔABC 任意頂點形成向量，則OA與ΔABC 的負法向量n的夾角就是點D 對ΔABC 的球缺角δ(0＜δ＜π)。

（4）四面體球缺角的特性1:由以上定義易知D 對ΔABC 的球缺角δ 唯一。

（5）四面體球缺角的特性2:在空間內當所有散亂點位于ΔABC 域內，一點與已知三角形構成Delaunay 四面體的條件是：這個點對ΔABC 的球缺角最大。

1.2 相關證明

（1）特性1 證明:根據定義3 可知，四面體球缺角即為球心與三角形所構成的圓錐的圓錐角的一半，所以得證。

（2）特性2 證明:圖1-b 為圖1-a 的正視圖，球O 過A、B、C、D 四點，r 為已知ΔABC 外接圓半徑，R 為球O 的半徑，h 為O 到ΔABC 的距離，H 為ΔABC 域內球缺高，則球缺的體積

若D-ABC 為Delaunay 四面體，則根據[定義1]只需證明球缺內不包含任何點。若證明[特性2]，只需證明球內任意點D1對ΔABC 的球缺角大于D 對ΔABC 的球缺角即可。

∵h=r×cotδ，R=r/sinδ，H=h+R

∴球缺體積V（δ）=（r/sinδ+r×cotδ）2×[3×（r/sinδ）-（r×cotδ+r/sinδ）] ×π/3

=（-cos3δ+3×cosδ+2）×πr3/（3×sin3δ）

由以上推導公式可知:

對于已經固定的ΔABC 外接圓半徑r 不變，即球缺體積為以δ 為自變量的函數，則可求出并化簡得到V-δ 函數的一階導數

V′（δ）=-（1+cosδ）2×πr3/sin4δ

∴由于(0＜δ＜π)，V′＜0

得證V 與δ 反相關

如圖1-b 易知在域內，D1點對ΔABC 的球缺體積必然小于點D 對ΔABC 的球缺體積。即D1點對ΔABC球冠角δ 更大。

當域內球缺小于半球時，δ 的范圍是（π/2，π），即cotδ 為負，證明過程相同。

所以［特性2］得證。

2 編程實現散亂點四面體劃分的步驟

2.1 散亂點數據讀入

數據文件要求：散亂點按統一格式存入inputA.txt 文件中，順序為“初始點號，x 坐標，y 坐標，z 坐標”。

2.2 程序運行思維圖

如圖2 所示，程序運行分為以下幾個步驟：

（1）三維空間點集獲取：讀取輸入文件中的所有點的坐標，讀取內容包括“初始點號，x 坐標，y 坐標，z 坐標”。

（2）初始化三角形：讀取文件中的第一、二、三點為初始三角形。

（3）點定位搜索：由初始三角形根據四面體球缺角的特性2 在散亂點中尋找與之構成Delaunay 四面體的點，并生成Delaunay 四面體。

（4）循環運行搜索：在生成每一個四面體的同時，有3 個新的三角形生成，分別以新生成的三角形為初始三角形繼續尋找對應的點。由此循環運行直至將所有的散亂點進行四面體劃分。

2.3 搜索方法

步驟（3）中搜索方法分為以下幾個部分：

（1）確定搜索區域:為減少程序的運算量，將所有散亂點分為域內、域外2 個部分（方法見3-(3)）,只對域內部分進行搜索。

（2）計算球心坐標：分別將域內的所有點按順序分別與初始三角形的3 個頂點組合，由子程序計算此空間四點的外接球的球心坐標O。

（3）計算保存向量：根據三角形頂點坐標求出三角形外接圓圓心坐標O1，向量OO1即為初始三角形負法向量。保存球心O 與任意三角形頂點形成的向量OA (根據[特性一]，O 與任意頂點形成的向量與OO1夾角大小相同)。

（4）計算OA 與OO1夾角

計算出OA 與OO1夾角，并保存，以便經歷一次搜索后找出最大的夾角，并確定為最大球缺角，該角對應的點與此初始三角形構成Delaunay 四面體。

3 程序的準確性和高效性

為提高程序運行的準確性和效率，本程序進行了以下幾方面的處理:

（1）順序記錄：在每次計算所得的新三角形點數據，都按照順時針方向進行記錄,保存的新生成的數據格式符合初始運算格式，使程序可以循環進行計算。

（2）查重：由于初始三角形不同，對于每次生成的新的四面體來說有可能與前面生成的四面體重復，因此，對新生成的四面體進行查重，假設生成的新四面體4 個點的序號分別為（a，b，c，d），將(a，b，c，d)分別與之前生成的各數組對比，如果出現與（a，b，c，d）數值構成相同的數組，則取消該四面體，不存入，由此循環下去，并不斷檢驗是否重復，直到令所有點都找到Delaunay 四面體則計算過程結束。

（3）域內域外判斷：對于新生成的三角形來說一部分點位于該三角形的域內，另一部分位于域外。本程序對所有點與初始三角形的位置關系進行了預判，域外的三角形不必要進行復雜的計算，因此提高了程序的運行。四面體體積判定法:根據三維行列式正負的意義，A（x1，y1，z1）、B(x2，y2，z2)、C（x3，y3，z3)為空間的三點，D（x0，y0，z0）與三點的行列式表示的是四面體的體積（當D 位于平面ABC 的域內時V 為正，否則V 為負）。

四面體體積

（4）避免計算誤差：本文采取了一系列方法避免計算過程中產生的誤差，其中以下2 個方法用于整個程序中相關的計算步驟。①為避免計算機在計算時出現除以0 的數學錯誤，本程序在編譯時避免了除法運算，全部使用正負號的判斷，避免了計算機儲存或計算中出現誤差的情況。②在步驟（2）、（3）中計算球心、外接圓圓心時，本程序采用克拉默法則對方程組進行求解，但是在程序編譯過程中發現，編程平臺自帶的IMSL 函數庫在計算行列式過程中DET（）函數使用lin_sol_lsq（）來解方程［7］，對于行列式值為0 或非常接近于0 的矩陣會失效，所以本程序中使用了重新編譯的行列式計算函數進行相關計算，消除了此類計算誤差。

4 時間復雜度分析

為了對本程序的計算效率進行比較，本文進行了程序運行時間復雜度的分析［8］。圖3 為不同點數的情況下程序運行時間的趨勢線情況。

根據趨勢線可以看出計算的點的數量y 與程序運行時間x 的關系擬合曲線為y=8E-06x2，即時間與點數總量呈二次關系，該計算方法有良好的應用特性。

5 結束語

球缺角定義的提出，使Delaunay 四面體的剖分有了更加量化的方法，區別于以往的定性判斷，量化的方法誤差更小，更容易控制和比較。雖然在本方法的實踐過程中仍有一些方面不完善，但是相信在不斷地改進后，將對空間散亂點集Delaunay 四面體剖分方法的研究做出貢獻，可用于各個領域的工程應用軟件的開發。

［1］關文革，武強，賈麗萍，等.約束數據域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.

［2］WU Q，XU H.An approach to computer modeling and visualization of geological faults in 3d［J］.International Journal of Geographical Information Systems，1993，7（6）:501-524.

［3］XU H，WU Q.Design & implementation of visualization for 3d samdwich geological bodies［J］. Computer Applictations，2001，21（12）:56-60.

［4］Lawson C L. Software for C1 surface interpolation，Mathematical Software III［M］. New York: Academic Press，1977: 161-194.

［5］陳學工，潘懋.空間散點集Delaunay 四面體剖分切割算法［J］.計算機輔助設計與圖形學報，2002（1）:93-95.CHEN X G，PAN M. Delaunay Triangulation Cutting Algorithm for A Set of Irregularly Located Spatial Points［J］. Journal of Computer Aided Design & Computer Graphics，2002（1）: 93-95.

［6］王建華，徐強尋，張銳.任意形狀三維物體的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 arbitrary boundary［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（5）:717-720.

［7］李麗.三維空間Delaunay 三角剖分算法的研究及應用［D］.大連：大連海事大學，2010.［8］彭國倫.Fortran95 程序設計［M］.北京：中國電力出版社，2010.