張 春 亢，趙 學(xué) 勝，鐘 新 科

(1.中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083； 2．中國(guó)科學(xué)院地理科學(xué)與資源研究所，北京 100101)

?

基于三角形驅(qū)動(dòng)的機(jī)載LiDAR數(shù)據(jù)柵格化算法

張 春 亢1，趙 學(xué) 勝1，鐘 新 科2

(1.中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083； 2．中國(guó)科學(xué)院地理科學(xué)與資源研究所，北京 100101)

針對(duì)基于TIN的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)插值為GRID過(guò)程中，大量TIN數(shù)據(jù)I/O操作導(dǎo)致的柵格化效率下降問(wèn)題，提出了基于三角形驅(qū)動(dòng)的點(diǎn)云柵格化流式算法:以基于直線正負(fù)區(qū)判別原理的TIN向GRID轉(zhuǎn)換新算法為基礎(chǔ)，通過(guò)遍歷三角形模擬流計(jì)算實(shí)現(xiàn)三角形生成、三角形插值為GRID以及內(nèi)存釋放，有效避免了TIN數(shù)據(jù)的I/O操作，并可以提高內(nèi)存利用率。試驗(yàn)表明：流式算法顯著提高了點(diǎn)云柵格化效率，為海量點(diǎn)云的柵格化及并行計(jì)算提供了一種算法支撐。

機(jī)載LiDAR點(diǎn)云；柵格化；I/O操作；流計(jì)算；三角形驅(qū)動(dòng)

作為一種高精度、高密度、快速獲取地面三維信息的有效手段，機(jī)載LiDAR已成為生成數(shù)字地面模型(DTM)的重要工具[1]。而作為DTM的兩種主要表達(dá)方式，規(guī)則格網(wǎng)(GRID)相對(duì)于不規(guī)則三角網(wǎng)(TIN)具有數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單、存儲(chǔ)量小、分析與計(jì)算方便等優(yōu)點(diǎn)[2]。將LiDAR點(diǎn)云柵格化，不但便于其管理與表達(dá)，也為其濾波與應(yīng)用創(chuàng)造了條件，但在柵格化大量LiDAR點(diǎn)云的過(guò)程中，因數(shù)據(jù)的I/O操作時(shí)間超過(guò)CPU計(jì)算時(shí)間而成為制約轉(zhuǎn)換效率的瓶頸[3-5]。在實(shí)現(xiàn)較大數(shù)據(jù)量的點(diǎn)云直接插值為GRID時(shí)，只需將點(diǎn)云以一定的數(shù)據(jù)結(jié)構(gòu)進(jìn)行劃分重組，即可提高I/O操作的效率，實(shí)現(xiàn)點(diǎn)云向GRID的高效轉(zhuǎn)換[4，5]。與直接插值方法相比，利用點(diǎn)云生成TIN，然后將TIN內(nèi)插為GRID具有更高的精度和效率[2]。傳統(tǒng)TIN插值為GRID通常采用對(duì)生成的TIN數(shù)據(jù)進(jìn)行分塊并建立索引，然后逐個(gè)判斷待插值格網(wǎng)節(jié)點(diǎn)在哪個(gè)三角形中，進(jìn)而計(jì)算格網(wǎng)點(diǎn)屬性值的插值策略[6]。因此，在實(shí)現(xiàn)基于TIN的點(diǎn)云數(shù)據(jù)插值為GRID時(shí)，需首先將點(diǎn)云剖分為TIN并進(jìn)行存儲(chǔ)，然后完成TIN向GRID轉(zhuǎn)換。這種基于格網(wǎng)節(jié)點(diǎn)驅(qū)動(dòng)的分步轉(zhuǎn)換方法雖然易于理解與實(shí)現(xiàn)，但轉(zhuǎn)換過(guò)程中涉及中間數(shù)據(jù)TIN的I/O操作。離散點(diǎn)的數(shù)量m(m≥3)與其所構(gòu)建的三角形數(shù)量n滿足條件[7]：m-2≤n≤2m-5。實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)m較大時(shí)，n逼近2m，并且TIN復(fù)雜的數(shù)據(jù)結(jié)構(gòu)決定其數(shù)據(jù)量要遠(yuǎn)大于離散點(diǎn)本身的數(shù)據(jù)量，所以點(diǎn)云數(shù)據(jù)比較大時(shí)，對(duì)計(jì)算機(jī)內(nèi)存也有較高要求。通過(guò)上述分析，本文引入流計(jì)算的思想，以TIN向GRID轉(zhuǎn)換新算法為基礎(chǔ)，對(duì)基于TIN的點(diǎn)云柵格化流程進(jìn)行了改造，以避免TIN數(shù)據(jù)的I/O操作，提高柵格化效率及內(nèi)存利用率。

1 基于三角形驅(qū)動(dòng)的流式算法

1.1 算法原理

流式算法的基本思想為：利用逆向思維，依據(jù)直線正負(fù)區(qū)判別原理，提出一種新的TIN向GRID轉(zhuǎn)換算法，使構(gòu)成TIN的每個(gè)三角形能夠獨(dú)立地實(shí)現(xiàn)其內(nèi)格網(wǎng)節(jié)點(diǎn)判斷及其屬性值的計(jì)算;同時(shí)引入流計(jì)算的思想，實(shí)現(xiàn)三角形的構(gòu)建、三角形插值為GRID和內(nèi)存的釋放在三角形驅(qū)動(dòng)下，以流水線的形式進(jìn)行組合，從而避免TIN數(shù)據(jù)的I/O操作并提高內(nèi)存的利用率。流式算法步驟如下：1)對(duì)于點(diǎn)云數(shù)據(jù)，利用Delaunay三角網(wǎng)構(gòu)建算法生成一個(gè)三角形T;2)利用本文提出的基于直線正負(fù)區(qū)判別原理的TIN向GRID轉(zhuǎn)換算法,對(duì)三角形T所包含的格網(wǎng)節(jié)點(diǎn)進(jìn)行判斷與插值;3)存儲(chǔ)完成插值的格網(wǎng)節(jié)點(diǎn)屬性值，釋放三角形T和格網(wǎng)節(jié)點(diǎn)所占用的內(nèi)存資源;4)重復(fù)步驟1-3，直至TIN構(gòu)建完成，實(shí)現(xiàn)所有三角形包含的格網(wǎng)節(jié)點(diǎn)的判斷與插值;5)對(duì)未包含在任何三角形內(nèi)的格網(wǎng)節(jié)點(diǎn)賦予一個(gè)特定屬性值，算法結(jié)束。

在流式算法的實(shí)現(xiàn)過(guò)程中，Delaunay三角網(wǎng)構(gòu)建算法需要選用三角網(wǎng)生長(zhǎng)算法或分治算法，以便在TIN構(gòu)建過(guò)程中不間斷地釋放符合Delaunay三角網(wǎng)準(zhǔn)則、不需再優(yōu)化的三角形，從而形成數(shù)據(jù)流。本文詳細(xì)闡述基于直線正負(fù)區(qū)判別原理的TIN向GRID轉(zhuǎn)換算法。

1.2 基于直線正負(fù)區(qū)判別原理的TIN向GRID轉(zhuǎn)換算法

存在一直線，假定直線兩邊分別為該直線的正區(qū)與負(fù)區(qū)，有如下公式：

(1)

其中：A=(y2-y1)/(x2-x1)，B=(x2y1-x1y2)/(x2-x1)，(x1,y1)、(x2,y2)為直線上兩點(diǎn)坐標(biāo)。對(duì)于任意一點(diǎn)P(x,y)，可利用下式判斷該點(diǎn)處于直線的正區(qū)或負(fù)區(qū)：

(2)

由式(1)、式(2)可知：平面內(nèi)兩點(diǎn)的表達(dá)式同號(hào)時(shí)，則兩點(diǎn)位于直線段的同側(cè)；表達(dá)式異號(hào)時(shí)，兩點(diǎn)位于直線段的異側(cè)。

依據(jù)上述直線正負(fù)區(qū)判別原理，TIN向GRID轉(zhuǎn)換算法描述如下：1)如圖1，對(duì)于構(gòu)成TIN的一個(gè)三角形ABC，根據(jù)其頂點(diǎn)坐標(biāo)，計(jì)算其頂點(diǎn)在x與y方向上坐標(biāo)的最大(小)值xmax(xmin)、ymax(ymin)；2)根據(jù)xmin、xmax、ymin、ymax與格網(wǎng)大小，計(jì)算格網(wǎng)的起始行Ymin和終止行Ymax以及起始列Xmin和終止列Xmax，則直線x=Xmin、x=Xmax、y=Ymin和y=Ymax相交構(gòu)成矩形EFGH，計(jì)算矩形內(nèi)格網(wǎng)節(jié)點(diǎn)的坐標(biāo)值；3)將三角形的3個(gè)頂點(diǎn)A、B、C與三角形的三邊AB、AC、BC組成三組點(diǎn)線數(shù)據(jù)A、BC，B、AC和C、AB，對(duì)于需要判斷的格網(wǎng)節(jié)點(diǎn)P，當(dāng)點(diǎn)P與三組點(diǎn)線數(shù)據(jù)中的點(diǎn)相對(duì)于該組中的線段運(yùn)用式(1)與式(2)計(jì)算的值都同號(hào)時(shí)，則點(diǎn)P在三角形ABC內(nèi)，否則，點(diǎn)P在三角形外；4)根據(jù)三角形頂點(diǎn)的坐標(biāo)及其屬性值內(nèi)插三角形內(nèi)格網(wǎng)節(jié)點(diǎn)的屬性值；5)重復(fù)步驟1-4，遍歷TIN中的三角形，實(shí)現(xiàn)TIN向GRID的轉(zhuǎn)換。

圖1 三角形插值為格網(wǎng)節(jié)點(diǎn)示意

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

2.1TIN向GRID轉(zhuǎn)換算法實(shí)驗(yàn)與分析

實(shí)驗(yàn)環(huán)境為Intel(R)Pentium(R)DCPU(3.40GHz)，1G內(nèi)存。圖2為本文算法的測(cè)試效率，其格網(wǎng)分辨率約1.6個(gè)節(jié)點(diǎn)/三角形。文獻(xiàn)[6]與本文算法效率對(duì)比見表1。分析可知：本文算法的時(shí)間復(fù)雜度為O(n)，n為構(gòu)成TIN的三角形個(gè)數(shù)；本文算法比傳統(tǒng)算法在效率上有明顯的提高，當(dāng)點(diǎn)云數(shù)量或格網(wǎng)分辨率增大時(shí)，效率優(yōu)勢(shì)更加突出。因此，在算法效率與時(shí)間復(fù)雜度方面，本文算法更適合海量點(diǎn)云的處理。

圖2 TIN向GRID轉(zhuǎn)換效率

點(diǎn)數(shù)(個(gè))三角形數(shù)(個(gè))行列數(shù)(行數(shù)×列數(shù))本文算法所用時(shí)間(s)文獻(xiàn)[6]算法所用時(shí)間(s)843816754800×6000．22．310986218043200×24003．372．5

2.2 流式算法實(shí)驗(yàn)與分析

表2為流式算法測(cè)試結(jié)果，其中TIN的生成利用文獻(xiàn)[8]提出的分治算法。作為對(duì)比，測(cè)試了基于格網(wǎng)節(jié)點(diǎn)驅(qū)動(dòng)的分步轉(zhuǎn)換算法的時(shí)間，包括TIN的生成、TIN的I/O操作與TIN插值為GRID三部分時(shí)間。由測(cè)試結(jié)果可知：1)TIN的I/O操作占用了分步轉(zhuǎn)換算法消耗的大部分時(shí)間，嚴(yán)重制約了轉(zhuǎn)換算法的效率，可見，避免TIN的I/O操作對(duì)于提高點(diǎn)云柵格化效率意義重大；2)因避免了TIN的I/O操作，流式算法在效率上比分步轉(zhuǎn)換算法有了顯著提高；3)相比于分步測(cè)試中TIN構(gòu)建與TIN插值為GRID的時(shí)間和，流式算法多消耗了部分時(shí)間，這些時(shí)間主要用于內(nèi)存空間的申請(qǐng)與釋放。另外，因文獻(xiàn)[8]是一種分治算法，其遞歸的特性使其對(duì)內(nèi)存空間要求較高，無(wú)法對(duì)更大數(shù)據(jù)量的機(jī)載LiDAR點(diǎn)云進(jìn)行處理。

表2 流式算法測(cè)試結(jié)果

并行計(jì)算是處理海量空間數(shù)據(jù)的有效方法[3,5]，串行算法并行化的一個(gè)難點(diǎn)在于任務(wù)分解。分步算法在進(jìn)行并行化時(shí)任務(wù)分解較困難，流式算法基于三角形驅(qū)動(dòng)的轉(zhuǎn)換方式使任務(wù)分解相對(duì)簡(jiǎn)單，點(diǎn)云的柵格化并行處理也易于實(shí)現(xiàn)。

3 結(jié)論

流計(jì)算是避免I/O操作的一種有效方式，本文針對(duì)基于TIN的點(diǎn)云柵格化方法存在的缺陷，對(duì)柵格化過(guò)程進(jìn)行了改進(jìn),即通過(guò)遍歷三角形模擬流計(jì)算，有效避免了TIN數(shù)據(jù)的I/O操作，提高了內(nèi)存利用效率，實(shí)現(xiàn)了大量點(diǎn)云的快速柵格化。下一步將著重研究兼顧算法時(shí)間和空間性能的 Delaunay三角網(wǎng)剖分算法，用于海量點(diǎn)云數(shù)據(jù)的柵格化流式處理及并行計(jì)算。

[1] KOBLER A,PFEIFER N,OGRINC P,et al.Repetitive interpolation:A robust algorithm for DEM generation from Aerial Laser Scanner data in forested terrain[J].Remote Sensing of Environment,2007,108(1):9-23.

[2] 李志林,朱慶.數(shù)字高程模型[M].武漢:武漢測(cè)繪科技大學(xué)出版社,2000.34-59.

[3] 宋效東,劉學(xué)軍,湯國(guó)安,等.DEM與地形分析的并行計(jì)算[J].地理與地理信息科學(xué),2012,28(4):1-7.

[4] AGARWAL P K,ARGE L,DANNER A.From point cloud to grid DEM:A scalable approach[A].The 12th International Symposium on Spatial Data Handing[C].2006.771-788.

[5] GUAN X F,WU H Y.Leveraging the power of multi-core platform for large-scale geospatial data processing:Exemplified by generating DEM from massive LiDAR point clouds[J].Computers & Geosciences,2010,36(10):1276-1282.

[6] 吳飛,吳凡.TIN向規(guī)則格網(wǎng)DEM轉(zhuǎn)換的快速算法[J].測(cè)繪科學(xué),2005,30(4):76-77.

[7] 章孝燦,黃智才,戴企成,等.GIS中基于拓?fù)浣Y(jié)構(gòu)和凸殼技術(shù)的快速TIN生成算法[J].計(jì)算機(jī)學(xué)報(bào),2002，25(11):1212-1218.

[8] 胡金星,潘懋,馬照亭,等.高效構(gòu)建Delaunay三角網(wǎng)數(shù)字地形模型算法研究[J].北京大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,39(5):736-741.

2014-06-16；

2014-08-19

高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金項(xiàng)目(20130023110001)

張春亢(1984-)，男，博士研究生，主要研究方向?yàn)長(zhǎng)iDAR數(shù)據(jù)處理及應(yīng)用。E-mail：chkang.chd@163.com

10.3969/j.issn.1672-0504.2015.01.026