Generating Raster DEM from Mass Points via TIN Streaming

1. Generating Raster DEMfrom Mass Pointsvia TIN Streaming LIDAR points raster DEM Martin IsenburgUC Berkeley Yuanxin LiuUNC Chapel Hill Jonathan ShewchukUC Berkeley Jack SnoeyinkUNC Chapel Hill Tim ThirionUNC Chapel Hill

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

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

4. Streaming TIN output StreamingComputationof Delaunay Triangulation ImmediateRasterizationof Streaming TIN output final DEM to disk small bufferin main memory store elevation rasters to temporary files(grouped by rows) DEMGenerationviaTINStreaming read points from disk

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

6. Related Work

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

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

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

10. 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

11. 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

12. 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

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

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

15. Spatial Finalization Tags • inject tags intothe point streamthat inform whengrid cells are freeof future points finalized region future pointsonly arrive inthese cells • use tags to certify triangles as beingpart of the finaltriangulation andoutput them

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

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

18. spfinalize -i points.raw -ospb | spdelaunay –ispb –osmb -ospb | tin2iso –ismb –o dem.bil Streaming Rasterization Pipeline

19. The Finalizer(spfinalize)

20. Spatial Finalization of Points  compute bounding box

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

22. 5 4 9 8 6 1 6 7 4 7 3 7 9 7 7 5 8 7 8 9 2 7 7 8 7 3 7 8 7 6 8 7 7 4 8 8 9 7 9 8 2 6 9 6 5 8 2 Spatial Finalization of Points 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

23. 1 4 9 7 6 0 3 5 4 7 3 7 9 7 7 5 8 7 8 9 2 7 7 8 7 3 7 8 7 6 8 7 7 4 8 8 9 7 9 8 2 6 9 6 5 8 2 Spatial Finalization of Points 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

24. 0 4 9 2 6 1 0 5 4 5 3 6 9 7 7 5 8 7 8 9 2 7 7 8 7 3 7 8 7 6 8 7 7 4 8 8 9 7 9 8 2 6 9 6 5 8 2 Spatial Finalization of Points 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

25. 0 3 9 0 4 1 0 0 4 3 3 3 9 7 7 5 8 7 8 9 2 7 7 8 7 3 7 8 7 6 8 7 7 4 8 8 9 7 9 8 2 6 9 6 5 8 2 Spatial Finalization of Points 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

26. 0 7 1 0 4 3 1 2 7 1 1 7 7 9 2 7 7 8 7 3 7 8 7 6 8 7 7 4 8 8 9 7 9 8 2 6 9 6 5 8 2 Spatial Finalization of Points 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

27. 0 2 1 0 4 1 6 3 2 5 1 7 6 4 5 8 6 7 4 8 8 8 7 9 8 2 6 9 6 5 8 2 Spatial Finalization of Points 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

28. 0 2 1 0 4 1 1 3 2 3 2 3 7 2 5 9 5 4 7 2 Spatial Finalization of Points 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

29. 0 2 1 0 4 Spatial Finalization of Points 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

30. 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

31. 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

32. finalization tags Unfinalized Region as Quadtree • as tags finalize cell after cell, maintain remaining cells in adaptive quadtree  used for efficientoverlap tests bythe triangulator

33. The Triangulator(spdelaunay)

34. Modified Incremental Delaunay • insert points • when reading afinalization tag • find certified triangles • output them to disk or pipe • deallocate data structure

35. Modified Incremental Delaunay • insert points • when reading afinalization tag • find certified triangles • output them to disk or pipe • deallocate data structure

36. Modified Incremental Delaunay • insert points • when reading afinalization tag • find certified triangles • output them to disk or pipe • deallocate data structure

37. The Rasterizer(tin2dem)

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

39. Results

40. Neuse-RiverBasin

41. Pipeline Details for the 20 ft DEM 10 MB60,388 triangles 6 min 5 MB512 x 512 grid 6 min 35 MB 2,249,268 points 43 min 8 min 270 MB 239 files 19 MB444 million rasters 46 MB128 rows per file

42. Comparison to (20 ft) [Agarwal et al. ‘06] • 53 hours • 3.4 GHz desktop • 10,000 RPM disks • scratch space: • quadtree: ~ 8 GB • rasters: ~ 2.5 GB • 67 minutes • 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: 1.1 hours • finding cell neighbors: 5.6 hours

43. Demos & Discussion

44. Demos • streaming generation of DEM • extraction of elevation contours

45. final DEM Streaming DEM Rasterization raw points

46. 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

47. 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?)

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

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