馬獻(xiàn)德，謝小峰

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

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

1 引言

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

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

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

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

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

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

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

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

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

2 M eyer方法概述

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

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

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

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

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

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

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

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

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

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

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

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

（3）采用公式（5）計(jì)算vj。

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

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

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

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

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

設(shè)

則有

注意vj'與vj并不相同，但對(duì)于同一個(gè)v，求和號(hào)中vj'的數(shù)目與vj的數(shù)目是相等的，均為m。因?yàn)?/p>

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

同理

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

則有

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

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

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

因而得改進(jìn)算法步驟：

（1）計(jì)算

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

（2）計(jì)算

由公式（9）可得：

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

（3）計(jì)算

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

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

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

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

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

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

（2）求和：

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

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

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

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

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

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

（3）計(jì)算

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

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

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

5.1 混合面積普通算法簡述

AΔijk一般可采用如下方法計(jì)算：

AΔijk，V的計(jì)算步驟如下：

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

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

（3）計(jì)算面積AΔijk，V：

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

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

5.2 混合面積改進(jìn)算法分析

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

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

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

5.3 混合面積計(jì)算方法的判斷

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

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

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

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

6 改進(jìn)算法的效果

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

（1）給定曲面形式、幾何特征和空間范圍，隨機(jī)生成滿足這些條件的離散點(diǎn)坐標(biāo)及其三角網(wǎng)格拓?fù)浣Y(jié)構(gòu)，并計(jì)算各點(diǎn)處平均曲率和Gauss曲率的真實(shí)值。

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

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

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

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

圖6 計(jì)算精度對(duì)比（相對(duì)誤差直方圖）

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

圖7 計(jì)算效率對(duì)比（計(jì)算耗時(shí)直方圖）

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

7 結(jié)束語

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

[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曲面的點(diǎn)云曲率估算[J].浙江大學(xué)學(xué)報(bào)：工學(xué)版，2005，39（6）：761-764.

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