Terrain Elevation Data Structure Operations

Wm Randolph Franklin
Michael B Gousie

Goal

  • Operations:

    • How the Parts Fit Together




    Theme

    What's New?

    ... since people have converted data for several decades
    now.
     

    Historical Context

    Research Program Components

    Conversion from DEM to TIN

    Example - Lake Champlain West DEM

    1. 1201x1201
    2. Can process to over 1,000,000 triangles.
    3. However, max error < 1 after 600,000 triangles.
    4. Here is original, then after 1024, 2048, 4096, 8192, 16384, 32768, 65536 points inserted.


    Extensions:

    Visibility

    Extensions:

    1. Unshaded (): almost certainly visible.
    2. Lightly shaded (): probably visible.
    3. Darkly shaded (): probably hidden.
    4. Black (): almost certainly hidden.

     Sample viewshed output
  • Following are two images:
      1. some sample terrain from South Korea, color-coded by elevation, and
      2. the visibility indices of every point, colored similarly.

    Sample Terrain; Visibility Indices

    Extensions:


    Siting Observers


    Programs viewshed and vix may be used to find a set of observers that jointly can see every point as follows.

    1. Use vix to list all the points by visibility index, and hence, to find the most visible point. Place the first observer, O1, there.
    2. Use viewshed to find the points that O1 cannot see.
    3. Filter the sorted list of points to delete points that O1 can see.
    4. Find the most visible point that O1 cannot see; that is the second observer, O2.
    5. Repeat until the set of observers can see every point.

    DEM Compression

    BPP  RMS Error, m 
    5.234  0.0 
    1.0  6.0 
    0.1  51.7 
    0.03  148 
    0.01  535 
    0.0  1096 

    Original

    Lossy Compression to 1 Bit/Point, RMSE=6.0 meters

    Lossy Compression to 0.1 Bit/Point, RMSE=51.7 m

    Lossy Compression to 0.03 Bit/Point, RMSE=148 m

    Lossy Compression to 0.01 Bit/Point, RMSE=535

    Interpolating from Contours to DEM

     

    Prior Art


    Figure: Interpolating 3 Nested Squares with Inverse Distance Weighting, and 3 Ways with Delaunay Triangles

    Interpolation & Approximation Introduction

    Desirable properties: The math is counterintuitive, and the desired properties mutually contradictory.
     
     
     

    Goodness Criteria


    Ideal:

    Next best:
     

    Partial differential equations (PDEs)

     Laplacian, or heat-flow equation,
    zxx+zyy = 0,

    where

    zxx = d2z/dx2
    solved by iteration on a grid,

    4zij = zi-1,j + zi+1,j + zi,j-1 + zi,j+1.

    Problem:  demonstrates the terracing artifact.

    The thin plate model minimizes total curvature.

    zxx2 + 2zxy2 + zyy2 = 0;

    20zij = 8(zi-1,j + zi+1,j + zi,j-1 + zi,j+1) - 2(zi-1,j-1+zi-1,j+1+zi+1,j-1+zi+1,j+1) - (zi-2,j + zi+2,j + zi,j-2 + zi,j+2)

    Intermediate Contours

    Our first new method repeatedly interpolates new contour lines between the original ones, somewhat similar to the medial axis. If two adjacent contour lines are: A, at elevation a, and B, at b, then we interpolate at elevation (a+b) thus.
    1. Pick a point on A.
    2. Find the closest point on B (approximating the gradient).
    3. The midpoint is on the intermediate contour.
    Our contrarian philosophy is that we don't find the elevation at certain points, but rather find points with a certain elevation.

    Figure: Crater Lake: Original Contours and Intermediate Contour Interpolation

    The Maximum Intermediate Contours (MIC) method iterates the process.

    Gradient Lines

    1. like lofting in CAGD.
    2. Calculate gradient lines.
    3. Refine their elevations.
    Example: 900x900 piece of the Crater Lake DEM.
     

    Figure 6. Crater Lake: Original Contours

    Figure 7. Thin Plate Approx

    Figure 8. Intermediate Contour Approx

    Figure 9. Gradient Line Approx

    Overdetermined Laplacian Solution

    1. The problem with many existing methods is that no information flows across the contours. This causes terracing.
    2. Make the known data points also be the average of their neighbors.
    3. assume that there are N2 unknowns; i.e., even the known points are unknown. There are now two types of equations. First, for all points:
      4. 4zij = zi-1,j + zi+1,j + zi,j-1 + zi,j+1
    4. Second, for the known points, we make an equation that they are equal to their known values:
      5. zij = hij
    5. Do a best-fit solution.

    Figure: Overdetermined Laplacian PDE, Weighting Different Equations Differently


    Table: Error versus Weight 
    R Max % Error Mean % Error
    0.1  0.27  0.01 
    1.0  5.5 0.6 
    3.0  12  2.7 

    Proposed Future Work

    Short Term

    Medium Term

    Long Term

      1. Investigate the general question of the proper representation of terrain elevation data.
      2. If using TINs, use triangular splines, not just piecewise multilinear flat triangles.
      3. Consider the idea of a conceptually deeper representation, based on the geomorphological forces that created the terrain. Devise a basis set of operators, such as uplift, downcut, etc. Deduce the operators that created the particular terrain under consideration, and store them.
      4. Slightly differently, represent the terrain by the features that people would use to describe it. This idea is many decades old. Can we get it to work? The problem is that you can't just say that there is a hill over there, you have to specify the hill in considerable detail, which would seem to take more space than simply listing all the elevations with a grid or TIN.
      5. Extend to geological data structures. Plan excavation strategy for open pit mine.
      6. The major creative force for terrestrial terrain is water erosion. This does not apply to the moon, Venus, and, to the same extent, to Mars. Therefore those surfaces might be statistically different. What impact, if any, would this have on the speed and output distribution of our algorithms and data structures?

    How the Parts Fit Together