2.3 Working with Me

Here are some ideas for projects relating to my research. They might lead to co-authored publishable papers.

2.3.1 Lossy Compression of Elevation Data

Generic lossy image compression algorithms, such as Said & Pearlman's, perform excellently on gridded terrain elevation data. We used 24 random USGS DEMs as test data. Lossless compression required 2.12 bits per point on average, with a range from 0.52 to 4.09. With lossy compression at a rate of 0.1 bits per point, the RMS elevation error averaged only 3 meters, and ranged from 0.25 to 15 meters. Even compressing down to 0.005 bits per point, or 902 bytes per cell, gave an indication of the surface's characteristics. This performance compares favorably to compressing with Triangulated Irregular Networks. If it is desired to partition the cell into smaller blocks before compressing, then even blocks work well for lossless compression. However, for lossy compression, specifying the same bit rate for each block is inadequate, since different blocks may have different elevation ranges. Finally, preliminary tests suggest that the visibility indices of the points are robust with respect to even a quite aggressive lossy compression.

See [#wrflossy##1#, #wrfelev95##1#], and generally ftp://ftp.cs.rpi.edu/pub/franklin/compress for more info on this. The term project would be to extend these ideas.

2.3.2 Triangulated Irregular Network

The Problem:

Compressing a square grid of raster terrain or greyscale data.

One Solution:

Partitioning the square into a net of planar triangles. At each step, find the triangle that is the poorest approximation to the points inside its projection, and split it into 3 triangles.

So What's New?

This method has been around a few years; e.g., I wrote a program to do it in 1973. There are various papers on it. However, there is still room for research.

More Info:

See ftp://ftp.cs.rpi.edu/pub/franklin/tin.tar.gz and related files.

Open Areas

• Test TINs on non geographic data, such as photos. There's lots of data available.

• There is much more geographic data available now, such as in the DEM format. Here a region is represented by a grid of elevations. Test many different cells, like all of them, and gather statistics. Are some cells harder to TIN?

• Do a serious statistical study of the performance. Study accuracy as a function of number of triangles. Does the criterion used to swap diagonals matter?

• Is there a scale-variance? If you use coarser data, how does the accuracy versus size graph change? There is clearly scale-variance in terrain, it is not a scale-invariant fractal.

• Don't tell ex-Sen Proxmire, but the USGS is shamefully exceeding its mandate by providing non-GEOlogical data, e.g., Mars and Venus. One project is to test whether Martian and Venusian (martial and venereal?) data is statistically similar to geological data.

• Make faster on large datasets, perhaps by selecting a random sample of points, TINning them, then inserting the rest.

• Parallelize it.

• Vary some criteria. What about curved surfaces? What if the output elevation can be discontinuous, so long as the max error is ok?

Files

• Look on the RCS in . For a simple run, reading the medical image file , just run there.
• there is DEM cell of elevations for part of the Adirondacks. You'd have to recompile with larger arrays to run it.
• is the lat 42-43N, long 73-72W. It has some header info before the actual elevations. (BTW, compressing it with took 75 CPU minues on an IPC; the compression was 86%).
• More DEM files are on .
• I have a CDROM with data from Mars. If more than one person wants to use it, I might load a few files onto the RCS.
• On the WWW, check out .

2.3.3 Terrain Visibility and LOS

This is based on work for the US Army Topographic Engineering Center, formerly the Engineer Topographic Labs, at Ft. Belvoir. The general problem was to defend an area against enemy air attack by planes or Scuds. Some terrain is more valuable and critical than other terrain. Our tools include, in order of range, cost, and scarcity, these weapons: Patriots, Hawks, and Stingers. Specifically we wish to site these defenses in the best places. There are mobility restrictions, e.g., we can't drive through a swamp, but perhaps can use a helicopter to lift something over it. There are line of sight restrictions; we can only defend against an attacker that we can see, either visually or with radar. First, determining the terrain or airspace visible from a proposed site, and, second, determining the sites with most combined visibility, are major problems. Figure 14 illustrates the problem. Mountain B is the single best point for viewing the most terrain. However, the smallest set of points that can see the whole terrain is , which does not include B.

  
Figure 14: Visibility of a Mountainous Terrain

Currently the siting is done manually, except for a program that runs 360 rays out from a given site to calculate approximate visibility. The problem is to do the siting more automatically.

Clark Ray did a PhD on this[#rayphd##1#]. However, there are many pieces suitable for smaller projects. Some possibilities are these.

• Efficient calculation of visibility from several points without restarting from scratch each time.

• Modification of the joint visibility of several sites when one changes w/o redoing everything.

• Sensitivity analysis: identify inaccessible peaks that would really help if we could get up them. Alternatively, suppose some peaks are more expensive, since we have to use a helicopter. Factor that in.

For more info, see [#frhinbv94##1#], available from me.

Files

I wrote a fast approximate LOS program, , which you can start from. It's in /dept/ecse/compgeom/los/. I'll move in some data files soon.

2.3.4 Image Edge Detection by Line of Sight

The visibility index work described above seems to identify edges in greyscale images. The method is as follows.

1. Consider the intensity of each pixel to be its elevation.
2. Apply the Line of Sight algorithm to calculate a visibility index for each pixel.
3. Transform the visibility indices back to grey-scale values.

There's no theoretical reason for this to work, except that it does. One point to note is that the LOS calculation is nonlinear, and has a step function. Either a pixel can see another pixel, or it can't. This contrasts to most usual image processing operators, which are linear. Since we have more power, it's reasonable that we might have more functionality.

The project is to explore this idea.

1. Try more images of different types.
2. Vary parameters of the LOS algorithm, such as radius of interest, and source and estination elevations above the terrain.
3. Transform the image first, such as by histogram equalization, or inversion. Since the LOS is nonlinear, these will affect the result.
4. Try iterating LOS and other operators.
5. Try rougher, but faster, LOS approximations.

2.3.5 Correlating Two Overlapping FEM Tesselations

In Finite Element Modeling (FEM), we may have two different tetrahedral tesselations of the same object. There may be some data stored with the tetrahedra of one tesselation that we wish to transfer to the other tesselation. E.g., for each output tetrahedron we might determine which input tetrahedra overlap it and construct a prorated value. The ideas used in my map overlay algorithm [#fcmopa90##1#] are relevant here. I describe this 3-D algorithm in more detail in [#fkvo3tp93##1#].

2.3.6 Fast Rational Approximations to Surface Illumination Functions

Consider a bicubic patch being illuminated so that each pixel reflects a different intensity. To display it we need this intensity as a function of x along each scan line. Even though this is a complicated function that contains an inversion, it is smooth. Therefore there should be a fast approximation. Rational, or Pade and Pade-Chebyshev functions appear better than purely polynomial approximations.

A Pade approximation is a formal transformation of the Taylor series into a rational form. A Pade-Chebyshev approximation is a formal transformation of the Chebyshev series. Although the Pade and Pade-Chebyshev approximations contain no more information than the approximations that they were derived from, they are often better. Thus they are a half-way step to full rational approximation.

A rational function is a ratio of a polynomial numerator and denominator. The homogeneous functions used somewhat in CAD are a restricted form of rational functions, which don't use the full power. Rational functions appear to be a powerful and underutilized tool in geometry. Rational approxes are harder to work with, partly since they are not linear and you cannot just linearly combine a set of basis functions. However, they often give much more economical approxes, expecially to functions with kinks, like |x|, or asymptotes, like , than do polynomials.

Here is some info on how much better a rational approx is that a polynomial approx for . would seem to be about the best function for polynomially approximating, so this result is interesting. For N d.f., the best polynomial approx (on some reasonable interval) has an error of . The is because Taylor is suboptimal. The best rational approx, which will have roughly coefficients in both numerator and denominator, has an error of . However this doesn't mean that you can use many fewer terms to get the same error since all these series converge so rapidly.

Before demonstrating the rational approximations by doing the 3D case, we would want to explore something simpler: approximating the tangent to a parametric cubic curve as a function of x. We would use Maple. For example, if the curve is , then we want a simple approximation to each component of the tangent as a function of x. It would be of a form something like this, where '5' is a variable.

Using a spline approximation loses the information that we are approximating a single mathematical function. Since a Pade approximation keeps this information, it might be much more efficient.

The biggest worry is that, although very efficient approximations might exist, it might be too hard to find them. A solution might be a two-step approximation:

• Use a first set of approximations to take the control points of the curve and produce the coefficients of the second approximation.
• Then use the second approximation to find the curve's tangent.

2.3.7 Implement the Chazelle-Edelsbrunner Optimal Line Intersection Algorithm

This algorithm [#chazelleedelsbrunner##1#] finds the K intersections of N line segments in time , which is optimal. It's complicated, but they describe it in great detail. Implement it and test it. Apart from the specific algorithm, this would be an excellent introduction to current computational geometry data structures.

This has been done by people at other universities, and there is a video describing the algorithm, but this is still a difficult enough project. Read the paper before picking this one.

2.3.8 Install Edelsbrunner's Fast Convex Hull Program

We have a Pascal program, written by Herbert Edelsbrunner and his students, that finds the convex hull of a set of points in 3D and 4D. (The convex hull is the smallest polyhedron or polytope that contains all the points.) Get it working here and test it.

2.3.9 Implement my Local Topological Polyhedron Volume Formula

Implement my local topological formulae for finding the volume, area, etc of a polyhedron. Since polyhedra are conventionally in a vertex, edge, face format, first you would have to convert a polyhedron to this format.

2.3.10 Voronoi Diagrams with Barriers

We are given a plane with sources and barriers, and wish to partition it into regions determined by the closest source to each point, with the distance measured around the barriers. Is there a fast algorithm? It seems that wave-marching might work.

2.3.11 Interpolating Irregular Points with Voronoi Diagrams

We want to find a Z value for a new point when given the Z values at an irregularly spaced set of old points. Chris Gold finds a Voronoi diagram on the old points, as shown in figure 15. To interpolate the new point, he temporarily inserts it, 16, and calculates how much area its Voronoi polygon steals from each of the old Voronoi polygons that it overlaps. Then the new Z is a weighted sum of the old Zs, with the weights being the stolen areas. vor1Voronoi Diagram on Four Pointswidth=3in vor2Adding a Fifth Point to the Voronoi Diagramwidth=3in

The project is first to use Maple to calculate the weights algebraically without needing to actually insert the point and calculate intersection lines. The next step is to calculate the gradiant at the new point. It appears the main difficulty will be determining which of many topologies occurs for the new point in each case.