大規模點云選擇及精簡

范 然， 金小剛

（浙江大學CAD&CG國家重點實驗室，浙江 杭州 310058）

逆向工程中最普遍的應用模式是利用基于光學原理的掃描設備測量零件或模具外表面形成點云數據，從中提取幾何特征進而重建多邊形或 NURBS曲面[1]。隨著掃描硬件精度的提升，原始掃描數據能夠較好逼近物理模型的曲面形態，但是，刪減處理始終是不可或缺的關鍵步驟。根據被測物體的擺放方式，掃描點云中不可避免地包含如墻面、支撐物等非目標背景數據，解決這一問題需設計能夠高效并靈活指定大規模點云中待刪除區域的選擇算法。此外，原始掃描點云通常由多次測量的單片數據拼合而成，接合處存在重疊冗余采樣，而且數據規模大、分布不均勻，普通計算機的計算及內存資源不易直接對其進行曲面重建。針對這一問題，點云精簡算法的目標是減少數據量并使采樣點依幾何特征自適應均衡分布。

由于常用計算機輸入方式是二維指針設備，二維套索UI是點云選擇最直觀的交互方式。研究領域通常使用簡單空間形體實現點云選擇[2]，并未給出支持套索UI的點云選擇算法。目前，在逆向工程應用領域中多使用Geomagic[3]、Rapidform[4]等商業軟件進行處理，相關算法并未公開，無法滿足研究領域構建低成本三維掃描系統的需求。因此本文提出了一種支持套索UI的大規模點云選擇算法，利用套索形狀的矩形覆蓋與輸入點云八叉樹編碼剔除大部分點在多邊形內判斷。不同于傳統的平面射線奇偶法及其改進算法，本文方法利用了掃描點云密集分布的特點，能夠整體判斷點云的局部區域。

點云精簡通常有兩種定義[5]：給定容許誤差閾值計算最小采樣數分布；給定目標采樣點數目搜索最小誤差分布。在逆向工程中后者相對合理，用戶通常希望點云精簡算法同時具備如下特性：1）任意指定目標采樣點數目；2）不改變原有采樣點位置以維持掃描精度；3）盡可能保留尖銳幾何特征，依曲面曲率自適應均衡分布。然而，傳統基于曲率自適應[6]以及基于QEM誤差[7]的點云精簡方法具有采樣點局部聚集以及不能保持尖銳邊特征與邊界的缺點。本文利用Poisson-disk采樣原理獲得簡化后點云均衡分布的效果，并以采樣點鄰域球布爾交運算來定義曲面上的圓盤半徑度量、擴展采樣邊界，獲得了尖銳邊特征及邊界保持的性質。

1 相關工作

點云選擇算法涉及如何判斷采樣點在視平面上的投影與套索多邊形之間的包含關系。Eric[8]總結了若干經典的點在多邊形內測試方法，包括射線奇偶、角度和、三角形扇等算法。在地理信息與圖形學領域，目標應用通常需實時測試大規模數據，經典方法無法滿足速度指標，后續研究采用基于區域分解的分治策略進行優化。Li[9]等分解多邊形為凸子域，將點在任意多邊形內測試轉換為更為簡單的點在凸多邊形內測試。Zalik[10]等將多邊形所在區域分解為均勻單元格并賦予內、外、邊界屬性，首先進行簡單的點在單元格內測試，最后僅對包含于邊界單元格的點執行完整測試。Yang[11]等利用近似最近點概念提高了均勻單元格方法的準確度。本文注意到現存方法僅注重分解或降低多邊形復雜性然后進行逐個單點測試，針對大規模掃描點云數據密集分布的特點，將輸入點云組織為八叉樹層次結構，以局部區域為單元判斷包含關系，剔除大部分點在多邊形內測試。

隨著基于點的繪制與造型技術的發展，已經存在許多算法簡化點云數據。根據精簡點云采樣分布形成的方式可將現有方法分為3類：迭代最優剔除、層次聚類、曲面重采樣。不同的方法分別側重減少原始點云與精簡后點云之間的距離、曲率自適應、分布質量等方面。迭代最優剔除方法的優點在于采樣前后的點云距離誤差小，缺點在于隨著點云規模的增加由于全局排序、屬性更新需要大量的內存與計算消耗，不適合逆向工程中海量點云高效精簡，且不容易保留尖銳邊特征與邊界。Pauly[7]等將三角網格的QEM簡化算法推廣到點云上。Linden[12]提出一種綜合的采樣點影響度，每次移除最小影響度的采樣點。Song[13]等首先識別出尖銳特征采樣點并予以保留，其余采樣點依次按照影響度排序迭代去除。層次聚類方法具有計算效率高的優點，其缺點是不容易控制采樣點的分布與誤差。Lee[6]等采用非均勻八叉樹對點云進行聚類并以法向變動控制八叉樹細分，效率較高，但精簡點云分布不均勻。Yu[14]等利用k-means算法對點云進行層次分割，采用局部協方差分析控制聚類過程，雖然分布質量提高但無法精確指定采樣點數目。Shi[15]等同樣利用k-means算法獲得了類似的效果。Song[16]等采用全局聚類方法具有更加優良的分布效果，但由于其全局優化特性計算效率較低。現有的曲面重采樣方法直接從分布特性角度出發，能夠獲得理論上最優采樣點分布，但由于通常需計算流形上的距離或維持采樣點之間的動態平衡，其所消耗的計算資源也是最大的。Moenning[17]等利用Fast Marching方法漸進計算內蘊Voronoi圖以實現最遠點采樣策略。Pauly[7]等人將粒子平衡模擬方法推廣到點云上。

上述3類方法分別體現了優化、數據分析、幾何結構分析的思想，本文的點云精簡方法屬于曲面重采樣，基于Poisson-disk采樣均勻分布采樣點、防止過度集中，利用鄰域球布爾交運算擴展可用采樣邊界，如圖9和圖10所示由于尖銳邊特征以及邊界數據不容易被周圍鄰域球覆蓋從而被密集采樣。與本文目標最接近的方法是Kang[18]等提出的均衡點云采樣方法，不同之處在于基于Poisson-disk采樣的方法能夠通過采樣半徑調整點云分布。

2 支持套索UI的大規模點云選擇算法

點云選擇的套索UI接口與圖像區域選擇類似，用戶旋轉點云模型至合適方位后草繪套索多邊形、圈選目標區域，如圖1(a)和圖1(b)所示，遴選出投影在套索多邊形之內的點云區域。最基本的實現方式是將每個采樣點投影在視平面上逐個進行點在多邊形內測試，然而對于百萬級甚至千萬級的掃描點云數據，該方法無法獲得理想的交互性能。

圖1 套索UI點云選擇

本文算法受到Zalik[10]等方法的啟發，該方法主要針對地理信息應用，不適合直接用于選擇掃描點云數據。首先，該方法只能處理簡單曲線，在點云區域選擇中，限制用戶輸入的套索形狀為簡單曲線并不合理。本文提出的算法構造了輸入多邊形的矩形區域覆蓋并利用改進的射線奇偶法賦予內、外、邊界屬性，圖2比較了Zalik[10]等方法與本文方法處理非簡單曲線的差異。

圖2 輸入非簡單曲線，從紙張點云模型中選擇數據

此外，本文注意到掃描點云數據通常是密集分布的，因此，將原始點云嵌入八叉樹層次結構，利用節點單元在視平面上投影與矩形覆蓋之間的幾何包含關系剔除大部分點在多邊形內測試。

2.1 套索多邊形矩形覆蓋

給定多邊形P，由邊列表{E1,E2,… ,En}構成，其中Ei連接頂點(xi,yi)和 (xi+1,yi+1)。設C為任意矩形，定義P的矩形覆蓋為{Ci}，要求矩形之間無交集且所有矩形的并集覆蓋了P所在區域，即：DC(P) = {Ci|∪iCi?Pand ?i≠jCi∩Cj=?}。本文利用四叉樹計算DC(P)，主要有兩個優點：能夠根據多邊形P的距離場自適應確定矩形數量；索引矩形計算量較小。設c為四叉樹節點單元，其位置由坐標函數x0(c),x1(c),y0(c),y1(c)確定，點(x,y)∈c：x0(c) ≤x＜x1(c)與y0(c) ≤y＜y1(c)。DC(P)由四叉樹的所有末端子節點構成。根據與P的相對幾何位置分別賦予矩形單元內部、外部、邊界屬性，其中邊界屬性是指P穿越該矩形，此外矩形單元的大小取決于相對于P的距離，即遠離P的矩形較大，邊界矩形較小。DC(P)由遞歸分裂四叉樹節點構建而成，根節點是覆蓋P的正方形。節點分裂由以下形狀誤差決定：

其中ci是四叉樹節點，vij是其四個頂點之一，dist(vij)表示到P的距離函數。由于距離函數在P內部中軸附近存在一階不連續使得線性逼近誤差較大、節點分裂過多，為減少包含測試計算量通常需以盡可能少的矩形分別覆蓋P的內、外、邊界區域。因此(1)中顯式地判斷了ci是否具有邊界屬性。DC(P)的復雜性取決于P以及分裂閾值Te，該閾值由用戶指定或者根據P的最小線段長度即P的分辨率確定，當err(ci)＞Te時分裂ci為4個子節點。

計算上述矩形覆蓋需判斷ic是否具有邊界屬性，如果P是凸多邊形則可以通過對ic四個頂點執行點在多邊形內測試得知其是否具有邊界屬性。對于一般的多邊形，如果ic的4條邊與P沒有交點，根據連續性推斷P沒有穿越ic即不具有邊界屬性。本文采用結合單元邊交點偵測的射線奇偶算法判斷單元的內、外、邊界屬性，如圖3所示該方法能夠正確地賦予邊界屬性。

對于ci的每個頂點分別發射水平與豎直方向且與ci的邊部分重合的射線，統計射線與P的交點數目并偵測交點是否落在ci的邊上。以從ci的左上角頂點 (x0(ci),y1(ci))發射水平射線rright為例，對于所有Ej∈P，判斷與rright的穿射關系。

圖3 基于四叉樹的矩形覆蓋，藍色、綠色、紅色分別表示外部、內部、邊界屬性

在DC(P)計算完成后，為了在后續測試中減少可能的求交運算，將組成P的所有線段嵌入矩形單元中。本文采用Lastra[19]等的射線穿透四叉樹訪問算法。設iE∈P，以t為參數可表示為：

其中(dx,dy)是單位方向向量。首先確定包含起始頂點(xi,yi)的末端節點cs，移除包含在cs中的Ei部分、更新ts，同時將Ei記錄在cs中；對于cs的父節點cp計算其4條邊所在直線與Ei方向射線r的交點參數：

并按照下式斷定Ei是否穿透cp：

對于滿足式(5)的cp，考察其所有子節點并移除Ej中被包含的部分，交替進行以上兩步直至到達終止頂點 (xi+1,yi+1)。

2.2 點在多邊形內測試

掃描點云具有密集分布的特點，可利用其空間一致性剔除大部分點在多邊形內測試。本文將點云嵌入在八叉樹G內，每個末端子節點不超過64個采樣點。從根節點開始，將節點立方體gi∈G沿視線方向投影為側影輪廓多邊形gi′，其邊界由同時鄰接可見面與不可見面的立方體邊投影而成，可見性由面法向與視線向量內積確定。判斷gi′與DC(P)矩形之間的幾何包含關系，分3種情況處理：

1）gi′包含在單個非邊界矩形內的，將gi內所有采樣點標識為與該矩形同樣的屬性；

2）gi′的所有頂點皆落在內部或者外部區域，偵測gi′的邊與P的交點，采用Feito[20]等方法，無需計算出交點，并且已將P的邊嵌入在DC(P)中因此只需做局部判斷，如果沒有交點將gi內所有采樣點標識為gi′所在區域的屬性；

3）gi′與具有邊界屬性的矩形有交集，遞歸處理gi′的子節點；

如果仍有末端節點gi′具有邊界屬性，將gi內的采樣點投影在視平面上，如果落于具有內、外屬性的矩形內，則賦予同樣的屬性。最后對落于邊界矩形內的采樣點執行基于射線奇偶的點在多邊形內測試。經過上述層次測試步驟只有很少的一部分采樣點需執行完整的點在多邊形內測試，因而能夠提升大規模點云區域選擇的效率。

帶預處理階段的包含測試策略其核心在于將針對復雜多邊形的測試轉化為針對一組簡單多邊形的測試，這一策略有利于快速判斷大量采樣點針對同一多邊形的位置。不同于之前的工作，本文基于點云密集分布、具有空間一致性的特性，將點云組織為八叉樹結構，利用投影八叉樹節點單元與多邊形矩形覆蓋之間的包含關系剔除了大部分點在多邊形內測試。

3 基于Possion-disk采樣的點云精簡

在點云精簡過程中顯式地控制采樣點分布對于尖銳邊特征保持、邊界保持、平坦與彎曲區域之間的均衡分布至關重要。其他基于減少簡化前后點云間距離、曲率自適應的精簡方法通常會造成采樣點局部聚集、且不能密集采樣尖銳邊特征以及邊界數據。后續應用如三角化、基于點的繪制、形狀約束動畫等的效果非常依賴于點云數據的分布質量。本文首先估計采樣點法向，進而利用Poisson-disk采樣控制采用點的分布，但由于Poisson-disk采樣的控制參數是采樣半徑因而不容易控制采樣點數目，本文對初步采樣結果依采樣點Voronoi鄰域面積進行排序，向稀疏區域添加以及從擁擠區域移除采樣點以達到目標精簡數目。

3.1 法向估計

原始點云來源于不同的掃描技術，可分為帶有法向與缺失法向兩種情況。對于不具備法向信息的數據本文采用Pauly[7]等的方法，對每一個采樣點pi利用局部協方差分析計算其法向ni，pi鄰域的協方差矩陣為：

其中Ni是采樣點pi的鄰域，包括k個最近采樣點，是位置均值。然后進行特征值分解并以其最小特征值向量作為該點的法向：

圖4 法向雙邊濾波

需要指出的是本文算法并不需要統一法向定向，法向僅用于確定采樣點的切平面。采用式(7)計算的法向對于光滑點云數據能夠較好地逼近被測曲面法向，但不能良好表示逆向工程中常見的具有尖銳邊特征的機械模型，而且增加采樣分辨率并不能解決這一問題[21]。因此，本文采用法向雙邊濾波[22]，使得尖銳邊特征兩側法向數據分區域平滑：

3.2 點云Possion-disk采樣

Poisson-disk采樣具有兩種性質，即采樣點符合隨機分布、采樣點之間的距離不小于2r，其中r是每個采樣點所屬圓盤的半徑。飛鏢投擲算法是最直接的生成方法，即隨機生成新的采樣點并判斷該點所屬圓盤是否被已經存在的采樣點圓盤覆蓋。投擲成功率隨著采樣點的增加而降低因而算法效率較低。Dunbar[23]等提出了基于邊界區域增長的平面Poisson-disk采樣算法，如圖5所示，采樣點及其所屬圓盤用紅色表示，可用采樣邊界用藍色表示，初始3個采樣點增長為6個采樣點。采樣過程主要包括兩步：在可用邊界上隨機采樣、更新可用采樣邊界。

圖5 平面Poisson-disk采樣，藍色表示可用采樣邊界

在曲面上進行Poisson-disk采樣需定義適合的距離度量以確定每個采樣點的影響域，Fu[24]等將Dunbar[23]的方法推廣到三角網格曲面上，利用測地距離計算網格曲面上的等距線從而擴展采樣邊界。但是對于掃描點云數據，拓撲信息缺失、分布不均勻使得很難直接在點云上計算測地距離。本文的關鍵想法是基于Klein[25]等提出定義點云曲面的球影響圖方法，對采樣點ip定義半徑為2r的鄰域球Bi，所有采樣點鄰域球的集合構成對點云曲面的覆蓋，將可用采樣邊界定義在采樣點切平面與鄰域球的交集上，可用采樣邊界的擴展通過鄰域球布爾交運算實現。如圖6所示，在加入新的采樣點后，Bi內可用采樣邊界消去了鄰域球Bj沿法向投影覆蓋的部分。

采樣邊界更新后，從當前所有采樣點鄰域球內可用采樣邊界的集合中隨機選取點p′，然后查找p′在點云中的最近點作為下一個采樣點。其中，對每個已有采樣點的采樣邊界賦予一概率值iq，隨長度減少以及采樣次數增加而減少。

圖6 鄰域球布爾交運算，可用采樣邊界用藍色表示，被消去的采樣邊界用紅色表示

本文算法同樣可用于點云重采樣，只需將最近點查找替換為移動最小二乘點云曲面投影計算即可。算法主要計算量在于鄰域球布爾交運算以及搜索最近點。采樣半徑r根據原始點云的的近似面積S按照下式計算[24]：

其中n是目標精簡數目，ρ是疏密控制參數一般取為0.7，S根據原始點云的球影響圖計算[25]。點云Poisson-disk采樣算法步驟如下：

Step 1根據原始點云近似面積估計采樣半徑，隨機選擇一個初始采樣點，生成初始可用采樣邊界。

Step 2根據每段可用采樣邊界的概率值，隨機選取新的備選點，查找該點在點云中的最近點作為新的采樣點。

Step 3利用鄰域球布爾交運算更新可用采樣邊界，降低當前更新的可用采樣邊界段的概率值。

Step 4迭代2，3步直到所有鄰域球內采樣邊界段最大概率值低于預先指定值Qτ。

初步采樣結果接近預先指定的采樣點數目，但由于Poisson-disk采樣的控制參數是采樣半徑因而不容易精確達到該數目，本文對初步采樣結果依采樣點Voronoi鄰域面積進行排序，迭代向稀疏區域添加或者從擁擠區域移除采樣點以達到目標精簡數目。

Poisson-disk采樣保證了相鄰采樣點之間的距離不小于2r從而防止了在高曲率區域產生局部聚集以及在平坦區域均衡分布。同時，如圖6所示，可用采樣邊界的更新依賴于鄰域球的交集在曲面切平面內的投影，因此，在尖銳特征附近，由于法向差異較大導致相對密集的采樣。

4 討論與比較

本文的實驗平臺是Intel?Core?2雙核處理器，1.86GHz，2G內存，NVIDIA GTS450顯卡，Windows XP?操作系統。

圖7顯示了點云選擇算法性能測試的真實掃描樣例數據，如圖8所示本文算法的速度快于射線奇偶法，對于規模較大的數據效果更加明顯。同時本文與商業軟件Geomagic[3]的點云選擇功能進行了比較，獲得了接近的用戶體驗。表1顯示了本文點云選擇算法的具體統計數據，主要包括矩形覆蓋計算時間、內外判斷時間以及單點包含測試剔除率。通過實驗可以看出，主要計算時間消耗在構建矩形覆蓋階段；采樣點數目與多邊形復雜性共同決定了計算時間，例如圖7(d)的時間消耗大于圖7(e)；大部分點在多邊形內測試被剔除，因而大大提升了交互速度。

圖7 點云選擇算法測試例子

圖8 點云選擇算法速度比較

如圖8所示，速度曲線圖中數據依頂點增長順序分別對應圖7中的(a)、(b)、(c)、(d)、(e)。紅色曲線表示射奇偶法，綠色曲線表示本文方法

表1 點云選擇算法統計數據表

本文與經典的曲率柵格以及均勻柵格[6]點云精簡方法進行了比較。首先為了測試對尖銳邊特征與邊界的保持，針對立方體模型以及紙張模型進行了實驗，其中立方體點云是由計算機生成的，紙張點云是真實的單幅掃描數據。如圖9與圖10所示，相較于曲率柵格與均勻柵格方法，本文算法加強了對尖銳邊特征以及邊界數據的采樣，并在平坦區域均衡分布采樣點。

圖11顯示了對真實掃描數據進行精簡的例子，結果表明本文的方法在保持尖銳邊特征與邊界數據的同時防止了采樣點局部聚集，獲得了更加均衡的分布，因而更加有利于高質量三角化，基于點的繪制以及形狀約束動畫[26]。

圖9 立方體模型點云精簡

圖10 紙張模型點云精簡

圖11 龍頭模型點云精簡

本文采用Cignoni[27]等方法定量計算了精簡點云重構曲面與原始點云之間的距離，本文算法與曲率柵格方法結果的平均距離誤差分別為0.0329 與0.1556。這是由于該數據背面是光滑曲面，曲率柵格方法未能均衡曲率采樣與均勻分布。如圖12所示，本文算法在平坦區域形成了良好的分布因而誤差比曲率柵格方法低，在曲率較大區域由于采樣點數量相對稀疏因而誤差比曲率柵格方法高。

圖13 龍頭模型點云精簡誤差分析

5 結 論

本文提出了面向逆向工程、用于消除原始掃描點云中存在無效背景數據、冗余采樣、分布不均勻等問題的點云選擇與精簡算法。點云選擇算法支持套索UI接口，通過套索多邊形的矩形覆蓋與點云八叉樹節點投影之間的包含關系剔除大部分點在多邊形內測試，達到了商業軟件的交互性能。點云精簡算法利用Poisson-disk采樣控制簡化后點云的分布，與經典曲率柵格方法的比較結果表明本文算法具有保持尖銳邊特征、保持邊界、分布均衡的優勢。

[1]Wang W. Reverse engineering technology of reinvention [M]. New York: CRC Press, 2011: 6-11.

[2]Weyrich T, Pauly M, Heinzle S, et al. Post-processing of scanned 3D surface data [C]//Proceedings of Eurographics Symposium on Point-Based Graphics.New York: ACM Press, 2004: 85-94.

[3]Geomagic studio. Geomagic studio user’s guide [M].Carolina: Geomagic Studio, 2008: 8-12.

[4]Inus Technology. RapidForm2004 user guide and tutorial [M]. Korea: Inus Technology Inc, 2004:59-62.

[5]張麗艷, 周儒容, 蔡煒斌, 等. 海量測量數據簡化技術研究[J]. 計算機輔助設計與圖形學學報, 2001,13(11): 1023-1119.

[6]Lee K H, Woo H, Suk T. Point data reduction using 3D grids [J]. The International Journal of Advanced Manufacturing Technology, 2001, 18(3): 201-210.

[7]Pauly M, Gross M, Kobbelt L P. Efficient simplification of point-sampled surfaces [C]//Proceedings of the Conference on Visualization,Washington D C: IEEE Press, 2002: 163-170.

[8]Haines Eric. Point in polygon strategies, graphics gems IV [M]. San Diego: Academic Press, 1994:24-26.

[9]Li Jing, Wang Wencheng, Wu Enhua. Point-in-polygon tests by convex decomposition [J]. Computers &Graphics, 2007, 31(4): 636-648.

[10]Zalik B, Kolingerova I. A cell-based point-in-polygon algorithm suitable for large sets of points [J].Computers & Geosciences, 2001, 27(10): 1135-1145.

[11]Yang S, Yong J H, Sun J, et al. A point-in-polygon method based on a quasi-closest point [J]. Computers& Geosciences, 2010, 36(2): 205-213.

[12]Linsen L. Point cloud representation. Universitat Karlsruhe [M]. CS technical report, 2001: 6-8.

[13]Song H, Feng H Y. A progressive point cloud simplification algorithm with preserved sharp edge data [J]. The International Journal of Advanced Manufacturing Technology, 2009, 45(5-6): 583-592.

[14]Yu Zhiwen, Wong Hansan, Peng Hong, et al. ASM:an adaptive simplification method for 3D point-based models [J]. Computer-Aided Design, 2010, 42(7):598-612.

[15]Shi Baoquan, Liang Jin, Liu Qing. Adaptive simplification of point cloud using k-means clustering [J].Computer-Aided Design, 2011, 43(8): 910-922.

[16]Song Hao, Feng Hsiyung. A global clustering approach to point cloud simplification with a specified data reduction ratio [J]. Computer-Aided Design, 2007,40(3): 281-292.

[17]Moenning C, Dodgson N A. Intrinsic point cloud simplification [C]//Proceedings of the 14th International Conference on Computer Graphics and Vision. Moscow, 2004: 1147-1154.

[18]Kang E C, Kim D B, Lee K H. Balanced feature-sensitive point sampling for 3D model generation [J]. The International Journal of Advanced Manufacturing Technology, 2007, 38(1-2): 130-142.

[19]Lastra M. An efficient parametric algorithm for octree traversal [J]. Journal of WSCG, 2000, 8(2): 212-219.

[20]Feito F R, Torres J C, Urena A. Orientation, simplicity,and inclusion test for planar polygons [J]. Computers& Graphics, 1995, 19(4): 595-600.

[21]Kobbelt L P, Botsch M, Schwanecke U, et al. Feature sensitive surface extraction from volume data [C]//Proceedings of the ACM SIGGRAPH. New York:ACM Press, 2001: 57-66.

[22]Zheng Y, Fu H, Au O K C, et al. Bilateral normal filtering for mesh denoising [J]. IEEE Transactions on Visualization and Computer Graphics, 2010, 17(10):1521-1530.

[23]Dunbar D, Humphreys G. A spatial data structure for fast poisson-disk sample generation [J]. ACM Transactions on Graphics, 2006, 25(3): 503-508.

[24]Fu Yan, Zhou Bingfeng. Direct sampling on surfaces for high quality remeshing [J]. Computer Aided Geometric Design, 2009, 26(6): 711-723.

[25]Klein J, Zachmann G. Point cloud surfaces using geometric proximity graphs [J]. Computers &Graphics, 2004, 28(6): 839-850.

[26]Zhao Hanli, Fan Ran, Wang Charlie C L, et al.Fireworks controller [J]. Journal of Visualization and Computer Animation, 2009, 20(2-3): 185-194.

[27]Cignoni P, Rocchini C, Scopigno R. Metro: measuring error on simplified surfaces [J]. Computer Graphics Forum, 1998, 17(2): 167-174.