夏 俊，李映華

(昆明理工大學(xué) 理學(xué)院，昆明 650500)

球面凸類圖形Delaunay三角剖分再分算法及其收斂性分析

夏 俊*，李映華

(昆明理工大學(xué) 理學(xué)院，昆明 650500)

在計算曲面Ricci Flow時，會因為三角網(wǎng)格中存在過小的角而出現(xiàn)不收斂的情況。針對這種不收斂的問題，提出一種提高最小角角度的球面凸類圖形Delaunay三角剖分再分算法。首先,給出球面凸類圖形Delaunay三角剖分再分算法。它的核心操作有兩個：1)如果某條Delaunay劣弧被“侵占”, 通過添加Delaunay劣弧中點分割Delaunay劣弧;2)如果存在“瘦”球面三角形, 通過添加球面三角形外接球面小圓圓心分解球面三角形。然后，利用局部特征尺度探索出所提算法的收斂條件并給出輸出頂點的一個上界公式。根據(jù)實驗輸出的網(wǎng)格驗證，所提算法網(wǎng)格生成的球面三角形沒有狹小的角，適合用來計算Ricci Flow。

球面；Delaunay三角剖分；劣弧；局部特征尺度；收斂

0 引言

Delaunay三角剖分在圖形學(xué)、地理學(xué)以及有限元分析等領(lǐng)域是一項極為重要的預(yù)處理技術(shù)。隨著計算機技術(shù)的發(fā)展以及分析對象的復(fù)雜性，Delaunay三角剖分在這些領(lǐng)域扮演著重要角色。例如，在有限元分析中，需要把研究的區(qū)域分解成三角形單位元。如果剖分的三角形角度過小，那么有限元算法的收斂性得不到保證。因為Delaunay三角剖分可以最大化最小角，所以Delaunay三角剖分是一種有效的剖分方法。

目前，平面上的Delaunay三角剖分的研究已經(jīng)相當(dāng)成熟。平面上的Delaunay三角剖分分為點集Delaunay三角剖分和約束Delaunay三角剖分。 對于平面點集Delaunay 三角剖分來說，Lee等[1]給出了點集Delaunay三角剖分的幾何性質(zhì)和算法。Dwyer[2]利用分治合并算法構(gòu)建了點集Delaunay三角剖分。近年來，對于點集Delaunay三角剖分的研究都是不斷降低時間復(fù)雜度，插入點混合定位算法[3]有效降低算法的時間復(fù)雜度。對于平面約束Delaunay三角剖分來說，Chew[4]提出了分治合并算法，Lee等[5]對平面直線圖形作約束Delaunay三角剖分。近幾年，關(guān)于平面約束Delaunay三角剖分的研究也都是不斷降低復(fù)雜度，如帶斷層約束的Delaunay三角剖分混合算法[6]，它在生成網(wǎng)格質(zhì)量和時間效率上都優(yōu)于傳統(tǒng)算法。隨著計算機技術(shù)的發(fā)展，約束Delaunay三角剖分已經(jīng)不能滿足生產(chǎn)需求。在有限元分析中，網(wǎng)格的質(zhì)量對有限元方程的收斂性有著直接的影響。Ruppert[7]通過添加“Steiner”點提出Delaunay Refinement算法以及Shewchuk[8]提出的Delaunay Refinement算法增大了輸出三角形的最小角角度，這些網(wǎng)格都適用于有限元方程。但是Ruppert’s Refinement算法復(fù)雜度的證明由Barbic等[9]完成。Rand[10]提出的“Off-Centers”不僅提高輸出三角形的質(zhì)量，還添加更少的輸出點。

近年來，出現(xiàn)了非歐空間的Delaunay三角剖分，如Caroli等[11]給出了球面點集Delaunay三角剖分和Zeng等[12]提出的雙曲曲面Delaunay Refinement三角剖分。在計算Ricci Flow方程時，一旦出現(xiàn)某個退化的球面三角形，就有可能出現(xiàn)不收斂的情況，文獻(xiàn)[13-14]研究就需要高質(zhì)量的球面三角形。因此，必須提高球面三角形的質(zhì)量。本文通過添加Delaunay劣弧中點以及球面三角形外接測地圓圓心，有效地提高整體球面三角形的質(zhì)量，避免狹小角的出現(xiàn)，有效解決Ricci Flow方程不收斂的問題。

1 球面凸類圖形Delaunay三角剖分再分算法

首先,給出幾個相關(guān)概念。

定義1 球面上兩點之間劣弧的長度,稱之為球面上兩點之間的距離, 簡稱距離。

定義2 球面n邊形A1A2…An總在它的每一邊所在大圓的半個球面內(nèi), 那么稱這個球面多邊形為球面凸n邊形。

定義3 球面凸n邊形以及包含在球面凸n邊形內(nèi)的點和劣弧稱為球面凸類圖形。

因為不同半徑的球面之間可以建立一一映射，所以不妨在單位球面上對輸入的球面凸類圖形進行Delaunay三角剖分再分，如果沒有特殊說明，下文中的球面指的是單位球面。

為了方便，假設(shè)球面凸類圖形中球面凸n邊形兩條相交劣弧夾角角度至少是90°，這個限制根據(jù)文獻(xiàn)[7]中的“Corner-Lopping” 可以解除，具體不再贅述；同時，記球面凸類圖形為X。另外本文中提到的角指的是球面角。

球面凸類圖形Delaunay三角剖分再分算法是在球面點集Delaunay三角剖分算法基礎(chǔ)上生成的，球面點集Delaunay三角剖分(Sphere Delaunay Triangulation(V))簡記為SDT(V)。球面點集Delaunay三角剖分指的是給定球面上一系列點，由球面點集Delaunay三角剖分算法產(chǎn)生的所有球面三角形最短邊對應(yīng)的外接球面小圓圓心角達(dá)到最大，也就是最大化最小圓心角。如果球面上不存在四點共圓，那么這個SDT(V)是唯一的。如果存在四點共圓，兩條對角線任取一條不影響最小圓心角的大小，不妨取對角線最短的作為邊。具體的算法實現(xiàn)見參考文獻(xiàn)[11]。 為了方便下文的分析，本文只需要對球面凸類圖形內(nèi)部進行Delaunay三角剖分再分，故球面凸類圖形外的球面三角形全部刪除。

把球面凸類圖形中的劣弧以及剖分產(chǎn)生的劣弧叫作Delaunay劣弧，用集合G表示。球面上任意一條劣弧AB記為AB。為了方便，也可用AB表示球面上A和B兩點之間的距離， 那么球面上AB的取值范圍為(0,π]。同時，把球面凸類圖形中的點以及剖分添加的點叫作頂點，用集合V表示。球面上任意一個點叫作點。對于這個算法，初始時，V就是所有輸入的頂點，G就是所有的輸入劣弧。

生成球面點集Delaunay三角剖分之后，利用本文算法繼續(xù)剖分，添加點的原因有兩個：保證G中的Delaunay劣弧是點集Delaunay三角剖分產(chǎn)生的測地線；整體提高球面三角形最短邊外接球面小圓圓心角。

對Delaunay劣弧g來說，以g為測地直徑的球面小圓稱為g的測地直徑圓。如果有一個頂點A在g的測地直徑圓內(nèi)，稱頂點A“侵占”Delaunay劣弧g。如圖1所示， 輸入的球面凸類圖形是圖中的實Delaunay劣弧，球面點集Delaunay三角剖分是圖中除Delaunay劣弧AB外，其他實Delaunay劣弧和虛劣弧組成的圖形。對于球面凸類圖形來說，這不是它的Delaunay三角剖分，因為AB不是點集Delaunay三角剖分中的測地線且頂點A“侵占”Delaunay劣弧CD。

本文算法的核心操作有兩個：第一，如果某條Delaunay劣弧被“侵占”，通過添加Delaunay劣弧中點分割Delaunay劣?。坏诙?，如果存在“瘦”球面三角形，通過添加球面三角形外接球面小圓圓心分解球面三角形，“瘦”球面三角形是球面三角形外接球面小圓中，球面三角形最短邊對應(yīng)的圓心角的角度小于φ(φ是預(yù)先給定的值)的球面三角形。

圖1 非法球面Delaunay三角剖分Fig. 1 Delaunay triangulation of illegal sphere

下面給出球面凸類圖形Delaunay三角剖分再分算法：

輸入 球面凸類圖形X和參數(shù)φ。

BEGIN

初始化 計算初始球面點集Delaunay三角剖分SDT(V)

while Delaunay劣弧g被“侵占”：

SpliltGeo(Delaunay劣弧g)

if “瘦”球面三角形t的外接球面小圓圓心P“侵占”Delaunay劣弧g1,g2,…,gk

fori=1 tok：

SplitGeo(Delaunay劣弧gi)

else SplitStri(球面三角形t)

end if

until 沒有Delaunay劣弧被“侵占”和沒有球面三角形最短邊所對應(yīng)的圓心角<φ

END

子算法SplitStri(球面三角形t)。

輸入 球面三角形t。

BEGIN

添加球面三角形t的外接球面小圓圓心到V中，更新SDT(V)

END

子算法SplitGeo(Delaunay劣弧g)。

輸入 Delaunay劣弧g。

BEGIN

添加Delaunay劣弧g的中點到V中，更新SDT(V)，將g從G中 移除，將g分成的兩個子段g1和g2添加到G中

END

2 算法收斂性分析

在本章證明本文算法在執(zhí)行有限的步驟之后是可以終止的，這與輸入球面凸類圖形的局部特征尺度有關(guān)，給出局部特征尺度的概念。

定義4 給定一個球面凸類圖形X，以球面上任意一點P為圓心的與X中兩個和P不相關(guān)的頂點或Delaunay劣弧相交的最小球面小圓的測地半徑叫作點P的局部特征尺度，記為lfs(P)。

在球面上，根據(jù)定義1知，lfs(P)的取值范圍為(0,π)。圖2是對這個定義的直觀解釋，lfs(A)、lfs(B)和lfs(C)分別對應(yīng)圖2中三個球面小圓的測地半徑，其中l(wèi)fs(C)有一部分頂點或Delaunay劣弧相交。對于給定的輸入球面凸類圖形X，lfs(X)是連續(xù)函數(shù)且滿足下列引理1。

引理1 局部特征尺度是Lipschitz函數(shù)，即任意給定一個球面凸類圖形，在球面上任取兩個點P和Q，有：

lfs(Q)≤lfs(P)+PQ

其中PQ是球面上兩點P和Q的距離。

證明 根據(jù)lfs()的定義，lfs(P)是以點P為圓心的與球面凸類圖形中兩個和P不相關(guān)的頂點或Delaunay劣弧相交的最小球面小圓的測地半徑，記這個球面小圓為D1。

1)若lfs(P)+PQ≥π，因為lfs()<π，顯然成立。

2)若lfs(P)+PQ<π，現(xiàn)在以點Q為圓心，r=lfs(P)+PQ為測地半徑的球面小圓記為D2，那么D1包含在D2中，因此：

lfs(Q)≤r

故得到：

lfs(Q)≤r=lfs(P)+PQ

每當(dāng)一個點添加之后，它不改變lfs()的性質(zhì)，lfs()只和輸入的球面凸類圖形有關(guān)，下面的引理2說明頂點的密度完全是由lfs()決定。

引理2 1)初始時，對于每個輸入的頂點P，P到它最近的頂點的距離至少是lfs(P)。

2)當(dāng)點P是球面“瘦”球面三角形的外接球面小圓圓心時，點P到它最近頂點的距離至少是lfs(P)/CT(點P屬于這個三角剖分中的頂點或者是“侵占”某條Delaunay劣弧)。

3)當(dāng)頂點P是某條Delaunay劣弧的中點時，點P到它最近的頂點的距離至少為lfs(P)/CG。

證明 當(dāng)頂點P是輸入的點時，根據(jù)lfs()的定義，頂點P到它最近的頂點的距離至少是lfs(P)。對于后來加進去的點，假設(shè)這個引理對之前所有的頂點都成立。

1)當(dāng)點P是球面“瘦”三角形t的外接球面小圓圓心時，那么點P到球面“瘦”三角形t的頂點最近。設(shè)點P到球面“瘦”三角形t的頂點A、頂點B和頂點C的距離等于r1，那么r1的取值范圍為(0,π/2)。如圖3所示，不妨假設(shè)弧AB是球面△ABC的最短邊，那么弧AB所對應(yīng)的圓心角最小，設(shè)為ψ。

不失一般性，首先考慮點A是在點B之后加進去的或者兩者都是輸入的頂點。當(dāng)點A是輸入的頂點時，則lfs(A)≤AB；當(dāng)點A是球面三角形外接球面小圓圓心時，設(shè)這個球面小圓的測地半徑為rA，根據(jù)這個引理有l(wèi)fs(A)≤CTrA≤CTAB；當(dāng)點A是某條Delaunay劣弧的中點時，在A點利用這個引理有l(wèi)fs(A)≤CGAB。

假設(shè)給定條件CG≥CT≥1，則上面幾種情況綜合起來為：

lfs(A)≤CGAB

如圖4所示，先把球面△PAB提出來，再由點P作AB的垂線，根據(jù)球面正弦定理有：

由此得到：

AB=2arcsin(sin(ψ/2)·sinr1)

令ρ(ψ,r1)=(AB)/(2r1)，因為:

lfs(P)≤lfs(A)+r1

則ρ(ψ,r1)在r1∈(0,π/2)上的最大值的上確界為sin(ψ/2)，得：

這里只需要取CT≥2sin(ψ/2)CG+1。

圖3 球面三角形及其圓心角Fig. 3 Spherical triangle and its central angle

圖4 球面三角形及其垂線Fig. 4 Spherical triangle and its vertical

2)如圖5所示，當(dāng)頂點P是Delaunay劣弧g的中點時，如果存在頂點A或某個球面“瘦”三角形外接球面小圓圓心A在以g為測地直徑的球面小圓內(nèi)時，Delaunay劣弧g需要添加中點并假設(shè)這個球面小圓的測地半徑為r2，易知，r2的取值范圍為(0,π/2)。

圖5 點A“侵占”劣弧BCFig. 5 Vertex A “encroaches upon” minor arc BC

a)點A位于Delaunay劣弧h上并且和Delaunay劣弧g不相交，之前假設(shè)所有的輸入球面凸類圖形的角度至少是90°，所以必然存在兩個不相交的Delaunay劣弧，一個包含點P，另一個包含點A，使得它們的距離比r2小。因此，得到：

lfs(P)≤r2

根據(jù)CG≥1知，這種情況是成立的。

b)點A是球面“瘦”三角形外接球面小圓圓心，但是這個點在以Delaunay劣弧g為測地直徑的球面小圓內(nèi)，所以這個點是不存在于三角剖分內(nèi)。那么存在一個以點A為圓心，測地半徑為rA的球面小圓CA使得圓內(nèi)沒有其他點，在A點應(yīng)用這個引理得到：

rA≥lfs(A)/CT

同時，點B和點C在CA之外，故當(dāng)rA最大時有：

cosrA=cosr2×cosr2

令ρ(r2)=rA/r2，根據(jù)引理1有：

lfs(P)≤lfs(A)+r2

根據(jù)上面的分析，得到如下不等式：

通過解這個方程組，得到：

所以ψ的一個上界為φ=41.4°。

引理2表明當(dāng)點P加進去之后，沒有其他頂點在以P為圓心、lfs(P)/CG為測地半徑的球面小圓內(nèi)。下面的這個定理表明當(dāng)點加進去之后，頂點P到其他頂點的距離存在下界。

定理1 給定輸出的球面凸類圖形的Delaunay三角剖分再分，對于任意一個頂點P，頂點P到它最近的頂點的距離至少是lfs(P)/(CG+1)。

證明 引理2解決了除P在Q之后加入的情況，利用這個引理2對Q進行處理，得到：

PQ≥lfs(Q)/CG

結(jié)合引理1得到：

PQ≥lfs(P)/(CG+1)

有了上面的定理之后,就能證明輸出頂點存在上界。

定理2 Delaunay三角剖分再分輸出球面三角網(wǎng)格中，頂點的數(shù)量至多是：

其中：B是球面凸類圖形的內(nèi)部；C1是具體的常數(shù)。

根據(jù)引理1，在DP中，lfs()的最大值是lfs(P)+rP，設(shè)DP的面積為SDP，故：

下面對rP分情況討論：

1)當(dāng)rP≤π/2時，因為SDP=4π sin2(rP/2)，所以:

綜合起來得：

2)當(dāng)π/2

因此，只需取C1≥(2CG+3)2π即可。

根據(jù)以上兩種情況知，只需取C1≥(2CG+3)2π。

由上面的分析知，當(dāng)輸入的參數(shù)φ≤41°時，本文算法在執(zhí)行過程中添加的點是有限的，也就是本文算法是收斂的。

3 算法實現(xiàn)及分析

球面凸類圖形的Delaunay三角剖分再分算法比歐氏空間的Delaunay Refinement Algorithm復(fù)雜。Delaunay Refinement Algorithm收斂條件由剖分的三角形最小角限定，本文算法收斂條件是由剖分三角形最短邊對應(yīng)的球面圓心角限定。這是因為球面三角形內(nèi)角和不是固定常數(shù)，導(dǎo)致無法通過球面三角形最小角限定收斂條件。在具體算法實現(xiàn)過程中，球面凸類圖形Delaunay三角剖分再分算法比Delaunay Refinement Algorithm更復(fù)雜。Delaunay Refinement Algorithm可以由Matlab直接實現(xiàn)并且可視化，但是球面凸類圖形Delaunay三角剖分再分算法只能通過程序語言輸出邊和點的數(shù)據(jù)，再由其他繪圖軟件可視化。

本文算法的核心操作和平面的情況是相同的，具體的證明過程見參考文獻(xiàn)[9],所以本文算法的時間復(fù)雜度和平面情況相同，都是O(n2)。本文算法的程序?qū)崿F(xiàn)如圖6所示。

圖6 程序?qū)崿F(xiàn)圖(最小圓心角為41.4°)Fig. 6 Program implementation diagram (minimum central angle is 41.4°)

本文算法還可以適當(dāng)提高最小圓心角的角度，可以在本文算法中插入一個限制球面三角形大小的語句，不妨以球面三角形外接測地圓半徑為限制。這種方法可以在球面上作近似均勻剖分，生成盡量同樣大小的三角形，并且可以加密剖分，滿足許多應(yīng)用的要求。

4 結(jié)語

近年來，有關(guān)Delaunay三角剖分研究重點仍然是在歐氏空間進行。對于非歐空間的Delaunay三角剖分研究得還不夠深入。在求解Ricci Flow方程時，Ricci Flow的收斂性是和剖分的三角形相關(guān)的，一旦出現(xiàn)某個角度過小的三角形，Ricci Flow就有可能出現(xiàn)不收斂情況。本文中的三角網(wǎng)格通過提高最短邊對應(yīng)的圓心角提高三角形最小角角度，有效解決不收斂的問題。本文算法中添加的是球面三角形外接測地圓的圓心，下一步可以研究添加“Off-center Steiner Vertices”，或許這種方法輸出的點更少，且可以突破本文收斂條件41.4°的限制，提高到更大的角度。

References)

[1] LEE D T, SCHACHTER B J. Two algorithms for constructing Delaunay triangulations [J]. International Journal of Computer and Information Sciences, 1980, 9(3): 219-242.

[2] DWYER R A. A faster divide-and-conquer algorithm for constructing Delaunay triangulations [J]. Algorithmica, 1987, 2(1/2/3/4): 137-151.

[3] 高莉.改進的Deluanay三角剖分算法研究[D].蘭州：蘭州交通大學(xué),2015:30-33.(GAO L. An improved Deluanay triangulation algorithm [D]. Lanzhou: Lanzhou Jiaotong University, 2015: 30-33.)

[4] CHEW L P. Constrained Delaunay triangulation [J]. Algorithmica, 1989, 4(1/2/3/4): 97-108.

[5] LEE D T, LIN A K. Generalized Delaunay triangulation for planar graphs [J]. Discrete & Computational Geometry, 1986, 1(3): 201-217.

[6] 張群會,解子毅.帶斷層約束的Delaunay三角剖分混合算法[J].西安科技大學(xué)學(xué)報,2014,34(1):52-56.(ZHANG Q H, XIE Z Y. Mixed algorithm of Delaunay triangular subdivision with fault constraint [J]. Journal of Xi’an University of Science and Technology, 2014, 34(1): 52-56.)

[7] RUPPERT J. A Delaunay refinement algorithm for quality 2-dimensional mesh generation [J]. Journal of Algorithms, 1995, 18(3): 548-585.

[8] SHEWCHUK J R. Delaunay refinement algorithms for triangular mesh generation [J]. Computational Geometry, 2002, 22(1/2/3): 21-74.

[9] BARBIC J, MILLER G. A quadratic running time example for ruppert’s refinement algorithm [EB/OL]. [2017- 04- 16]. http://www-bcf.usc.edu/～jbarbic/BarbicMiller-RupperQuadraticExampleTechReport.pdf.

[10] RAND A. Where and how Chew’s second Delaunay refinement algorithm works [C/OL]// Proceedings of the 2013 23rd Annual Canadian Conference on Computational Geometry, [2017- 04- 16]. http://www.w.cccg.ca/proceedings/2011/papers/paper91.pdf.

[11] CAROLI M, DE CASTRO P M M, LORIOT S, et al. Robust and efficient Delaunay triangulations of points on or close to a sphere [C]// Proceedings of the 2010 9th International Conference on Experimental Algorithms. Berlin: Springer, 2010: 462-473.

[12] ZENG W, SHI R, GU X F. Global surface remeshing using symmetric Delaunay triangulation in uniformization spaces [C]// Proceedings of the 2011 Eighth International Symposium on Voronoi Diagrams in Science and Engineering. Piscataway, NJ: IEEE, 2011: 160-169.

[13] GU X F, GUO R, LUO F, et al. A discrete uniformization theorem for polyhedral surfaces II [EB/OL]. [2017- 04- 16]. http://math.oregonstate.edu/～guoren/docs/u10-10.pdf.

[14] GU X F, GUO R, LUO F, et al. A discrete uniformization theorem for polyhedral surfaces I [EB/OL]. [2017- 04- 16]. https://arxiv.org/pdf/1401.4594.pdf.

[15] PREPARATA F P, SHAMOS M I. Computational Geometry — an Introduction [M]. New York: Springer, 1985: 241-276.

This work is partially supported by the Natural Science Foundation of Kunming University of Science and Technology (KKSY201507066).

XIAJun, born in 1991, M. S. candidate. His research interests include computational geometry.

LIYinghua, born in 1978, Ph. D., lecturer. His research interests include computational geometry.

Delaunaytriangulationsubdivisionalgorithmofsphericalconvexgraphanditsconvergenceanalysis

XIA Jun*, LI Yinghua

(SchoolofScience,KunmingUniversityofScienceandTechnology,KunmingYunnan650500,China)

When calculating curved Ricci Flow, non-convergence emerges due to the existence of undersized angles in triangular meshes. Concerning the problem of non-convergence, a Delaunay triangulation subdivision algorithm of spherical convex graph of enhancing the minimum angle was proposed. First of all, the Delaunay triangulation subdivision algorithm of spherical convex graph was given. The proposed algorithm had two key operations: 1) if a Delaunay minor arc was “encroached upon”, a midpoint of the Delaunay minor arc was added to segment the Delaunay minor arc; 2) if there was a “skinny” spherical triangle, it was disassembled by adding the center of minor circle of its circumscribed sphere. Then, the convergence criteria of the proposed algorithm was explored on local feature scale and an upper-bound formula of the output vertex was given. The grids based on the output of experiment show that the spherical triangle generated by the grids of the proposed algorithm has no narrow angle, so it is suitable for calculating Ricci Flow.

sphere; Delaunay triangulation; minor arc; local feature scale; convergence

2017- 06- 29;

2017- 09- 04。

昆明理工大學(xué)自然科學(xué)基金資助項目(KKSY201507066)。

夏俊(1991—)，男，安徽合肥人，碩士研究生，主要研究方向：計算幾何； 李映華(1978—)，男，廣東梅州人，講師,博士，主要研究方向：計算幾何。

1001- 9081(2017)12- 3558- 05

10.11772/j.issn.1001- 9081.2017.12.3558

(*通信作者電子郵箱yinghuali@kmust.edu.cn)

TP301.6

A