三角網(wǎng)格曲面上離散曲率改進算法

馬獻德，謝小峰

1.南京電子技術(shù)研究所一部，南京210039

2.南京祿口國際機場動力技術(shù)部，南京211113

1 引言

曲率是曲面上一點附近的重要局部幾何性質(zhì)，對于有精確解析形式的曲面模型，其上任一點處的曲率可以根據(jù)其坐標，由經(jīng)典的微分幾何方法的解析公式直接計算[1]。對于沒有精確解析形式，而只有表面離散點坐標數(shù)據(jù)及離散點之間拓撲關(guān)系的曲面模型，則只能完全借助于待求曲率的點附近各點的分布以及相互之間的拓撲關(guān)系，采用數(shù)值方法求解。求解曲面上一點處的曲率是獲取該曲面幾何信息的重要步驟，在離散幾何、計算幾何、計算機圖形學、虛擬現(xiàn)實技術(shù)等許多領(lǐng)域廣泛應(yīng)用。

由離散曲面上一點的鄰近局部點計算離散曲率，一般先把離散點劃分為網(wǎng)格曲面模型，根據(jù)待求點各鄰近網(wǎng)格的幾何信息計算待求點處的曲率。

Chen等[2]在1992年提出，將經(jīng)典微分幾何和曲面論中描述法曲率n、主曲率k1和k2與主方向e1和e2之間關(guān)系的Euler公式改寫為參數(shù)形式，并以圓弧擬合法和最小二乘法求解法曲率參數(shù)，從而計算主曲率。此方法容易實現(xiàn)，但圓弧擬合對真實的離散曲面而言過于粗糙，精度不高。

Taubin等[3]在1995年根據(jù)曲面微分幾何特性，以n、e1和e2積分構(gòu)造一個矩陣B，并證明B的特征向量就是n、e1和e2，特征值中就包含了k1、k2。通過待求曲率的頂點的鄰近三角面加權(quán)求和來估算矩陣B，進而計算主曲率。

M eyer等[4]在2003年采用離散化Laplace-Beltram i算子并積分來直接計算平均曲率kh，采用離散化Gauss-Bonnet定理直接計算Gauss曲率kg，并進而求解主曲率。

Dong等[5]在2005年提出基于弧長參數(shù)微分計算法曲率的方法，即點p1處的一個法曲率可由點p1及其附近一點p2的法向量及兩點各自坐標估算。在此基礎(chǔ)上借助上述Chen等提出的方法求解主曲率和主方向，克服了Chen法中圓弧擬合有時導致很大誤差的問題。

此外，還有M oreton等[6]提出的擬合法，Chen等[7]提出的重心加權(quán)法，Page等[8]提出的法向量投票法，柯映林等[9]提出的通過Shepard曲面插值點線性組合法等。

其中，M eyer的方法具有簡明的幾何意義，簡單易行，計算結(jié)果能達到一定的精度，且計算量較小。文獻[10]對多種算法進行了對比分析，認為M eyer方法的計算效果最優(yōu)。

經(jīng)分析發(fā)現(xiàn)，此方法的主要計算步驟可以通過轉(zhuǎn)換得以簡化，在計算精度和效率上可得到改進。

2 M eyer方法概述

M eyer提出的計算平均曲率kh和高斯曲率kg的方法[4]為：

其中，N1(pi)是pi的所有直接鄰近點（與pi共邊的點）的集合，n為pi點處的曲面單位法向量。vij=pj-pi，其余依此類推。符號·表示向量內(nèi)積運算。各參數(shù)的定義如圖1所示。

圖1 離散曲率計算參數(shù)示意圖

AMixed稱為混合面積，其定義如下：其中，AΔijk是Δpipjpk的面積，而AΔijk，V為Δpipjpk當Δpipjpk為銳角三角形時，圖2中四邊形ohkpihj的面積（點o為Δpipjpk的外心）。

圖2 銳角三角形的Voronoi面積(AΔijk，V)示意圖

上述公式需要多次計算三角函數(shù)值、反三角函數(shù)值，由于涉及到三角形的邊長、面積的計算，還需要多次開平方，計算中容易引入截斷誤差和舍入誤差，且計算效率較低。以下試圖精簡計算步驟，提高計算精度和計算效率。

假設(shè)pi共有m個直接鄰近點，從而有m個鄰近三角面（即有一個頂點為pi的三角面），pi及其所有直接鄰近點的直角坐標，pi處的曲面單位法向量n均已知。研究目標是根據(jù)這些已知條件，以最精簡的計算步驟完成以M eyer方法計算pi的離散平均曲率和離散Gauss曲率的過程。設(shè)

把v稱為平均曲率構(gòu)造向量，把θΣ稱為Gauss曲率構(gòu)造角?？梢妚、θΣ、AΔijk，V的計算是算法中的主要步驟，決定了整個計算過程的計算量。以下分別尋找這三個參數(shù)的改進算法。

3 對平均曲率構(gòu)造向量的改進算法設(shè)計

3.1 平均曲率構(gòu)造向量的普通算法簡述

一般通過計算三角形三邊長來計算公式（5）中的兩個余切值，則計算vj的步驟如下：

（2）采用三角形的性質(zhì)來計算三角形內(nèi)角的余切值。對于cotβij也按同樣方法計算：

（3）采用公式（5）計算vj。

以上步驟共17m次實數(shù)乘法，2m次實數(shù)除法，4m次實數(shù)開平方（加減法計算量很小，實數(shù)乘除2的整數(shù)次冪可由位運算實現(xiàn)，故不考慮其計算量）。

3.2 平均曲率構(gòu)造向量幾何意義

首先分析平均曲率構(gòu)造向量v的幾何意義。把平均曲率計算公式中涉及的各個角度變換到一個三角網(wǎng)格內(nèi)，如圖3所示。

圖3 Meyer法中各角度和向量變換示意圖

比較圖1和圖3可知，v可改寫為：

設(shè)

則有

注意vj'與vj并不相同，但對于同一個v，求和號中vj'的數(shù)目與vj的數(shù)目是相等的，均為m。因為

其中，vkj×vki表示向量vkj與vki的外積運算，并依此類推。所以

同理

如圖4所示，在平面pipjpk內(nèi)過點pi作直線pjpk的垂線pih，垂足為h。延長線段pih到h′，使得

則有

由此可見，v實際上是pi點的所有鄰近三角面片中，pi點所在三角形的對邊的正交向量的加權(quán)和，權(quán)值就是pi點對邊的長度。這就是平均曲率構(gòu)造向量v的幾何意義。

圖4 vj'的幾何意義示意圖

3.3 平均曲率構(gòu)造向量的改進算法分析

因而得改進算法步驟：

（1）計算

每2個相鄰三角形有1條公共邊，因而在N1(pi)內(nèi)a2和b2的總計算量為3m次實數(shù)乘法。

（2）計算

由公式（9）可得：

該步驟實際上附帶得到了Δpipjpk的面積，可在后續(xù)的計算過程中使用。

（3）計算

以上步驟共有22m次實數(shù)乘法，m次實數(shù)除法，m次實數(shù)開平方。

4 對Gauss曲率構(gòu)造角的改進算法設(shè)計

4.1 Gauss曲率構(gòu)造角的普通算法簡述

通常，Gauss曲率構(gòu)造角按如下步驟計算：

（1）對pi的每個直接鄰近三角形（即有一個頂點為pi的三角形），根據(jù)下式計算角θ：

其中a和b在平均曲率構(gòu)造向量的計算過程中已算出，不必重復計算。

（2）求和：

以上步驟共有4m次實數(shù)乘法，m次實數(shù)除法，m次實數(shù)反余弦運算。

4.2 Gauss曲率構(gòu)造角的改進算法分析

利用了正切函數(shù)及兩角和的正切公式，所得的改進算法步驟如下：

（1）對于pi的每個直接鄰近三角形，計算cotθj的值（或θj=π/2）。根據(jù)公式（10）有：

（2）計算cotθΣ的值（或得到θΣ=(2r+1)π/2，r為非負整數(shù)），依次計算：

其中t的取值依次為1，2，…，m。注意須先處理分母為零的特殊情況，同時根據(jù)參與計算的兩個余切值和所得結(jié)果的正負號，確定所得結(jié)果對應(yīng)的角度所在的區(qū)間。最后得到cotθΣ=x。

（3）計算

其中，m0為非負整數(shù)，由上述得到的θΣ所在的區(qū)間而確定。一般而言，m0的值不大，此處的m0π可用少數(shù)幾次加法實現(xiàn)。

以上步驟共有2m次實數(shù)乘法，m次實數(shù)除法，1次實數(shù)反余切運算。

5 對混合面積的改進算法設(shè)計

5.1 混合面積普通算法簡述

AΔijk一般可采用如下方法計算：

AΔijk，V的計算步驟如下：

（1）設(shè)Δpipjpk外接圓半徑為R，則：

（2）計算圖2中線段ohj和ohk的長度。

（3）計算面積AΔijk，V：

AΔijk，V=(a|ohj|+b|ohk|)/4

注意到在平均曲率構(gòu)造向量的普通算法中，實際上已經(jīng)計算出a2、b2、c2和。因而以上步驟共有5m次實數(shù)乘法，2m次實數(shù)開平方（假設(shè)所有三角面都是銳角三角形，下同）。

5.2 混合面積改進算法分析

如圖5所示，連結(jié)pio并延長到點pm，使得pkpm||ohj，連結(jié)pkpm與pjpm。推導可得

故共需m次實數(shù)乘法。

圖5 Voronoi面積改進算法構(gòu)造示意圖

5.3 混合面積計算方法的判斷

前文已提到，計算混合面積時，對不同的三角形，其計算方法不同。根據(jù)前述計算，有

其中ui、uj和uk已在前述計算步驟中得到，因此對每個三角面片，決定采用哪種方法計算其混合面積時，并不需要額外的計算量。

綜上所述，當所求點處的直接臨近三角形全為銳角三角形時，普通算法共需26m次實數(shù)乘法，3m次實數(shù)除法，6m次實數(shù)開平方，m次實數(shù)反三角函數(shù)；而改進算法共需25m次實數(shù)乘法，2m次實數(shù)除法，m次實數(shù)開平方，1次實數(shù)反三角函數(shù)。

可見，改進算法中除了乘法次數(shù)與普通算法基本持平外，計算效率很低、耗時很長的除法、開平方和反三角函數(shù)的次數(shù)減少，因此提高了計算效率，同時更重要的是由于減少了計算量，避免了不必要的截斷誤差和舍入誤差，從而有望減小計算誤差，提高計算精度。

6 改進算法的效果

為驗證所提算法的有效性，構(gòu)造離散的橢球面、橢圓柱面、橢圓拋物面、雙曲拋物面、雙葉雙曲面、單葉雙曲面等6種離散標準二次曲面，分別按普通算法和所提改進算法計算對比。步驟如下：

（1）給定曲面形式、幾何特征和空間范圍，隨機生成滿足這些條件的離散點坐標及其三角網(wǎng)格拓撲結(jié)構(gòu)，并計算各點處平均曲率和Gauss曲率的真實值。

（2）計算各點處的法向量真實值，在上述算法中的離散法向量均直接采用真實值，這樣可避免采用離散算法計算法向量帶來的誤差。

（3）在同等計算條件下，分別采用普通方法和所提改進方法來實現(xiàn)M eyer算法，得到兩種方法的相對誤差（剔除因真實曲率為0而無法計算相對誤差的點）；并記錄兩種算法所耗費的計算時間，得到計算精度和計算效率的評價。

采用的計算軟硬件平臺為：Intel?CoreTM2 Duo T5870@2.00 GHz CPU，4.00 GB內(nèi)存，Window s 7操作系統(tǒng)，Visual C++2008開發(fā)平臺。

曲率相對誤差的統(tǒng)計直方圖如圖6所示。普通方法的平均曲率相對誤差的均方根值為7.722%，Gauss曲率相對誤差的均方根值為7.712%；改進方法平均曲率相對誤差的均方根值為1.066%，Gauss曲率相對誤差的均方根值為1.058%。

圖6 計算精度對比（相對誤差直方圖）

對一個離散點的平均曲率與Gauss曲率的計算耗時的統(tǒng)計直方圖，如圖7所示。普通方法的耗時均方根值為1.222m s，改進方法耗時均方根值為1.036m s。

圖7 計算效率對比（計算耗時直方圖）

由計算結(jié)果可見，所提的改進方法大大提高了M eyer算法對離散平均曲率和Gauss曲率的計算精度，同時還提高了其計算效率。

7 結(jié)束語

通過分析M eyer所提的離散曲面一點處平均曲率和Gauss曲率的算法，提出了平均曲率構(gòu)造向量和Gauss曲率構(gòu)造角的概念，分析了其幾何意義，并據(jù)此提出了提高計算精度和計算效率的改進算法。對6種典型離散標準二次曲面的仿真計算結(jié)果，表明了本文改進算法可大大提高曲率計算精度，同時還可提高計算效率。

[1]do Carmo M.Differential geometry of curves and surfaces[M].Englewood Cliffs，NJ：Prentice-Hall，1976.

[2]Chen Xin，Schmitt，F(xiàn).Intrinsic surface properties from surface triangulation[C]//Proceedings of the European Conference on Computer Version.Heidelberg：Springer-Verlag，1992，588：739-743.

[3]Taubin G.Estimating the tensor of curvature of a surface from a polyhedral approximation[C]//Proceedings of the Fifth International Conference on Computer Vision.Piscataway，NJ：IEEE，1995：902-907.

[4]Meyer M，Desbrun M，Schr?der P，et al.Discrete differential geometry operators for triangulated 2-Manifolds[J].Visualization and Mathematics，2003，3：35-57.

[5]Dong Chenshi，Wang Guozhao.Curvatures estimation on triangular mesh[J].Journal of Zhejiang University Science，2005，6A：128-136.

[6]Moreton H P，Sequin C H.Functional optimization for fair surface design[C]//Computer Graphics Proceedings，Annual Conference Series（ACM SIGGRAPH），Chicago，Illinois，1992：167-176.

[7]Chen S，Wu J.Estimating normal vectors and curvatures by centroid weights[J].Computer Aided Geometric Design，2004，21（5）：447-458.

[8]Page D L，Sun Y，Koschan A F，et al.Normal vector voting：crease detection and curvature estimation on large，noisy meshes[J].Graphical Models，2002，64（3/4）：199-229.

[9]柯映林，陳曦.基于4D Shepard曲面的點云曲率估算[J].浙江大學學報：工學版，2005，39（6）：761-764.

[10]方惠蘭，王國瑾.三角網(wǎng)格曲面上離散曲率估算方法的比較與分析[J].計算機輔助設(shè)計與圖形學學報，2005，17（11）：2500-2507.