官亞勤 趙學勝 王鵬飛 李大朋

摘要：針對傳統(tǒng)點云簡化算法效率低且處理點數(shù)少的缺陷，結(jié)合快速成型領域的切片原理顧及特征計算復雜度低的特點，設計并實現(xiàn)了適合千萬級海量激光雷達（LiDAR）點云的并行切片簡化算法。該算法根據(jù)切片原理對點云模型分層并按照角度排序，利用NVIDA的統(tǒng)一計算設備架構(gòu)（CUDA）和可編程圖形處理器（GPU）高度并行的性能優(yōu)勢，使用GPU多線程高效并行地執(zhí)行單層切片點云簡化，提高了算法效率。最后，應用3組不同數(shù)量級點云模型分別進行簡化對比實驗。實驗結(jié)果表明：在保持模型特征與壓縮比不變的情況下，所提算法效率高出傳統(tǒng)基于CPU的串行切片算法1～2個量級。

關(guān)鍵詞：

海量點云；簡化；切片法；計算設備架構(gòu)；圖形處理器；并行計算

中圖分類號： TP391.413 文獻標志碼：A

0引言

隨著大規(guī)模精細三維模型獲取技術(shù)的不斷發(fā)展，三維激光掃描技術(shù)憑借其數(shù)據(jù)獲取速度快、精度高、覆蓋廣的特點，成為高精度三維模型數(shù)據(jù)獲取的主流方式之一，獲取的點云數(shù)據(jù)量也呈幾何級數(shù)增長，因此，如何對海量散亂點云數(shù)據(jù)進行簡化已成為計算機圖形學、快速成型、三維測繪、地理信息系統(tǒng)、數(shù)字城市、軍事仿真、游戲娛樂等點云模型應用領域的重要研究課題之一。

傳統(tǒng)的點云簡化方法主要分為兩個大類：第一類是顧及特征的簡化[1-3]，此類算法需要依據(jù)單點的K鄰近點集擬合曲面，并構(gòu)建曲面的法向量和曲率等相關(guān)特征度量因子判定單點是否為特征點，從而實現(xiàn)點云簡化。這些算法能夠保持模型特征，但是涉及K鄰近點集等復雜計算操作，耗時多，僅適用于小數(shù)據(jù)量的點云簡化。第二類是規(guī)則采樣簡化算法[4-6]，此類算法依據(jù)一定規(guī)則對原始點云進行采樣，然后以采樣點作為特征點保留，剔除其他點實現(xiàn)點云簡化。這類算法簡化效率高，但是不能有效地保持模型特征，由于采樣標準單一，在特征變化明顯的尖銳彎曲處會導致局部細節(jié)過度光順丟失細節(jié)。由此可見，傳統(tǒng)算法的主要問題是點云簡化過程中計算復雜與模型特征保持不能兼顧。

近年來通用計算圖形處理器（General Purpose Graphics Processing Unit， GPGPU）的快速興起，尤其是NVIDIA公司2006年推出的圖形處理器（Graphics Processing Unit， GPU）并行計算框架——統(tǒng)一計算設備架構(gòu)（Compute Unified Device Architecture， CUDA）[7]憑借其高性價比、低通信開銷、卓越的并行計算能力，讓海量化或者計算復雜度高的三維點云模型數(shù)據(jù)快速處理成為可能。文獻[8]使用GPGPU實現(xiàn)基于邊緣點的激光雷達（Light Detection And Ranging， LiDAR）點云濾波算法，文獻[9]提出一種基于CUDA的雙邊濾波的點云濾波算法。兩者都將K鄰近點、曲面擬合、法向量以及曲率等計算復雜度高的步驟，利用CUDA編寫不同的kernel并行化，從而加速點云簡化，但是，由于在單個線程中完成類似求解單點K鄰近點的計算需要消耗太多的全局內(nèi)存等GPU資源，這嚴重制約此類算法處理點云的規(guī)模，僅適用于數(shù)十萬級的小規(guī)模點云處理。

本文借助眾核GPU通用計算高性能并行的特點，結(jié)合快速成型領域切片點云簡化算法[10-12]顧及特征的優(yōu)勢，實現(xiàn)了基于CUDA的顧及模型特征且適合千萬級點云的并行簡化算法，并從該算法的時間效率方面的優(yōu)勢闡述了CUDA用于海量數(shù)據(jù)處理的優(yōu)勢和潛力。

1切片點云算法原理

基于CUDA的切片算法實現(xiàn)原理如下：首先，在CPU端根據(jù)點云的幾何分布特征對點云進行分層并降維投影至相應的投影平面上；然后，依次對不同投影面上單層點云中的每一點和相應投影平面坐標原點所連直線與投影面某一坐標軸固定方向的夾角大小進行升序排序；最后，使用本文提出的利用CUDA在GPU端對每層排序后的切片點云依據(jù)弦高差法并行計算各點的弦高差值和各層切片的弦高差均值作為閾值，通過比較各點弦高差與閾值的關(guān)系，確定該點是否為特征點從而完成該層切片簡化。

1.1特征點的提取原理

利用弦高差法來確定切片中各點是否為特征點，原理詳見文獻[13]，其中弦高距離由式（1）求得：

di=|Axi+Byi+C|/A2i+B2i+C2i（1）

如圖1所示，pi為目標待判定點，直線pi-1pi+1所構(gòu)成直線方程為l：Ax+By+C=0，由計算幾何的原理可求得pi到直線l的垂直距離為di。

其中：mj表示第j層點云的總數(shù)；di表示第j層點云中的第i個點的弦高差值。

1.2基于CUDA的切片算法實現(xiàn)

由上述原理可以看出：依次計算單點的弦高差、單層切片的閾值σ以及單點弦高與閾值的比較等操作均是計算密集的操作，不涉及對原始切片數(shù)據(jù)的寫操作，不會因為數(shù)據(jù)的復寫而引發(fā)數(shù)據(jù)的二義性，具備良好的數(shù)據(jù)獨立性，因此本文借助NVIDIA推出的GPGPU平臺CUDA的單程序多數(shù)據(jù)（Single Program Multiple Data， SPMD）特性[7]，利用GPU中大規(guī)模并行處理器的并行計算能力，使用相互獨立的線程并發(fā)執(zhí)行這些計算，實現(xiàn)基于數(shù)據(jù)并行性的點云簡化，具體算法如下。

步驟1將依據(jù)角度排序后的單層點云切片從CPU的cpuvector中傳入至GPU的gpuvector中。

步驟2結(jié)合CUDA，設定核函數(shù)的線程塊數(shù)量blockDim.x和線程塊中線程數(shù)量threadDim.x，啟動核函數(shù)BowstringCaculate_kernel并行計算出每一點的弦高差值并存入GPU顯存的數(shù)組Height中。

步驟3在GPU中用并行歸約算法求出該層切片的弦高差均值σ，作為該層切片的特征判定閾值。

步驟4啟動核函數(shù)isFeaturePoint_kernel，根據(jù)BowstringCaculate_kernel返回的各點弦高Height與閾值σ，確定該層切片中的特征點，并將計算結(jié)果存入數(shù)組isFearturePoint中。

步驟5將isFearturePoint數(shù)組中的元素利用CUDA的核函數(shù)cudaMempy傳回至CPU端，在CPU端根據(jù)對應索引位置的值決定cpuvector相應位置上的點保留與否，從而得到簡化后的點云SimplyfiedVector。

步驟6回到步驟1繼續(xù)對其他切片層的點云簡化。算法中的流程如圖2所示。

其中，執(zhí)行單層切片點云弦高差計算的核函數(shù)BowstringCaculate_kernel的偽代碼如下。

算法BowstringCaculate_kernel。

有序號的程序——————————Shift+Alt+Y

程序前

輸入：LayerPoints為排序后的單層切片點云的坐標數(shù)組；Size表示該層切片的點云數(shù)量。

輸出：Height為用于保存每個線程計算的弦高差值的數(shù)組。

1）

Dim index ← threadInx.x+blockIdx.x*blockDim.x；

2）

Dim *leftPoint，*rightPoin；

3）

if（index

4）

Loop index from 0 to Size Do

5）

leftPoint ← LayerPoints+index-1；

6）

rightPoint ← LayerPoints+index-1；

7）

if（index==0）Then

8）

leftPoint ← LayerPoints+Size-1；

9）

End if

10）

if（index==Size-1）

11）

rightPoint ← LayerPoints+index-（Size-1）；

12）

End if

13）

Height[index] ← calculateHeight（leftPoint；

14）

LayerPoint[index]，rightPoint）；

15）

index ← index+blockDim.x*GridDim.x；

16）

End Loop

17）

End if

程序后

其中，calculateHeight為計算單點的弦高函數(shù)，其實現(xiàn)原理如式（1）和圖1所示。

2實驗分析

實驗平臺配置如下：Windows 7操作系統(tǒng)，CPU為Intel Core i53470 @ 3.20GHz 3.60GHz，內(nèi)存為4.0GB（3.28GB可用），顯卡為NVIDIA GeForce GT 640。在該平臺下使用C++語言結(jié)合Visual Studio2010和NVIDIA的CUDA6.0框架實現(xiàn)本文算法。

算法實驗的模型數(shù)據(jù)如圖3所示：龍模型共437645個點，大佛模型共有1753052個點，高觀音模型共11807207個點。

本文使用上述3個不同數(shù)量級的點云模型，對本文提出的基于CUDA的切片點云簡化算法（簡稱GPU切片算法）與基于CPU的切片點云簡化算法[13]（以下簡稱CPU切片法）進行對比實驗。

其中：弦高差閾值由式（1）求出，切片方向均為z軸，依次選取切片數(shù)量laynum為10，25，50，75，100作5組對比實驗，結(jié)果如表1所示。

通過表1，可依次求得上述3個不同數(shù)量級模型的算法耗時與切片層數(shù)的關(guān)系曲線如圖4所示。

其中圖4中的折線圖（a）、（b）、（c）依次表示龍模型、觀音模型以及高觀音模型的算法耗時與切片層數(shù)的對應關(guān)系，（d）表示龍模型、觀音模型以及高觀音模型對應的CPU切片法與GPU切片法的加速比（其中加速比等于CPU算法執(zhí)行時間除以GPU算法執(zhí)行時間）與切片層數(shù)的對應關(guān)系。

從表1呈現(xiàn)的數(shù)據(jù)以及圖4的折線圖（a）、（b）、（c）呈現(xiàn)的線條走勢可以看出：在兩種算法對應的壓縮率基本一致基礎上，本文提出的GPU切片算法的效率比傳統(tǒng)CPU算法高出10～30倍，約1～2個量級；但是，隨著切片層數(shù)的增加，本文提出的GPU算法耗時有一定程度的增加。因為在CUDA架構(gòu)的切片算法中（流程如圖2所示），隨著切片層數(shù)的增加，CPU和GPU之間進行I/O交互的次數(shù)也隨之增加，最終導致算法的執(zhí)行時間有一定程度的增加。從圖4（d）可以看出：龍模型和觀音模型的加速比，均隨著切片層數(shù)增加有一定程度的減小。而數(shù)據(jù)量多達千萬個點的高觀音模型，其加速比曲線變化相對平穩(wěn)。這是由于GPU更適合于密集型的計算，當數(shù)據(jù)量（計算量）較小時，算法在GPU上的執(zhí)行時間無法隱藏訪問和數(shù)據(jù)傳輸?shù)难舆t，而隨著數(shù)據(jù)量的增大，這些延遲逐漸被隱藏，因此加速比逐漸增大。而當數(shù)據(jù)量增大一定的程度，GPU近乎滿負荷的工作，所有的訪問和數(shù)據(jù)傳輸?shù)难舆t都已被很好地隱藏，加速比也趨于穩(wěn)定如高觀音模型的加速比曲線所示。

此外，以模型特征最為明顯的觀音模型為例，依次選取該模型切片數(shù)為25，50，75的底座前側(cè)簡化局部視圖，如圖5所示。發(fā)現(xiàn)當切片層數(shù)為25時，由于壓縮率粒度太大導致底座衣服褶皺與蓮花形等多處被過度平滑，特征丟失太嚴重；切片層數(shù)為75時，由于壓縮粒度太小導致殘留的冗余點較多；而切片層數(shù)為50時，底座衣服褶皺與蓮花形的特征細節(jié)完整保持，而且冗余的點較少，簡化效果是較為理想，因此，對不同的模型選取合適的切片層數(shù)對模型簡化至關(guān)重要。

3結(jié)語

本文利用CUDA的高性能并行計算優(yōu)勢，對傳統(tǒng)基于CPU的串行切片點云簡化算法進行了改進，將傳統(tǒng)算法的核心步驟：單點的弦高差計算與特征點判定算法邏輯并行化，通過對不同數(shù)量級的3個點云模型的簡化實驗，得出以下結(jié)論：

1）對于同一模型，GPU算法盡管隨著切片層數(shù)的增加，耗時由于數(shù)據(jù)交互次數(shù)增加有一定程度的小幅震蕩，但均遠少于傳統(tǒng)CPU算法。

2）對于不同數(shù)量級的模型，加速比曲線隨著點云數(shù)量的增加而逐漸穩(wěn)定，且加速效果更優(yōu)，驗證了本文算法應對海量點云簡化的優(yōu)勢和潛力。

3）模型的特征保留完整性與切片層數(shù)無直接關(guān)系，僅僅與模型表面特征有關(guān)。

下一步主要工作是在CPU端使用多線程并行技術(shù)，提高點云切片分層排序的速度；應用GPU架構(gòu)中不同訪問性能的內(nèi)存模型和基于任務并行性的流水線模型對算法進行優(yōu)化。

參考文獻：

[1]

周煜，張萬冰，杜發(fā)榮，等.散亂點云數(shù)據(jù)的曲率精簡算法[J].北京理工大學學報，2010，30（7）：785-790.（ZHOU Y， ZHANG W B， DU F R， et al. Algorithm for reduction of scattered point cloud data based on curvature [J].Transactions of Beijing Institute of Technology， 2010， 30（7）： 785-790.）

[2]

李義琛，龐明勇.基于二次誤差度量的點云簡化[J].小型微型計算機系統(tǒng)，2012，33（11）：2538-2543.（LI Y C， PANG M Y. Decimating point cloud based on quadric error metric [J]. Journal of Chinese Computer Systems， 2012， 33（11）： 2538-2543.）

[3]

王先澤，李忠科，張曉娟，等.特征保持的基于緊支徑向基函數(shù)的點云簡化[J].計算機工程與設計，2013，34（1）：201-206.（WANG X Z， LI Z K， ZHANG X J， et al. Feature preserving simplification of point cloud based on CSRBF [J]. Computer Engineering and Design， 2013， 34（1）： 201-206.）

[4]

胡志勝，于敬武，述濤.一種結(jié)合了柵格化和特征判斷的點云壓縮方法[J].遼寧工程技術(shù)大學（自然科學版），2015，34（8）：958-963.（HU Z S， YU J W， SHU T.A point cloud compression approach combined with rasterizing and feature estimate [J]. Journal of Liaoning Technical University （Natural Science）， 2015， 34（8）： 958-963.）

[5]

邢正全，鄧喀中，薛繼群.基于柵格劃分和法向量估計得點云數(shù)據(jù)壓縮[J].測繪通報，2012（7）：50-54.（XING Z Q， DENG K Z， XUE J Q. Point cloud data compression based on grid division and normal vector estimation [J]. Bulletin of Surveying and Mapping， 2012（7）： 50-54.）

[6]

邵正偉，席平.基于八叉樹編碼的點云數(shù)據(jù)精簡方法[J].工程圖學學報，2010，31（4）：73-77.（SHAO Z W， XI P. Data reduction for point cloud using octree coding [J]. Journal of Engineering Graphics， 2010， 31（4）： 73-77.）

[7]

趙開勇，汪朝輝.大規(guī)模并行處理器編程實戰(zhàn)[M].2版.北京：清華大學出版社，2013：36-40.（ZHAO K Y， WANG C H. Programming Massively Parallel Processors： a Handson Approach [M]. 2nd ed. Beijing： Tsinghua University Press， 2013： 36-40.）

[8]

崔放，徐宏根，王宗躍.基于GPGPU的并行LiDAR點云濾波算法[J].華中師范大學學報（自然科學版版），2014，48（3）：431-436.（CUI F， XU H G， WANG Z Y. A point cloud filtering algorithm based on GPGPU parallel computing [J].Journal of Huazhong Normal University （Natural Science）， 2014， 48（3）： 431-436.）

[9]

唐杰，徐波，宮中樑，等.一種基于CUDA的三維點云快速光順算法[J].系統(tǒng)仿真學報，2012，24（8）：1633-1638.（TANG J， XU B， GONG Z L， et al. Fast fairing of 3D point cloud using CUDA [J].Journal of System Simulation， 2012， 24（8）： 1633-1638.）

[10]

PARK H T， CHANG M H， PARK S C. A slicing algorithm of point cloud for rapid prototyping [C]// Proceedings of the 2007 Summer Computer Simulation Conference. San Diego， CA： Society for Computer Simulation International， 2007： Article No. 24.

[11]

SHIN H， PARK S. Direct slicing of a point set model for rapid prototyping [J]. ComputerAided Design and Applications， 2004， 1（1/2/3/4）： 109-115.

[12]

PERCOCO G， GALANTUCCI L. Localgenetic slicing of point clouds for rapid prototyping [J]. Rapid Prototyping Journal， 2008， 14（3）： 161-166.

[13]

方芳，程效軍.海量散亂點云快速壓縮算法[J].武漢大學學報（信息科學版），2013，38（11）：1353-1357.（FANG F， CHENG X J. A fast reduction method for massive scattered point cloud based on slicing [J]. Geomatics and Information Science of Wuhan University， 2013， 38（11）： 1353-1357.）

[14]

徐工，程效軍.基于小波技術(shù)的散亂點云自適應壓縮算法[J].同濟大學學報（自然科學版），2013，41（11）：1738-1743.（XU G， CHENG X J. Adaptive reduction algorithm of scattered point clouds based on wavelet technology [J]. Journal of Tongji University （Natural Science）， 2013， 41（11）： 1738-1743.）