梁增凱,孫殿柱,薄志成,李延瑞,林偉

(1.山東理工大學機械工程學院,255049,山東淄博;2.西安交通大學機械工程學院,710049,西安)

曲率是曲面信息度量的重要指標,可反映曲面任一點處局部幾何性質,在數(shù)字化設計與制造、醫(yī)學輔助設計、文物保護等領域具有廣泛的應用價值。三維掃描設備的發(fā)展使得獲取反映物體細節(jié)特征的密集點云數(shù)據(jù)成為現(xiàn)實,但所得數(shù)據(jù)通常只包含樣點位置坐標且拓撲信息缺失,因此準確估計散亂點云曲率對不同視角的點云配準、特征識別以及曲面重建等后續(xù)處理具有重要意義[1]。

現(xiàn)有點云曲率估計方法主要分為基于樣點k鄰域點集擬合的方法和基于三角網(wǎng)格的方法。Hoppe等采用主元分析法擬合采樣點的k鄰域點集,通過曲面變分求取其近似曲率[2]。Gumhold等選用移動最小二乘法對樣點鄰域進行二次曲面擬合,構造逼近原始曲面的局部參數(shù)曲面,基于所得參數(shù)曲面估計樣點曲率[3-4]。Ernest等通過坐標變換構造拋物面,并基于加權擬合有效降低了曲率估計值的誤差[5]。然而,基于樣點k鄰域點集擬合的方法在高復雜度的模型外表面和尖銳特征區(qū)域難以尋找合適的擬合曲面形式,限制了曲率估計的適用范圍。基于三角網(wǎng)格的方法通常先對散亂點云進行曲面重建,依據(jù)重建所得網(wǎng)格結構建立數(shù)據(jù)間拓撲關系,進而準確估計樣點曲率。其中Taubin在三角網(wǎng)格頂點處建立近似矩陣,通過求解矩陣特征值估計頂點曲率[7];Meyer等基于Voronoi元與有限元定義離散的微分幾何算子,利用微分幾何中的Green公式進行離散化求解,使計算結果準確性進一步提升[8-9]?；谌蔷W(wǎng)格的方法參考樣點局部屬性,能對不同區(qū)域樣點曲率進行準確估計,應用范圍更為廣泛。

在實際掃描原表面獲取點云過程中,因測量距離及掃描角度的變化,即便對于光滑曲面,仍會存在采樣誤差。特別當點云模型坐標出現(xiàn)法向偏移誤差時,曲面G2連續(xù)性[10]難以延續(xù),從而導致現(xiàn)有基于重建網(wǎng)格曲面的曲率估計方法計算結果存在偏差且過渡不夠平滑。為解決這一問題,本文將目標樣點的鄰域點集作為局部樣本進行曲面重建,獲得插值于采樣點集且與原表面拓撲同構的局部網(wǎng)格曲面,基于網(wǎng)格曲面中一階鄰域面估計樣點曲率,進而依據(jù)網(wǎng)格曲面上測地距離對樣點曲率估計結果進行修正。在曲率估計的過程中,基于一階鄰域面形心權重計算所需的樣點法向,提高曲率估計的精度與穩(wěn)健性。實驗結果表明,該算法具有較高的計算效率,能夠有效抑制曲面采樣誤差的干擾,使得樣點曲率估計結果過渡更為平滑。

1 曲面局部拓撲重建

為獲得插值于采樣點集的高質量Delaunay三角網(wǎng)格曲面,本文依據(jù)局部平坦程度,采用二維Delaunay網(wǎng)格剖分與三維Delaunay網(wǎng)格過濾相結合的策略重建局部樣本。對于采樣點集S中目標樣點p,將p的k鄰域點集作為局部重建樣本,通過適度擴大鄰域搜索范圍對該樣本進行輔助點添加。設kη為獲取λ(p)時所需的鄰域點集數(shù)量,kξ為添加輔助點后鄰域點集數(shù)量,局部重建算法的具體流程如下。

(1)搜索p的kη鄰域,獲得局部重建樣本λ(p)。

(2)搜索p的kξ鄰域點集,獲得包含輔助點的局部重建樣本φ(p)。在φ(p)中,若?pi,pi∈φ(p)且pi?λ(p),則將pi標記為輔助點。

(3)對φ(p)中樣點進行共面檢測,若共面,則對φ(p)進行二維Delaunay網(wǎng)格剖分,輸出剖分的局部重建網(wǎng)格,程序結束。

(4)對φ(p)進行三維Delaunay網(wǎng)格剖分,獲得四面體集合D(φ(p)),并求出對應的Voronoi圖V(φ(p))。

(5)對D(φ(p))中任一面片Ti進行夾角檢測,若檢測未通過,則將Ti從D(φ(p))中刪除,否則保留。

(6)對D(φ(p))進行流形提取[11],輸出局部重建網(wǎng)格,程序結束。

在上述步驟(5)中,對Ti進行夾角檢測的具體流程如下:

(1)提取面片Ti的頂點{v1,v2,v3};

(2)檢測每個頂點的標記,若三個頂點中存在輔助點,則該面片未通過檢測,程序結束;

(3)獲取面片Ti對偶Voronoi邊的端點b1、b2;

(4)對于頂點v1,計算Voronoi單元中與其距離最遠的Voronoi頂點,并將由v1指向該頂點的向量記為Vv1;

(5)計算v1b1和v2b2與Vv1的夾角θ1和θ2,其中θ1=∠(v1b1,Vv1),θ2=∠(v1b2,Vv1)。對于給定閾值θp,若θ1與θ2均小于θp,或θ1與θ2均大于π-θp,則該面片未通過檢測。對v2和v3進行相同操作。

局部樣本λ(p)的重建效果如圖1所示。

圖1 曲面局部重建結果

2 基于形心權重的樣點法向估計

基于三角網(wǎng)格估計樣點曲率時通常需要計算樣點法向,且樣點法向估計的準確度對計算曲率的準確度有很大影響。由局部曲面重建獲得的Delaunay網(wǎng)格曲面包含著樣點之間的拓撲關系,可基于網(wǎng)格頂點的鄰接信息計算樣點法向。對網(wǎng)格曲面,頂點的一階鄰域如圖2所示?？蓪⒁浑A鄰域面法向的加權和表示為頂點法向,即樣點法向。在分析網(wǎng)格頂點一階鄰域面幾何特性的基礎上,提出一種一階鄰域形心權重,使樣點法向計算更為穩(wěn)定與精確。網(wǎng)格曲面上的頂點法向計算公式可表示為

(1)

式中:m為頂點一階鄰域面的數(shù)量;ck為鄰域面fk的權重值;nfk表示鄰域面fk的法向,可表示為

k=1,2,…,m;vm+1=v1

(2)

圖2 網(wǎng)格頂點的一階鄰域

在三角網(wǎng)格曲面中三角面片以形狀為等邊三角形為質量最優(yōu),在同等條件下,三角面片形狀越規(guī)整,它對頂點法向的貢獻越大。為對三角面片形狀規(guī)整程度進行量化,以三角形外接圓圓心與內(nèi)切圓圓心的歐氏距離與外接圓半徑的比值作為形狀評價因子

(3)

式中:Ro、ro分別為三角形的外接圓圓心和內(nèi)切圓圓心;R為外接圓半徑。在[0,1]區(qū)間內(nèi),隨著λ的增大,表示三角面片形狀規(guī)整程度越高。當λ等于0時,三角面片的形狀為退化三角形,即三點共線,可認為該面片對頂點法向沒有影響;λ等于1時,三角面片形狀等邊三角形,此時三角面片規(guī)整度最高。

在一階鄰域面形狀評價因子相同的條件下,鄰域面的尺寸也對頂點法向有一定的影響,尺寸越小的鄰域面幾何性質與頂點處相似性越高,其法向對頂點法向的貢獻越大。因質心可以很好地反映三角面片尺寸和形狀,所以以鄰域面的質心到頂點的距離來刻畫鄰域面尺寸這一影響因素,記為

γ=d(pc,vi)

(4)

式中:d(·)表示兩點之間的歐氏距離;pc表示鄰域面fk的質心;vi表示三角網(wǎng)格頂點。

由形狀和尺寸可確定一個鄰域面,故頂點一階鄰域面fk權重值可表示為

式中:λk為fk的形狀評價因子;γk為頂點vi到fk質心的距離。故樣點法向即頂點法向可表示為

(5)

基于網(wǎng)格頂點一階鄰域面形心權重的樣點法向計算方法考慮了鄰域面的形狀和尺寸對頂點局部幾何性質的綜合影響,所得法向更為精確。

3 樣點曲率穩(wěn)健估計

曲率估計過程包含初步計算和加權修正兩個階段。在獲得局部Delaunay網(wǎng)格曲面后,基于網(wǎng)格中樣點一階鄰域面估計樣點曲率,通過加權計算對其修正,提高樣點曲率估計的穩(wěn)健性。

對于局部網(wǎng)格中樣點pi,其平均曲率和高斯曲率的估計公式[8]分別為

(6)

(7)

式中:αij、βij分別為連接樣點pi和pj的對角;vij為pi指向pj的向量;N(i)為樣點pi的一階鄰域點集;npi為頂點pi的法向,按式(5)計算所得;θj表示鄰域三角形中以為頂點的角度值;Amax為網(wǎng)格曲面中pi所在Voronoi單元面積;?M表示局部網(wǎng)格邊界樣點集合。局部網(wǎng)格中樣點曲率估計示意圖如圖3所示。

圖3 局部網(wǎng)格中樣點曲率估計示意圖

獲取樣點曲率估計結果后,計算局部網(wǎng)格曲面中鄰域點與目標樣點間測地距離,基于所得測地距離確定不同鄰域點權重系數(shù),對目標樣點曲率進行加權修正,最終實現(xiàn)樣點曲率的平滑過渡。對于目標樣點p,其平均曲率修正結果可表示為

(8)

式中:di表示在局部網(wǎng)格所構成圖結構中鄰域樣點pi至目標樣點p的最短路徑長度,可近似表示該兩點間的測地距離;h為目標樣點p至所有鄰域樣點中最短路徑長度的最大值;G(x)為核函數(shù),表達式為

(9)

在式(8)中,將kH(p)替換為kG(p),即可獲得高斯曲率估計結果。

在局部重建過程中,為散亂點集建立合理的索引結構可顯著提高局部樣本獲取效率。目前有多種用于提高目標樣點k近鄰查詢效率的空間索引,如空間八叉樹[12]、k-d樹[13]、R樹[14]以及R*樹[15]等。鑒于k-d樹優(yōu)越的空間索引性能和便于實現(xiàn)的特點,本文將其作為散亂點云的空間索引,實現(xiàn)鄰域點集的快速查找。

綜上所述,對S中任一樣點p,其曲率估計完整流程如下:

(1)基于k-d樹索引查詢p的鄰域點集并進行局部重建,獲得局部重建網(wǎng)格D(φ(p))。

(2)基于樣點一階鄰域形心權重,按式(4)計算D(φ(p))中任一樣點法向。

(3)基于式(5)與式(6)估計D(φ(p))中任一樣點的曲率。

(4)采用Dijkstra算法[16]計算目標樣點p至鄰域點集中任一樣點的近似測地距離,并利用式(7)求解p點曲率。

4 實驗與分析

為驗證本文所提算法的有效性,在硬件配置為HP xw8600 Workstation(2.5 GHz,4.0 GB內(nèi)存),操作系統(tǒng)為GNU/Linux的測試環(huán)境中,分別利用文獻[3]算法、文獻[5]算法、文獻[8]算法與本文算法對如圖4所示的不同點云模型進行測試,測試結果如圖5所示。

(a)安雅(點數(shù)為1.02×106) (b)彌勒佛(點數(shù)為1.49×106) (c)維納斯(點數(shù)為2.50×106)

(d)風扇盤(點數(shù)為0.21×106) (e)引擎罩(點數(shù)為0.62×106) (f)龍(點數(shù)為0.00×106)

將不同算法所得高斯曲率進行歸一化處理,采用顏色索引表示曲率變化,可得模型高斯曲率云彩圖,如圖5所示。對于圖4e所示引擎罩模型,文獻[8]算法曲率估計結果與其他3種算法計算相比存在一定誤差,這是因為文獻[8]算法易受采樣誤差的影響,而其他3種算法對采樣誤差均具有一定的抑制作用。通過對比可知,本文算法與文獻[5]算法所得曲率過渡較為光滑,尤其在曲率變化較大特征處,本文算法所得曲率過渡最為光滑。

圖5 引擎罩點云數(shù)據(jù)不同算法高斯曲率估計結果

平均曲率是模型的重要幾何信息,可作為型面特征度量標準應用于點云簡化過程。為驗證4種算法對平均曲率估計的準確性,以不同算法所得平均曲率為輸入數(shù)據(jù),對圖4d所示的非均勻采樣的風扇盤模型點云數(shù)據(jù)進行簡化。簡化后點云數(shù)據(jù)規(guī)模為原始數(shù)據(jù)的30%。圖6a所示為基于文獻[3]算法所得曲率的精簡結果,精簡后的點云分布較為均勻,特別在尖銳特征區(qū)域,未能保留足夠多的點云數(shù)據(jù)。如圖6b所示,基于文獻[8]算法所得曲率進行精簡的結果在尖銳特征處樣點數(shù)據(jù)較基于文獻[3]算法的精簡結果有所增加。如圖6c和圖6d所示,文獻[5]算法和本文算法精簡所得點云數(shù)據(jù)在特征區(qū)域分布較為密集,并且通過對比可知,本文算法精簡結果隨曲率變化的趨勢最為明顯。

(a)文獻[3]算法 (b)文獻[8]算法

(c)文獻[5]算法 (d)本文算法

為對本文算法曲率估計結果進行定量分析,構造半徑為10 mm的球面的散亂點云數(shù)據(jù)(點數(shù)為104),以不同算法計算球面樣點的高斯曲率和平均曲率。重復10次實驗,分別統(tǒng)計各算法所得平均曲率高斯曲率與標準曲率的相對誤差均值ηH、ηG和曲率標準差σH、σG。ηH、ηG可反映各算法曲率計算的精度,σH、σG可反映各算法曲率計算的穩(wěn)健性,σH、σG值越小,表明曲率估計穩(wěn)健性越好。如表1所示,文獻[8]和文獻[5]算法平均曲率計算精度最高,對于無采樣偏差的點云數(shù)據(jù),文獻[5]算法采用的是文獻[8]所提的Voronoi算法,因此,具有較高的計算精度,本文算法與文獻[8]算法計算精度相當,高斯曲率計算精度高于其他算法且穩(wěn)健性最好。

表1 不同算法所得球面樣點曲率結果比較

為分析各算法對具有采樣誤差的點云數(shù)據(jù)的曲率估計效果,對球面點云數(shù)據(jù)添加均值為0、方差為0～1.0 mm的高斯噪聲,分別統(tǒng)計各算法曲率計算結果的ηH、ηG和σH、σG,結果如圖7所示。從圖7中可以看出:隨著噪聲程度的增大,本文算法的曲率計算精度和穩(wěn)健性均高于其他算法;相比于計算精度較高的Voronoi算法,本文算法在噪聲為0時,計算精度與之相當,隨著噪聲程度的增大,本文算法優(yōu)勢明顯,計算精度和穩(wěn)健性可保持在Voronoi算法的1～2倍。

為驗證本文算法的運行效率,分別采用文獻[3]算法、文獻[5]算法、文獻[8]算法和本文算法對如圖4所示的6組不同采樣數(shù)據(jù)進行曲率估計,實驗結果如圖8所示。由圖8可見,本文算法的運行效率低于基于樣點一階鄰域面進行曲率估計的文獻[8]算法,但明顯高于基于局部擬合的文獻[3]算法。文獻[5]算法通過識別曲率估算異常區(qū)域并對其進行加權拋物面擬合估計離散曲率,其運行效率略低于本文算法。

(a)4種算法的ηH比較

(b)4種算法的ηG比較

(c)4種算法的σH比較

(d)4種算法的σG比較

圖8 不同算法曲率估計時間對比

本文算法所需主要參數(shù)有kη、kζ和θp。kη為局部樣本的樣點數(shù)量,其值大小決定著樣點曲率估計質量?；诟咚购撕瘮?shù)將局部樣本中樣點預估計曲率的加權和作為樣點曲率的估計結果,kη值越大,參與加權計算的樣點數(shù)量越大,樣點曲率過渡越光滑,對于有噪聲以及采樣不均勻等缺陷的點云,其樣點曲率計算越穩(wěn)定,但過大的kη值會導致樣點曲率過于平滑,且計算效率低下;kη值也不宜過小,kη值過小會導致曲面局部樣本不夠完備,無法準確反映曲面局部形狀,導致樣點曲率估計不準確,且無法穩(wěn)健處理有缺陷點云的樣點曲率估計問題。通過大量實驗表明:對于采樣均勻的點云,kη取15時可獲得較為理想的實驗效果。為更為穩(wěn)健地處理非均勻采樣點云的曲率估計問題,在以樣點的15近鄰點集為局部樣本的基礎上對其進行增益優(yōu)化[17],獲得可準確反映曲面局部形狀的局部樣本。kζ為添加輔助點后樣點鄰域點集樣點數(shù)量,通過添加輔助點使得重建所得網(wǎng)格曲面中局部樣本邊緣樣點可以有完整的一階鄰域面,輔助點的數(shù)量不宜過多,否則會導致局部重建效率低下。根據(jù)增益優(yōu)化后局部樣本樣點數(shù)量N和樣點局部采樣密度確定kζ的取值,計算公式如下

(10)

5 結 論

通過對局部樣本進行Delaunay網(wǎng)格剖分并濾除多余面片,可為樣點構建鄰域同構曲面,基于所得網(wǎng)格曲面并結合測地距離進行加權修正估計樣點曲率,可得到如下結果:

(1)局部重建所得Delaunay網(wǎng)格曲面插值于采樣點集,并且當采樣密度符合要求時,拓撲同構于原表面,因此在復雜表面及尖銳特征區(qū)域均可實現(xiàn)理想的估計效果;

(2)通過對局部樣本獲取時所需參數(shù)kη、kζ合理取值,以一階鄰域形心權重計算樣點法向進行曲率預估計,并基于測地距離對曲率估計結果進行平滑修正,能夠有效抑制曲面采樣誤差對曲率估計結果的影響,提高了算法的穩(wěn)健性;

(3)相比于Meyer提出的Voronoi算法,本文算法對采樣精度較高的點云數(shù)據(jù)可保證與其相當?shù)挠嬎憔?對存在噪聲的點云數(shù)據(jù)計算精度和穩(wěn)健性均可提高1～2倍,基于k-d樹構建空間索引獲取局部樣本并進行輔助點添加,使本文算法具有較高的計算效率,適用于大規(guī)模點云數(shù)據(jù)的處理。