官亞勤 趙學(xué)勝 王鵬飛 李大朋

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

關(guān)鍵詞：

海量點(diǎn)云；簡(jiǎn)化；切片法；計(jì)算設(shè)備架構(gòu)；圖形處理器；并行計(jì)算

中圖分類號(hào)： TP391.413 文獻(xiàn)標(biāo)志碼：A

0引言

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

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

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

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

1切片點(diǎn)云算法原理

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

1.1特征點(diǎn)的提取原理

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

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

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

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

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

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

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

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

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

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

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

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

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

算法BowstringCaculate_kernel。

有序號(hào)的程序——————————Shift+Alt+Y

程序前

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

輸出：Height為用于保存每個(gè)線程計(jì)算的弦高差值的數(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為計(jì)算單點(diǎn)的弦高函數(shù)，其實(shí)現(xiàn)原理如式（1）和圖1所示。

2實(shí)驗(yàn)分析

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

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

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

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

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

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

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

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

3結(jié)語(yǔ)

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

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

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

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

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

參考文獻(xiàn)：

[1]

周煜，張萬(wàn)冰，杜發(fā)榮，等.散亂點(diǎn)云數(shù)據(jù)的曲率精簡(jiǎn)算法[J].北京理工大學(xué)學(xué)報(bào)，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]

李義琛，龐明勇.基于二次誤差度量的點(diǎn)云簡(jiǎn)化[J].小型微型計(jì)算機(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ù)的點(diǎn)云簡(jiǎn)化[J].計(jì)算機(jī)工程與設(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é)合了柵格化和特征判斷的點(diǎn)云壓縮方法[J].遼寧工程技術(shù)大學(xué)（自然科學(xué)版），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]

邢正全，鄧喀中，薛繼群.基于柵格劃分和法向量估計(jì)得點(diǎn)云數(shù)據(jù)壓縮[J].測(cè)繪通報(bào)，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]

邵正偉，席平.基于八叉樹編碼的點(diǎn)云數(shù)據(jù)精簡(jiǎn)方法[J].工程圖學(xué)學(xué)報(bào)，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ī)模并行處理器編程實(shí)戰(zhàn)[M].2版.北京：清華大學(xué)出版社，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點(diǎn)云濾波算法[J].華中師范大學(xué)學(xué)報(bào)（自然科學(xué)版版），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的三維點(diǎn)云快速光順?biāo)惴╗J].系統(tǒng)仿真學(xué)報(bào)，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]

方芳，程效軍.海量散亂點(diǎn)云快速壓縮算法[J].武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），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ù)的散亂點(diǎn)云自適應(yīng)壓縮算法[J].同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版），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.）