Generating Raster DEM from Mass Points via TIN

Generating Raster DEM from Mass Points via TIN Streaming LIDAR points raster DEM Martin Isenburg Yuanxin Liu Jonathan Shewchuk UC Berkeley UNC Chapel Hill UC Berkeley Jack Snoeyink Tim Thirion UNC Chapel Hill

Generating DEMs from LIDAR resampling LIDAR points triangulating raster DEM j i temporary TIN interpolating

Billions of Points courtesy of www. ncfloodmaps. com Neuse-River Basin Elevation map 0. 5 billion points 11 GB Contour map

DEM Generation via TIN Streaming read points from disk Streaming Computation of Delaunay Triangulation small buffer in main memory Immediate Streaming Rasterization TIN of Streaming output TIN store elevation rasters to temporary files (grouped by rows) output final DEM to disk

Example Result 500, 141, 313 Points 11 GB (binary, xyz, doubles) LIDAR points 50, 394 × 30, 500 DEM 3 GB (binary, BIL, 16 bit, 20 ft) raster DEM on a household laptop with two harddisks in 67 minutes 64 MB of main memory 270 MB temporary disk space

Related Work

Rasterizing Billions of Points • popular GIS packages [Agarwal et al. ‘ 06] – Arc. GIS 9. 1 (limit: ~ 25 million points, 1 GB) – QTModeller 4 (limit: ~ 50 million points, 1 GB) – GRASS (limit: ~ 25 million points, 1 GB) • break data into (overlapping) tiles

Interpolation Methods • • • Inverse Distance Weighting (IDW) Natural Neighbors Kriging Splines TIN – linear – quintic • …

External Memory Quadtree From Point Cloud to Grid DEM: A Scalable Approach [Agarwal et al. ‘ 06] • rearrange points on disk into quadtree • process cell – find its neighbors – interpolate all (i, j) – write rasters to temporary file • sort temporary file on (i, j) & write DEM

Our Approach: TIN Streaming • reorganize computations … not points – compute TIN in places where all points have already arrived – raster, output & deallocate – keep parts that miss neighbors • enhance points with spatial finalization

Our Approach: TIN Streaming • reorganize computations … not points – compute TIN in places where all points have already arrived – raster, output & deallocate – keep parts that miss neighbors • enhance points with spatial finalization

Our Approach: TIN Streaming • reorganize computations … not points – compute TIN in places where all points have already arrived – raster, output & deallocate – keep parts that miss neighbors • enhance points with spatial finalization

Delaunay Triangulation • standard algorithm for TIN construction • good properties for interpolating points • a triangulation in which every triangle has an empty circumscribing circle

Incremental Point Insertion [Lawson ’ 77] [Bowyer ’ 81] [Watson ’ 81] • locate triangle enclosing the point • find and remove all triangles with nonempty circumcircles • triangulate by connecting new point

Spatial Finalization Tags • inject tags into the point stream that inform when grid cells are free of future points • use tags to certify triangles as being part of the final triangulation and output them finalized region future points only arrive in these cells

Streaming TIN Generation input from disk another small buffer in main memory output to disk output (pipe) streaming TIN to rasterizer finalization tags Streaming Computation of Delaunay Triangulations [Isenburg et al. ‘ 06] • triangulate point stream seamlessly • output completed parts & reuse memory

Rasterization of Streaming TIN input from disk input (pipe) streaming TIN from triangulator yet another small buffer output final DEM to disk • points from single UTM time zone • triangles with long edges are killed • output to a single row-ordered DEM

spfinalize -i points. raw -ospb | spdelaunay –ispb –osmb -ospb | tin 2 iso –ismb –o dem. bil Streaming Rasterization Pipeline

The Finalizer (spfinalize)

Spatial Finalization of Points compute bounding box

Spatial Finalization of Points 1 7 4 2 4 5 3 6 7 1 7 8 1 7 7 2 4 5 8 6 1 7 2 9 7 9 1 6 9 6 compute 4 9 5 6 7 5 7 8 8 6 1 8 8 7 8 9 3 7 7 7 9 8 7 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2

Spatial Finalization of Points 1 7 4 5 3 6 7 8 7 7 2 5 8 7 9 9 6 compute 4 9 6 7 7 8 8 6 8 8 7 8 9 3 7 7 7 9 8 7 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2 output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 0 5 4 5 3 3 7 8 7 7 2 1 7 7 9 9 6 compute 4 9 6 7 7 8 8 6 8 8 7 8 9 3 7 7 7 9 8 7 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2 output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 1 5 4 5 3 0 6 8 7 7 2 0 2 5 9 9 6 compute 4 9 6 7 7 8 8 6 8 8 7 8 9 3 7 7 7 9 8 7 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2 output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 1 0 4 5 3 0 3 8 7 7 2 0 0 3 9 9 6 compute 3 9 4 7 7 8 8 6 8 8 7 8 9 3 7 7 7 9 8 7 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2 output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 0 1 0 4 1 1 1 2 3 7 9 7 7 6 8 9 8 2 6 6 compute 7 7 8 9 3 7 7 7 9 8 7 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2 output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 0 compute 2 1 0 4 1 6 5 6 4 8 8 8 2 6 6 9 1 3 5 7 9 8 6 7 8 7 5 bounding box create finalization grid – count number of points per cell 2 4 8 2 output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 0 2 1 0 4 compute bounding box create finalization grid – count number of points per cell output 1 1 3 3 2 7 2 5 5 9 7 4 2 finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Spatial Finalization of Points 0 1 0 4 2 compute bounding box create finalization grid – count number of points per cell output finalized points – buffer per grid cell – if full, output points in one randomized chunk followed by finalization tag

Exploit the Existing Coherence • document with finalization tags • enhance by releasing chunks of points benefits of global sort, but cheaper raw point stream finalized point stream

Exploit the Existing Coherence • document with finalization tags • enhance by releasing chunks of points benefits of global sort, but cheaper raw point stream finalized point stream

Unfinalized Region as Quadtree • as tags finalize cell after cell, maintain remaining cells in adaptive quadtree used for efficient overlap tests by the triangulator finalization tags

The Triangulator (spdelaunay)

Modified Incremental Delaunay • insert points • when reading a finalization tag – find certified triangles – output them to disk or pipe – deallocate data structure

Modified Incremental Delaunay • insert points • when reading a finalization tag – find certified triangles – output them to disk or pipe – deallocate data structure

Modified Incremental Delaunay • insert points • when reading a finalization tag – find certified triangles – output them to disk or pipe – deallocate data structure

The Rasterizer (tin 2 dem)

Rasterizing the Streaming TIN write raster directly into final DEM on disk write raster to temporary files on disk (use small buffer, compress when writing buffer)

Results

Neuse. River Basin

Pipeline Details for the 20 ft DEM 10 MB 60, 388 triangles 6 min 5 MB 512 x 512 grid 43 min 35 MB 2, 249, 268 points 8 min 270 MB 239 files 46 MB 128 rows per file 19 MB 444 million rasters

Comparison to [Agarwal et al. ‘ 06] • 53 hours (20 ft) • 67 minutes – 3. 4 GHz desktop – 10, 000 RPM disks – scratch space: • quadtree: ~ 8 GB • rasters: ~ 2. 5 GB – 2. 1 GHz laptop – 5, 400 RPM disks – scratch space: • compressed rasters: 270 MB • 86 % of CPU time for expensive RST – remaining 7. 5 hours • constructing quadtree: • finding cell neighbors: 1. 1 hours 5. 6 hours

Demos & Discussion

Demos • streaming generation of DEM • extraction of elevation contours

Streaming DEM Rasterization raw points final DEM

Main Features • order of magnitude faster than method that uses external memory • avoid temporary query structure on disk – data-driven, no gathering of neighbors • interleave input & computation & output – 100% CPU usage while interpolating • no global sort – don’t fight the data, use existing coherence

Future Work • breaklines / barriers (e. g. rivers, roads) – use Constraint Delaunay Triangulation • other interpolation methods – natural neighbors – inverse distance weighting • stream other tasks – thinning of points – compression of raw LIDAR data (useful? ) – hydrological enforcement (streamable? )

Acknowledgements • data – Kevin Yi, Duke • support – NSF grants 0429901 + 0430065 "Collaborative Research: Fundamentals and Algorithms for Streaming Meshes. " – NGA award HM 1582 -05 -2 -0003 – Alfred P. Sloan Research Fellowship

Thank You executables, video, slides, paper, and soon also source code: http: //www. cs. unc. edu/~isenburg/