ISPRS, Vol.34, Part 2W2, "Dynamic and Multi-Dimensional GIS”, Bangkok, May 23-25, 2001 
Terrain Modelling from Contours 
Christopher Gold and Maciej Dakowicz 
Department of Land Surveying and Geo-Informatics 
Hong Kong Polytechnic University, Hung Horn, Kowloon, Hong Kong 
Tel: (852)2766-5955; Fax: (852)2330-2994 
christophergold@voronoi.com; lsmaciek@polyu.edu.hk 
KEYWORDS: Terrain Modelling; Contours; Skeleton; Interpolation; Delaunay Triangulation; Voronoi Diagram. 
ABSTRACT 
Good quality terrain models are becoming more and more important, as applications such as runoff modelling are being developed that 
demand better surface orientation information than is available from traditional interpolation techniques. A consequence is that poor- 
quality elevation grids must be massaged before they provide useable runoff models. 
Rather than using direct data acquisition, this project concentrated on using available contour data, for two reasons. Firstly, despite 
modern techniques, contour maps are still the most available form of elevation information. Secondly, manual contour tracing has 
imposed a subjective interpretation of the form of the landscape that is lost with automation, yet which is of considerable value. 
The maximum slope of the terrain is perpendicular to the contour, and this permits us to visualize the relationships between pairs of 
contours. With care this may be modelled by triangulation methods, as the spatial relationships can be preserved, although standard 
grid interpolation methods based on "n nearest neighbours" often have problems. However, whenever we have relationships between 
portions of the same elevation contour, such as in peaks, pits or valley heads, our interpretation based on triangle slope is insufficient - 
we get "flat triangles". In this case we need to re-examine our spatial model. 
The usual triangulation approach, the Delaunay triangulation, is effective because it is locally stable - a property based on its dual, the 
Voronoi diagram. These two spatial structures have been much studied by workers in the field of computational geometry - largely in 
terms of efficient calculation, but also in terms of their properties. In particular, recent work on the automatic reconstruction of curves 
from point samples, and the generation of medial axis transforms (skeletons) has greatly helped in the visualization of the relationships 
between sets of boundaries, and families of curves. This provides us with tools to enrich our original contour data for "flat triangles". The 
insertion of skeleton points in these cases guarantees the elimination of all flat triangles. Additional assumptions about the local 
uniformity of slopes, either along or across valleys and other features, give us enough information to assign elevation values to these 
skeleton points. If required, appropriate interpolation techniques may generate an elevation grid for visualization purposes that 
preserves reasonable slopes at all points on the model - even at the data points themselves - and that are faithful to the input data. In 
addition, the algorithms used are only moderately more complex than the underlying Voronoi diagram or Delaunay triangulation. The 
result provides us with a surprisingly realistic model of the surface - that is, one that conforms well to our subjective interpretation of 
what a real landscape should look like. 
In addition, comparisons of the methods used in a variety of 
weighted-average techniques throw a lot of light on the key 
components of a good interpolation method, using three- 
dimensional visualization tools to identify what should be “good" 
results - with particular emphasis being placed on reasonable 
slope values, and slope continuity. This last is often of more 
importance than the elevation itself, as many issues of runoff, 
slope stability and vegetation are dependent on slope and 
aspect - but unfortunately most interpolation methods can not 
claim satisfactory results for these parameters. 
GEOMETRIC PRELIMINARIES. 
The methods discussed here depend on a few fundamental 
geometrical constructs that are now fairly well known - the 
Voronoi diagram and its dual, the Delaunay triangulation, as 
shown in Fig. 1. The first is often used to partition a map into 
regions closest to each generating point; the second is usually 
used as the basis for triangulating a set of data points, as it is 
guaranteed to be locally stable. It may easily be constructed 
using its “empty circumcircle” property - this circle is centred at 
the Voronoi node associated with each triangle. As will be seen 
later, these nodes, and circles, are associated with the skeleton, 
or “medial axis transform”. 
GENERATION OF RIDGE AND VALLEY LINES. 
Amenta, Bern and Eppstein (1998) examined the case where a 
set of points sampled from a curve, or polygon boundary, were 
triangulated, and then attempted to reconstruct the curve. They 
showed that this “crust” was formed from the triangle edges that 
did not cross the skeleton, and that if the sampling of the cun/e 
was less than 0.25 of the distance to the skeleton the crust was 
guaranteed to be correct. Their algorithm consisted of inserting 
all the Voronoi vertices into the diagram. Gold (1999) and Gold 
and Snoeyink (2001) simplified the approach, showing that, in 
each Voronoi/Delaunay edge pair, one edge could be assigned 
INTRODUCTION 
This paper concerns the generation of interpolated surfaces from 
contours. Whiie this topic has been studied by many people 
(including the first author) for over 20 years, this work is 
interesting for a variety of reasons. Firstly, contour data remains 
the most readily available data source. Secondly, valid theorems 
for the sampling density along the contour lines have only just 
been discovered. 
Figure 1: Delaunay triangulation and Voronoi diagram 
Thirdly, the same publications provide simple methods for 
generating the medial axis transform, or skeleton, which 
definitively solves the “flat triangle” problem which often occurs 
when triangulating contour data, by inserting additional points 
from this skeleton. Fourthly, the problem of assigning elevation 
values to these additional ridge or valley points can be resolved, 
using the geometric properties of this skeleton, in ways that may 
be associated with the geomorphological form of the landscape.