GeoStar annual report — June 2006

!!!Frequency of Use of Various Drill Types

We do not propose this specific set of drills as a polished technique to use, yet, but as an example of the experiments that we are conducting in order better to understand this idea.

Degree-1 Scoop

This version differs in the following ways. Instead of a set of overlapping drills, selected in a greedy manner, here there is only one drill, a 7x7 square, and the footprints of its applications partition the dataset. The face of each scoop is a tilted plane whose slope is chosen to minimize the elevation error.

Therefore, a 7x7 Scoop size represents 49 elevations using only 3 coefficients. 7 is not a magic number, and indeed we have experimented with various sizes. However, 7 is good enough for Level 2 DTED cells Large errors are rare and the mean error is very low, often less than 2 meters. The regularity brings simplicity to this representation.

Here are some test results on the `W111N31 DTED II cell, shown in Figure 5, with about 3600'^2^' posts with elevations from 811 to 2882, shown below. (To simplify the programming in these experiments, a few rows and columns were shaved off the sides of the cell to make the size an exact multiple of the tile size.)

Figure 6 shows the effect of using fewer and larger tiles to approximate the cell. Even with 7x7 tiles, about 90% of the posts are represented to an absolute error under 4, or 0.2% of the elevation range.

Another test was performed with 5x5 drills, optimized with Singular Value Decomposition (SVD), on the Lake Champlain West cell, with a data range of 1576. The maximum error was 101, but that was rare. Over 50% of elevation posts were exactly represented. (This was helped by Lake Champlain being large and flat.) The input data size was 1442401 numbers, and the output data size 172800 numbers, for a compression ratio of 8.35:1. The contour images in Figure 7 contrast the original and reconstructed terrain. We see no visible difference.

For another example, we used 7x7 drills to represent a 100x100 detail of the `W111N31 cell, with an elevation range of 512. The reduction in data size is a factor of 16.3:1. Figure 8 shows the error histogram.

Next, we used observer viewsheds to evaluate this representation, as follows.

1. Pick 81 observers,
2. Assume that each observer is elevated 10 units above the terrain elevation at its location.
3. Assume that each observer can see out to a radius of interest that is 300 posts.
4. Assume that potential targets are also 10 units above their local terrain.
5. For each observer, compute its viewshed, or what posts it can see, on the original dataset.
6. Redo the computations on terrain reconstructed with the alternate representation.
7. Measure the area of the symmetric difference between the two sets of viewsheds, as a fraction of the cell's area.

!!!Scoop Compression - Error Tradeoff

The tile3 representation approximates the surface with a regular grid of 3x3 tiles. tile5 and tile7 are similar. The tile75 representation uses first 7x7 and the 5x5 tiles.

Figure 8 and Figure 9 compare the joint viewsheds for the original terrain and the tile7 representation. The visible difference is small.

Sculpting

This is an extension of scooping being studied by Westort. She adapted a generic DTM sculpting tool definition and implementation from Java and C++ to Matlab, including blade, path, and target model matrices. With the help of Tracy and Luk, she significantly improved the efficiency of the code for large datasets. She evaluated the use of Arc/Map drainage and slope functionality for adoption as paths and blades, respectively. Together with Steve Martin, she is developing a stand-alone hydrological modeling solution for path and blade generation for DTM compression and reconstruction.

This subproject has the largest potential payoff, for several reasons. It will be able to represent discontinuities (cliffs), and its nonlinearity is powerful. Dr Westort is creating an initial testbed in Matlab to help us understand the idea. This initial testbed will be of the nature of a terrain reconstruction assistant. The user will specify the terrain, the X-Y paths that the scoop (also called a blade) will follow, and the blade's circular cross-section. The system will run the blade along the path at the lowest elevation that does not cut into the terrain. It will show the error at each post and report the average and max error. We anticipate demonstrating a prototype at the review meeting in San Diego.

This testbed will provide a facility to gain experience with different blade shapes, including variable shapes. That experience will allow us to freeze some decisions, such as whether to use variable shapes, or not. The former choice will require fewer cutting paths, but with the latter choice, each path will take less space to represent. An analogy to splines would be, what is the appropriate polynomial degree to use?

Cutting paths are appropriate for canyon-like terrain. However, reconstructing, say, a mesa, might be better done with other operators. Designing such operators will be the next step, after we can reconstruct eroded terrain.

Postprocessing the reconstructed terrain with a smoothing operator might improve the accuracy at the cost of adding only one bit to the size of the representation ("apply the builtin smoother"). Smoothing will effective if the terrain satisfies some smoothness criterion, which contradicts our design principle that terrain is often discontinuous. Nevertheless, for completeness, this should be checked.

Finally, after obtaining a sufficient understanding of this idea, we might think about formulating and proving some theorems about it.

ODETLAP (Overdetermined Laplacian PDE)

This surface representation and reconstruction technique was pioneered by Franklin, and is being investigated in this project primarily by Xie. It extends the traditional Laplacian (heat flow) partial differential equation (PDE), which interplates unknown points in a matrix of known and unknown elevation posts, as follows. A system of linear equations is constructed with one equation for each unknown elevation, except those in the first or last row or column of the matrix:

4z'_ij_' = z'_i-1,j_' + z'_i+1,j_' + z'_i,j-1_' + z'_i,j+1_'

Choosing appropriate equations for the unknown border points is tricky, for some uninteresting technical reasons. The resulting sparse system easily can be solved by, e.g., Matlab, for 1201x1201 matrices. We call this technique DETLAP for exactly `DETermined `LAPlacian system. ODETLAP extends DETLAP as follows.

1. One instance of the above equation is created for every (nonborder) point, whether known or unknown.
2. Each known point also induces a second equation: z'_ij_' = h'_ij_'
3. Each Border point induces a slightly different equation or two.

The resulting system of linear equations is overdetermined, i.e., inconsistent, so no exact solution is generally possible, but only a least square approximation, as follows. If the system is

Ax=b

we instead solve

Ax=b+e

so as to minimize the error e. The solution is

x = (A'^t^'A)'^-1^'A'^t^'b

The two types of equations may be given different relative weights, depending on the relative importance of accuracy versus smoothness. Matlab easily processes 400x400 matrices of elevation posts. More specialized techniqes can handle larger systems. See ¹, which uses a Paige-Saunders conjugate gradient solver with a Laplacian central difference approximation solver for the initial estimate.

ODETLAP has many advantages, which are not, to our knowledge, combined in any other surface approximation technique.

1. It will infer a local maximum, i.e., a mountain top, inside a set of nested contours. Therefore, it does not generate artificial buttes all the time. This property can be surprising because points interpolated with a Laplacian PDE must fall inside the range of boundary value points. However, we're not doing an exact Laplacian PDE, but an extension of it. Therefore, we can generate values outside the range of boundary value points.
2. The generated surface doesn't droop between contours. For reasonable weighting parameters, the contours are invisible in the resulting surface.
3. It will utilize isolated data points, if available.
4. If will interpolate broken contours. This is not a property of methods that fire out lines until they hit a contour.
5. The resulting surface is at least approximately conformal. That is, if a conformal mapping is applied to the input points the resulting surface is that conformal map of the surface computed from the untransformed points. Therefore, sets of nested kidney-bean contours are handled reasonably. Again, this property is not shared by all surface approximation methods.
6. ODETLAP will conflate inconsistent data, with user-defined weights, producing one merged surface. This allows one to combine a large, low-resolution dataset with one or more small high-resolution datasets covering only parts of the large dataset.

Figure 11 shows a stress test of ODETLAP. The test data is a set of four nested square contours. A surface is fit to them with ODETLAP run three times with different values of the smoothness parameter, R. Larger values for R mean a smoother but less accurate surface. The surfaces are plotted with an oblique projection so that the silhouette edge is obvious. The lower right surface shows no evidence at all of the original square contours, and interpolates a smooth moutain top in the center. The maximum absolute elevation error is 12% of the elevation range, and the mean only 2.7%. These errors would be much smaller for smoother contours.

The impact of ODETLAP for this project is that surprisingly few points are needed to generate a good surface.

We tested ODETLAP on a mountainous 400x400 section of the Lake Champlain W DEM, with elevation range 1378m, by selecting every Kth post in both directions and using ODETLAP to fit a surface (with R=0.1). The following table show that average absolute errors. Data set size % is the ratio of posts used to reconstruct the ODETLAP surface compared to the original number. It is 1/K'^2^'.

||!ODETLAP Compression - Error Tradeoff

That is, with ODETLAP we reduced the data size by a factor of 100, at a cost of an average error of 0.46%

Choices for Input Points

Xie has have been mainly working on ODETLAP in the past year. To be more specific, he is trying different sampling methods on some 400x400 elevation matrices and then using ODETLAP to reconstruct the surface.

We have tried three different sampling methods.

1. The first technique, the regular method is just regular sampling, i.e., choosing 1 point for each K points in both horizontal and vertical directions. The compression ratio is 1/K'^2^'.
2. The second technique, the TIN method, is to use output points from the TIN program, which iteratively selects points furthest from the current surface.
3. The third technique, the improved TIN method, operates as follows.
   1. Run TIN to produce a predetermined number of important points.
   2. Use ODETLAP to reconstruct the surface from those points. Note that the surface will pass near, but not exactly through, the points.
   3. Find the worst points in that surface, that is the points whose reconstruction is farthest from their correct value.
   4. Add those points to the set of points from TIN, and rerun ODETLAP to compute a better surface.

Here is a comparison of the three methods on three 400x400 test cells, with ODETLAP smoothness parameter R=0.3. Both the regular and TIN methods used 1000 points. The refined TIN method added the 100 worst points from the TIN method. For each test, we measured the average absolute error and the maximum absolute error. The column headings are as follows.

||!Abbreviations Used in Following Table

||!ODETLAP Variations Experimental Results

The TIN method produced a somewhat worse average error, but a much better maximum error. That is, the TIN method produced a much better conditioned surface. The refined method was slightly better then the TIN method under both the average and maximum metrics. Modifying the refined method to add even more bad points would presumably make it even better. That is an area for future research.

Figures 12 to 14 compare the regular and TIN methods, each with 100 points, on the `W111N3110 dataset. Note how the TIN method captures small features better. Note also how well even 100 points represent the original surface.

Figures 15 to 17 show an example of approximating a more complex surface with fewer points. Even only 36 TIN points capture the broad essence of the surface.

We have also done some preliminary experiments on a new sampling method, namely, the visibility method. This method fits a surface to the points of greatest visibility.

Filling in Missing Data Holes

This topic was pursued because of sponsor interest expressed at the review meetings.

Terrain (elevation) data sometimes has holes or regions of missing data, perhaps because of the production technology. This section compares some algorithms for filling in a circle of radius r=100 of missing data. This missing region is large enough to stress-test the various algorithms, most of which are unacceptably bad. Any algorithm that succeeds here can also handle smaller holes.

The two best algorithms, which produce realistic results, are a thin plate PDE, and an overdetermined Laplacian PDE. The thin plate PDE is faster, while the overdetermined Laplacian PDE may also be used to smooth known data.

The test data and environment are as follows.

1. The terrain format is an array of elevations.
2. The test data sets are drawn from the 1201x1201 Lake Champlain West USGS level 1 DEM, seen in Figure 18. It combines part of the Adirondack Mountains, the highest in NYS, with Lake Champlain and its lowlands. North is to the right.

1. The region of missing data, which needs to be interpolated, is always a circle of radius r=100.
2. Each test missing circle was embedded in a 400x400 square of data. That is relevant only for the two interpolation algorithms that use data beyond the boundary of the missing data. Those are the two best algorithms, overdetermined Laplacian PDE and thin plate PDE. They actually use data within 2 posts of the circle of missing data. The poorer algorithms use only data on the circle.

Figure 19 nicely summarizes the results. The data outside the black circle (r=100) is real; the data inside was interpolated by our overdetermined Laplacian PDE. Although the thin plate PDE produces similar results, our method has other advantages that that lacks, such as the ability to conflate inconsistent data.

1. Note the realistic appearance of the contour lines.
2. Note the interpolated local mountain top near the bottom just inside the circle. That was caused by the large slope of the known terrain just outside the circle. This interpolated feature is remarkable since few interpolation algorithms can produce a value that is outside the range of the known values.

More detailed results are available at a web site with a javascript program allowing the user to choose from various preprocessed datasets, ². That presents the results of 6 interpolation algorithms on 20 different sections of Lake Champlain West. Each cell of the table shows six algorithms on one dataset:

R=.1 overdetsol()

This is our overdetermined Laplacian PDE. R=.1 means that the equations for all points being equal to the average of their neighbors are weighted 0.1 as much as the equations for the known points being equal to their known values. That results in no visible smoothing of the known data.

detsol()

This is the Laplacian PDE for the unknown points.

thinplate()

This is the thin plate PDE for the unknown points.

Matlab nearest

This is the Matlab builtin nearest point interpolation algorithm, chosen because it is well known.

Matlab linear

This is the Matlab builtin linear interpolation algorithm.

Matlab cubic

This is the Matlab builtin cubic interpolation algorithm.

A fourth Matlab interpolation algorithm is not shown since it failed. Figure 20 compares the different techniques.

Planned future work on this topic is as follows.

1. Using the builtin Matlab overdetermined linear equation solver for our overdetermined Laplacian algorithm is very slow. We are currently researching faster solution techniques.
2. This algorithm can interpolate general classes of data, such as contour lines, possibly broken, and isolated points, perhaps derived from a TIN. We plan to study this in more detail.

Multiobserver Siting

This is prototype software originally written by Franklin, and extended Christian Vogt as a masters thesis. It is a set of C++ programs that exchange data in files and are controlled by Linux shell scripts and makefiles. The toolkit works, albeit being somewhat complicated and messy. It can quickly compute high resolution viewsheds as demonstrated in Figure 21.

Dan Tracy is using the siting toolkit to evaluate our compression. The following figures demonstrate the toolkit's output. The test dataset is the Lake Champlain West DEM. A particular radius of interest and observer and target height are picked and the 60 observers are sited, so that their viewsheds will jointly maximally cover the cell. See Figure 22. Figure 23 adds an intervisibility requirement. That means that the observers must form a connected graph, if there is a graph edge between pairs of observers that can see each other. This allows messages to be passed from observer to observer, which is useful if they are radio base stations.

Of course, the area of the joint viewshed is smaller since the individual viewsheds must overlap, as quantified in Figure 24.

RPI has issued a subcontract to Environmental Systems Research Institute, Inc. to productize our siting toolkit it as an `ArcGIS module. A copy of the contract was sent to DARPA/DSO to check that the government will be able to use this module for free.

Effect of Reducing Resolution

We earlier conducted initial studies of the effect of lowering the horizontal or vertical resolution of the dataset before performing the multiobserver siting. We observed that lowering the horizontal resolution lowers observer siting quality. However, lowering the vertical resolution does not have as large as effect.

The visibility, computed on the lower resolution data, is too high. In other words, when an observer's viewshed is computed on a lower resolution matrix of elevation posts, it may have a larger area (after scaling) than when computed on the original, high resolution, data. That means that it is important to compute the viewsheds on the best available data.

Careful consideration was given to the proper procedure for evaluating the effect of reduced resolution; the chosen procedure is given in Figure 25.

Our point is this. When evaluating a set of observers, it does not matter where they are, but how much they can see. There may be many sets of observers, each as good as the other. Therefore, when performing a multiple observer siting on the reduced resolution dataset, we compute their joint viewshed, and compare its area to to that of the joint viewshed of multiple observers sited on the original, best quality, data. Also, when computing the viewshed area of the observers sited on the lo-res data, we do the computation on the hi-res data, since that is more accurate.

Figures 26 and 27 show some experiments with reducing the horizontal and vertical resolution. The visibility index is the fraction of the cell covered by the joint viewshed of a fixed number of observers sited with our testbed. The hi res line shows the result of using the original data. The lo-res line shows the siting and viewshed computation formed on the lo-res data, for three different reductions in horizontal resolution. The transferred res line shows the effect of siting the observers on the lo-res data, but then computing their joint viewshed on the hi-res data. It is this line that should be compared to the hi res line. The extent to which it is smaller shows the deleterious effect of reducing the data resolution.

The difference between the lo res and transferred res lines shows the error caused by evaluating viewsheds on lo-res data. That the lo res line is higher shows that the error is biassed up. There may be no deep significance to this; it may be an artifact of how we interpolate elevations between adjacent posts.

Figure 28 shows how imprecise viewshed computations can be. The test data is the Lake Champlain West cell; the observer is near the lower left corner, on Mt Marcy. A particular observer and target height are chosen, and then the visibility of each point in the cell is computed by running lines of sight from the observer. When the line runs between two posts with elevations A and B, and some elevation must be computed there, one of four rules is used.

1. min(A,B)
2. mean(A,B)
3. linear interpolation between A and B depending on their relative distances from the line of sight.
4. max(A,B)

The choice of interpolator affects the visibility of one half of all the points in the cell. The black region is hidden for all the interpolators. The dark grey region is hidden for three interpolators, the light grey region for one, and the saturated color region is visible for all of the interpolators. In other words, the visibility of (at least!) one half of the cell is actually unknown.

Smugglers' Path Planning

A path planning program was created by Tracy. Given a cost metric for the terrain, a start point, and an end point, the program computes the path between the two points that minimizes the line integral over the cost metric by implementing the A* algorithm. In particular, with a strictly binary go/no go cost metric (e.g. avoid the area visible by the planted observers,) the Euclidean distance of the path is minimized. We call this smugglers' path planning.

Future work may entail computing other error metrics related to the visibility and path planning.

Using Smugglers' Path Planning to Evaluate Alternate Representations

The purpose for researching alternate terrain representations is to use the terrain for some application, such as path planning. Since the alternate representation is good only so far as it supports such applications, we have started to use them to evaluate our representations, in addition to using elevation error.

The assumption is that the cost of a path between some source and goal depends on the how visible that path is to a set of observers. If our client is a smuggler, then low visibility is desirable. Conversely, if our client is using a radio and the observers are base stations, then high visibility is good.

Initially, we are computing paths that stay out of the joint viewshed of a set of observers optimally sited to see as much as possible of the terrain. Finding no suitable implemenation on the web, we implemented this ourselves. (The main difficulty is that the cost matrix is large and sparse. A second problem is that it is impossible explicitly to compute and store the distance between every pair of mutually visible points; remember that we are not working with toy-sized datasets.) Here is one test of our techniques.

Figure 29 shows the test dataset, `W111N31, 3595x2595 posts with 12,924,025 coefficients (elevations). The elevation range is 2071.

The alternate representation is to use 7x7 scoops. The reduced the number of coefficients to 791,267, that is, by a factor of 16:1. The mean absolute elevation error was 1.7, or 0.1%. There is no visible difference in the image of this representation.

Next we optimally sited 324 observers, with radius of interest 100, to maximize their joint viewshed, shown in Figure 30.

Then we picked a source and goal that would have a complicated path between them, computed the shortest path that avoided the joint viewshed, and plotted the terrain, viewsheds, and path together in Figure 31.

Our first evaluation of the alternate representation went as follows. We transferred the 324 observers back to the original terrain and recomputed their viewsheds. Then we overlaid the path on them and counted how many path pixels were inside the more accurate viewsheds. The path was 4767 points long; only 14, or 0.3% were visible, as shown in the Figure 32.

Coding the Representation's Coefficients

The representations discussed in this report all reduce the terrain to a set of coefficients. The raw dataset is a matrix of elevations. The degree-1 scoops are triples of coefficients of planes, and so on. As a final step, these coefficients must be coded, or reduced to a minimal number of bits. There is a wide range of coding research to draw on, such as JPEG when the coefficient matrix looks like an image. However, we hope to do better. This is future research.

Comparison of Actual Achievements with Goals

This project's progress is aligned with the proposal except for the following change. Because of sponsor feedback, we are researching path planning, instead of hydrology.

Cost Overruns

None.

Interactions Between RPI and Sponsors

In addition to written reports, there have been several face-to-face meetings during this contract.

1. March 2005: kickoff meeting in Savannah GA, including presentation by WR Franklin.
2. June 2005: site visit by Paul Salamonowicz and Ed Bosch of NGA to RPI.
3. Nov 2005: Program review meeting in San Diego CA, including presentation by WR Franklin.
4. April 2006: Program review meeting in Snoqualmie WA, including presentation by WR Franklin.

In addition, a password-protected website has been established, ³, to communicate results in a less formal manner.

Publications

Papers presenting results from this project have already been accepted at the following scientific meetings.

1. Two Novel Surface Representation Techniques, W. Randolph Franklin, Metin Inanc, and Zhongyi Xie, Autocarto 2006, June 2006
2. Tradeoffs When Multiple Observer Siting on Large Terrain Cells, W. Randolph Franklin and Christian Vogt, 12th International Symposium on Spatial Data Handling, July 2006
3. Alternative Sculpting Hypotheses for Terrain Data Compression, C. Westort, W. R. Franklin, Proceedings of SPIE Vol. 6313 Advanced Signal Processing Algorithms, Architectures, and Implementations XVI, Tuesday-Wednesday 15-16 August 2006, San Diego.
4. Compressing Terrain Datasets Using Segmentation, W. R. Franklin, M. Inanc, Proceedings of SPIE Vol. 6313 Advanced Signal Processing Algorithms, Architectures, and Implementations XVI, Tuesday-Wednesday 15-16 August 2006, San Diego.

More papers are in preparation for other meetings and then journals.

Press and Blog Mentions

With the prior approval of DARPA, RPI announced the contract with a press release that was mentioned in the following RPI places.

1. RPI Press release, 10/31/2005
2. RPI School of Engineering Press Release
3. Rensselaer School of Engineering News
4. Improving Terrain Maps, Rensselaer Alumni Magazine Winter 2005-06

That press release was picked up by these news sources.

1. Better terrain maps of Earth... and beyond, in Roland Piquepaille's Technology Trends, 5 nov 2005
2. ZDNet, 11/5/2005
3. Surfwax Government News
4. Defence Talk, 11/1/2005
5. Advanced Imaging (no longer online)
6. ACM Technews 7(862), 11/2/2005
7. Interview on WGY AM-810 radio, 11/3/2005