[comp.graphics] Contour to grid algorithm

yuleat@prlhp1.prl.philips.co.uk (yuleat) (10/25/88)

Can anyone help me with the following problem?

I have data for a relief map in terms of height contours.
The data stored for each contour is :
    1. A height above sea-level
    2. A linked list of points, which are joined by
       straight line segments to make up the contour.

The area mapped by this data is rectangular and all component
points of the contour must be within this area.  As a result,
some contours do not join up, as they start/finish off the map.
All the contour line segments are guarenteed not to touch/cross.

The problem I have is that I wish to convert this relief data
stored in contour form into a form where a uniform rectangular
grid is used.  Each point of the grid will have a spot height.

I have considered a number of "simple" algorithms, but they
all tend to fail in particular cases (mostly near the corners
of the rectangular map area, when contours are sparse).

How do I work out what the height should be at each point?

                Thanks in advance,

        Andy Yule (yule@prl.philips.co.uk)

petersen@uicsrd.csrd.uiuc.edu (11/01/88)

In general I know of no "simple" algorithms that will do a good job,
the basic method is to compute a surface to approximate the contours
and then sample that surface at the desired grid spacing.

The system I worked with used a DTM (digital triangulation model)
to represent the surface, so we never went the final step of re-sampling
the surface at grid points.  However since we already had the routines
to create the DTM from random points, the problem was half solved.

To create a reasonable DTM from a contour map you need to worry about
the flat artifacts that can occur when the triangulation procedure
connects 3 adjacent points on a contour line.  In general this is the
wrong thing to do (maxima, and minima excepted).  Most methods for DTM
construction have a property that we refered to as the "Great Circle
Test", this being that for any three points in the plane fit a circle
to those 3 points, and if any other point lies within that circle then
those 3 point will NOT be connected into a triangle.  Thus the crux of
the method we used was to scan each contour line in sequence and to
compute the circle for each 3 adjacent points, and find the point on
another contour line that was closest to the circle.  If that point
was inside the circle then we are ok.  Otherwise do a linear
interpolation from the new point to the middle of the three points.

This should give a flavor for what is needed to use this method.  It is
resonably accurate, we were able to get 95%-99% volume comparisions,
buy creating a contour map from a DTM, creating a new DTM from the
contour map, and then comparing the two surfaces (depending on the
spacing of the contour lines).  If anyone knows of a simpler of more
elegant method, I would like to know about it.

-Paul Petersen

University of Illinois, Urbana-Champaign
Center for Supercomputing Research and Development

    UUCP:    {uunet,convex}!uiucuxc!uicsrd!petersen
    ARPANET: petersen%uicsrd@uxc.cso.uiuc.edu
    CSNET:   petersen%uicsrd@uiuc.csnet
    BITNET:  petersen@uicsrd.csrd.uiuc.edu

kg0r+@andrew.cmu.edu (Kenneth Gober) (11/04/88)

Dr. Dobb's Journal of Software Tools had an article on this in the November 1987
Graphics Issue.  Some of the references cited include:

Ammeraal, L.  Programming Principles in Computer Graphics.  Chichester, England:
Wiley, 1986
Pavlidis, Theo.  Algorithms for Graphics and Image Processing.  Rockville, Md.:
Computer Science Press, 1982.

I didn't read the article very carefully; I don't know if it can deal with
contours that leave the edge of the map.  You may have to enlarge the map,
artificially connect them, do the conversion, then clip the borders to get your
original area back.

Kenneth Gober
kg0r@andrew.cmu.edu

aramini@apollo.COM (Michael Aramini) (11/04/88)

If the contours were generated from a grid by a standard algorithm, it may be
possible to partially invert the algorithm to at least partially reconstruct
the original grid.

If you don't know the positions of the initial grid points, you might be able to
infer at least some of them, and thus the grid spacing.  Note that if the
standard contour algorithm had been used, all of the contour verticies will fall
grid edges.  Thus a given contour vertex will either have the x coordiate of a
column of grid points or the y coordinate of a row of grid points.  You can't
tell from an individual contour vertex which of its coordinates matches that of
grid points.  However, you can construct a histogram of x coordinates of contour
verticies and a histogram of y coordinate of contour verticies.  Peaks in each
histogram will correspond to coordinate values of grid points.  There won't
necessarily be a peak at every grid point coordinate value, but you can at least
infer the spacing for the peak positions.

Once you know where the grid points are, you can derive what some of the grid
point heights must be.  Note that:

    1. Whenever two or more contour verticies fall on the grid segment between 2
       adjacent grid points, you can compute the height of both ofthose grid
       points by linear extrapolation.  

    2. If only one contour vertex falls between 2 adjacent grid points, and you
       know the height of the one of those grid points, you can compute the
       height of the other by linear extrapolation.

By first applying rule 1 wherever possible, and then by repeatedly applying rule
2 wherever possible, you should be able to compute the heights of a large number
of grid points.

For the rest of the grid points, you can determine their approximate heights by
looking at the heights of the contour lines between which each grid point falls.
Note that in regions of the grid where you must resort to this must be fairly
level (since they are regions where contour lines are fairly sparse), so I would
expect that low order curve fitting methods will probably give good results.

-Michael

yuleat@prlhp1.prl.philips.co.uk (yuleat) (11/09/88)

About a fortnight ago I posted the following problem:

> I have data for a relief map in terms of height contours.
> The data stored for each contour is :
>     1. A height above sea-level
>     2. A linked list of points, which are joined by
>        straight line segments to make up the contour.
> 
> The area mapped by this data is rectangular and all component
> points of the contour must be within this area.  As a result,
> some contours do not join up, as they start/finish off the map.
> All the contour line segments are guarenteed not to touch/cross.
> 
> The problem I have is that I wish to convert this relief data
> stored in contour form into a form where a uniform rectangular
> grid is used.  Each point of the grid will have a spot height.

As I've now had quite a lot of mail on this subject
(for which thanks to those concerned), I shall give
a quick summary.

The proposed solutions fell into two main categories:
    i)  Line based. Interpolation (anything from linear
        to use of splines) along the grid lines.
    ii) Fitting a surface to the contours and 
        determining the point values from that.
The problem in all cases cocerns not knowing the edge
of the map values (and as you can see, below, I have
cheated a bit!)

What I didn't mention in my original posting was the fact
that I really have two goals. 
    i)    To get a prototype system working a.s.a.p.
    ii)   A long term view to get the "best" solution.
As a result I have implemented a simple system that seems
to work and will be investigating improvements/new methods
in the future.

I will briefly describe the algorithm I have used for 
your interest/comments:
    i)    Determine/Get user to supply the spot height of
            the four corners of the "map".
    ii)   Determine any intersections of contour lines with
            the edges of the map.
    iii)  Use the above data to determine values for all the
            edge values. This is done by interpolating using
            the "know" points (i.e. corners + contour
            intersections from i & ii).
With this I can guarantee that every row & column has two
known values, one at each end. In addition, I then :
    iv)   Determine all intersections between contour lines
            and individual mesh columns & rows. This data I
            store in ordered linked lists; one list for each
            row & each column, containing all the "fixed"
            points along that line (including the edge
            values determined by step iii).
    v)    Calculate the values at each mesh point by averaging
            the values obtained by interpolating along both
            the relevant column & row lists.

I hope this makes some kind of sense (it's a bit difficult to
explain without a diagram!).  At present I am using simple 
linear interpolation for steps iii & v but obviously more
sophisticated methods could be used.

Andy Yule (yuleat@prl.philips.co.uk)