Fast computation of Delaunay triangulations

Fast computation of Delaunay triangulations Y. C O R R E C Centre d'l~lectronique de l'Armement, 35170 Bruz, France E. C H A P U I S *

Chantiers de l'Atlantique, 44600 Saint-Nazaire, France

We propound a method designed for the triangulation (in Delaunay's sense) of large sets of scattered points in R 2. This method is derived from the well known Watson's algorithm. After the insertion of each new point, in a sequential way, we use Lawson's circle criterion in the local retrlangulation process. Preliminary sorting of data, a 'good' numbering of triangles and use of a tree structure in local retriangulation, results in a 0 (N logN) computing time, as shown in examples. Key words: interpolation, triangulation, surfaces, computational geometry. 1. INTRODUCTION The aim of this study is the development of surface approximation software for digital terrain modelling: given a large set (10.000 to 1.000.000 points) of scattered data (xi, Yi, zi, i = 1...IV) extracted from a geographical data base (for example contour lines), construct a bivariate function f (x, y) which interpolates the data values, i.e. f (xi, Yt) = zi,

i=l...N. Among the interpolation methods tested in ref. 2, only local methods induce reasonable computing times, and most of them are based on a triangulation of the domain. Moreover, it must be pointed out that triangulation process may be involved in other applications, such as synthetic image generation. We first recall briefly some results on Delaunay triangulations, then we present the statement of our method, and we compare our results with those reported in refs. 3, 5, 6 and 8 for N < 2000, with those issued from our own implementation of algorithms4 for N < 500, and ref. 8, for N < 5000. Finally we state our results for data sets containing up to 100.000 points and the corresponding conclusions about the implementation of interpolation algorithms.

of S, is the geometric dual of the Dirichlet (or Voronoi, or Thiessen) tessellation P constructed on S. This tessellation divides the plane into polygonal regions, called tiles. Each tile Pj contains all the points of the plane closer to the tile generating point Nj than to other generating points N k (in the sense of euclidean distance):

Pj = ( N ( x , y ) ~ 2 ; V

k~/

d(N, Ni)<~cl(N, ND)

The edges of these tiles are the perpendicular bisectors of the segment lines that join the vertices of the triangulation, i.e. of the edges of I~launay triangles. Figure 1 shows an example of Dirichlet tessellation and its dual triangulation. From interpolation point of view, Delaunay triangulation is optimal, for its triangles are as much equiangular as possible ('locally equiangular '7) which is quite interesting for finite element method grids. Following Watson, s the basis of our algorithm is the insertion of the points of S, one by one, in a dynamic Delaunay triangulation ~i built on (N/,/= 1 . . . i ) , using the circumscribed circle criterion (see refs. 1, 4-6 and 8). Its equivalence with Thiessen's and max rain angle criterions has been established in refs. 5 and 6, as the equivalence between local and global triangulation optimisation, for strictly convex domains. Figure 2 shows application of this criterion to a quadrilateral.

1

/

2. DELAUNAY TRIANGULATION Given a set S of arbitrarily distributed distinct points (Nj: (xj,y]), ] = 1... N) in ~2. A Delaunay triangulation & Accepted December 1986. Discussion closesJune 1987. * This work was performed while E. Chapuis was at the Centre d'Electronique de l'Armement.

0141-1195/87/020077-07 $2.00 © 1987 Computational MechanicsPublications

s t

~

Figure 1. Dirichlet tessellation ( - - - ) triangulation ( )

and Delaunay

Adv. Eng. Software, 1987, Vol. 9, No. 2

77

12t----

:3 16

I non optimal triangulation

Figure 2.

optimal triangulation

Circle criterion

We will not treat the degeneracy problem, approached in refs. 4-6, as its occurrence probability in actual, physical, data can be considered as zero. Anyway, we believe (as in ref. 1) that a simple perturbation technique is sufficient to solve this problem. Moreover, the S domains we deal with are subregions of a much more wider one, so that the boundary problem, abundantly treated in literature 3-6, is discarded in this study: the boundary of & will be very naturally the convex hull of S.

9

Figure3.

lnsertion polygon [

) o f the new point (e)

16

3. TRIANGULATION ALGORITHM The method we propose is an improvement of the one developed in ref. 8 by Watson, and used in refs. 6 and 1.

3.1. Watson's method Given S, the process is initiated with a first fictitious triangle large enough to enclose all the data points. Then the iterative process works by inserting the i + lth point of S in triangulation gi (built up on the first i points of S), which results in triangulation &i.l. This operation consists of two steps:

5

3

B

4

9

5

2

(1) Search among the 2i + 1 triangles of gi those which circumscribed circle contains the point to be inserted

(Ni÷O. From the list of the vertices of these peculiar triangles (noted T~k), extract the so-called insertion polygon defined as by the union of the Tg (as shown in Fig. 3). (2) Destroy the edges interior to the insertion polygon, and connect its vertices to the new point, achieving the local retriangulation stage. Figure 4 shows triangulation ~i+, thus obtained. It must be pointed out that in ref. 8, search for the Tg of the insertion polygon is made among all the triangles of gi, which leads to a 0 (N l+a) time (0.5 < a ~< 1) for the algorithm published by Watson (as band width is 0(x/N)): this is inconsistent with Watson's claim of 0(NtS), but verified by our tests.

3.2. Proposed method: link and sort 3.2.1. Link: definition of a tree structure on the triangulation (insertion tree). We will call strong neighbour of a triangle Tk any triangle having exactly one edge in common with it, and weak neighbour, any triangle having one vertex only in common. Then suppose that we dispose at this stage of the list of each triangle's strong neighbours. The insertion of node Ni. , in triangulation gi is performed in four steps:

78

Adv. Eng. Software, 1987, Vol. 9, No. 2

Figure 4. Local retriangulation inside the insertion polygon: new edges ( )

(1) Search among ~i's triangles the one (T~) containing Ni+ l (triangle no. 7 on Fig. 3). (2) To find the insertion polygon, construct the following tree: - the root is T~ - at the next level, build the branches (and nodes) corresponding to the current triangle's strong neighbours (except its father, already stored upper in the tree). Stop the process as soon as the explored neighbour's circumscribed circle does not contain the point to be inserted (Ni+0, then, when completed, the insertion tree's leaves are triangles sharing at least one edge with the insertion polygon, but are exterior to it. And the other triangles' union (the other nodes of the tree) is the insertion polygon. Figure 5 shows insertion tree associated with Fig. 3. (3) Local retriangulation is identical with Watson's method step 2.

(4) Update the neighbouring list of all the triangles involved in the tree. The mean number of operations required is then half the one in Watson's method, for random data points.

Then, to speed up the search of the tree's root, i.e. Tt~, we propose the following improvement of the algorithm. 3.2.2. Sort. Basically, the first good thing to do is to sort the data points by ascending values of the co-ordinate (x or y) corresponding to the largest variation

D x = max(x/) -- min(x/) or

Dy = max(y/) --rain(y/) ('main' co-ordinate).

Figure 5.

(

)

Insertion tree and insertion polygon's boundary

Figure Z

Figure 6. Partial g2so triangulation of Fig. 7"s &soo. Fourteen shaded triangles along the moving front are not in their final configuration

500 points. Delaunay triangulation gsoo

Figure8. 50 points. Delaunay triangulation. Triangles numbering follows ascending abscissae order

Adv. Eng. Software, 1987, Vol. 9, No. 2

79

This preliminary step, widely suggested in literature, substantially reduces the later search calculations. But we also modify the retriangulation stage in the Watson algorithm (step 2) in the following way: - Sort the triangles included in the insertion polygon (i.e. the insertion tree's nodes except its leaves) by ascending numbers. Sort the insertion polygon edges by ascending values of their midpoint's main co-ordinate. These three sorts will insure us that, if ~i's triangles are already numbered by ascending order of their barycenter's main co-ordinate, so will be &i.l's triangles. So that computing mean time required to find T~ in &i, by examining triangles in reverse order from the last one is 0 (log N) (instead of the 0 (~/N) crude search time. 4.

TESTS AND PERFORMANCE EVALUATION

and 0(N 1'6) + 0(N 1'2) + 0(N logN) + 0(N) for the number of operations required by the triangulation process (the (N logN) coming from the initial sort). These two authors come then to the conclusion that 0(N logN) is a lower bound for the expected number of operations (which is debatable as a complete sort of the data is not necessarily an obligation). Renka 6 quotes a result from Shamos giving 0(N logN) as a lower bound for a triangulation with memory requirements of 32 N

WATSON

CpU

WATSON

LAWSON

/~AWSON

[8]

tested

[5]

i n [6]

N 1`5

N].£

NT M or NlogN

Hl.i 2

N min

200

N max

1200

Hemory

16N

{0 5000

RENKA me

)

N].16

RE{~A

~p,I~LIRE

meth.2

tested

HI,9

N2

[3] N 1.7

pROPOSED t,~THOD N log N

25

IO00

I000

I000

I0

I000

I0

500

2000

2000

ZOO0

500

2952

100 000

lSN

gN

ISN

20N

23N

t

Theoretical estimates of Lawson s and Renka 6 are respectively

CFU

16N

Figure 9.

0(N ''3) + 0(N logN) + 0(N)

'

I

'

'

'

'

'

,

I

,

Test results comparison

,

,

,

,

[

T R [ A ~ U G A T I ~ N

/i /

tO ~

/ /

// tO'

/ /

10 °

/

/

/

,/

/

/

,/

10t

NOI'IBRE OE P O I N T S

lO t

10 =

- I I_~FRAI~HE ~- • blAq'SON COC 1 7 0 - 7 3 0 0

• klATSON U N I V A C :

NOTRE

CCS^/L^ 17

t180

I'~THOOE

Figure 10

80

lO 3

Adv. Eng. Software, 1987, Vol. 9, No. 2

EG&

lIAR5

r~cE ~ - LlU~I3~IIE

1986 ~{IVI~

In the table of Fig. 9, we present a comparison between several methods, their asymtotic CPU estimates, their memory requirements, and the size of the treated examples (n.b. there is a little inconsistency between columns 1 and 2, i.e. Watson's publication 8 and our implementation of Watson's software given in ref. 8. We believe that the number of examples given in ref. 8 was not enough for a good regression on the exponent of N: cf. Fig. 1). Figure 10 shows CPU time evolution as N increases, for the three methods we tested on UNIVAC 118013 (Lafranche, Watson and ours), and the results obtained on CDC 170-730 by Watson.8 The difference between Watson's results quoted in ref. 8 and ours, using his software, appears clearly here. Figure 11 shows CPU]Iog(N) as N increases, for the two versions we have implemented of our method: the 'all in core' first version, limited to N = 9000, and the 'virtual memory' second version, theoretically unlimited, in fact tested up to N = 100000. Their CPU times are in proportion of 1 to 1.35 or 1.71 (optimised or automatic swapping on UNIVAC). Figure 12 shows CPU/(NlogN) as N increases, and constitutes an experimental verification of the assertion

i

CPU

i

i

,

i

,

TREAI'~I~UI..,,&rION/I.,OG[

i

,

I

,

,

,

,

that our method works in O(NlogN) time (0.54x10 -3 N logN for the first version). It must be noticed that preliminary data sort time is not included in these statistics, but, as shown in Fig. 13 and in ref. 9, the CPU time required is 0(0.7 x 10-SN logN), which is negligible (1%) compared to triangulation time (and that is the same convergence order...).

5. CONCLUSION The introduction of some changes, described here as link and sort procedures, in the local retriangulation process involved in Watson's algorithm, leads to a substantial reduction of CPU time required to construct Delaunay triangulation on large scattered data sets, thus reaching the 0(NlogN) theoretical lower bound so often put forward in hterature, but apparently, at least at our knowledge, never attained. This allows us to consider surface interpolation with finite elements for a reasonable cost, comparable with more classical methods but profiting by finite elements' flexibility and accuracy, which is especially interesting in the treatment of the large data sets involved in digital terrain modelling.

,

, !

v

~

B

i

w i

i

i

|

i , / ~

N )

/ r

101

/

./

10o

10"

10..a

OE POINTS i~

f

,

q

. ~,

I I0

I SANS r'IErIOIRE VIRTUEL,I.,E - : .A'qEC r'IEI101FE ¥[RTUEI.LE

~ z

,

'~

,

'F

r

F,I

~ 'I.0='

,

F

,

~i

, ~,

I

i=

~"

,

'F,F

10"

10

t7 CI~&

C C S A / L ^1 ~ 5

mRS "r~E N

- I,.lllIl~,ll~ll[

I

2 g~lrlC

Figure 11

Adv. Eng. Software, 1987, Vol. 9, No. 2

81

=01. l ' i, '

i,

,

1

'

ff "ON'6 "lOA 'Zg6 [ 'axa°UJ°S 'lu3 "nPV Cll

i

*0~. I ' il

'

i,

,

~I

g[ dartS!d ,01.

, _....I..~

S,l.N]Od

6~Oi

/ /"

~]-

lI

~WJ

Vq/iS33 i

)06 :"

O 0 0 ;B

,

,

O00Z •

SJ,N ] O d

30

l

,

,

i

,

•

•

,

,

,

I

.

.

.

.

.

]ui

,

|

000~

0009

I

0007 '

""

"

'

I

i

,

,

,

O00F., .

.

.

.

•

'

'

"

"

na'~

,

O00C

I

0

000 1

I

'

'

I

'

'

'

7~000"0

'

31:lftJON 9~000"0

8~000"0

09000"0

Z9000"O

~9000=0

99(XX)'O

89C~G°0

OLO00"O

ZL(X~O'O

9LO00"O !

£ 9861

N 33V1~ ~ li

8LO00"O

I

V~I/VSO3

8LIX)O'O it I ' N t i ' l ' i " l l N ) / N r i l l V " l l " l ~ N ¥ l l = L l

. . . .

I

,

,

,

,

i

I

I"

i

i

i

I

J

i

,

I

I

I

I

i

I

i

I

fld'l

•

•

•

i

08000"0

Figure 14.

Digital terrain model. 6 Ion x 9 kin. Isometric view (Alpes, near Nyons)

T h e a l g o r i t h m d e s c r i b e d in this p a p e r is p a r t o f a s o f t w a r e d e v e l o p e d at C E L A R for digital t e r r a i n m o d e l l i n g . This D T M s o f t w a r e is c u r r e n t l y being used to i n t e r p o l a t e e l e v a t i o n d a t a ( f r o m c a r t o g r a p h i c o r s t e r e o p h o t o g r a m m e t r i c origin) to a u n i f o r m grid f o r m i s c e l l a n e o u s a p p l i c a t i o n s . T h e size o f d a t a files p r o c e s s e d ranges f r o m 50000 t o 300000 points. A n e x a m p l e is given in F i g 14, with an i s o m e t r i c view o f a 6 x 9 k m D T M p r o d u c e d in t h a t way. This w o r k m a d e a p p a r e n t s o m e r e l a t i o n s between the c r i t e r i o n d e f i n i n g the t r i a n g u l a t i o n , the i n t e r p o l a t i o n process used o v e r triangles a n d the physical n a t u r e o f d a t a . But the choice o f " g o o d " elements, a n d their efficient use, will be d e a l t with in a f o r t h c o m i n g p a p e r .

REFERENCES 1 Cavendish, J. C. et al. An approach to automatic three dimensional f'mite element mesh generation, Int. J. Num. Mech. Eng. 1985, 21,329

2 3 4 5 6 7 8 9 10

Franke, R. Scattered data interpolation: test of some methods, Math. o f Comp. 1982, 38,157 Hermeline, F. Triangulation automatique d'un poly~dre en dimension N, RAIR O - Analyse Numdrique 1982, 16 (3), 211 Lafranche, Y. Application de l'interpolation polynomiale au d6pouillement graphique de valeurs de R 3, ThJse 3d cycle, 1984, Rennes Lawson, C. L. Software for Ct surface interpolation, in Mathematical Software 11I, Academic Press, New York, 1977 Renka, R. J. A storage efficient method for construction of a thiessen triangulation, ORNL CSD-I01, Oak Ridge National Laboratory, 1982 Sibson, R. Locally equiangular triangulation, Comp. J. 1978, 21 (3), 243 Watson, D. F. ACORD = Automatic Contouring of Raw Data, Computer & Geoscience 1981, 8 (1), 97 Am61ioration de la m6thode de Hibbard et du Quick sort, note CCSA/AN no. 35, Centre d'Electronique de l'Armement Bruz Correc, Y. C-6n6ration automatique d'une triangulation sur des points donn6s, note CCSA/LA no. 49, Centre d'Electxonique de l'Armement ~ Bruz, 1980

Adv. Eng. Software, 1987, Vol. 9, No. 2

83