王先澤，李忠科，張曉娟，呂培軍，王 勇

（1.第二炮兵工程大學(xué)401教研室，陜西 西安710025；2.北京大學(xué) 口腔醫(yī)學(xué)院，北京100871）

0 引 言

對點(diǎn)云數(shù)據(jù)進(jìn)行簡化是數(shù)據(jù)預(yù)處理中一個非常重要的技術(shù)，它有助于提高數(shù)據(jù)的存儲、傳輸、計算以及顯示的效率。點(diǎn)云的簡化大體上可以分為基于網(wǎng)格的簡化和基于點(diǎn)的簡化?；诰W(wǎng)格的簡化方法需要先將點(diǎn)云網(wǎng)格化且在簡化的過程中需要維持拓?fù)涞囊恢滦?，因此比較復(fù)雜。而基于點(diǎn)的簡化方法則不需要維護(hù)這些拓?fù)湫畔ⅲ鄬碚f更有優(yōu)勢。Pauly等［1］首次將幾種基于網(wǎng)格的聚類簡化方法應(yīng)用到點(diǎn)集曲面上，該方法速度很快，但是產(chǎn)生的誤差較大，且不能像基于網(wǎng)格的簡化方法那樣預(yù)先用一個全局誤差來控制簡化過程，Brodsky［2］等提出的聚類簡化法也有類似問題。鄒萬紅［3］等提出了一種大規(guī)模點(diǎn)云模型的簡化算法，該算法可以分為兩步：首先根據(jù)幾何近似原理將點(diǎn)模型分割成與平面近似的多個分片，然后在分片的基礎(chǔ)上進(jìn)行層次空間聚類簡化，該方法的效果比較好，但是復(fù)雜度較高。王仁芳等［4］提出了一種基于相似性的曲率自適應(yīng)點(diǎn)模型簡化算法，將點(diǎn)模型劃分為強(qiáng)邊性和非強(qiáng)邊性兩部分，首先對非強(qiáng)邊性部分進(jìn)行表面區(qū)域幾何特征相似性聚類，然后對各類簇分別進(jìn)行劃分和簡化，該算法有效地保持了特征邊界部分和曲面的細(xì)節(jié)，且能夠生成高質(zhì)量的簡化點(diǎn)集曲面。

由于隱式曲面表示在曲面造型和和三維可視化中的良好特性［5］，徑向基函數(shù)［6］在點(diǎn)云簡化之中已經(jīng)得到廣泛應(yīng)用，但是它的計算效率非常低，只能處理小的點(diǎn)云數(shù)據(jù)。而緊支徑向基函數(shù)（CSRBF）［7－9］能夠用于大量點(diǎn)云數(shù)據(jù)的重構(gòu)且在隱式表示方面有和徑向基函數(shù)相似的特性，本文在此基礎(chǔ)上，提出一種特征保持的基于緊支徑向基函數(shù)的點(diǎn)云簡化方法。該方法通過建立與點(diǎn)曲率和支撐半徑相關(guān)的評估函數(shù)來評價點(diǎn)的重要性，評估函數(shù)值會隨著支撐半徑內(nèi)點(diǎn)的數(shù)量而變化，然后迭代地更新并刪除函數(shù)值最小的點(diǎn)。實驗結(jié)果表明，該方法能夠很好地控制簡化后點(diǎn)云的數(shù)量，產(chǎn)生的誤差較小，且能夠很好的保持點(diǎn)云的特征。

1 緊支徑向基函數(shù)

1.1 概 念

式 中：λi——第i點(diǎn)影響的權(quán)重，φσ（r）＝φ（r／σ），φσ（r）——緊支徑向基函數(shù)。σ——支撐半徑，也就是點(diǎn)pi所能影響點(diǎn)的范圍，這是一個與采樣密度相關(guān)的參數(shù)。

1.2 確定支撐半徑

采用CSRBF進(jìn)行點(diǎn)云簡化時，支撐半徑的大小決定了落在某確定散亂點(diǎn)作用范圍內(nèi)的點(diǎn)的數(shù)量，直接影響點(diǎn)云簡化以及后續(xù)曲面重建的效果。對于稀疏的點(diǎn)云，如果支撐半徑σ選取得太小，則散亂點(diǎn)之間的相互作用減小，達(dá)不到點(diǎn)云簡化的目的；對于密集的點(diǎn)云，如果支撐半徑σ選取得太大，則在支撐半徑中的點(diǎn)的數(shù)量會增加，這就會增加算法的求解時間；另外，對于一些密度不規(guī)則的點(diǎn)云，支撐半徑的選取會直接影響到重建模型的質(zhì)量。

本文參考文獻(xiàn)［10］中的方法來確定算法中緊支徑向基函數(shù)的支撐半徑大小，大體上可以分為三步：①使用基于八叉樹的方法對點(diǎn)云p的包圍盒進(jìn)行分割；②計算葉子的平均對角線長度diagonal；③設(shè)置σ值為diagonal的3／4，即σ＝0.75×diagonal。

1.3 緊支徑向基函數(shù)選擇

緊支徑向基函數(shù)應(yīng)用非常廣泛，迄今為止，為了解決不同問題以及不同的連續(xù)性劃分，人們提出了許多種類的基函數(shù)。Wendland為了解決徑向基函數(shù)的最小度的對稱正定問題，提出的緊支徑向基函數(shù)都有如下形式

針對不同的維數(shù)d以及不同的連續(xù)性Ck，Wendland提出了一些緊支徑向基函數(shù)，表1列出了三維情況下的幾個函數(shù)。

表1 Wendland提出的緊支徑向基函數(shù)

文獻(xiàn)［10］中用上述4種緊支徑向基函數(shù)對不同采樣密度的一些模型進(jìn)行了重構(gòu)，并且將重構(gòu)后的模型與原始點(diǎn)云進(jìn)行了誤差分析。實驗結(jié)果表明使用C2連續(xù)的緊支徑向基函數(shù)進(jìn)行重構(gòu)得到的曲面質(zhì)量比其它3種所得到的更好。因此，在本章中，我們使用Wendland的C2連續(xù)的緊支徑向基函數(shù)作為算法中使用的函數(shù)。

2 特征保持的點(diǎn)云簡化算法

本文提出了一種基于緊支徑向基函數(shù)的特征保持的三維散亂點(diǎn)云簡化算法，算法基本步驟描述如下：

步驟1 假設(shè)點(diǎn)p為點(diǎn)云中任意一點(diǎn)，則根據(jù)改進(jìn)的k＿近鄰算法求點(diǎn)云中距離點(diǎn)p中最近的k個頂點(diǎn)即求點(diǎn)p的k鄰域；

步驟2 采用移動最小二乘法（moving least squares，MLS）在局部坐標(biāo)系內(nèi)擬合點(diǎn)p的鄰域，然后通過擬合得到的局部曲面擬合多項式計算該點(diǎn)的高斯曲率；

步驟3 根據(jù)2.3節(jié)的描述選擇合適的緊支徑向基函數(shù)，并在此基礎(chǔ)上建立評估函數(shù)Eva（i），用以評估頂點(diǎn)的重要性，然后將步驟2中求得的高斯曲率值代入，計算得到評估函數(shù)Eva（i）的值；

步驟4 重復(fù)步驟1到步驟3，直到所有點(diǎn)的高斯曲率與評估函數(shù)值都求出，然后將Eva（i）值最小的點(diǎn)刪除；

步驟5 判斷刪除后點(diǎn)的數(shù)目是否達(dá)到預(yù)先的要求，如果達(dá)到，則執(zhí)行步驟7，否則，更新與刪除點(diǎn)相關(guān)的點(diǎn)的評估函數(shù)值，然后刪除所有函數(shù)值中值最小的點(diǎn)；

步驟6 迭代執(zhí)行步驟5；

步驟7 算法結(jié)束。

2.1 改進(jìn)的k＿近鄰算法

當(dāng)點(diǎn)云為非均勻采樣時，傳統(tǒng)的k＿近鄰算法來具有一定的局限性。如圖1（a）所示，傳統(tǒng)的k＿近鄰算法求得的k鄰域只包含了該點(diǎn)周圍一半的環(huán)境，不能充分描述周圍的形勢，從而可能使求得的曲率發(fā)生偏差?；诖耍疚奶岢鲆环N改進(jìn)的k＿近鄰算法，該算法通過設(shè)置一個最大角限制條件來保證點(diǎn)p的k鄰域圍繞著該點(diǎn)是球形分布的，從而在最大程度上保證所求得的k鄰域能夠反應(yīng)點(diǎn)所處的環(huán)境。

圖1 選擇點(diǎn)P的鄰域所需的角度限制

使用改進(jìn)的k＿近鄰算法計算輸入點(diǎn)p的k鄰域NNp，大體上可以分為三步：①找到k個鄰域點(diǎn)加入NNp中，通常8≤k≤18在大多數(shù)情況下能滿足需要；②將這k個鄰域點(diǎn)映射到一個合適的局部支撐平面H上，該平面可以通過求局部法向量來確定；③確保映射到支撐平面H上的兩個相鄰鄰域點(diǎn)與中心點(diǎn)之間形成的角度不大于預(yù)先設(shè)定的門限值，如果不滿足條件，則繼續(xù)增加鄰域點(diǎn)到NNp中，直到找到滿足條件的替代鄰域點(diǎn)為止。

首先是使用傳統(tǒng)的k＿近鄰算法計算點(diǎn)p的k鄰域，并加入到NNp中。為了估計以cNNp為中心的鄰域點(diǎn)的局部法向量，可以先計算如式（1）所示的半正定權(quán)重協(xié)方差矩陣C＝｛cvw｝∈R3×3

式中：θj——以si為中心的二次B樣條B（t）：θj＝，d——?dú)W式距離，r——支撐半徑，表示以sisi為中心的包圍球的半徑，包圍球包含對輸入點(diǎn)pj使用k＿近鄰算法求得的鄰域點(diǎn)；k、cNNp——鄰域點(diǎn)的個數(shù)、中心，k＝｜NNp｜。

求得協(xié)方差矩陣C后，對C進(jìn)行特征分析可以得到k鄰域NNp的主成分，假設(shè)λ1，λ2，λ3分別為C的3個實特征值，且有λ1≤λ2≤λ3，而v1，v2，v3分別為它們所對應(yīng)的一組正交特征向量，則特征值λ1對應(yīng)的特征向量v1就是局部法向量。估計出局部法向量就可以求得局部支撐平面H，它是一個垂直于局部法向量且過點(diǎn)p的平面，然后把k鄰域NNp中的點(diǎn)映射到平面上。

將NNp中的點(diǎn)映射到局部支撐平面上之后，下一步是確保映射到支撐平面H上的兩個相鄰鄰域點(diǎn)與中心點(diǎn)之間形成的角度不大于預(yù)先設(shè)定的門限值max＿angle。在支撐平面上，對于點(diǎn)p，我們考慮它所有的鄰點(diǎn)qi以及按逆時針方向搜索的qi的相鄰點(diǎn)qi＋1，需要它們形成的角∠qipqi＋1不大于max＿angle，比如π／2，如圖1（b）所示。如果在這些鄰域點(diǎn)和中心點(diǎn)p之間形成的角度中沒有大于門限值的，則這些鄰域點(diǎn)就是使用改進(jìn)的k＿近鄰算法求得的該點(diǎn)的k鄰域NNp，算法結(jié)束；否則，就需要增大搜索的包圍球半徑radii，從而增加鄰域，然后直接映射到支撐平面上，繼續(xù)判斷角度，更新NNp，直到都滿足門限條件的替代鄰域點(diǎn)被找到。

2.2 高斯曲率計算

使用改進(jìn)的k＿近鄰算法得到點(diǎn)的k鄰域后，就可以通過MLS方法擬合pi點(diǎn)的k鄰域NNp（pi）＝｛pj｝，j＝0，…，k中的點(diǎn)，得到pi點(diǎn)的局部擬合多項式gi，然后計算pi點(diǎn)的曲率。計算局部曲面擬合多項式可以分為兩步：①尋找支撐平面，建立局部坐標(biāo)系；②計算多項式。

首先，建立NNp（pi）的協(xié)方差矩陣，B＝其 中，是NNp（pi）的重心，B是一個半正定的3×3的對稱矩陣。同2.1節(jié)中的方法類似，對矩陣B進(jìn)行特征分析，v2，v3，oi共同定義逼近NNp（pi）的支撐平面Hi，v1即為oi的法向，記為ni。在Hi平面上建立局部坐標(biāo)系（oi，u，v，ni），oi為原點(diǎn)，ni為Z軸。然后在（oi，u，v，ni）坐標(biāo)系內(nèi)用二次多項式gi擬合NNp（pi），使pj點(diǎn)與其在擬合曲面上投影點(diǎn)的距離平方和最小。設(shè)gi表達(dá)式為

則通過對多項式的各個系數(shù)求一階導(dǎo)數(shù)可以得到6個線性方程，然后利用高斯消元法計算得到gi的各個系數(shù)。求得gi的表達(dá)式后，參照文獻(xiàn)［11］中的方法，可得曲面的高斯曲率kg為

Dragon和Armadillo模型的高斯曲率如圖2所示。其中深色部分為曲率較大的區(qū)域，淺色部分則為曲率較小的區(qū)域。

2.3 建立評估函數(shù)

本文提出的簡化算法的核心思想就是建立一個判斷點(diǎn)云中的點(diǎn)重要性的標(biāo)準(zhǔn)，通過這個標(biāo)準(zhǔn)來計算點(diǎn)的重要性，而簡化就是迭代地刪除最不重要的點(diǎn)的過程。反映一個點(diǎn)在點(diǎn)云中的重要性有兩個方面：一是該點(diǎn)的曲率，高曲率區(qū)域是特征區(qū)域，能夠反映模型的特征；二是該點(diǎn)周圍點(diǎn)的分布情況，分布密度越大，表示反映該區(qū)域特征的點(diǎn)越多，則該點(diǎn)的重要性就會降低，反之亦然。因此，我們建立評估函數(shù)就從這兩個方面考慮。點(diǎn)的曲率用上節(jié)中的高斯曲率來表示，而點(diǎn)的分布情況則可以通過緊支徑向基函數(shù)反映出來。

圖2 Dragon和Armadillo模型的高斯曲率

假設(shè)p＝為輸入點(diǎn)集，我們定義評估函數(shù)為

其中，φσ（r）＝φ（r／σ），φ（r）＝（1－r）4＋（4r＋1）為前面介紹的緊支徑向基函數(shù)。σ為支撐半徑，可以通過點(diǎn)云p的密度估計出來；Ni表示pi對鄰域貢獻(xiàn)的權(quán)重；kgi表示式（2）中計算得原始點(diǎn)云中pi點(diǎn)的高斯曲率，由于曲率值可能為負(fù)，所以在這里取它絕對值｜kgi｜；ε＞0為一個小的正實數(shù)，用來避免式（3）中分母為0。由于使用Eva是來評估點(diǎn)的重要性，所以ε的絕對取值并不影響比較的結(jié)果。在我們的實驗中，選取ε＝1，當(dāng)某點(diǎn)的σ半徑內(nèi)不包含鄰點(diǎn)時，Eva值由曲率決定。

當(dāng)刪除一個點(diǎn)后，所有支撐半徑σ內(nèi)包含被刪除點(diǎn)的點(diǎn)的σ半徑內(nèi)的點(diǎn)的個數(shù)就會減少，則也會減小，而高斯曲率｜kgi｜值不變，從而導(dǎo)致這些點(diǎn)的Eva（i）值增加。更新了評估函數(shù)Eva（i）的值之后，再一起比較，迭代刪除值最小的點(diǎn)，直到達(dá)到預(yù)先的簡化要求為止。

2.4 算法實現(xiàn)

本文所提簡化算法的偽代碼如下所示：

3 算法試驗及分析

由于在評估函數(shù)中涉及到點(diǎn)的曲率，而曲率計算與k＿近鄰的個數(shù)的選取密切相關(guān)。如果k鄰近的個數(shù)選取過大，不能保證單凹或者單凸，曲率估算就會出現(xiàn)誤差；如果k鄰近的個數(shù)取得過小，則不能被二次曲面逼近，且容易受噪聲影響，同樣不能正確估算誤差。在本算法中，取k＝10取得了不錯的效果。

為了驗證本文算法的有效性，在Windows平臺上（CPU主頻2.0G，內(nèi)存1G），基于VC＋＋2008和OpenGL實現(xiàn)了本文算法，并將實驗結(jié)果與文獻(xiàn)［2］中所提方法的結(jié)果進(jìn)行了比較，Armadillo模型的部分效果對比如圖3所示。由于本文算法使用迭代刪除最不重要的點(diǎn)的方法，能夠很好地控制簡化后點(diǎn)的數(shù)目，從而可以保證在相同簡化率的條件下與其它算法進(jìn)行比較。從圖3中可以看出，本文算法在一些特征點(diǎn)較多的部位保留了較多的點(diǎn)，而在其它一些部位則保留點(diǎn)比較少，相對于文獻(xiàn)［2］中方法來說，能夠更好地保持模型的特征。另外，本文也使用尖峰信噪比PSNR［10］（peak signal to noise ratio）進(jìn)行誤差分析，PSNR的定義如下：PSNR［dB］＝其中，peak為模型包圍盒的對角線長度，d為點(diǎn)到模型的距離。表2中描述了一些模型的誤差與運(yùn)行時間對比。從表2中可以看出，本文算法運(yùn)行效率比較低，但是尖峰信噪比高，簡化之后重構(gòu)出來的模型質(zhì)量好。

圖3 對于Armadillo模型本文算法和文獻(xiàn)［2］算法簡化效果比較

表2 簡化后PSNR與運(yùn)行時間對比

4 結(jié)束語

本文主要解決散亂點(diǎn)云數(shù)據(jù)的簡化問題，主要思想是采用迭代簡化的策略，通過建立一個與點(diǎn)的曲率和基函數(shù)支撐半徑內(nèi)點(diǎn)云密度相關(guān)的評估函數(shù)來評估點(diǎn)的重要性，然后更新并迭代刪除評估函數(shù)值最小的點(diǎn)。實驗結(jié)果表明，該方法能夠精確地控制簡化后點(diǎn)云的數(shù)量，產(chǎn)生的誤差比較小，且能夠較好地保持點(diǎn)云的特征，適合對簡化后質(zhì)量要求高的點(diǎn)云。

［1］Pauly M，Gross M，Kobbelt L P.Efficient simplification of point－sampled surfaces［C］／／Boston：Proceedings of IEEE Visualization，2002：163－170.

［2］Brodsky D，Watson B.Model simplification through refinement［C］／／Montreal：Proceedings of Graphics Interface，2000：221－228.

［3］ZOU Wanhong.Techniques for shape modeling from large scale point clouds［D］.Zhejiang：Zhejiang University，2007（in Chinese）.［鄒萬紅.大規(guī)模點(diǎn)云模型幾何造型技術(shù)研究［D］.浙江：浙江大學(xué)，2007.］

［4］WANG Renfang，ZHANG Sanyuan，YE Xiuzi.Similaritybased simplification of point－sampled surfaces［J］.Journal of Zhejiang University，2009，43（3）：448－454（in Chinese）.［王仁芳，張三元，葉修梓.基于相似性的點(diǎn)模型簡化算法［J］.浙江大學(xué)學(xué)報，2009，43（3）：448－454.］

［5］Craig Schroeder，Wen Zheng，Ronald Fedkiw.Implicits sur－face tension formulation with a Lagrangian surface mesh on an Eulerian simulation grid［J］.Journal of Computational Physics，2011（1）：1－27.

［6］Bjorn Andersson，Stefan Jakobsson，Andreas Mark，et al.Modeling surface tension in SPH by interface reconstruction using radial basis functions［C］／／Manchester：The 5th International SPHERIC Workshop，2010.

［7］TAN Yehao，JIANG Zhifang，DU Xiaoliang，et al.Visualization of multi－dimensional data with interpolation based on compactly supported radial basis functions［J］.Computer Engineering and Applications，2010，46（9）：220－223（in Chinese）.［譚業(yè)浩，蔣志方，杜曉亮，等.緊支徑向基函數(shù)插值實現(xiàn)多維數(shù)據(jù)可視化［J］.計算機(jī)工程與應(yīng)用，2010，46（9）：220－223.］

［8］Wu Jinming，Lai Yisheng，Zhang Xiaolei.Radial basis functions for shape preserving planar interpolating curves［J］.Journal of Information ＆ Computational Science，2010，7（7）：1453－1458.

［9］TIAN Jianlei，LIU Xumin，GUAN Yong.Implicit surface algorithm for holding features［J］.Computer Engineering and Applications，2011，47（1）：208－210（in Chinese）.［田建磊，劉旭敏，關(guān)永.一種保特征的隱式曲面算法［J］.計算機(jī)工程與應(yīng)用，2011，47（1）：208－210.］

［10］Masaki Kitago，Gopi M.Efficient and prioritized point subsampling for CSRBF compression［C］／／Boston：Eurographics Symposium on Point－Based Graphics，2006：121－128.

［11］PANG Xufang，PANG Mingyong，XIAO Chunxia.An algorithm for extracting and enhancing valley－ridge features from point sets［J］.Acta Automatica Sinica，2010，36（8）：1073－1083（in Chinese）.［龐旭芳，龐明勇，肖春霞.點(diǎn)云模型谷脊特征的提取與增強(qiáng)算法［J］.自動化學(xué)報，2010，36（8）：1073－1083.］