ETL-0472

AD-A 188 108

Contour-to-Grid Interpolation with "41

~p

Nonlinear Finite Elements:

A Feasibility Study

Betty Mandel D 3.0

Javier Bernal S0T !)

Christoph Witzgall

U.S. Army Corps of Engineers

Engineer Topographic Laboratories

Fort Belvoir, VA 22060-5546

Center for Applied Mathematics

National Bureau of Standards

Gaithersburg, MD 20899

September 1987

APPROVED FOR PUBLIC RELEASE; DISTRIBUTION IS UNLIMITED.

~~ 87 12
UNCLASSIFIED

SECURITY CLASSIFICATION OF THIS PAGE

IForm Apporoved

REPORT DOCUMENTATION PAGE oMB No. 0704-0ed

lI& REPORT SECURITY CLASSIFICATION lb. RESTRICTIVE MARKINGS

UNCLASSIFIED

2a. SECURITY CLASSIFICA';ON AUTHORITY 3. DISTRIBU ION /AVAILABILITY OF REPORT

2b. DECLASSIFICATION /DOWNGRADING SCHEDULE Approved for public release

distribution is unlimited.

4. PERFORMING ORGANIZATION REPORT NUMBER(S) 5. MONITORING ORGANIZATION REPORT NUMBER(S)

ETL-0472 ETL-0472

6a. NAME OF PERFORMING ORGANIZATION 6t. OFFICE SYMBOL 7a. NAME OF MONITORING ORGANIZATION

U.S. Army Engineer Topographic (if applicable) Center for Applied Mathematics

Laboratories INatinal Bureau of Stndardqs

6c. ADDRESS teity, State, and ZIP Code) 7b. ACDRESS (City, State, and ZIP Code)

Fort Belvoir, VA 22060-5546 Gaithersburg, MD 20899

1a. NAME OF FUNDING/SPONSORING 8b. OFFICE SYMBOL 9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER

ORGANIZATION (if applicable)

St ADDRESS (Chi, State, and ZIP Code) 10. SOURCE OF FUNDING NUMBERS

PROGRAM PROJECT TASK WORK UNIT

IELEMENT NO. INO. NO. 1ACCESSION 110.

11. TITLE (Include Security Classification)

Contour-to-Grid Interpolation with Nonlinear Finite Elements: A Feasibility Study

12. PERSONAL AUTHOR(S)

Betty Mandel, Christoph Witzgall, Javier Bernal

13a. TYPE OF REPORT 13b. TIME COVERED 14. DATE OF REPORT (Year,0Month,Day) ". PAGE COUNT

Technical FROM Jan86 TO Se September 1987 71

16. SUPPLEMENTARY NOTATION

17. COSATI CODES ". S;UBJECThTERMS (Contiue on reverse if necessary and identify by block -umber)

FIELD GROUP SUB-GROUP cartography, Cloug-T ocher element, contour line, Delaunay diagram, digitized data, finite element methQJ, line generalization, surface

interpolation, synthetic surface, terrain modeling, Thiessen diagram,

19. ABSTRACT (Continue on reverse if necessary and identify by block number), erance an, triangulation, oronoi diagram.

The U.S. Army Engineer Topographic Laboratories and the National Bureau of Standards investigated the

computational feasibility of developing a contour-to-grid algorithm by: (1) sampling a large set of digitized cartographic contour data with a tolerance band technique, (2) its subsequent triangulation with a Voronoi (Delaunay,

Thiessen) method and (3) the construction of an internolating smooth synthetic surface from whi'ch grids can be

generated at any given interval. This report presents a discussion on the results obtained.

20. DISTRIBUTION /AVAILABILITY OF ABSTRACT 21. ABSTRACT SECURITY CLASSIFICATION

[YUNCLASSIFIED/UNLIMITED [0 SAME AS rPT. 0[ OTIC USERS UNCLASSIFIED

22a. NAME OF RESPONSIBLE INDIVIDUAL 22b. TELEPHONE (Include Area Co&p) 22c. OFFICE SYMBOL

E. James Books 202-355-3039 -FETI.-IM-T

DD Form 1473, JUN 86 Previous editions are obsolete. SECURITY CLASSIFICATION OF THIS PAGE
PREFACE

Our interest in this area was stimulated by Dr. Stephen Grotzinger and Mr.

Douglas Caldwell. We are grateful to Mr. Arthur Noma of the Defense Mapping

Agency Hydrographic/Topographic Center (DMAHTC) for providing the data for these

demonstration efforts, and for his continuing interest and support. We

appreciate Mr. Stephen Rogers' assistance in the prepatation of plots needed to

represent and update the data. Finally, we thank Mr. Steve Hasenfus, Mr. Larry

Cook and Mr. Barry Holecheck for being there whenever an "emergency" occurred.

or 0

ý~es oi ý i or

NTIS CkA&I

DTIC TAU

by .. ....... ....... ID D',,.t :bJ,'jc• : / ,,

i

ILA
CONTENTS

1 INTRODUCTION 1

2 DIGITAL DATA STRUCTURE 4

2.1 Digital Graphic Recorder Data 4

2.2 Editing the Input Data 7

3 TEINNING BY TOLERANCE BAhD ALGORITHMS 9

3.1 Definitions and Overview 9

3.2 A One-Pass Version of the Reumann and Witkam Method 11

3.3 Determining Contour Tangents 19

4 VORONOI TRIANGULATION OF LARGE SETS OF CONTOUR POINTS 21

4.1 Definitions and Algorithms 21

4.2 Implementation 28

4.3 Results of a Compotational Experiment 33

5 SURFACE GENERATION WITH CLOUGH-TOCRER ELEIENTS 35

5.1 The Clough-Tocher Element 35

5.2 Surface Generation over a Triangulated Region 49

5.3 Finding the Right Triangle 53

6 RESULTS AND CONCLUSIONS 58

6.1 Setting Up the Experiment 58

6.2 Results 62

6.3 Conclusions 65

REFERENCES 67

ii
1. INTRODUCTION

This report concludes the second phase of a three-phase collaborative effort

between the U.S. Army ingineer Topographic Laboratoriesc7WI4)and the National

Bureau of Standardso(NBS). This phase is the object of Interagency Agreement

E8786K041 cover the time period from January to September 1986. The first

phase of 7s effort took place from May to December 1985 under Interagency

r4ement E8735K137. The third phase of this effort is in the planning stage.

•The main objective of this effort is to develop and implement a contour-togrid algorithm, that is, an algorithm capable of converting digitized contour

data into Digital Terrain Elevation Data (DTED). The approach is to create a

synthetic terrain surface based on digitized contour information and then to

calculate any grid from that surface. Five tasks are to be performed:

02)Edlt a digitized contour data set (input)

oZ)Thin (reduce, sample) the input data,I

c )Triangulate the selected points,

S4)Generate a smooth synthetic surface '

p5)enerate the desired grid.)

•The first phase, which was completed in December 1985, addressed the first three

tasks., A describtion of-the york done and the results obtained appear in

Witzgall, Bernal and Mandel (85). )The second phase, to which this report

refers, elaborates on tha previous results and includes work in all five areas.

The effort demonstrates the feasibility of rformini such tasks for large data

sets, such as a Digital Graphic Recorder 44)GR) data tape of roughly 500,000

points. &---

Details about the data structure used, the problems encountered, and the

necessary transformations will be discussed in the second section of the report.
In the third section, we will describe a tolerance band algorithm for thinning

the data set to a size suitable for triangulation and subsequent interpolation.

This tolerance band algorithm is a one-pass version of the method by Reumann and

Witkam (74) in analogy to techniques proposed by Williams (78,81). It was used

to create sets of sample contour points, ranging from roughly 40,000 to 70,000

points. An algorithm previously developed by NBS was then employed to determine

the Voronoi triangulations of these sets of points in a plane as illustrated in

Figure 1. This process required, however, the development of a decomposition

algorithm for data sets consisting of more than 50,000 points. This and other

necessary adaptations are discussed in the fourth section of this report.

Section Five will include a detailed dircussion of the methodology employed to

generate the surface, including the algorithm developed to produce rectangular

grids for any given unit sizes. The testing procedure, along with a description

of the structural layout used while performing the different tasks appears in

the sixth section. This section will also contain a discussion of additional

issues that are expected to be relevant to further work in this area.

Plots were obtained for the contour lines in their original form as well as in

generalized form corresponding to 40,000 and 70,000 point samples using a Gerber

Plotter. The Voronoi triangulation of a 40,000 point sample was also plotted.

The surface generation software which war prepared as part of this effort, was

applied to a a 40,000 point sample. The resulting surface was tested using the

original data, independently of the sample data set.

All computations were initially carried out on a VAX 11/750 system and later

(July 1986) transported to a VAX 11/780 system at ETL. They were found to be in

the range of extensive, but routine computations of the kind that can be

expected as part of the normal load of such systems.

Triangulation based surface modelling has been considered for some time, e.g.

Peucker and Douglas (75). For other rel ted work see Davis, Downing, and

Zoraster (82), and Grotzinger, Danielson, Caldwell, and Mandel (84).

2

I
MUrKMy ,,u v Mu wwuvvu U UU ffMUU I.MTUM Wwwwqvvxww~ UMUUI( MWM U U LW~,K-wvt- U W U7 Uv ýUV V W TKWVRWV. -WYYT~r

4,

.4.

Figue Copute-Geeratd Potio of Plt o theVornoi riaai.lti.

Figre . ortioneofea Contputr-Geinerted Plor Mutnouontain the o Tria.lo
2. DIGITAL DATA STRUCTURE

The first task performed in this effort consisted of obtaining and editing a

digital contour data set making it compatible with the software developed for

the othe'r tasks. All the necessary computations were performed originally in a

VAX 11/750 system and then transported to a VAX 11/780 with a VMS operating

system. Details of the modifications made to the particular Aata set usad

follow.

2.1 DIGITAL GPHRIC RECODER DATA

The process of digitiain$ cartographic information has been used for various

purposes over the past decades. In particular, the Digital Graphic Recorder

(DCR) has been employed to trace the lines of a map in order to generate

sequences of coordinate pairs representing these lines. This digital data set

io then stored on magnetic tape and is available in this form for computer

processing. Specifically, we are interested in the problem of determining a

computer internal terrain surface representation from this information.

An example of a data tape generated by a DGR was made available by the Defense

Mapping Agency Hydrographic/Topographic Center (DMAHTC) together with a graphic

plot derived from this information. The information had been generated from a

1:24,000 map of the Mustang Mountain area in Fort Huachuca, Arizona, -- an area

which exhibits a large range of elevations and slopes. The digitizer resolution

used, 0.0O. inches, corresponds to a horizontal distance of 20 feet.

Contour related data were extracted from this cape and the information was

reformatted and installed in a VAX system. Two basic problems were encountered

in this effort concerning the format of the data and the size of the files.

First of all, the records in a DGR magnetic tape contain 192 CDC-1700 words.

Each word contains 18 bits of whtch only the high order 16 bits are used. We

were to install this data set in a VAX with a VMS operating system. The words

4
I

in thls system are 32 bits long and the bytes are interchanged in each word. To

solve this problem a program called DGR2TAPR was written based on an example

provided by DKUTC (85). This program interchanges the high order with the low

order bytes iu each word and then it regroups the bits into 16-bit words by

extracting the high order 16 bits from each 18-bit group.

The DGR data tape used Included digitized contours, ridges, drains and

neatlines. In the original file layout (D*AWTC, 85) that follows, seguents of

such traced lines are referred to as "scan lines":

VIMST RJC0RDt HADIR

WORD CONTENTS

6-8$$$$ 1-- Sheet Number

21-30 Dat(

314 Scale (1:24000) il1-404 Series

45-50n

52-60 100M

61-68 registration Marks in Uinary (in 0.01 inches) 69-192 Coments

DATA RECORDS

WORD CONTENTS

1 Scan Line Fl

2 Scan Line Nuaamer (I.D.) 3 ypea of Data

4 X

5 Y1 (in 0.01 inches)

6 Z1

3n+l Xn

Yn

3:3 Zn

5
For each line or line segment that it digitized, the following information in

recorded: a "Scan Line Starting Flag" (octal 200000) that indicates a new

segment; a "Scan Line Number" that identifies the segment; a "Data Typa Index"

that classifies the segments as contours, neatlines, and ridges or drains; and

one or more data points at illustrated in the Data Records Layout above. YRach

record may contain one or more traced lines, or depending on the size of the

segment, it may contain only a portion of the segment. Note that a segment is

any portion of a traced line, covering possibly the whole lint, but more often

dividing the traced line into several non-contiguous segments.

The data file consisted of 9023 blocks of 192 words each. Special constraints

on the amount of disk space available and data accessibility for the type of

manipulation intended guided the design of the new file layout. It was

necessary to eliminate all the data not useful at this point of our

investigation. The first record or header record was eliminated and from the

data records only the contour line segments were kept. The resulting data

format is described below.

DATA KC*ORDS

WORD CONTENTS

1 Segment Flag 2 Se ment Number (I.D.)

3 Z-Value (Elevation of Contour Line) 4 Xl

5 Yl

2n+2 Xn

2n+3 Yn

6
2.2 IbITING THll INPUT DATA

An examination of the data showed that the coQtour lines extended beyond the

neatlines. Leaving the contours unedited eaused some undesirable triangles to

appear in the triangulation. These triangles are very long and narrow and can

be the source of large errors when a surface is generated. They are caused by

the fact that the VORONOI triangulation produces a conve. hull. A program

called CLOSCONT was written to solve this problem and to join "neighboring"

segments. This program edits the contours of the input data to include only

information that lies within a rectangular area that can be chosen arbitrarily.

The contour points falling beyond this boundary are deleted, so that the

boundary includes the last point of each contour which was intersecting it. The

resulting boundary contains points that are closer together, eliminating the

undesired triangles. For the effort reported here, the boundaries were chosen

so that the number of points to be deleted was kept to a minimum. CLOSCONT

then searchee the digital contour segments that remain for adjoining segments

(neighbors) and joins them.

The nature of the digitized data is such that all points are in a grid with a

poia separation of 0.01 inches. Thus, neighboring points in a line are

separated by either a vertical space, a horizontal space, or a diagonal, with a

maximum separation between points of 0.01414 inches. An illustration of a

digitized line appears in Figure 2. Integer numbers are used to identify each

grid point although the distance between points is given in hundredths of an

inch. Therefore, there is a finite number of points rep- esenting a line. It

was also observed thaL there are cuppressed contours in steep areas where the

contour density is high. In these steep areas a small change in the horizontal

location of a point has greater impact on its elevation computation i:nd the

resolution of the digitizer becomes more important. The effects caused by this

fact are discussed in Section 6.1. Recall that the resolution of the digitizer,

0.01 inches, corresponds to a horizontal ground distance of 20 feet.

7
""0O-4

00

%0----0

% % k0 

o --- 0

oso

0.01 ~o-0o-o - 0-• 

,o-- -o, 0---o-- -o"

,6%%

/% "o----

I

" -0 '0 -- 00 --- 0 /

Fiue2 oto n ofaDiiiedLn

8

"" - - - -o.o . .9

9-• -- 0S

0i ur.01or i n f ig ti e L n

08
3. THINNMG BY TOLERANCE BAND ALCiRITHMS

For the purpose of reducing or thinning digitized line data, the second task of

this effort, we selected a tolerance banc method. Thia approach pernits the

linking of the sampling process to the accuracy of the appr3ximation, selecting

more points in areas of high curvature than in those of lower ones. We also

calculate the contour tangents at most selected points except at the ends of

segments. A program called THIN that performs all these computations was

developed.

3.1 DEFINITIONS AND OVERVIEW

In general, the task of tolerance band methods is to select a suitable sample of

points from an initial line . The points in the selected sample are frequently

called "critical points". Connecting successive critical points by straight

lines leads to a "ieneralized line" often desired for car. Aphic purposes

(Figure 3).

The critical points are selected in such a fashion that the line points between

two successive critical points are contained in a rectangular strip, the

"tolerance band", whose bandwidth is determined from a given "tolerance" S. The

deviation of the initial line from the general line can thus be bounded in terms

of E.

Several varieties of tolerance band methods have been proposed by Lang (69),

Douglas and Poiker (73), Reumann and Witkam (74) and Williams (78,81). For

details the reader is asked to consult the comprehensive survey by Zoraster,

Davis, and Hugus (84). For this project, we designed and used a one-pass

algorithm for the method of Reumann and Witkam (74) in the spirit of Williams

(78). It will 1_ described below.

9
Figure 3. Covering of an Initial Line by Tolerance Bands .Solid Dots

Represent Critical Points. Their Heavy Straight-Line Connections

Define the Generalized Line.

10
3 e 2 A .ONE-PA3S VERSION OF THE RKUNN AND WITKAM METHOD

By

(3.2.1) Pi - (xi,Yj) ..ni~, n > !

we denote the sequence of points in the initial line. In what follows, we will

describe an algorithm th-'t iiaplements the method of Reumann and Witkam (74).

Inherent in this method is the convention that the first initial point P1 and

the last initial point Pn are automatically selected as critical points. In

addition we impoee the condition that there should be no "doubling back" by the

initial points within a tolerance band, that is, '-hat

(3.2.2) distances from the first (critical) point to subsequent points i

the tolerance band do not decrease when considered in sequence.

Furthermore we require that

(3.2.3) two critical points have at least distance E from each other,

unless they represent the first and last initial points.

With these conventions, the algorithm proceeds as follows. The point P1 is

automatically selected as the first critical point. As to subsequent points Pi,

we consider the following conditions:

(3.2.4) (i) i < n

(ii) distance (Pl, Pi) <

(iii) distance (P1 , Pi) < distance (Pl, Pi+l)

(iv) there exists a tolerance band of the second kind anchored

at P1 which contains the points Pk' 1 < k < i+l.

The point Pi is then noncritical if, in terms of the conditions (3.2.4),

(3.2.5) (i) and ((ii) or ((iii) and (iv))).

11
The algorithm checks P2, P3 ,-. successively until a first critical point Pj is

encountered. If j<n, that critical point will play the role of Pl, that is, it

will anchor the search for a third critical point, and so on. In Figure 4 we

give a pseudo-code description of the algorithm.

The major part of the algorithv concerns the verification of the tolerance band

condition. Return to the two points P1 and P2 , and assume that n>2 and that P 2

satisfies the £-separation ccndition (3.2.4.ii). Then there exist unique

left-most and right-most tolerance bands of width 2 anchored at P1 which also

contain P2 (Figure 5). The flanks of these tolerance bands are tangent to the

S-circle around P1 . In particular, the right flank of the right-most tolerance

band touches that circle at point R, whereas the left flank of the left-most

tolerance band touches at point L. We are particularly interested in the area

which lies to the right of the line from R to L, and which is bounded by rays

extending from the points R to L in the direction of the right-most and

left-most flanks, respectively (Figure 6). If the point P 3 satisfies the

no-doubling-back condition (3.2.4.iii), then the subsequent condition(3.2.4.iv)

holds clearly if and only if P3 lies in this area.

In general, suppose that Pj is the last critical point to have been determined,

and that Pip i<n, as well as the points between it and Pip are contained in at

least one suitable tolerance band. Then, there are left-most and right-most

tolerance bands -- in extreme cases they may coincide -- which contain those

points also. Suppose further that Pi+l is not closer to P than Pi is. Then,

there exists a tolerance band containing Pi+1, and all previous points back to

Pip if and only if Pi+l lies in the area Indicated in Figure 6.

In this latter case, there exist again a left-most and a right-most tolerance

band, each also covering the extended set J Pj, Pj+1,"', iD Pi+1 ." These

extremal positions are determined by the location of Pi+l with respect to the

original extreme tolerance bands. To explain the situation, it is advantageous

to introduce the notion of a "tolerance strip" associated with a tolerance band.

By this we mean the infinite straight continuation of the

12
N

tolerance band in its direction. The tolerance strip thus io bounded by the

origin-end of the tolerance band and the two infinite rays that extend the

flanks of the band. We call those rays again the "flanks" of the tolerance

strip. They are divided symetrically by the "center line". in addition, we

call the area described in Figure 6 the "area of flexibility" of Pi with respect

to the anchor point P This area of flexibility turns out to be the convex

hull of the two tolerance strips associated with the right-most and left-most

tolerance bands, respectively.

The point Pi+l now lies in (1) both tolerance strips, (2) the right-most

tolerance strip only, (3) the left-most tolerance strip only, or (4) none. Two

of these four cases are illustrated in Figures 7 and 8. In case (1), the

original extreme tolerance strips remain unchanged. However, the tolerance

bands cut from these strips are now, in general, longer since Pi+l is to lie on

their destination ends. In case (1), only the right-most tolerance strip

remains the same, as the direction of the left-most strip must now be changed to

cover Pi+l" The smallest adjustment that will achieve this is a rotation of the

left-most strip until its right flank meets Pij,. This then characterizes the

new position of the left-most tolerance strip. Case (3) is like case (2),

except that the roles of the right-most and left-most tolerance strips are

reversed. In case (4) finally, both tolerance strips need to be adjusted. The

adjustment is such that the right flank of the new left-most strip and the left

flank of the new right-most strip both intersect at Pi÷1"

After a critical point Pj, J<n, has been established, and Pj+I has at least

distince 8 from Pj, the left-most and right-most tolerance strips are first

established as if they were arising from case (4) above. This then establishes

the area of flexibility. Each time a subsequent point is found to lie in this

area, the tolerance strips are adjusted and the area of flexibility narrowed, If

a subsequent point falls outside the area of flexibility, then its predecessor

will be selected as a critical point. This also will be the case if doubling

back occurs, or if there are no points left (compare the pseudo-code description

in Figure 4).

13
define sets

LS a LEFT-MOST TOLERANCE STRIP

RS - RIGHT-MOST TOLERANCE STRIP

AF - AREA OF FLEXIBILITY

define variables

N - NUMBER OF INITIAL POINTS (INPUT)

P(I) - I-TH INITIAL POINT (INriJT)

EPS - TOLEIRANCE PARAMETER (INPUT)

H - NUMBER OF CRITICAL POINTS (OUTPUT)

if N1<2 abort

while J<N do

LS:-EMPTY; RS:-E14PTY; I:-J; CRITICAL POINT:-FALSE

while not CRITICAL POINT do

if IP(J)-P(I)I>EPS then

if P(IM LS then ADJUST LS

if P(I) R S then ADJUST RS

DETERMINE AF

if I-N or IP(J)-P(I)I>IP(J)-P(I+l)1 or P(I+1) 4Al then

CRITICAL POINT: -TRUE

K:-K+1; J:-I; C(K):-P(J)

M: -K

Figure 4. Pseudo-Code Description of the One-Pass Reumann and Witkam Algorithm.

14
Cag2

15/

IA

of

#0

P

i s

Pk r
hI

L

Figure 6. Area of Flexibility Defined by the Two Extreme Tolerance Bands

Anchored at P1 and Containing P2.

16

I
left-most

tolerance-strip

Adjusted

left-most

tolerance-strip

Original

right-most

tolerance-strip

Pk • Pi Pl.1

Figure 7. CASE (3): The Subsequent Point P t+I Lies in the Rigsht-Most

Tolerance Strip of the Initial Points JPJ, Pj+l,..', P*•

but not in the Left-Most.

17
I Adjusted

Original left-most

left-most tolerancetolerance-strip strip

Adjusted

right-most

tolerance-strip

Original

right-most

tolerance-strip

Pi Pk 

Figure 8. CASE (4): The Subsequent Point PI+I Does Not Lie in any of the Two

Extremal Tolerance Strips of the Initial Points 4Pj, PJ+1,..., Piý"

18
3.3 9DTIh = COUTOUR ?AMUT

Consider three subsequent contour points:

(xlyl), (x 2 -y 2 ), (x 3 ,y 3 ).

We define the tangent to the contour at the second point as the tangent to the

circle through these three points (Figure 9). A vector normal to the tangent at

(x2,y2) then !j given by

Y12 Y32 X12 x32 (3.3.1) )-

r 1 2 r 3 2 r 1 2 r 3 2

where

x1 2 - x1 - x 2 , x3 2 -x3 - x2 ,

Y1 2 Yl - Y2 6 Y3 2 " Y3 - Y2D

r 1 2 m x12 + Y1 2 1 r 3 2 m z32 + Y32"

This formula is used to calculate an apriori estimate of the normal to the

contour tangent at each point that lies in the interior of a contour segment.

Note that the terrain gradient is parallel to this normal. The same will hold

for the gradient of the subsequently generated synthetic surface. This will

ensure a more faithful representation of the terrain.

At the endpoints of contour lines, the tangent was left undetermined. Also if

the :hree points through which the circle is to be passed form an acute angle at

the second point, then the above determination of the contour tangent has little

meaning. In most of the later runs, the tangent was therefore left undetermined

at such points.

19
~p4~*p

Figure 9. Determining the Tangent at a Contour Point.

20
4. *ORONOIT lIANGULATIOUOF LARGE SETS OF CONTOUR POINTS

After a suitable sample of the contour points has been selected, the map area is

partitioned into triangles whose corners are sample points. We speak of a

"triangulation" of the map area. As will be described in the next section, a

synthetic surface will be generated by defining a surface patch or "element" in

each triangle of the triangulation. For this purpose, a triangulation method

that ca.)tures the proximity relationships in the sample set is desired. One

such method is the Voronoi triangulation, also called Delaunay triangulation,

which is obtained from the Voronoi diagram by dualization in a fashion described

below.

4.1 DIRFIITIONS AND ALGORITUIS

Consider a finite set S of points in the plane. For any point P in S the

"Voronoi polygon of P relative to S" is the set of all points in the plane, such

that P is as close to any point in this set as is any other point in S. The

Voronoi polygon of a point P in S is the intersection of the half-planes which

contain P and which are determined by the perpendicular bisectors of the line

segments connecting P and the other points in S. Thus, the Voronoi polygon of a

point P is a convex polygon, possibly unbounded, which contains P in its

interior. Given an arbitrary point X in the plane and the Voronoi polygons

associated with a set S, then one and only one of the following statements is

true:

(1) X lies in the interior of one and only one Voronoi polygon.

(2) X lies in the interior of an edge shared by two Voronoi polygons.

(3) X is a vertex of three or more Voronoi polygons.

It follows that the Vorcn,)i polygons of the points in S cover the plane without

overlapping, that is, without common interior points. The union of their edges

forms a diagram, the "Voronoi diagram for S", which partitions the plane into

21
the Voronoi polygons. A point that lies in the plane and satisfies condition

(3) above is said to be a "vertex" of the Voronoi diagram for S. It is called a

"degenerate vertex" whenever it is a vertex.of more than three Voronoi polygons.

Figure 10 illustrates a Voronoi diagram that has a degenerate vertex.

The "dual Voronol diagram" for a finite set S of points in the plane is the

diagram obtained by connecting with straight-line segments those pairs of points

in S whose Voronoi polygons relative to S have an edge in common. Figure I1

shows how such a dual diagram is obtained from a Voronoi diagram. The dual

diagram defines a collection of non-overlapping convex polygons which cover the

cornvex hull of S. Since each edge of the Voronoi diagram is a line segment

whose end-points are vertices of the Voronoi diagram, it follows that there is a

one-to-one correspondence between the vertices of the Voronoi diagram and the

polygons determined by the dual diagram. In fact, it follows from the

definition of a Voronoi polygon that a vertex of the Voronoi diagram is

equidistant from the points in S that are vertices of the corresponding polygon

determired by the dual diagram, and it is closer to these points than it is to

any other point in S. In general, most vertices of the Voronoi diagram are

non-degenerate, so that most of the polygons defined by the dual diagram are

triangles. In the presence cf degenerate vertices, the corresponding polygons

determined by the dual diagram can be partitioned into non-overlapping triangles

by introducing suitable diagonals. Thus a "Voronoi triangulation" results from

the dual diagram for tthe set S. We note that, while the dual Voronoi diagram is

uniquely determined, there are several compatible Voronoi triangulations in the

presence of degeneracies. Figures 12 and 13, respectively, illustrate the dual

diagram and a Voronoi triangulation that results from it.

J. Bernal and S. E. Howe of the National Bureau of Standards (NBS) have

generalized, extended and combined algorithms by Bentley, Weide and Yao (80),

and Bou-rer (81), to obtain an algorithm which constructs the Voronoi diagram

and, therefore, a Voronoi triangulation, in "linear expected time" for a set S

of points distributed uniformly in the interior of a rectangle in the plane.

This material is being readied for publication under the title "Expected 0(N)

and 0(N4/3) Algorithms for Constructing Voronoi Diagrams in Two and Three

22

~ WAWA K L~A bJ ~PJA~ JL~JX~~1F I
0 0

0 

0

0

0 0

0

0 0

0

0 

0

0

0

Figure 10. Voronoi Diagram of Point Set Indicated by Small Circles.

23

---- ------------
Figure Voronol 11. Diagram SppentdbisDul

qb go 2

.24
Figure 12. The Dual Voronoi Diagram or Delaunay Diagram Representing Neighbor

Relations Between Points.

25
Figure 13. Voronoi Triangulation Obtained by Adding a Diagonal to

a Non-Triangular Cell in the Dual Diagram in Figure 12.

26
Dimensions" by J. Bernal and S. E. Howe. A draft copy is available upon

request.

In the resulting algorithm, the set S which is to be triangulated is enclosed

into a rectangle. The algorithm consists of three steps. The first step

divides the rectangle into approximately N (number of points in S) equally sized

square cells, and assigns each point in S its proper cell. On the av~erage there

will be one point per cell, although there may be empty cells as well as cells

with more than one point. Each cell is of the form

4(x, y): x < x < x+b, y < y < y-Ib]-,

for some x, y, b, and a point in S is assigned that cell if it lies in that

area. Cells within a distance of two cells from the boundary of the rectangle

are called "~outer cells" and all others, "inner cells."

Thke second step c~onstructs the Voronoi polygons of the points belonging to inner

cells. Given a point P in an inner cell, a search for other points in S is

conducted through each of the layers of cells surrounding P. This search

procedure, called a "spiral search", starts with the -cell that contains P, and

then proceeds in outward direction to each of the layers of cells surrounding

this cell. The Voronoi polygon of P is progressively built by intersecting the

half planes which contain P and which are determined by the perpendicular

bisectors of the line segments connecting P and the points in S found through

the search. A geometrical test is available, which permits to ascertain whether

the Voronoi polygon has achieved its final form. In most cases the Voronoi

polygon of P is obtained after examining only a small number of cells and

points.

The third step, finally, builds the Voronoi polygons of points in the outer

cells by applying a modified version of Bowyer's insertion algorithm to this set

vi points.

27

I. NU
4.2 INPLEMMUATION

J. Bernal and S. E. Hove have implemented the above algorithm on a Control Data

Cyber 205 at NBS. The im'plementation consists of about 5,000 FORTRAN

statements. It requires as input a tolerance 6, and a list of the x and y

coordinates of the points in the set for which a Voronoi diagram is desired.

This list of points must be free of duplication, that is, the distance between

any two points must always be above the tolerance E. The execution of the

package requires approximately 34N words of memory, where N is the number of

points to be triangulated. Versions of this package were successfully transferred first to a VAX 11/750 and then to a VAX 11/780 at ETm. This required

adaptations which are described below.

One of the main objectives of our work was to demonstrate the feasibility of

triangulating data sets of about 40,000 to 70,000 points. Two difficulties

arose in the course of the demonstration. The first difficulty was that of

furnishing the package with the ability to deal with data sets whose members are

not all necessarily distinct. The second difficulty had to do with memory

restrictions that would not allow the execution of the package for more than

50,000 points at a time.

The "check for duplication" posed a difficulty because it was not possible to

consider every pair of points in view of their large number. In the case of a

40,000 point set this would have amounted to the examination of about

800,000,000 pairs. Fortunately, it was discovered that a portion of the Voronoi

j triangulation package already provided a useful tool for the solution of this

problem, namely the cell structure set up by the implementation of the first

step of the algorithm. Thua, a procedure was developed which takes advantage of i

this structure to check for the duplication of points. Given a point P in a set

S, use a spiral search through each of the layers of cells surrounding P to

search for other points in S. Eliminate any point found through the search

whoae distance from P does not exceed the tolerance E.

28
Terminate the search as soon as every cell that intersects the closed circular

disk with center P and radius has been searched. In general, the disk defined

above is contained in the cell that contains P, so that this is the only cell

that has to be searched. This procedure was incorporated into the package, and

was able to identify and remove duplications for data sets of 25,0n0 points in

less than 3X seconds of CPU time.

The restrictions on available memory mentioned above required the development of

a "decomposition procedure" which allows the separate triangulation of a finite

number of subsets of a data set in such u way that the correct total triangulation results. This procedure will now be described.

Given a positive integer k we belect numbers x0 , xk, yo, Yls such that x0 < xk,

yo < yl, and the rectangle

R - {(x, y): xo < x < xk, yo y Y ,

contains the data set. We select numbers

xoL, X1L, xI, XlR, x2L, x2, X2R,..," xk-1,L, xk-1,- Yk-1,R, XkR,

such that

Xo<Xl<... <Xk-4<Xk and

Xo'xOL, XlL<Xl<XlR, x2L<x2<x2R,..., xk..k-l,L<Xk-l<xk-l,R, xkR-xk"

For each i, 1-1, 2,..., k, we define rectangles Ri' and Ri by:

Ri' yi(x, Y): xi_,L xjSxi,R, YO Y<y y1 k,

Ri - 4(x, y): xil < x < xi, Y0 < y < yj"}

k k

It follows that Ri C Ri' for each i, i-, 2,..., k, and R Ui Ri i•U 1Ri'.

29
We also assume that the points (xO, y0)9 (xO, yl), (x 1 , yO), (xl, yl),...,

(xk_l, yo), (xk_l, yl), (xk, yo), (xk, yl), that Is, the corner points of the

rectangles Ri igi,..., k, belong to the data set (Figure 14). Furthermore,

we define Si to be the set of points in the data set that belong to Ri' for

each i, i-l..., k. Assuming that k and the numbers x0L, XlL x1l, XlR,...,

1k-iL, xk~i, 1k-i,., xkR have been properly selected, the decomposition

procedure consists of obtaining separate triangulations Ti, 1.l,..., k, for the

sets Si, ig1,..., k, respectively (Figure 15). The correct triangulation of

the entire data set is then given by:

k

U~ 4 t 6 Ti : Ri intersects the interior of t '.

In order to properly select k and the numbers xOL, xlL, xl, XlR ... , Xk-l,L,

xk-_i xk-_,R, xkR, a separate procedure was developed. In what follows, we

define, for a given triangle t in a Voronol triangulation, x(t) and y(t) to be

the x,yocoordinates of the vertex in the Voronoi diagram that corresponds to t.

Accordingly, we define d(t) to be the distance from (x(t), y(t)) to any one of

the vertices of t. Since (x(t), y(t)) is equidistant from the vertices of t,

d(t) is well defined. In the following procedure, m denotes the maximum number

u points that can be triangulated with a single run of the package:

Step 1. Let k - 1 and obtain Rl, Rl' and S1. Let j - 1.

Step 2. If the number of points in Sj does not exceed a go to step 3.

Else increase k to the next positive integer for which xOL, XlL, xl,

" lR,..., xk-l,L, xk-l, xk-1,R, xkR, can be defined with: xO<xl<...<xk_l<xk;

xlL<xl<xlR•., xk-l,L<xk.l<xk-i,R; XO=xOL; xkR=xk; and the corresponding

Ri, Rif, Sj ... , k, can be obtained with the number of points in each

Si i ,..., ..iot exceeding m. Let j - 1.

Step 3. Obtain the Voronol triangulation T for S If j is equal to 1

let XL=XOL" Else define xL by:

xL - min x(t) - d(t): t * Tj, 2 vertices of t lie in Rj_ 1 ,

1 vertex of t lies in the interior of Rj J.

30

Il
x %00

C~C4 N

Kx

x %00

600

K x

Figure 14. Decomposition of a Rectangular Point Set for Voronol Triangulation.

31
(X0,YJ IStl) 42"41 

(xL~yO) X-O I

Figure 15. Voronoi Triangulations of Three Sections Superimposed Over the

Voronoi Triangulation of the Entire Set.

32j
If j is equal to k, let xR-xkR. Eloe define xR by:

xKR - max x(t) + d(t): t * Tj, 2 vertices of t lie in RJ+l,

1 vertex of t lies in the interior of Rj J.

If ljl$L does not exceed xL and xR does not exceed x ij go to step 4.

Else,

if xjil,L exceeds xL let xjlL-XL, and if xR exceeds xjR let xjR-xR.

Obtain R , Rj', Sj and go to step 2.

Step 4. If j equals k, stop. Else let j - j + l and go to step 3.

This procedure was incorporated into the triangulation package and is currently

operational.

4.3 RESULTS OF A CO1PUTATIONAL WKMPRIMUT

A Voronoi triangulation package that includes the adaptations described above

was implemented at ETL and NBS. The feasibility of performing a triangulation

with large sets of data was demonstrated. For a specific experiment conducted

on the VAX 11/750, we divided the triangles into three classes:

Class 1: Those whose vertices lie exactly on one contour line of the map.

Class 2: Those whose vertices lie exactly on two contour lines of the map.

Class 3: And those whose vertices lie exactly on three contour lines.

Each class may contain triangles that are intersected by contour lines other

than those containing their vertices, and any given triangle belongs to one and

only one class.-

A first data set contained 39,645 points and was decomposed into two subsets. A

second data set contained 70,249 points and was decomposed into four subsets.

The tables below illustrate some of the results obtained when triangulating

these data sets with the Voronoi triangulation package. The times, given in CPU

seconds, indicate the rate at which the package ran per point.

33
TAME FM 39,645

Element of I Class 1 Z Class 2 2 Ciao@ 3 Number of CPU Time

Decomposition Triangles Triangles Triangles Triangles Sec/Point

1 17.20 49.96 32.84 41,875 0.044175

2 17.05 48.09 34.86 36,451 0.046810

TANBL 0 1 " 70,249

Element of % Class 1 2 Class 2 2 Class 3 Number of CPU Time

Decomposition Triangles Triangles Triangles Triangles Sec/Point

1 22.58 59.58 17.84 41,134 0.048542

2 22.65 58.65 18.70 32,824 0.043350

3 18.74 59.53 21.73 33,715 0.044793

4 26.00 55.76 18.24 31,755 0.0521S9

34
S - SURFACR GIUtRKLON WITE CLOUGW-TOMcUI M WS

In this section we discuss the construction of a surface function passing

through irregularly spaced given points for most of which both elevations and

contour tangents are specified. The planar projections of all given points are

triangulated, that is, the map area is tiled with triangles whose vertices are

these points. The description of the surface is in terms of these triangles: in

order to find the surface elevation for an arbitrary point in the map region, a

triangle containing this point must be found. An evaluation formula or

"element" is then invoked using triangle-specific parameters. The "Clough-Tocher

element" employed in our work requires the determination of suitable tangent

planes at the given points. For this purpose, "local" as well as more expensive

"global" methods are available. A global sethod based on energy minimization

"has been implemented and tested for computational feasibility.

5.1 THE CLOUGR-TOCHER ELEMENT

Triangulation-based surface interpolation is a classical computational problem.

Various versions of the Finite Element Method (see Zienkiewicz (71), Birkhoff

and Mansfield (74)) are usually employed in its solution. The "linear element"

represents linear interpolation by the plane through the vertices of the

triangle at hand. It yields a surface of continuous elevation. However, this

surface is not smooth since "creases", that is, tangential discontinuities,

occur along the boundaries of the triangles.

Nonlinear elements are needed for smooth surface interpolation. A major

advantage of smooth surfaces is that their corresponding surface functions are

uniquely differentiable at each point of their domain; in other words, there are

unique gradients. The greater flexibility and information content of nonlinear

elements also allows for a more precise representation of the original data than

that provided by linear elements, given triangulations of comparable densities.

35
- al , - - . w -' -,Lrt rff l f•u•lS -mn • "

The "Clough-Tocher" element (see Clough and Tocher (65), Lawson (72, 76, 77)) is

a particularly attractive tool for smooth surface interpolation. It requires

that at each vertex i of the triangle over which it is defined, the elevition ti

and the partial derivatives zix and ziy be given. The Clough-Tocher element

then is described by a function z - f(xy) on the given planar triangle. It

represents a "surface patch" above this triangle (Figure 17, see also Figure

16). The surface patch meets the prescribed elevation as well as the prescribed

derivatives at each vertex. It is fully defined by these quantities, in other

words, by the three elevations and the three tangential planes (gradients).

The following considerations concern the construction of an entire surface from

such surface patches. In order to ensure continuity of elevation between

adjacent triangles, the Clough-Tocher element satisfies the following

(5.1.1) Cubic Boundary Condition: Along each triangle edge, the

Clough-Tocher element agrees with a cubic (degree 3 or less)

polynomial in terms of a variable linearly traversing the edge.

A cubic polynomial in one variable is completely determined by two elevations

and two derivatives (in the direction of the edge). Since the tangential planes

at the vertices agree, so do the derivatives in the direction of the edge.

Therefore, the Clough-Tocher elements of two adjacent triangles determine the

same cubic polynomial on the edge they share. It follows that not only the

elevations, but also the derivatives in the direction of the edge agree along

that common boundary.

In order to ensure smoothness across the triangle boundaries, the Clough-Tocher

element satisfies the following

(5.1.2) Linear Derivative Condition: Along each triangle edge, the A

derivative taken in the direction perpendicular to the edge

varies linearly between the values it assumes at the ends of the

edge.

36

N-
(X 2,Y2)

Figure 16. Linear Element.

37
f:Z2xv Z2y

Zig, Zly :::.

Zi

Figure 17. Clough-Tocher Element.

38
The edge-perpendicular derivatives along an edge are thus uniquely determined by

their location on the edge and the derivatives at the vertices at the ends of

the edge. Since the latter agree for two adjacent triangles, the edgeperpendicular derivatives along their common edge agree also. It was seen that,

as a consequence of the cubic boundary condition, the edge-parallel derivatives

agree. Hence, adjacent Clough-Tocher elements share tangential planes

everywhere along their common boundary.

Functions over triangles are best expressed in terms of their "barycentric

coordinates", also called "triangle coordinates." These are three real numbers

).Il x2, X31

such that

(5.1.3) X1 + X2 + X3 = 1

X1x + X2 x2 + X3 x3 = x

XIyI + X2 Y2 + X3 Y3 - y,

where x, y are the planar coordinates of the point in question. The barycentric

coordinates are functions of these planar coordinates. To express these

functional relationships, we use Zienkiewicz notation!

xij xi - xi, Yij = Yi - Yj 1 i,j - 1, 2, 3.

We then have

(5.1.4) ),l X I(xy) . (Y 2 3 .x + x3 2 -Y + x2Y3 - y 2 x 3 )/D

X2 X2 (x,y) = (y 3 1 .x + xl3.y + x3 y1 - y3xl )/D

X3 = X3 (x,y) ' (Y12,x + x 2 l.y + xlY2 - YlX2)/D,

39
where the denominator is given by the determinant

1 1 1

D = x, x 2 x3

Y1 Y2 Y3

We also note that the partial derivatives of the barycentric coordinates with

zespect to x, y, are given by

(5.1.5) klx m Y2 3 /D, k2x = Y3 1 /D, k3x -Y2/D,

Xly = x3 2 /D, '2 y = xl 3 /D, X3y - x2 1 /D.

The advantage of barycentric coordinates lies in the symmetric way the vertices

of the triangle are treated. Aleo, their signs indicate immediately whether the

point (xy) lies inside or outside the triangle: a negative barycentric

coordinate indicates that the point lies outside the triangle. If all

barycentric coordinates are positive, then the point lies in the interior of the

triangle. On the boundary at least one barycentric coordinate vanishes.

Vertices are characterized by single nonzeru barycentric coordinates of value 1:

(1, 0, 0) , (0, 1, 0) , (0, 0, 1).

The barycenter of the triangle or "centroid" is given by

( 1/3, 1/3, 1/3 ).

It defines what we call a "barycentric partition" (Figure 18) of the triangle

(5.1.6) T - B1 u B2 u B3 u B.,

where

M IlI
CI

Figure 18. Barycentric Partition of Triangle.

41

3 LAw 3B
B, - 40L1,k 2,x 3 ): kl <')2, kl < k3 1

B2 - i(kLk2,X3): k2 < k3, k2 kl k

B3 - 4(XL,)i 2,)k3 ): )3 < 'L-, )3 L2•

BO - {( 1 ,X2, ) 3 ): 61 - 2 X 3 1 l/31.

In each of the major triangle regions Bi, i - 1,2,3, the corresponding

barycentric coordinate is dominated by the remaining ones:

ki -min4 Xl, k2, X3 I.

Each function representing a Clough-Tocher element is a cubic polynomial of the

barycentric coordinates in each of the major regions Bi of the barycentric

partition of the triangle. The function is continous and smooth at the

boundaries of these regions of the barycenter. We call such a functton

"piecewise cubic" with respect to the barycentric partition.

In his seminal work, Lawson (76) introduces three "correction functions"

(5.1.8) Pi - Yi("1,) 2 ,) 3 ), i - 1,2,3,

with which to describe the Clough-Tocher element. They are piecewise cubic with

respect to the barycentric partition, and given by

klk2k3 + 5/6)L13 - 1/2k12 for (k1,k2,3) E- B1

-1/6X2

3 + 1/2k22 >-3 for ( 1 ,k2 ,X3 ) E B2

.91~

-1/6)33 + 1/2X3

2X2 for (Xl,Xv2,3) 6 B3

1/81 for (Xi,)2,)3) * B0

42
11-)3+ 5/6kL23 - 1/2)L22 for (X1)L ,)2 3)1k G-B

1/633+ 1/2k3 

2 x for (ki % 2,L3) B33

92

-16,+ 1/2)1, k3 for (01 ,X2,X3) G B1

1/81 for (X1,XL2,k3) 6 O

)L)23+ 5/6k~3 3- 1/2)13 2 for ()1,12,k3) B33

-16L3+ 1/2k L2)L2) for (k1,X2 ,)x3) 4 B,

93-

-1/6)L2 3 + 1/2k2 

2)x for ()Ll,k 2 ,X3 ) 6 B2

11/81 for (X1,X2, 3) S BO.

The motivation for the choice of the correction functions is given in Lawson

(76). The partial derivatives of these correction functions are:

(5.1.9)

(xX35/)L2_Lly2 + 3X)y 3 l + lkl2y12 in Bi

(1222x23yl+ 1/2Xh22Y12 in B

Dpx-

(-/X k)2Y2+ 1/2X~3 .Y31 in B3

1l/18Y23 in B

Nk+/),2-lx2+ Xlx3+ XXxlin B

(-1/2k2 

2+k2)L3)xl3 + 1/2X2 

2x2 l in B2

D y " (-1/2 )`3 

2+X3XL2 )x 2 l + 1/2X13 2X13 in B3

-1/ 18x32 in B

43
r (>3kl+5/2L2 2-X2)y3l + Xl~l 2k'3y23 in B2

(-/k +3ly2+ 1/2)L3

2 Y23 in B3

D-?

(-1/2k,2 +klk3)y23 + l/2kj2 Y12 in B

-1 /18Y31 in B

kxl526 _2l3+ XlLxl+ k2Lx in B

(-/k +)lxl+ 1/2kb3 

2x32 in B3

D - 2 -/ k 2 k ) 3 x 2 + 1 / 2 k l 2 x 2 , i n B

-1/18X13 in B0

(xx+/23213Y2+ "2k3y23 + '-3Xl1y3l in B3

(-/)l+1)2y3+ 1/2)12Y31 in B1

(-/)22+2lyl+ 1/2XL2

2 Y23 in B2

.-1 /8y12 in Bo

r(X 1)k2+5/2XL3 

2- 3)x2 l + +~~3 XL3)'lx13 in B3

(-'~2+ 1.2)x32 + 1/) 1 x13 in B

(122kx + 1/2k2 

2x32 i B2

in Bo.

44
Now let

(5.1.10) Mji - zixXji + Ziyyji, ij - 1, 2, 3, iIj.

These six quantities represent the directional derivatives at vertex i with

respect to the direction vector (xji,Yji), which represents a directed edge of

the triangle. It will be convenient to use the additional abreviations

(5.1.11) Qij - l/2(Mji + Mij), Cij - l/ 2 (Mji - Mij) - Zji"

Note that for a linear function z - f(x,y), Mji - -M j - zji, and for a

quadratic function, M t - Mij - 2 zji. As a consequence, the coefficients Qij

and Cii vanish for linear functions, and the coefficients Cij for quadratic

functions. We also denote by

(5.1.12) Li, i - 1, 2, 3,

the euclidean length of the edge opposite to vertex i. The Clough-Tocher

element is now of the form:

(5.1.13) z - Xlz1 + X2z2 + X3z3

+ Q2 3X2 X3 + Q31X3'l + Q12Xl"2

+ C2 3Vl + C3 1V2 + Cl2V31

where the functions Vi, i - 1,2,3, are given by

V1 " 1 2 x3 (X2-X 3 ) + [3(L 2+L3 )(L 2-L 3 )/Ll 2 ]?1 P2 +

V2 )X 3XI(X 3 -k 1 ) + [3(L 3+L1 )(L 3 -L 1 )/L 2

2 ]1?2 ?3 + ?l

V3 - XlX2(kI-X 2 ) + [3(Ll+L2 )(L 1 -L 2 )/L 3

2 ]? 3 - Yl + ?2"

Note that the expression

xlZ1 + k2z2 + X3z3

45
represents the linear elemert, namely the plane passing through the elevations

given at the vertices of the triangle. The portion

Q23ý-213 + Q3)

1 31 1 + Q1

is a quadratic function in the entire triangle. Note that quadratic functions

satisfy the Linear Derivative Condition (5.1.1). The remaining portion is a

piecewise cubic function with respect to the barycentric partition. This

function also satisfies the Linear Derivative Condition as a result of the

choice of correction functions Pi"

It is easy to derive various expressions for the gradient of the Clough-Tocher

function. The gradient components are the partial derivatives with respect to

x, y, and may be written as

(5.1.14) D.zx M zlY2 3 + z 2y31 + z3Y12

+ Q23(X 3y31+X2 y 1 2 ) + Q3 1 (X 1yl 2+X3y 2 3 ) + Q12(L2Y23+'lY31)

+ C2 3D-Vlx + C3 1D.V 2 x + C12D.V3x

D-zy - z 1 x3 2 + z 2 x1 3 + z3x21

+ Q2 3 (0 3 x1 3 +k2x 21 ) + Q3 1 (klx21 +)h3x3 2 ) + Q1 2 (O2x 3 2+)lX1 3 )

+ C2 3D-Vly + C3 1 D-V 2 y + C1 2D.V 3 y,

where

DnVlx W X3 (2X 2 -k 3 )Y3 1 - X2(2k3-`2)Y12

+ [3(L 2 +L3 )(L 2 -L 3 )iL 1

2 ]D.? 1 x - V. 2 + 3x

D.V 2 x - X1 (2)L3 -) 1 )y 1 2 - A3(2),--3)Y23

+ [3(L 3+Ll)(L3 -L1)/L 2 2]D.?2x - D.p 3 x + V-?jx

V.V3x - k2( 2 X1 -x 2 )y 2 3 - kl(2k2-ý-)Y31

+ [3(Ll+L2 )(Ll-L 2 )/L 3

2 ]D.P3x - D.lxp +

46
DVy- X3(2)2-)3)xl3 - 221kxl

+ [3(L 2+L3)(L2-L3 )/Ll2 ]D.pl7 - D. ?2y + D ~

DVy- )L1 (2) 3 -) 1 )x 2 l - 32lk)2

+ [3(L 3+Ll)(L 3-Ll)/L2 

2 ]DiP 2y - D*9y+ ol

D.V3 y - )h2 (D) 1-XL2 )x 32 - l22Xx3

+ (3(L1+L2)(L1-L2 )/L3 

2]D*93 - Degly + D.92y.

Amore transparent formula for the gradient can be derived. For )kj - 1,

i-1,2,3, that is, at the three corners of the triangle, 9i - Piy " 0. It thus

follows from (5.1.14) that

Dezlx m 12 + zY1+ zl2+ (Q12+'Cl2 )y3 l + (Q3l-C3 l)yl2

D -~ z1y23 + z31+ zy2+ (Q23+C23)y1 2 + (l-l)2

D -~ ~2 + zY1+ z12+ (Q31+C3l)y23 + (2-2)~

Dozl, zlx 32 + zx3+ zxl+ (Ql2+'Cl2)xl3 + (Q3 1 -C3 1 )x2 l

D~z2, z1x32 + z"3+ zxl+ (Q23+C23)'21 + (l-l)3

D -~ l3 + zx3+ zxl+ (Q3 1 +'C3 1 )x3 2 + Q3C3x ,

where zix ziy, 1-1,203, are the prescribed gradient components at the vertices

of the triangle. It follows that

kiD.zix1+ 2D'z2x'*N3D z3x zly 23+z2y3 l+z3yl2

-C23 (X3y3l-x2yl2 )-C31 (X'lyl 2 X)3y2 3)-Cl 2(X2y23)1ly3 l)

47
XlD.5ly+k2D.5 2 y+X3D.e3y - glX32+22x13+z3x2l

+ Q23(-3x13+k2x21)+Q31()ý1X21+k3x32)+Q12()L2x32+klXl3)

- C2 3 (k 3x 1 3-) 2 x21 )-C 3 1 (OlX21 -) 3x3 2 )-Cj 2 (X2x 3 2 -)lx 1 3 ).

Formula (5.1.14) can nov be rewritten as

(5.1.15) z llx + k2z2x + "3 '3x + C2 3Alx + C3 1A2 x + C1 2A3x

sy M kXlly + k2'2y + k 323y + C2 3 Aly + C3 1A2 y + C1 2 A3 y,

where

Aix a V1X + (1 3Y3 1- 2Y12 )/D, Aly M V ly + (O3x1 3 -) 2 x21 )/D

A2x M V2x + ()lY1 2 -• 3 y 2 3 )/D, A2 y - V2 y + ('lx2l-"3x32)/D

A3U - V3U + (X2y 2 3-)ly 31 )/D, A3y - V3y + (1 2 x3 2- 1xl 3 )/D.

To sum up, given three points with their elevations and gradients, that is,

given 15 quantities

xi, Yi, zi, zix, Ziy, i - 1, 2, 3,

the Clough-Tocher element defines a surface that assumes those prescribed

elevations and derivatives (gradients). In order to calculate the elevation at

an arbitrary specified point, we

"o compute the auxiliary quantities Qij' Cij and Li from the above

15 quantities

"o evaluate the correction functions 91 for the barycentric coordinates I

kI' X21 X3 with respect to the three given points in the plane

" determine the values Vi and enter them into the Clough-Tocher formula

(5.1.13). The gradient components are computed analogously using

formula (5.1.15).

48I
5.*2 SURFACE GENERATION OVR A TRIJANGULATED REGION

In order to pass a smooth surface through given elevations at irregularly spaced

points, one may partition the region in which the surface is to be defined into

triangles, specify gradients along with the evevations at the vertices of these

triangles, and generate the resulting Clough-Tocher patch in every triangle. The

Clough-Tocher patches are designed in such a fashion that they fit together

smoothly along common boundary edges. Ho"iever, there is still a missing link.

In order to fully define our synthetic surface, we have yet to provide the

gradients (tangential planes) at the vertices of the triangulation. There are

several approaches to finding suitable gradient values. One is to examine

neighboring elevations and to estimate a suitable position of the tangential

plane at the point in question by using local interpolation or least squares

regression. The success of spline techniques suggests a different approach. A

spline can be interpreted as an idealized mechanical structure consisting of

"1.thin beams" which, when forced through specified points, assume a position in

which a surrogate elastic energy is minimized. Since oscillatory behavior is

associated with high elasti, energy, minimizing elastic energy tends to minimize

oscillations. In this work, we extend this approach to two dimensions. We

consider a mechnical structure consisting of thin beams of equal "thickness"

along the edges. They are joined together at vertices by small "thin plates"

which represent tangential planes forced to be met by the adjacent thin beams.

The sum of the surrogate energies of all thin beams is then minimized. This is

achieved by varying the positions of the thin plates to which the adjacent thin

beams must be tangential while passing through prescribed elevations and

satisfying other side conditions that may have been specified for selected

vertices. Once the idealized mechanical structure has found its optimal, that

is, energy-minimal position, the gradient at each vertex is defined by the tilt

of its thin plate. The resulting surface is supported by thin beams much as the

fabric of an umbrella is spanned by its ribs.

Mathematically, the surrogate elastic energy is a positive definite quadratic

form in those parameters that are permitted to vary. Setting the partial

49
derivatives of the energy with respect to these variables to zero yields the

optimality conditions in the form of a large system of equations.

The classic Gauss-Seidel method for solving a system of linear equations is an

iterative procedure where each iteration consists of passing in sequence through

all the equations of the system. Each equation is used for determining a new

value for a particular unknown variable while keeping the other unknowns fixed.

The method thus requires an initial value for each unknown. To start, the first

equation is transformed into a linear equation of only the first unknown by

substituting into this equation the initial values of all remaining variables.

This equation then yields an improved value for the first unknown variable. This

value, along with initial values for the third and subsequent variables, enters

the second equation, which then yields an equation for the second variable

alone, and so on.

We modify che Gauss-Seidel procedure slightly. To this end we observe that each

variable of the system of linear equations is "associated" with a particular

vertex. The associated variables at each vertex satisfy "local optimality

conditions". These optimality conditions have a structural interpretation: They

express the conditions for the structural parameters represented by the

associated variables to assume minimum energy values, supposing that the

structural parameters at all other vertices remain fixed. The local optimality

conditions again take the form of linear equations in the variables associated

~with the vertex. In a typical case, the thin plate at the vertex is tilted into

the best position it can assume, given the tilts and elevations at neighboring

vertices. The local optimality conditions then define this locally optimal

tilt. It can be shown that:

(5.2.1) THEOREM: The linear system of equations for optimizing the

position of the idealized mechanical structure is equivalent to

the combination of all local optimality conditions.

Our variation of the Gauss-Seidel method now is to pass through all vertices in

sequeftce, solving at each vertex the local optimality conditions. It was a

50
major concern at the start of this project whether this procedure was

computationally feasible. We found. that the time needed for a Gauss-Seidel

iteration of the kind described above was well within an acceptable time frame.

As discussed in Section 6, this marks a main accomplishment of this feasibility

study.

In what follows, we will list the local optimality conditions for various sets

of associated variables. We distinguish several "types ol vertices" according

to their kinds of associated variables. A particularly important case is the

one in which the tangent to the contour curve is given along with its elevation.

In this case, the direction of the gradient is given -- it is perpendicular to

the contour tangent -- and the only variable to be determined is the length or,

rather, a positive or negative "gradient multiplier" 1. It represents the only

variable associated with a vertex of this type.

We use a three letter code to characterize vertex types. The first letter of

the code refers to elevation: it is 'E' if the elevation is given, and IN'

otherwise. In the latter case, the elevation is an associated variable. The

second and third letters refer to gradient components zx and zy, respectively.

Letter 'X' in the second position indicates that the x-component of the gradient

is given, and 'N' in this position indicates that it is not. Analogously, letter

IV' in the third position indicates that the y-component of the gradient is

given, and IN' indicates that it is not. Finally, the combination 'RN' is found

in positions two and three, if the gradient direction, but not the gradient

itself, is specified. For example, 'ERN' signals the case in which both

elevation and gradient direction are prescribed, leaving the gradient mul~tiplier

as the only associated variable. Any occurrence of 'N' indicates an unknown

associated variable. In particular 'NNN' is used if all three quantities,

elevation and gradient components, are to be determined.

The surrogate energy of a thin beam of length L is given by

sinL

(5.2.2) E -J z"(s) 2ds.

51
It is a well-known result of the theory of thin beams that the above expression

is minimized by cubic functions of the distance s along the projection of the

beams into the plane. These cubic functions are uniquely determined by the

elevations and slopes at the end points ("Hermite interpolation"). The

following local optimality conditions are derived using these facts.

In describing local optimality conditions at vertex i, we let the vertices j run

through the "star" of vertex i, namely the following set of vertices J:

(5.2.3) S(i) - { j : j is connected by an edge to vertex i].

The abbreviations (5.1.10) and

(5.2.4) Lji - xji2 + Yi2

are used to denote the edge-directional derivatives alnd the distances between

vertices, respectively. In addition, we will need the quantities

(5.2.5) Pji = Pixxji + piyYji'

where

Pix' Piy' Pix2 + Piy2 1

are the components of the gradient direction at vertex i normalized to length 1.

The resulting optimality conditions for all types of vertices, except type

'EXY', are displayed in Figure 19. 'EXY' represents the fully specified case in

which there are no associated variables to be determined.

To sum up, the vertices of the triangulation are divided prior to surface

generation into types depending on which surface parameters are given and which

have to be determined. These types are described by the letter codes

(5.2.6) 'ERN', 'ENN', 'EXN', 'ENY', 'EXY',

'NRN', 'NNN', 'NXN', 'NNY', 'NXY'.

52
For each of these ten types -- except the fully specified case 'EXY' -- there is

an update formula by which the surface parameters at the triangulation vertices

are calculated from the current values of the surface parameters at neighboring

vertices. Sequentially updating every vertex in this fashion constitutes a

Gauss-Seidel iteration. Such iterations are repeated until the changes in the

variables remain within a given tolerance or a specified limit on the number of

iterations is reached.

5.*3 FINDING THE RIGHT TRIANGLE

Suppose a surface is specified in terms of triangular Clough-Tocher patches. In

order to evaluate the elevation at a given point in that surface, a triangle

containing that point needs to be found. In what follows we describe a method

for finding such a triangle. This method is intended for applications in which

not just one point, but sequences of such points are given, most of which are

moreover close to each other. The sets of sequential contour points as they

arise from digitized contour information are a case in point. A closely spaced

regular grid is another. In these cases, a given point in the sequence will lie

with high probability in the same triangle or in a triangle directly adjacent to

the triangle of the previous point.

The method we use for finding a triangle containing a given point relies for its

efficiency on the above observation. Before searching for the triaragle of a

given point, the method requires that an arbitrary st~arting point be specified

for which a triangle containing it is known. We then move from this starting I point straight towards the given point until we reach the boundary of the

starting triangle or the given point itself, whichever happens first. If the

boundary is reached at a non-vertex point, that is, somewhere in the interior of I. a boundary edge, then the unique adjacent triangle can be readily identified

using our triangulation data structure, and we continue moving in that triangle

as far as possible or necessary. In the unlikely case that a vertex is

encountered, all triangles adjacent to this vertex are examined in sequence

until one is found in which progress can be made towards the given point. The

processes are repeated until the given point is reached. The last triangle in

53
Cue ERN." elevation xj and gradient direction (Ai., A,1 ) specified, gradient

multiplier "yj to be determined.

y) x 2[3z., + Mij

Case EXN: elevation zi and gradient component Aj, specified, gradient component z3⁄4, to be determined.

4. x 2E (. E (,) - 2;,,,. + Mi,] , (Li3 , (Ljj)3

Case ENY: elevation A and gradient component 3⁄4s specified, gradient component 4. to be determined.

xi, 2- x (Li) )3tz31 - 2y31,,3⁄4 + M,,]

Case ENN: elevation x specified, gradient (zi,, 3⁄4) to be determined.

z,,, x 2 + zi, x 2 z'3, i

zjx !Eýi+ 3⁄4j, x 2 lgd!) Y" [Sj +M

Case NRN: gradient direction (j ,, Pjv) specified, elevation A and iradient

multiplier -y, to be determined.

(L,)3 + ,•i x

(4)L

j (L,,)3 ' (L,,)3 ,(L,,.)3 +3f

Figure 19. Optimality Conditions for yertices by Type.

54
CaseNXN: gradient component 4,, specified, elevation z, and gradient component zi, to be determined.

Xj x 6F1 ~- +q x 3(L) - 3zxiis. + 3M,,1

Fzji3 +-Li3 (j)

j L,) (L,,)3 ,(L,,) 3

Case NNN: gaino cozmpoetters specified, elevation zj and gradient ,, c~ om

pnn ,tobe detrmined.

I + 1 xX3F 3 x -j =6z - 3y 6zw+ 3Mj]

j (Lji)3 (L,,)3 x2 (Li)3

Case NXY:, gopradientersa. specified, elevation and torbedientexi;ANed.o

Zj x Yiji

F Ligure 19. (Cotined

Z-j 3 . x2 ziyji= E Z" 3zj55j

I-3ji3xv j Tjj3(j)
the chain of triangle. obtained in the course of this procedure will contain the

given point (Figure 20).

It is clear that this procedure will work best for a sequence of points in 'which

the previous point can serve as a close starting point f or the task of locating

the given point in a triangle. For points in a regular grid arranged by

sequential rows the following procedure is used. Locate the first point of the

first row from an arbitrary starting point. Make a note of the first row and

its triangle for further reference. Use it also as starting point for the

second point in the first row. Then use the second point and its triangle as

starters for locating the third point, and so on, until the end of the row is

reached. Then retrieve the first point and its triangle and use them as

starters for locating the first point in the second row. This point and its

triangle are again kept for further reference, while the second row is traversed

in the same manner as the first row. This process is repeated for the remaining

rows.

56

.......... j
Specified point whose

triangle is to be

determined

Starting point

"with known triangle

Figure 20. Line Search for a Triangle Containing a Specified Point.

57
6. RESULTS AND CONCLUSIONS

In order to test the computational feasibility of the approach proposed in the

previous sections, a pilot implementation was applied to the Mustang Mountain

data set (see Section 2). CPU times were recorded and experience about memory

requirements was gained. As a preliminary testing procedure the residuals~

obtained when trying to recover the contour elevations were gathered and

analyzed. After describing the set-up of the experiment, we report its results

and the observed computational effort.

6A~ SETTING UP THE EXPERIMENT

The pilot implementation consists of several independent modules whose output

files serve as input files to subsequent modules. The interrelationship of

these modules and their interconnecting files is schematically described in

Figure 21.

The original input file has the format (see Section 2) of Digital Graphic

Recorder Data (DGR) and refers in our case to the Mustang Mountain, Fort

Iluachuca, area. These data are input into the module EDIT, which extracts and

edits contour information as descri.bed in Section 2 of this report. The

resulting edited digital contour file is then sampled and contour tangents are

determined in module THIN. The sampling method is described in Section 3 of

this report. The resulting sample is first fed to the module VORONQI (see

Section 4), which determines a Voronoi triangulation of the sample points.

Accordingly, the main output o' this module is a "triangle table" that lists the

vertices and neighbors of each triangle. In addition, duplicate sample points

are identified and points that need to be added, such as map corners, are

recorded. This information is utilized by the UPDATE module which creates the

final "data base". This data base is needed along with the triangle table to

generate the surface in module SURFACE. This module also provides options for

evaluating the elevation of either random points or points in a regular

58j
Digital Graphics Records Da*la(DGR)

Edited Contour Data THIN

(Pin Sapl & DrvtesVORONOI

Figre 1. ystm ayot.

(-Verte5ata Base

( for testing only)

L. ~SU RFAC E

(Digital Terrain Elevation Data (DTEDo

Figure 21. System Layout:.

59
rectangular grid. The design of the testing procedure is based on the following

observation. When the sample points enter the triangulation process they are

considered as points in a plane. No longer do they hold any direct relationship

to other points in the contour as such. Therefore, these original contour

points are in essence independent of the sampled ones. The testing procedure

consists of evaluating the surface at the original contour points and comparing

the results to the given elevations, yielding a set of residuals. These

residuals provide an indication of the ability of our algorithm to "recover" the

original data.

The full original data set contains features such as lake shores, lake

hatching., dams, and peak elevations, which the algorithm in its present state

of development does not yet handle in an accurate fashion. Indeed, the primary

goal of the effort reported here is to establish computational feasibility. Also

some limitations of the testing procedure itself need to be pointed out.

First, the digitized data themselves carry a "digitization error" so that it is

not necessarily desirable to reproduce the digitized data precisely. Indeed, if

the recovered contours represent a "smoothing" of the digitized lines, they may

be more representative of the true surface than the given digitized contour

points.

Second, in flat terrain, the nonsampled digitized contour points tend to be

close to the boundary of triangles, so that the behavior of the surface in their

interior is monitored to a less extent than in steep terrain, where more contour

lines cut through the interior of triangles. In order to retain interior

monitoring, the sample was deliberately chosen somewhat smaller than indicated

for the purpose of improved accuracy, but still large enough to testI

computational feasibility.

Third, the vertical deviations measured by the residuals tend to be dispropor-I

tionally large in steep terrain. Ideally, the 3-dimensional distances from the

true data points to the generated surface should be evaluated in order to

determine the accuracy of the latter. In flat terrain the vertical

60
deviation, that is, the residual, provides a good approximation to the true

distance. However, this is not the case for steep terrain. For this reason we

introduce an

residual

(6.1.1) "adjusted residual" -

4 M + 1

based on a measure a of "steepness" of terrain. The length of the gradient

(5.1.16) appears to be a natural measure of steepness% However, the example of

a mountain peak, where the slope is by definition zero, shows that the slope at

a single point is not a good indicator of steepness. For the purposes of this

report, we determine the triangle containing the point in question and then

chose the vertex slope of largest magnitude, taking into account that the unit

length in the plane is 20. feet (Section 2.1). The adjusted residual would

represent the 3-dimensional surface precisely if the surface were linear, that

is, a tilted plane in a suitable neighborhood of the data point. In general,

however, it is still an approximation, but a better one than the unadjusted

residual. For horizontal terrain, m - 0 and the adjusted residual equals the

original one. In all other cases, the adjusted residual is smaller. In our

experiment, we collected statistics on both types of residuals.

In Mandel, Witzgall, and Bernal (86), the results of a first test run were

reported. For this run, the unadjusted residuals for the full set of digitized

contour points were collected, including lake hatchings and dams, even though

the algorithm is not yet equipped to handle nonsmooth terrain, as pointed out

above. Nevertheless, 95% of the residuals were between +12.5 feet, indicating

t~at at least 90% did not deviate more than halfway to the next contour line.

Furthermore, an analysis of the biggest residuals led to the discovery of

several contour lines whose altitudes had been apparently miscoded. For the

purpose of the more extensive experiments reported here, the original data set

was purged of lake hatchings and the altitudes were recoded f or the above

contour lines. In addition, one contour line representing a dam was removed

from consideration. In what follows, the results are based on this "sanitized"

data set.

61
6.2 RESULTS

The experiment reported here comprises three runs, all for the same "sanitized"

Mustang mountain data set of about 38,000 sample points giving rise to about

75,000 triangles. For the first and main run, a histogram of the (unadjusted)

residuals is displayed in Figure 22. In addition, the standard statistical

quantities such as expected value (-average), standard deviation, maximum and

average absolute value are reported, the latter two also for the adjusted

residuals (6.1.1). The second run did not utilize the contour tangent

information and the third run used linear rather than Clough-Tocher

interpolation on the given set of triangles.

The timing of the procedure is broken down by major steps, including the

generation of a 901 X 901. rectangular grid, which we expect to be of a size

relevant to prospective applications. The calculations were timed and carried

out on a VAX 11/780 system at the Engineer Topographic Laboratories. The

observed CPU times for the different steps follow

o Step 1: Editing Digitized Input Data ......... 25 min

o Step 2: mhining and Tangent Determination .... 9 min

o Step 3: Voronoi Triangulation and Update ..... 21 min

o Step 4: Surface Generation ................... 20 min

o Step 5: Grid Determination (901 X 901) ....... 26 min

Total CPU Time ..................... 101 min

The timing of Step 5 represents an improvement over the one reported in Mandel,

Witzgall, and Bernal (86).

As discussed before, the testing procedure consisted of calculating and

analysing the restiduals of the elevation of the original contour points. The

histogram of these residuals is displayed in Figure 22. other statistics

calculated include:

62
300HISTOGRAM OF RESIDUALS

200

. . t , % * 5 3 1 7 ..

6% 5%2.4%.O 0* 

21

-2 1 100 -5.0 5 0152

RESIDUALS (FT)

Figure 22. Histogram of Residuals.

63
"o Expected Value (average residual) ........... 0.294 FT

"o Standard Deviation.......................... 5.934 FT

"o Average Absolute Error ...................... 3595 FT

"o Maximum Absolute Error ...................... 99.610 FT

"o Percent Absolute Values < 12.5 FT ........... 95.141 2

"o Number Absolute Values > 99.5 FT ........... 1

"o Number of

Positive Residuals (overestimates) ... 222 417

Negative Residuals(underestimates) ... 205 365

Zero Residuals ....................... 38 678

"o Maximum Absolute Adjusted Residual .......... 65.588 FT

"o Average Absolute Adjusted Residual .......... 3.12. FT

"o Percent Absolute Adjusted Values > 40 FT .... 0.020 X

The adjusted residuals are seen to be much smaller than the unadjusted ones.

This indicates that big values of the (unadjusted) residuals are largely

confined to steep terrain. Results for the second run are displayed below. It

is seen that not to use tangential information results in a definite

deterioration of accuracy:

"o Expected Value (average residual) ........... 0.209 FT

"o Standard Deviation .......................... 6.034 FT

"o Average Absolute Error ...................... 3.657 FT

"o Maximum Absolute Error ...................... 160.619 FT

"o Percent Absolute Values < 12.5 FT ........... 95.044 %

"o Number Absolute Values > 99.5 FT ........... 15

"o Number of

Positive Residuals (overestimates) ... 221 420

Negative Residuals(underestimates) ... 206 419

Zero Residuals ......................... 38 621

64

- - - - - - - - - - - - - - - -
As expected, the third run, featuring linear interpolation, did not achieve the

accuracy of the first run. It comes as somewhat LtZ a surprise, however, that

the loss of accuracy is not more pronounced. An analysis of the results shows

that this is due to the extraordinary large number of zero residuals. This

phenomenon, in turn, is an artefact of the digitization: digitized contour lines

contain many groups of successive points that lie on a common line. If two of

such points are vertices of the same triangle, then they and all intermediate

points reproduce elevation under linear interpolation. As we pointed out

earlier, the precise reproduction of digitizcd points is not necessarily

desirable because the latter carry digitization errror. Below are the

atatistics for the linear run:

o Expected Value (average residual) ............. 0.326 FT Ic Standard Deviation .......oo..........o...... 5.856 FT

o Average Absolute Error .................... 3.557 FT

o Maximum Absolute Error.......... o.......... 106.557 FT

oPercent Absolute Values < 12.5 FT ......... 95.277%

oNumber Absolute Values > 99.5 FT ............... 4

o Number 3f

Positive Residuals (overestimates) ... 190 778

Negative Residuals(underestimates) ... 175 542

Zero Residuals ................... ....'100 140

6.3 CONCLUSIONS

Two particular concerns were our capability of obtaining the Voronoi

triangulation of a sufficiently large set of sample points and of solving the

large linear system for the gradients at the sample points, using commonly

available computer resources. With respect to these concerns we found that

samples of up to 70,000 points could be triangulated within a reasonable time

frame. In fact, CPU time has been less of a problem than memory space, a

limitation that was overcome by developing a decomposition method. The

65
Gauss-Seidel method described in Section 5.2 was found to be faster than

anticipated, both with respect to the CPU time required by a single iteration

(approximately 150 CPU seconds/iteration) and the speed of convergence (ten

iterations). Alternative methods will perform as well or faster.

To sum up, we conclude that the computational effort of using nonlinear

techniques for the generation of a smooth synthetic surface and subsequent

regular grid is substantial, but within the bounds for routine calculations on a

computer of medium size such as the VAX 11/780. Thus we were succesful in

achieving the major goals of this feasibility study.I

The residuals obtained when trying to recover the original digitized contour

points are by and large comparable to the resolution of the digitized data.

However, there are instances of very large residuals which may well be

unacceptable for subsequent applications. Such discrepancies were expected

because our method at this point still lacks the capability to handle

cartographic features at which the actual terrain surface is not smooth, that

is, it exhibits discontinuities of slope. Such features include lake shores,

river banks, as well as some ridge and drainage lines. Modeling by a smooth

surface may not be sufficiently accurate at such locations.

In order to achieve the full accuracy of our approach, the following measures 

"o Cancel the smoothness requirements along lake shores, river banks, as

well as along certain ridge and drainage lines.

"o Smooth digitized contour lines and determine tangent directions prior

to thinning.

"o Investigate adaptive tolerance selection schemes for higher density

sampling in ro'ugh terrain.

"o Compare local and global methods for specifying tangent planes.

66
REFERENCES

Bentley, J.L., Weide B.W. and A.C. Yao, 1980, "Optimal Expected Time Algorithms

for Closest Point Problems", ACM Transactions on Mathematical Software,

Vol. 6, No. 4, pp. 563-580.

Birkhoff, G. and L. Mansfield, 1974, "Compatible Triangular Finite Elements",

Journal of Mathematical Analysis and Applications, Vol. 47, pp. 531-553.

Bowyer, A., 1981, "Computing Dirichlet Tesselations", The Computer Journal, Vol.

24, No. 2, pp. 162-166.

Clough, R.W. and J.L. Tocher, 1965, "Finite Element Stiffness Matrices for

Analysis of Plates in Bending", Proceedings of the Conference on Matrix

Methods in Structural Mechanics, Air Force Institute of Technology,

Wright-Patterson A.F.B., Ohio.

Davis, D.M., Downing, J.A., and S. Zoraster, Algorithms for Digital Terrain

Data Modeling, USAETL Technical Report ETL-0302, July 1982.

Defense Mapping Agency Hydrographic/Topographic Center, Informal Communications,

1985.

Douglas, D.H. and T.K. Poiker, 1973, "Algorithms for the Reduction of the Number

of Points Required to Represent a Digitized Line or Its Character", The

Canadian Cartographer, Vol. 10, pp. 112-123.

Grotzinger, S., Danielson, B., Caldwell, D., and B. Mandel, 1984, The

Contour-to-Grid Problem. An Evaluation of Simplicially Based Local

Surface Estimators (SWIRL, CONTOCES, and TRENDVOL), Interim Paper, U.S.

Army Engineer Topographic Laboratories, Fort Belvoir, Virginia.

Lang, T., 1969, "Rules for the Robot Draughtsmen", The Geographical Magazine,

Vol. 42, pp. 50-51.

67
Lawson, C.L., 1972, Generation of a Triangular Grid with Application to Contour

Plotting, Technical Memorandum No. 299, Jet Propulsion Laboratory,

Pasadena, California.

Lawson, C.L., 1976, C 1 -Compatible Interpolation Over a Triangle, Technical

Memorandum No. 33-770, Jet Propulsion Laboratory, Pasadena, California.

Lawson, C.L., 1977, "Software for C1 Surface Interpolation", in Mathematical

Software III, J.R. Rice, editor, Academic Press, New York.

Mandel, B., Witzgall, C. and J. Bernal, 1'36, "Non-Linear Contour-to-Grid

Digital interpolation", to appear, Proceedings of the M,C&G-Conference,

Oktober 1986.

Reumann, K. and A. Witkam, 1974, "Optimizing Curve Segmentation in Computer

Graphics", International Computing Symposium, Amsterdam, Holland, pp.

467-472.

Peucker (Poiker), T.K. and D.H. Douglas, 1975, "Detection of Surface-Specific

Points by Local Parallel Processing of Discrete Terrain Elevation Data",

Comput. Graphics Image Process. Vol. 4 ,pp. 375-387.

Williams, C.M., 1978, "An Efficient Algorithm for the Piecewise Approximation of

Planar Data", Comp,-ter Graphics and Image Processing, Vol. 8, pp.

286-293.

Williams, C.M., 1981, "Bounded Straio,•-Line Approximation of Digitized Planar

Curves and Lines", Computer Graphics and Image Processing. Vol. 16, pp.

370-381.

Witzgall, C., Bernal, J. and B. Mandel, 1985, On Sampling and Triangulating

Large Digitized Contour Data Sets, Interim Paper, U.S. Army Engineer

Topographic Laboratories, Fort Belvoir, Virginia.

68
Zienkiewicz, O.C., 1971, The Finite Element Method in Structural and Continuum

Mechanics, 2nd edition, McGraw Hill, New York.

Zoraster, S., Davis, D. and M. Hugus, 1984, Manual and Automated Line

Generalization and Feature Displacement, Final Report for Contract

DAAK7O-82-C-0149, U.S. Army Engineer Topographic Laboratories, Fort

Belvoir,Virginia, Report No. ETL-0359, 121 Pages.

69