Terrain Elevation Data Structure Operations

Wm Randolph Franklin

Rensselaer Polytechnic Institute

http://wrfranklin.org/

Goal

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

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 more complicated thin plate model minimizes total curvature, similar to fitting a thin sheet of metal to fixed points, while minimizing the energy of bending. Gaussian curvature, commonly used in geometry, is inappropriate here, since it is zero on a ``developable'' surface, such as a bent sheet of paper. Instead we use the (scaled) divergence, or

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

    for the curvature at a point.

    The partial differential equation is:

    zxx2 + 2zxy2 + zyy2 = 0;

    while the corresponding iterative equation on a grid is:

    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)
    For example, see this:

    The following shows a vertical slice thru the interpolating surfaces. It runs thru opposite corners, the better to show any artifacts. The triangulated line shows the original points. The other lines show the cubic, inverse distance, linear, and nearest neighbor interpolations. We see that none of the interpolations is acceptable since all show where the original contours are.

    Figure: Figure: Vertical Slice Thru Four Interpolating Surfaces

    Desirable properties include:

    1. local control
    2. variation minimization
    3. interpolation
    4. conformal
    5. Don't want to see the contours in the result.
    6. Peak inside top contour and valley outside outer contour.
    The new methods described here are based on the PhD work of Gousie (1998).

    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, filling ever-finer contours. Alternatively, we can find intermediate contours one or more times, then do a thin-plate, Hermite-spline the peaks, inverse-distance weight other small gaps, and Gaussian-smooth the surface.

    Gradient Lines

    This method answers the problem of the generated surfaces terracing between contour lines by importing the idea of lofting from CAGD. It proceeds as follows. Generate a first version of a surface by any method. Find its gradient lines. Note that, on each gradient line, we know only the elevations where it crosses the contours. Interpolate its elevations in-between. This process smooths the surface, while keeping the interpolation. Springs may be added to any method, producing an approximation.

    The following figures show three approximation methods applied to a 900x900 piece of the Crater Lake DEM. The first shows the original contour lines, which are seen to be separated by many pixels, which makes the process harder. The next three show the classic thin plate, and our intermediate contour and gradient line techniques. We see that both our techniques are much smoother than the thin plate, and that the gradient method is best. Our longer papers quantify this, and show more examples.
     

    Figure 6. Crater Lake: Original Contours

    Figure 7. Thin Plate Approx

    Figure 8. Intermediate Contour Approx

    Figure 9. Gradient Line Approx

    Overdetermined Laplacian Solution

    The problem with many existing methods is that no information flows across the contours. This causes terracing. Our solution is to make the known data points also be the average of their neighbors. The method is to 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:
    4zij = zi-1,j + zi+1,j + zi,j-1 + zi,j+1

    Second, for the known points, we make an equation that they are equal to their known values:

    zij = hij

    Since there are now more equations than unknowns, none of the equations will be satisfied exactly. That is, no points will be exactly the average of their neighbors, and the known points will not be exactly their known values. If there are K points whose elevations we know, then there are N2+K equations for the N2 unknowns.

    We then do a least-squares solution in Matlab, to give an approximate solution, rather than an interpolation. With a least-squares solution, scaling an equation up makes it more important. Let R be the relative weight of the average-value equations compared to the others. A lower R will cause a more accurate surface, while a higher R will cause a smoother surface. A little inaccuracy goes a long way towards smoothing the surface.
     
     

    Here is a test case constructed to stress any contour interpolation algorithm. It consists of four concentric squares, at a distance of five from each other. The squares' sharp corners facilitate artifacts, while the large gap between adjacent contours allows interpolated surfaces to droop or terrace. Our intuition would desire that a square pyramid with triangular sides result. However, few algorithms are likely to produce such a surface with discontinuities in the slope.

    The next figure shows the surfaces that do result. R is the relative weight of the average-value equations. The four subfigures in order are the original contours, and the approximated surfaces with R=0.1, 1.0, and 3.0.

    Figure: Overdetermined Laplacian PDE, Weighting Different Equations Differently

    This figure shows the same four surfaces cut by a vertical plane cutting from corner to opposite corner, the better to show any artifacts. We see that with R=3, the location of the original contours is completely hidden, which is desirable.

    Figure: Vertical Slice of the Surfaces

    As the following table shows, the average error at the known points in this case is under 3%. That is, by allowing the elevations of known contour points to vary by an average of 3%, we can fit a quite smooth surface to the nested squares. Considering the accuracy standards of most hypsographic data, that may be acceptable.
     

    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 

    Note that our overdetermined solution concept is quite different from putting springs on the data points. That is also an approximation, but it does not try to make the data points to be the average of their neighbors. Absent that, it is impossible to make the surface not show the original contours.

    Various extensions of the overdetermined solution technique are possible. We might calculate the error on each data point, then for the less-accurate points, increase their weights and re-solve. This reduces the max absolute error, but increases the mean. The resulting surface is slightly less smooth. Alternatively, we might guess that the points with large errors represent breaks in the surface slope. Then an reasonable response might be to reduce their weights and re-solve.

    We have tested the overdetermined solution technique on grids of up to 257x257 points, and are working on larger cases.

    3D Connected Components


     
     

    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