AlllDE TDEfllfi

DF COMMERCE

tandards and Technology

M«««^'^^''-^ '®^" °^ Standards)

Qdod .U5753 NO.-

NIST Technical Note 1252

On Constructing Delaunay Triangulations

for Sets Constrained by Line Segments

Javier Bernal

QC

100

.U5753

No. 1252

1988

C.2

QjCIDO

NIST Technical Note 1252

On Constructing Delaunay Triangulations

for Sets Constrained by Line Segments

Javier Bemal

Mathematical Analysis Division

Center for Computing and Applied Mathematics

National Engineering Laboratory

National Institute of Standards and Technology

(formerly National Bureau of Standards)

Gaithersburg, MD 20899

Issued September 1988

NOTE: As of 23 August 1988, the National Bureau of

Standards (NBS) became the National Institute of

Standards and Technology (NIST) when President

Reagan signed into law the Omnibus Trade and

Competitiveness Act

^*?-4TESO^*^

U.S. Departmoit of Commerce

C. William Verity, Secretary

National Institute of Standards and Technology

(formerly National Bureau of Standards)

Ernest Ambler, Director
National Institute of Standards U.S. Government Printing Office For sale by the Superintendent

and Technology Washington: 1988 of Documents

Technical Note 1252 U.S. Government Printing Office

33 pages (Sept. 1988) Washington, DC 20402

CODEN: NBTNAE
CONTENTS

Page

Abstract 1

1

.

Introduction 2

2. The UPDTRI algorithm 5

3

.

The EDGSTR algorithm 10

4. The SEGMNT algorithm 18

5. Proof of complexity 21

References 29

111

On constructing Delaunay triangulations for

sets constrained by line segments

Javier Bernal

National Bureau of Standards, Gaithersburg, MD 20899, U. S. A.

Abstract. In this paper, we propose a simple algorithm for constructing a Delaunay triangulation for a finite set of points in the plane constrained by a finite collection of line segments.

This algorithm constructs first a Delaunay triangulation for the set, and then generates from it a

sequence of triangulations as each line segment is incorporated into the previously obtained triangulation. An expected time analysis shows the algorithm to be linear if the number of line segments

is kept constant.

Key words: algorithm; computational complexity,- computational geometry; constrained Delaunay triangulation*, expected time analysis; multiply-connected polygonal regions', Voronoi diagram

AMS(MOS) subject classifications. 68U05
1. Introduction

Let 5 be a finite set of points in the plane. By a triangulation for S we mean a finite collection

of triangles in the plane having pair-wise disjoint interiors, each of which intersects S exactly

at its vertices, and the union of which is the convex hull of S. Given a triangulation T for

5, we say that T is a Delaunay triangulation for S if for each triangle in T there does not

exist a point of S inside the circumcircle of the triangle. Delaunay triangulations have been

thoroughly studied in [8,11,12].

We say that a region in the plane is a multiply- connected polygonal region if it is a bounded and closed set whose interior is nonempty and connected, and whose boundary

is the union of a finite nuniber of closed line segments. Here and in the sequel, given a multiply-connected polygonal region i?, we denote by INT{R), BND{R), and VTX{R),

respectively, the interior of i?, the boundary of i?, and the set of vertices of the boundary

of R. In addition, given points P and Q in the plane, P ^ Q, we denote by [P,0], (P^Q),

[P, Q), and (P, Q], respectively, the closed, open, and half-open line segments that P and Q

determine.

De Floriani, Falcidieno, and Pienovi [5] have studied the concept of a Delaunay triangulation for a finite set of points in the plane constrained by the boundary of a multiply-connected

polygonal region. Let i? be a multiply-connected polygonal region, and let 5 be a finite subset of R that contains VTX{R). De Floriani, et al. define a triangulation for S constrained

by BND{R) as a finite collection of triangles having pair-wise disjoint interiors, each of

which intersects S exactly at its vertices, and the union of which is R. Given a triangulation

T for S constrained by BND{R), De Floriani, et al. call T a Delaunay triangulation for

S constrained by BND{R) if for each triangle in T there does not exist a point of S inside

the circumcircle of the triangle such that the line segments that connect the point and the

vertices of the triangle are contained in R.

Let T be a triangulation for a set S constrained by the boundary of a multiply-connected

polygonal region R. Given a triangle t in T, we denote by ADJ{t) the set of points in S

that are vertices of triangles in T adjacent to t, and say that t satisfies the circle criterion

in T if there does not exist a point in ADJ{t) inside the circumcircle of t. De Floriani, et

al. state that using arguments similar to those developed in [6,8], it can be shown that T

is a Delaunay triangulation for S constrained by BND{R) if each triangle in T satisfies the

circle criterion in T.

In this paper, we study a more general type of triangulation. Let 5 be a finite set of

points in the plane, and let E he a, finite collection, possibly empty, of line segments that

intersect only at points in 5 and whose endpoints belong to S. We say that a triangulation
T for 5 is a triangulation for S constrained by E if for each e \n E and each f in T, e does

not intersect INT{t). Given T, a triangulation for S constrained by E, we say that T is a Delaunay triangulation for S constrained by E if for each t in T there does not exist a point

P of 5 inside the circumcircle of t such that no e in £^ intersects the interior of the convex

hullof iU{P}.

Delaunay triangulations constrained by line segments have been studied by Lee [7] and

Chew [3,4] for line segments that do not cross. Lee and Chew have also shown how to

construct them in 0{N (logTV)^) and 0{N logN) worst-case complexity, respectively, where

A^ is the number of points in S. In this paper, we present an algorithm, based on a different

approach, that computes this type of triangulation in expected linear time for fixed E. It consists of two steps. In the first step, a Delaunay triangulation for S is constructed by

applying to S one of several existing algorithms [1,2,9,10]. In the second step, a sequence

of triangulations for S, the last of which is the desired triangulation, is generated from the

initial triangulation as each line segment in E is incorporated into the previously obtained

triangulation.

Given a simple polygon R, and points Pq and Pi in VTX{R) such that [Pq, Pi] is a linear

component of BND{R), we say that R is star-shaped relative to [Pq, Pi] if for each point P

in P \ [Po,Pi], there exists a point Q in (Po,Pi) such that {P^Q) is contained in INT{R).

We say that R is edge star-shaped if it is star-shaped relative to at least one of the linear

components of BND{R). Figure 1 illustrates an edge star-shaped simple polygon.

Given a region R in the plane, we say that P is a generalized multiply- connected polygonal

region if it is a bounded and closed set whose interior is nonempty, and whose boundary is the union of a finite number of closed line segments. Given a generalized multiply-connected

polygonal region P, for some unique positive integer A;, let Pj, z = 1, . . . , A:, be the unique

multiply-connected polygonal regions that are the closures of the connected components

of the interior of P. Let 5 be a finite subset of P that contains VTX{Ri) for each ^,

i = 1, . . . , fc, and for each z, z = 1, . . . , A;, let S^ be the intersection of S and Pj-. Given for

each z, z = 1,.. .,^% a triangulation Ti for Si constrained by BND{Ri), we call the union

of the triangulations Tj, z = 1, . . . , A;, a triangulation for S constrained by BND{R), where

BND[R) denotes the boundary of P. We call it a Delaunay triangulation for S constrained

by BND{R), if for each z, z = 1, . .

. , A;, Ti is a Delaunay triangulation for Si constrained by

BND{R,).

A crucial aspect of the algorithm we propose involves edge star-shaped simple polygons

and generalized multiply-connected polygonal regions: the incorporation of each line segment of E into the previously obtained triangulation for S. To do this, the algorithm first identifies those triangles in the triangulation whose interiors intersect the line segment. A
Figure 1: An edge star-shaped simple polygon. For each point P in the complement of

[Po^-Pi] with respect to the polygon, there exists a point Q in [Pq.Pi) such that {P,Q) is contained in the interior of the polygon.
finite collection of edge star-shaped simple polygons with pair-wise disjoint interiors is then

produced as the union of the identified triangles is intersected with each closed half-plane

that the line segment determines. Next, for each polygon in this collection, a monotonically

increasing (ordered by inclusion) finite sequence of generalized multiply-connected polygonal

regions that converge to the polygon is produced; a Delaunay triangulation for the intersection of S with each of these regions constrained by the boundary of the region is computed;

and in the limit a Delaunay triangulation for the intersection of S with the polygon constrained by the boundary of the polygon is obtained. Finally, a new triangulation for S

results as the triangulations obtained in the limit as above are appended to the previously

obtained triangulation for 5, and the triangles identified as having interiors that intersect

the line segment in E are eliminated from this triangulation.

Three algorithms, which we call SEGMNT, EDGSTR, and UPDTRI, are outlined and

justified in this paper. Algorithm SEGMNT, which makes use of algorithm EDGSTR, constructs a Delaunay triangulation for a finite set of points in the plane constrained by a finite

collection of line segments. Algorithm EDGSTR, which makes use of algorithm UPDTRI,

computes a Delaunay triangulation for a specialized subset of the boundary of an edge

star-shaped simple polygon constrained by the boundary of the polygon. Finally, algorithm

UPDTRI constructs a Delaunay triangulation for a set constrained by the boundary of a

generalized multiply-connected polygonal region from a specialized triangulation for the set constrained by the same boundary. In Sections 2, 3, and 4, we outline and justify algorithms

UPDTRI, EDGSTR, and SEGMNT, respectively. In Section 5, we show SEGMNT to be

of expected linear complexity if the cardinahty of the collection of line segments is kept

constant. Finally, in Section 6, we offer a summary of the results presented in this paper.

2. The UPDTRI algorithm

Consider a generalized multiply-connected polygonal region i?, and for some positive integer

k, let Ri, i — 1,. . . ,fc, be the multiply-connected polygonal regions that are the closures of

the connected components of the interior of R. Let S be a finite subset of i? that contains

VTX{Ri) for each i, i — 1, . .

. , k. For each i, f = 1, . .

. , ^, let Si be the intersection of S and

i?,, and let T,- be a triangulation for Si constrained by BND{Ri) such that for an integer /,

1 < / < A;, a triangle t* in T/, and a vertex P* of T, the following conditions hold:

1. If for some triangle t' in T/, t' ^ t* , P* is inside the circumcircle of t\ and the hne

segments that connect P* and the vertices of t' are contained in i?/, then one of these

line segments intersects INT{t*).
2. For each triangle t' in Ti^ t' ^ t* ^ there does not exist a point of 5/ different from P*

inside the circumcircle of t' such that the hne segments that connect the point and the

vertices of t' are contained in /?;.

3. For each 2,1 < i < k, i ^ I, Ti is a. Delaunay triangulation for S^ constrained by

BND{Ri).

With T defined as the union of the triangulations Tj, i — !,...,/:, we present and justify

algorithm UPDTRI, which produces from T, t*, and P*, a Delaunay triangulation T* for S

constrained by BND[R). Here and in the sequel, given a triangle t with vertices Pi, P2, P3,

in the order in which they appear in BND{t) in a counterclockwise direction around INT{t),

we denote t by AP1P2P3 and say that AP1P2P3 identifies t. Based on this convention, we let P' and P" be points in 5, for which AP*P'P" identifies t*. The outline of UPDTRI follows.

procedure UPDTRI(r, T\P% P', P")

begin

J- • J1 '^ adj • -^ ) '^ CUT • -' 5 flag := 1; while {flag = 1) do

begin

if (there does not exist P such that APPcurPadj

identifies a triangle in T* or (there exists P such that

APPcurPadj identifies a triangle in T* and

P is not inside the circumcircle of AP*PadjPcur)) then

if (Pad, ^ P') then

begin

P := point such that AP*PPadj identifies a triangle in T*\

p p.p. p -'^ cur • -• adj) -^ adj • — -* end

else flag :=

else

begin

P := point such that APPcurPadj identifies a triangle in T*;

T* := (T* \ {AP*PadjPcur. APPcurPadj}P

{AP*P,djP.^P*PPcur]\
Padj '= P

end

end

end

Justification o/UPDTRI. That the execution of UPDTRI produces a triangulation T* from

T follows from the outline of the algorithm. Thus, assume that UPDTRI has been executed.

For each z, z = 1,...,/:, let T* be the set of triangles in T* whose vertices are all in Si.

Since i? is a generalized multiply-connected polygonal region, it follows from the outline of

UPDTRI that for each i, z = 1, . . . , A;, i ^ /, T* equals T^, so that by the third condition

stated above, each triangle in T* satisfies the circle criterion in T* . Also, from the outline

of UPDTRI, by the first and second condition stated above, each triangle in T* fl Ti satisfies

the circle criterion in T* . Thus, if T* ^ Ti, it suffices to show that each triangle in T* \ Ti

satisfies the circle criterion in T* . We assume this is not the case. Therefore, since P* is a vertex of every triangle in T* \ Ti^ points Q\ Q'\ and Q'" exist in 5*/ such that AP*Q'Q"

and AP*Q"Q"' identify triangles in T; \ 7) with Q'" inside the circumcircle of AP*Q'Q".

Clearly, it follows that a triangle can be identified with AQ'Q"Q"' and P* is not inside the

circumcircle of this triangle.

Let Q' and Q'" be points in Si such that AQ'Q"Q"' identifies the first triangle in T; efiminated

by the algorithm with Q" as a vertex.

If Q' equals Q' and Q'" equals Q'" then P* is not inside the circumcircle of AQ'Q"Q"\ a contradiction to the definition of AQ'Q"Q"'

If Q' does not equal Q' then from the outline of UPDTRI and the definitions of Q' and Q",

there exists a point Q* in Si separated from Q" by the line through P* and Q' such that

AQ*Q'Q" identifies the first triangle in T/ eliminated by the algorithm with Q' as a vertex.

Figure 2 illustrates an example of this situation.

Clearly, [Q\Q"'] and [Q",Q"'] are contained in Ri, and we show that so is [Q*,Q"']. Let Ci

and C2 be the quadrilaterals with vertices P*, Q', Q", Q"\ and P*, Q*, Q', Q", respectively.

Since the fine through P' and P" separates P* from Q\ Q", Q'", and Q* , and Ci and

C2 are convex, it follows that Ci U C2 is convex so that [Q*, Q'"] is contained in Ci U C^-

Thus, [Q*,(5'"] is contained in P/, and by the second condition stated above Q'" is not

inside the circumcircle of AQ*Q'Q". In addition, since triangles can be identified with

AQ*Q'Q\ AQ'Q"Q'\ and AP*Q'Q'\ Q* and Q'" lie in the half-plane that contains P* and

that is determined by the line through Q' and Q". Therefore, as shown by the example in

Figure 3, the circumcircle of AQ'Q"Q"' contains the part of the circumcircle of AQ*Q'Q"

that is contained in this half-plane. Thus, by the definition of AQ*Q'Q", P* is inside the
i

Figure 2: AQ'Q"Q"' and AQ^Q'Q" when Q' does not equal Q'. The hne through P* and Q'

separates Q* and Q".
Figure 3: The circumcircles of AQ'Q"Q"' and AQ*Q'Q".
circumcircle of i^Q'Q"Q"'^ a contradiction.

Similarly, if Q'" does not equal Q'", a contradiction is obtained, and the justification of

UPDTRI is complete.

3. The EDGSTR algorithm

Let R be an edge star-shaped simple polygon, and let e be a linear component of BND{R)

such that R is star-shaped relative to e. Let 5 be a finite subset of BND{R) such that

S contains VTX{R) and does not intersect the relative interior of e. Let r be the number

of points in S, and let Pi, . .

. , P,. be the points in S in the order in which they appear in

BND{R) in a counterclockwise direction around INT{R) with [Pi,Pr] equal to e. Define a one-to-one function F from {!,..., r} onto S by setting F{i) equal to P, for each i,

i = l,...,r.

In what follows, we present and justify algorithm EDGSTR which computes a Delaunay

triangulation T for S constrained by BND{R). Here, given points Qi, Q2, Q3 in the plane,

Q2 ¥" Qs ^-iid Q2 ^ Qi, we let m{Q2Q3,Q2Qi) represent the size in radians of the angle

produced by a counterclockwise rotation around Q2 from ray Q2Q3 to ray Q2Q\- The

outline of EDGSTR follows.

procedure EDGSTR(r, P, r)

begin

1. T := 0; G(l) := P(l); G{2) := P(2); J := 2;

2. for / := 3 until r do

begin

3. J:^J + 1-G{J):=F{I);

4. Qi := G{J - 2); Q2 := G( J - 1); Q3 := G{J)-

5. while {m{Q2Q3,Q2Qi) < tt and J ^ 2) do

begin

6. T:=TU{AQ,Q2Q3};

7. UPDTRi(P,r,g3,Qi,Q2);

8. J -.^ J - 1; G{J) := Q3;

9. if (J 7"^ 2) then

begin

10. Q,:^G{J-2); Q2 := G{J - 1)

end

end

10
end

end

Justification o/EDGSTR. We show that the execution of EDGSTR produces a collection of

triangles T that is a Delaunay triangulation for S constrained by BND{R). First, we must

prove the following observation:

Given integers j, k, I, l<j<k<l<r, if R contains [Pj^Pk] and [Pk^Pi], m{PkPi, PkPj)

is less than tt, and [Pj^Pk] and [Pk^Pi] intersect S exactly at Pj and Pk, and P^ and Pi,

respectively, then APjP^Pi identifies a triangle that is contained in R and that intersects S

exactly at Pj, P^, and Pi. In addition, given a finite collection of triangles T with pair-wise

disjoint interiors, each of which is contained in R and intersects S exactly at its vertices, if {Pj,Pk) and {Pk,Pi) do not intersect the interior of each triangle in T, and INT{APjPkPi)

does not intersect the interior of each triangle in T with P^ as a vertex, then INT{APjPkPi)

does not intersect the interior of any triangle in T

Figure 4 illustrates an example of i?, P^, P^-, P/, Pi, P^, as described above.

In what follows, we prove the first part of the observation. The second part will follow from

it. That APjPkPi identifies a triangle follows from the fact that m{PkPi^PkPj) is less than tt.

It follows from the ordering of Pj, 2 = 1, . . . ,r, in BND{R), that in order to prove the rest of the first part of the observation, it suffices to show that for each integer i, 1 < i < i or

I < i < r^ Pi does not belong to C = APjPkPi \

([Pj, Pk] U [P/c, Pi])- We assume the contrary

and let H be the set of integers i at which m{PiPi, PiPk) attains its minimum value, with

1 ^ ^ ^ i or / < z < r, and Pi in C. Clearly, each i in H must be different from j and /, so that either l<i<jovl<i<r for each i in H.

li I < i < r for some i in H, we show that for each point Q in (Pi, P^), (P/, Q) is not entirely

contained in INT{R), a contradiction to the fact that R is star-shaped relative to [Pi,Pr].

Thus, let i* be an integer in H such that I < i* < r. It follows that the ray Pi Pi* intersects

[Pj, Pk). Thus, there exists an integer n, j < n < k, and a point P* in BND[R), such that

PiPi' intersects [P„,P„^.i) exactly at P*, and [P*,Pi*] is contained in R. Let R' be the

simple polygon whose vertices are contained in {P*} U {Pj, n < i < i*} and whose boundary

contains P*, Pj, i = n + 1, . . . , z*, in this order in a counterclockwise direction around its interior. Figure 5 illustrates an example of R', Pj, Pk, Pi, Pj>, P-n, Pn+i, P*, A, Pr- Let

Q be a point in (Pi,Pr) and assume that {Pi,Q) is contained in INT{R). It follows from

the definition of R' that [P/, Q] must intersect INT{R'). However, from the definitions of

i*, P*, and the ordering of Pi, i= 1, . . . ,r, in BND(R), Q is not contained in INT{R'), so that [Pi,Q] must intersect BND{R') at a point other than P/. Every point in BND{R')

11
Figure 4: The polygon R and the points P^, Pk, Pi, Pi, Pr- AP^PkPi C R and

AP,PkP,nS = {P,,P,,Pi].

12
Figure 5: The polygon R' (shaded region) and the points P,, P^, P/, P,. , P^, Pn+i, P*, -Pi, Pr. The ray PiPt^ intersects [P„,P„+i) at P*

13
that is not in (P*, Pi») belongs to BND{R) and differs from Q. Thus, [P/, Q] must intersect

[P*^Pi-). But P*, Pj., and P; are coHnear with P,-. in {P*,Pi), so that (P/, Q) must contain

Pi*. Since Q is an arbitrary point in (Pi,P7.), we have the desired contradiction.

Similarly, if 1 < z < j for each i in H, it can be shown that for each point Q in (Pi, Pr),

{Pj,Q) is not entirely contained in INT{R). This is also a contradiction, and it completes

the proof of the first part of the observation.

Let /, J, G, and T be as defined in the outhne of EDGSTR. Let K be a variable which is set to zero at the start of the execution of EDGSTR and which is incremented by one at the

start of each execution of line 6 of EDGSTR. It follows from the outlines of EDGSTR and

UPDTRI, the above observation, the definition of R, and an induction process on K, that

at the start of each execution of line 6, the criteria listed below are met. Figure 6 illustrates an example of R, F, G, T, /, J that satisfies these criteria. The criteria follow:

1. {G(j), J = 1,...,J} C {P(z), ?: = 1,...,/}, G(l) and ^(J) equal P(l) and P(/),

respectively, and given integers k, /, k\ /', 1 < A;, /</,!< k', /' < J, such that F{k)

and F{1) equal G{k') and G(/'), respectively, then A; < / if and only if k' < I'.

2. If J > 3 then m{G{j)G{j + l),G{j)G{j - 1)) > tt for each j,

i - 2 J -2

3. The triangles in T have pair-wise disjoint interiors.

4. Each triangle in T is contained in R and intersects

{F{i), i = 1, . .

. , /} exactly at its vertices.

5. Given integers k, l^ 1 < k < I < I such that [F{k)^ F{1)] is an edge of a triangle in

T, if / equals ^ + 1 or if for some integer j, 1 < j < J — 1, F{k) and F{1) equal G[j)

and G{j + 1), respectively, then [F{k),F{l)] is an edge of exactly one triangle in T.

Otherwise, it is an edge of two triangles in T.

6. Given an integer i, 1 < z < /— 1, [P(^),P(^ + 1)] is not an edge of a triangle in T if and

only if for some integer j, 1 < j < J — 1, F{i) and F{i + 1) equal G{j) and G{j + 1),

respectively.

7. Given an integer j,l<j<J— I, [G{j)^G{j + 1)] is not an edge of a triangle in

T if and only if for some integer z, 1 < i < / — 1, G{j) and G{j -f 1) equal F{i) and

F{i +1), respectively.

14
G(1)=F(1)

F(2)

F(7)

F(9)

G(2)= \

F(4)

G(3)=F(5) \

/^F(12)

F(11)

G(4)=F(8)

G(5) = /FdOt^

Figure 6: A polygon R, points F{i), z = 1, . . . ,r, for r equal to 12, and the corresponding

triangulation T and points G{j), j = 1,...,J, at the start of the execution of line 6 of

EDGSTR for / and J equal to 10 and 5, respectively.

15
8. Given t in T, there exist integers j, A:, /, ij, Z2, is, 1 < J < -^ — 1,

1 < A; < ij < 22 < ^3 < / < /, such that G'(j') and G(j + 1) equal F{k) and -F(/), respectively, and AF{ii)F{i2)F{i3) identifies t.

9. Given integers i, A:, /, 1 < j < J— 1, 1 < /: < / < /, such that G{j) and (^(j + 1) equal

F{k) and i*'(/), respectively, if A: + 1 < / then (i) the union of the triangles in T of the

form AF{ii)F{i2)F{i3), k < ii < 12 < is < /, is a simple polygon that is star-shaped

relative to [G{j),G{j + 1)], and (ii) the closed half-plane determined by a clockwise

rotation of tt radians of the ray G{j)G{j -\- 1) around G{j) contains F{i) for each i, k <i< I.

10. Given an integer k, 1 < k < I, such that G{J — 1) equals F{k), the closed half-plane

determined by a clockwise rotation of tt radians of the ray G{J — 2)G{J — 1) around

G{J — 2) contains F{i) for each i, k < i < L

11. T is a Delaunay triangulation for a subset of {F{i), i = 1, . .

. , /} constrained by the

boundary of the generalized multiply-connected polygonal region that the union of the

triangles in T defines.

Let jr be the initial value of J the first time line 6 is executed for I equal to r. From the

first criterion above and lines 3, 4 and 8 of EDGSTR, at the start of the execution of line 6

for / equal to r and J equal to jr, we must have that G{1), G{J — 1), and G{J) equal Pi,

Pr-i, and P^, respectively. This, together with the second criterion (for the same / and J)

and the ordering of P„ 2 = 1, . .

. , r, in BND{R), implies that m{G{j)G{J), G{j)G{j - 1))

is less than tt for each j, j = 2, . .

. , J — 1. Figure 7 illustrates this situation for the same R

and F of Figure 6. Thus, from the remaining criteria (for the same / and J), it can be seen from lines 6-10 that at the end of the execution of EDGSTR for / equal to r, the following

conditions are met:

1. The triangles in T have pair- wise disjoint interiors.

2. Each triangle in T is contained in R.

3. For each z, i = 1, . .

. , r — 1, there exists a triangle in T with [P{, Pi+i] as an edge.

4. There exists a triangle in T with [Pi, P^] as an edge.

5. Given an edge of a triangle in T, if it is not contained in BND{R) then it is an edge

of two triangles in T

16
G(1)«F(1)-Pi

F(2)

F(7)

F(9)

>G(6) =

G(Jr)

=F(12)=Pi2

G(5)=F(11)=P^^

Figure 7: A polygon i?, points F{i), i — 1, .. . ,r, for r equal to 12, and the corresponding

triangulation T and points G(ji'), j = 1,..., J, at the start of the execution of hne 6 of

EDGSTR for I and J equal to r and _;V, respectively.

17
6. T is a Delaunay triangulation for a subset of S constrained by the boundary of the

union of the triangles in T.

Since these conditions imply that the execution of EDGSTR produces a Delaunay triangulation T for S constrained by BND{R), the justification of EDGSTR is complete.

4. The SEGMNT algorithm

Let 5 be a finite set of points in the plane. Given a positive integer 71/, let ej-, ^ = 1, . . . , M,

be line segments in the plane that intersect only at points in S and whose endpoints belong to

S. Define a one-to-one function D from {!,..., M} onto E = {ci^ i = 1, . . . , M} by setting

D{i) equal to e^ for each z, i — 1, . .

. , M.

In this section, we present and justify algorithm SEGMNT which computes a Delaunay

triangulation T* for S constrained by E. The outline of SEGMNT follows.

procedure SEGMNT(r*, 5, i9, M)

begin

1. obtain a Delaunay triangulation T* for 5;

2. for / := 1 until M do

begin

3. T; := 0;

4. Qq, Qi := endpoints oi D{I);

—

= closed half-plane to 'left' of QoQi] = closed half-plane to 'left' of Q1Q2',

8. flagl := 1; 9. if (each triangle in T* with Qo as a vertex is contained in Hi)

then flagl := else J := 0; 10. while {flagl = 1) do

begin

11. J:=J-M;

12. g", Q'" := points in S for which AQj_iQ"Q"'

identifies a triangle in T* with Q" e Hj and Q'" ^ Hj]

13. if ( J = 1 and Q" e INT{Hj)) then

r;:=T;u{Agj_igv"};

14. Fj{l) := Qj_i; Ej{2) -.^ Q";

5. Q,

6. i/i

7. H,

18
15. L - 2; Q' := Q'": flag2 := 1; 16. while {flag2 = 1) do

begin

17. Q'" := point in S for which AQ'Q"Q"'

identifies a triangle in T*;

18. if ( J = 1 and {Q" e INT{Hj) or Q'" e INT{Hj)))

then r; := T/ U {AQ'Q"Q'"h

19. if (g'" / Hj) then Q' := Q'" else

begin

20. X:=X + 1; Fj{L) := Q'"-

21. if (g'V Qj) then Q" := g'"

else if ( J = 1 and T/ = 0) then flagl, flag2 :=

else

begin

22. rj := L;

23. EDGSTR(rj, Fj,rj);

24. if (J = 1) then //a52 :=

else

begin

25. T* := (r*\r;)uTi ur2;

26. //a^l, //a^2 := end

end

end

end

end

end

end

Justification o/SEGMNT. We show that the execution of SEGMNT produce a collection of

triangles T* that is a Delaunay triangulation for S constrained by E. Here, given a subset

E' of E, a triangulation T for S constrained by E\ and triangles t and f in T, we say that

f is adjacent to t in T relative to E' if the intersection of t and t' is a line segment that is not contained in any of the line segments in E'. In addition, we denote by ADJ{t,E') the

set of points in S that are vertices of triangles in T adjacent to t in T relative to E', and

say that t satisfies the circle criterion in T relative to E' if there does not exist a point of

19
ADJ{t,E') inside the circumcircle of t. As in [5], we note that using arguments similar to

those developed in [6,8], it can be shown that if each triangle in T satisfies the circle criterion

in T relative to E' then T is a Delaunay triangulation for S constrained by E'. Let /, J, T*, r;, Hj, Tj, Fj, rj be as defined in the outline of SEGMNT. We show that

at the end of each execution of line 25 of SEGMNT, T* is a Delaunay triangulation for S

constrained by {e,- : I < i < I}. Clearly, for / equal to M this is the desired conclusion.

From line 1, at the start of the execution of line 3 for / equal to 1, T* is a Delaunay

triangulation for S constrained by {cj : 1 < i < I — I}. Thus, we assume by induction on

/ that given k, 1 < k < M, at the start of the execution of line 3 for / equal to k, T* is a Delaunay triangulation for S constrained by {ci : 1 < i < / — 1}, and show that at the end

of the execution of line 25 for / equal to A:, T* is a Delaunay triangulation for S constrained

by {e,- : 1 < i < /}.

For each J, J = 1,2, at the end of each execution of hne 22, let Rj be defined as the

intersection of Hj with the union of the triangles in Tj. From the definitions of T/, Hj,

Fj, and rj, for each J, J = 1,2, if Rj is not empty for / equal to k, we may assume

without any loss of generality that Rj is a simple polygon that is star-shaped relative to

[Fj{l), Fj{rj)] = D{I), {Fj{L), L = 1, . .

. , rj} equals S D Rj and contains VTX{Rj), the

points Fj{L), L = l,...,rj, appear in BND{Rj) in a counterclockwise direction around

INT{Rj), and D{I) intersects S exactly at its vertices. Thus, for I equal to k and for

each J, J = 1,2, if Rj is not empty EDGSTR can be executed in line 23 for Fj and rj

to produce a Delaunay triangulation Tj for S f] Rj constrained by BND{Rj), and at the

end of the execution of line 25 for the same /, T* is a triangulation for S. In addition,

since D[I) does not intersect the interior of each triangle in T*, it follows from the induction

hypothesis and the definitions of E, Tj, and T2, that T* is a triangulation for S constrained

by {e,- : 1 < ^ < /}. Finally, let i be a triangle in T* at the end of the execution of line 25

for / equal to k, and assume Rj for the same / is not empty for each J, J = 1,2. If t belongs to T* \ (Ti U T2) then from the induction hypothesis it follows that t satisfies the

circle criterion in T* relative to {cj- : I < i < I}. Else, if t belongs to Ti then all triangles

adjacent to t in T* relative to {e,- : I < i < 1} must belong to T* \T2, and since Ti is a Delaunay triangulation for S (1 Ri constrained by BND{Ri), we must have that t satisfies

the circle criterion in T* relative to {cj- : I < i < I]- Similarly, if t belongs to T2 the same

conclusion is obtained. Thus, by induction on /, we have shown that at the end of each

execution of line 25, T* is a Delaunay triangulation for S constrained by {ci : I < i < /},

and the justification of SEGMNT is complete.

20
5. Proof of complexity

In this section, we assume that 5, E, Af, e^, i = 1,. . . , M, are as defined in the previous

section, and that the points in S have been chosen independently from a uniform distribution on a square in the plane. Letting TV represent the cardinality of 5", we then prove that for

fixed M the expected time required by SEGMNT to produce a Delaunay triangulation for

S constrained by E is 0{N).

First, we introduce additional terminology. Let B denote the square from which S has

been chosen, and let d denote the Euclidean length of each of its sides. Let K be defined

as the floor of A'^^'^, i. e. the largest integer less than or equal to N^'"^, and assume that

B has been divided into K^ equal-sized square cells. Let L{N) be defined as the floor of

logA'^ where log denotes the natural logarithm, and let B' denote the square obtained by

surrounding B with L[N) + 1 layers of cells of the type into which B has been divided. Let

S' denote the set of points that are centroids of cells in the L{N) + 1 new layers. Given e in

E, let |e| denote the Euclidean length of e. Given a finite set X of points in the plane and

points P and P' in X^ P ^ P', call P' a Delaunay neighbor of P relative to X if in some

Delaunay triangulation for X there exists a triangle with P and P' as vertices. Given G, a

finite nonempty set in the plane, and Gi, G2, nonempty subsets of G, let V{G, Gi, G2) denote

the set of points in G2 that are Delaunay neighbors of points in G\ relative to G. Given e in

£', let ^(e) represent the union of the sets V[S U {-P}, {-P}, S) for P in e, let V'{e) represent

the union of the sets V{S U 5' U {-P}, {P}, S) for P in e, and let 5(e) represent the union of

the sets V{S U 5' U {P.}, 5', S) for P in e. Finally let 5 represent the set V{S U 5', 5', S).

Next, based on the terminology just defined, we list and prove some observations, the

last of which establishes our assertion about the complexity of SEGMNT. Here, /, J, Fj,

and rj are as defined in the outline of SEGMNT.

1. For each e in £', 5(e) C 5.

2. For each e in £, V{e) C V\e) U 5(e) C V'{e) U 5.

3. For each e in i?, 5 n e C V{e).

4. Given points Qi and Q2 in 5, if Qi is a Delaunay neighbor of Q2 relative to 5 then Qi

and Q2 belong to V{S U {P}, {P},S) for each point P in (Qi, (?2)-

5. Given an integer ^, 1 < i < M, for / equal to i the points Fj{L), L = l,...,rj,

J = 1,2, are contained in ui^T^V{ei) C (^^^^'(e,)) U 5.

21
6. For each e in E, the expected number of points in V'[e) is bounded above by {\e\/d) •

7. The expected number of points in S is 0{N^^^).

8. For fixed M, the expected time required by SEGMNT to produce a Delaunay triangulation for 5* constrained by E is bounded above by

M i 0{N) + i:(((E \^i\ld) + 1)' • 0{N)) = (1 + 0{M')) . 0{N) = 0{N).

Proof of first observation. Let e be a line segment in £^, and let Q be a point in ^(e). By the

definition of 5(6), Q is in 5 and for some point P in e there exists a point Q' in S' such that

Q is a Delaunay neighbor of Q' relative to 5 U 5' U [P]. Thus, we must have that either

(i) Q and Q' are adjacent points on the convex hull of 5 U 5' U {P} or (ii) the line segment

[Q, Q'] is an edge of two triangles in some Delaunay triangulation T' for 5 U 5' U {P}. If

(i) occurs then Q and Q' are also adjacent points on the convex hull of 5 U 5' so that Q is a Delaunay neighbor of Q' relative to 5 U S" and, therefore, Q is a point in S. Else, if (ii) occurs then for some point Q" in 5 U 5' there exists a triangle t in T' such that Q, Q\ and

Q" are vertices of t. By the definition of T' there does not exist a point of 5 U 5' inside

the circumcircle of t. Thus, i is a triangle in some Delaunay triangulation for 5 U 5' and,

therefore, Q belongs to S.

Proof of second observation. Through slight modifications of a proof in [2], it can be shown

that if X and Y are finite nonempty sets in the plane and X' is a nonempty subset of X

then

v{x, x\ X) c v{^x u y, X', X) u v{x u y, y, x).

Thus, given e in £" and P in e, we must have that

y(5 u {P}, {P}, s) = v^s u {p}, {P},s u {P})

C V{S U 5' U {P}, {P}, 5 U {P}) U V[S U 5' U {P}, S\ S U {P})

C V{S US'U {P}, {P}, S) U V{S US'U {P}, 5', S) U {P},

and since P does not belong to V{S U {P},{P}, S) it follows that

y(5 U {P}, {P}, 5) C V{S US'U {P}, {P}, 5) U V{S US'U {P}, 5', 5),

22
which together with the first observation implies the desired result.

Proof of third observation. Let e be a line segment in E, and let P' be a point in 5 Pi e. The

assertion follows, since for P in e sufficiently close to P\ P' is in V[S U {P}, {P}, S).

Proof offourth observation. Let Qi and Q2 be points in 5, and assume that Qi is a Delaunay

neighbor of Q2 relative to S. Let P be a point in (Qi, Q2)i let T' be a Delaunay triangulation

for 5u{P}, and let T'{Qi) be the collection of triangles in T' with Qi as a vertex. If Qi does

not belong to V{S U {P}.i{P},S) then P is neither inside nor on the circumcircle of each

triangle in T'(Qi). Thus, (i) the points that are vertices of triangles in T'{Qi) belong to 5,

and (ii) Q2 is neither inside nor on the circumcircle of each triangle in T'{Qi). But this is a

contradiction, since from (i) and the definition of T' a Delaunay triangulation for S exists in

which the collection of triangles with Qi as a vertex equals T'[Q\) so that from (ii) Qx can

not be a Delaunay neighbor of Q2 relative to S. Thus, Qi belongs to V{S U {P}, {P}, S)., and similarly, so does Q2.

Proof of fifth observation. Let T* be as defined in the outline of SEGMNT. Given an integer

^, 1 ^ ^ ^ M1 let T*^ equal T* at the start of the execution of line 3 of SEGMNT for / equal

to i, and let 5' equal the set {Fj[L)., P = 1,.. .,rj, J = 1,2} for / equal to i. In what

follows, we show by induction on i that 5' C \J\_-^V[ei) for each z, 1 < z < M. The rest of

the observation will follow from the second observation above.

Given i, 1 < i < M, let Q be a point in S\ It follows that either (1) Q is a point in 5 fl e^, or (2) for some point Q' in S and some triangle t in T*\ Q and Q' are distinct vertices of t, and {QiQ') intersects the relative interior of Cj. If (1) occurs then by the third observation

above Q is in V(ej). Else, if (2) occurs, we assume inductively that if Pi and P2 are distinct

vertices of a triangle in T*' then either (a) i > 1 and Pi is a Delaunay neighbor of P2 relative

to 5, or (b) i > 1, and Pi and P2 belong to U]ZlV{ei). It then follows from (a), (b), and the

fourth observation above that Q is in U}^-xV{ei), and from the definition of T*\ for i < M

the induction assumption holds when i is replaced by i + 1

Proof of sixth observation. As in [1], given a point P in the plane, a finite subset X of the

plane, and a collection O of octants around P as shown in Figure 8, we say that P is closed

in X relative to O if at least one point in X different from P is contained in each of the

octants in O.

Let e be a Hne segment in E., let /<'(e)-represent the floor of {\e\l d)- N^^'^ ., and without regard

to any specific order, let Qq and QK{e) denote the endpoints of e. For each A:, A: = 0, ... , /'^(e),

let Pk denote the point

Qo^{klK{e))-{QK(e)-Qo).

and let O^ denote the collection of octants around P^^ for which the interior of each octant

23
Figure 8: A point P in the plane and a collection of octants around it.

24
does not intersect e. Figure 9 illustrates an example of Pk, O^, k = 1, . .

. , K{e).

Let pi be the probability that for each k, k = 0, . . . , K{e), Pk is closed in 5 U 5' relative to

Ok, and let ii be the expected number of points in V'{e) when for each fc, A; = 0, . . . ,/i'(e),

Pk is closed in 5 U 5' relative to Ok. Let p2 and t2 be similarly defined, respectively, for the

event in which for some k, < k < K{e), Pk is not closed in S U S' relative to OkGiven k, < k < A'(e), Pk is not in any of the outermost L{N) layers of cells of B'.

Thus, it follows from the selection of S and 5', and results in [1], that if Pk is closed in

S L) S' relative to Ok then the expected number of points in V{S U 5' U {Pk} , {Pk} , S U S')

is 0(1). In addition, we note that for each A:, k = 1, . .

. , K{e), the length of [Pk-i, Pk] equals

\e\/K{e) w djN^/^ (the length of each side of a cell), and e equals \Jf,=\[Pk-\tPk\- Hence, it follows from the geometry of the octants around the points Pjt, A; = 0, . .

. , K{e), as illustrated

in Figure 9, that if for each k, k = 0,. . . ,K{e), Pk is closed in SU S' relative to Ok, then the

expected number of points in Up^eV{S U 5' U {P}, {P}, S U S') is K{e) • 0(1). Therefore,

h = K{e)-0{l) = {\e\/d)-0{N'^').

On the other hand, if for some k, < k < K{e), Pk is not closed in S U S' relative to

Ok, then using an argument from [1], the probabiUty of this occurring is bounded above by

8.(/i:(e) + l)-e-0(^(^)'). Hence,

P2<8-(/i:(e) + l).e-^(^(^)^).

Finally, <2 is clearly bounded above by the number of points in S. Therefore,

h<N.

Thus, if t is the expected number of points in V(e), then

t = Pi-h+P2-t2 = 1 • {\e\/d) 0{N^I'') + 8 • {K{e) + 1) • e'^^^^^)') • N

= (|e|/cf)-0(iV^/2) + 8.(/^(e) + l).0(l) = i\e\/d)-0{N'/') + {\e\/d)-0{N'^') = {\e\ld)-0{N''').

Proof of seventh observation. Let S" represent the set of points in S' contained in the first layer of cells in B' that surrounds B. From the definition of 5', points in S" are the only

points in S' that can be Delaunay neighbors of points in S relative to 5 U S' . Since points

25
Figure 9: The points and collections of octants P^, Ojt, A: = 0, . .

. , K{e).

26
in S" are not contained in any of the outermost L{N) layers of cells of i?', it follows from

results in [1] that the expected number of points in 5 U 5' that are Delaunay neighbors of a given point in S" relative to 5 U 5' is 0(1). Therefore, since the number of points in S" is 0{N^^^), we must have that the expected number of points in

S = V{SU 5", 5', S) = V{SU S', S", S)CV{SU S\ S", S U S')

is 0{N'/^) 0(1) - 0(7Vi/2).

Proof of eighth observation. Let Tj be as defined in the outline of SEGMNT. Given z,

1 < ^ < M, it follows from the fifth observation above that for / equal to i the points Fj{L),

L = 1, . . . ,rj, J = 1, 2, are contained in Ul_-^V{ei) C {L}]_-^V'{ei)) U S. In addition, from the

outline of EDGSTR, we obtain that for I equal to i the worst-case time required by the

execution of EDGSTR to produce Tj, J = 1,2, in line 23 of SEGMNT is 0(r5). Therefore,

from the sixth and seventh observations above, we must have that for / equal to i the

expected time required by the execution of EDGSTR to produce Tj, J = 1,2, is bounded

above by

{i±{{\ei\/d).OiN'/'))) + 0{N'/')y = {{{i:\ei\/d) + l).OiN'/')r

1=1 1=1 = {d:\ei\/d)+ir.o{N).

1=1 The result now follows by assuming that the Delaunay triangulation for S obtained in line 1 of SEGMNT is constructed by applying to S an expected 0(7V) algorithm [1,2,9,10].

6. Summary

We have presented a simple algorithm for computing a Delaunay triangulation for a finite set 5 of points in the plane constrained by a finite collection E of line segments. The algorithm

is based on the simple idea of constructing first a Delaunay triangulation for 5, and then

'repairing' this triangulation as each line segment in E is incorporated into the previously

obtained triangulation. As it turns out, the 'repairing' of the triangulation is an easy process,

since it involves computing the Delaunay triangulation of subsets of edge star-shaped simple

polygons constrained by the boundaries of the polygons. Given that A'^ is the number of

points in S, and M is the number of fine segments in E, we have also shown that if the

points in S are chosen independently from a uniform distribution on a square, then the

27
expected time required by the algorithm to produce the desired triangulation is bounded

above by 0{M^ N). Actually, depending on the lengths of the line segments and their

locations with respect to the boundary of the square, we expect this time to be 0{M N) or better for N large relative to M, since in this case most of the triangles that are eliminated

during the 'repairing' of the triangulation belong to the Delaunay triangulation for S. At

any rate, for fixed M, the expected time required by the algorithm to produce a Delaunay

triangulation for S constrained by E is 0{N). Finally, we report that an implementation of

the algorithm has just been completed at the National Bureau of Standards.

Acknowledgments

The author gratefully acknowledges Professor James F. Lawrence and Mr. Donald J. Orser

for reading the manuscript and offering their criticism and helpful suggestions. He also

wishes to thank Dr. Christoph Witzgall for his encouragement during the preparation of the

manuscript.

28
References

[1] J. L. Bentley, B. W. Weide and A. C. Yao, Optimal expected-time algorithms for closest

point problems, ACM Trans. Math. Software 6, No. 4, 1980, 563-580.

[2] J. Bernal, Expected 0{N) and 0{N'^^^) algorithms for constructing Voronoi diagrams in

two and three dimensions, NBSIR 87-3679, National Bureau of Standards, Washington,

D. C, 1987.

[3] L. P. Chew, There is a planar graph almost as good as the complete graph, in Proc. 2'^'^ ann. Symp. on Computational Geometry, Yorktown Heights, 1986, 169-177.

[4] L. P. Chew, Constrained Delaunay triangulations, in Proc. 3'"'' ann. Symp. on Computational Geometry, Waterloo, Ontario, 1987, 215-222.

[5] L. De Floriani, B. Falcidieno and C. Pienovi, Delaunay-based representation of surfaces

defined over arbitrarily shaped domains, Comput. Vision Graphics Image Process. 32,

1985, 127-140.

[6] C. L. Lawson, Software for C^ surface interpolation. Math. Software III, (J. R. Rice,

Ed.), pp. 161-194, Academic Press, New York, 1977.

[7] D. T. Lee, Proximity and reachability in the plane. Tech. Rep. R-831, Coordinated

Science Laboratory, University of Illinois, 1978.

[8] D. T. Lee and B. J. Schachter, Two algorithms for constructing a Delaunay triangulation. Int. J. Comput. Inform. Sci. 9, No. 3, 1980, 219-242.

[9] A. Maus, Delaunay triangulation and the convex hull of A'^ points in expected linear

time, BIT 24, 1984, 151-163.

[10] T. Ohya, M. Iri and K. Murota, A fast Voronoi-diagram algorithm with quaternary tree bucketing, Inform. Proc. Letters 18, No. 4, 1984, 227-231.

[11] M. L Shamos, Computational geometry, Ph.D. dissertation, Yale University, 1978.

[12] M. L Shamos and D. Hoey, Closest-point problems, in Proc. 16*'^ Ann. IEEE Symp. on

the Foundations of Comp. Sci., 1975, pp. 151-162.

29

NBS-n4A (REV, z-ec) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NlST/TN-1252

2. Performing Organ. Report NoJ 3. Publication Date

September 1988

4. TITLE AND SUBTITLE

On Constructing Delaunay Triangulations of Sets Constrained by Line Segments

5. AUTHOR(S)

Javier Bernal

6. PERFORMING ORGANIZATION (If joint or otijer than NBS, see instructions)

NATIONAL BUREAU OF STANDARDS

U.S. DEPARTMENT OF COMMERCE

GAITHERSBURG, MD 20899

7. Contract/Grant No.

8. Type of Report & Period Covered

Final

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State. ZIP)

Same as item 6

10. SUPPLEMENTARY NOTES

I I Document describes a computer program; SF-185, PIPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual sumn)ary of most significant information. If document includes a significant

bibliography or literature survey, mention it here) In this paper, we propose a simple algorithm for constructing a Delaunay triangulation for a finite set of points in the plane constrained by a finite collection of line segments. This algorithm constructs first a Delaunay triangulation

for the set, and then generates from it a sequence of triangulations as each line

segment is incorporated into the previously obtained triangulation. An expected

time analysis shows the algorithm to be linear if the number of line segments is

kept constant.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

algorithm; computational complexity; computational geometry; constrained Delaunay

triangulation; expected time analysis; multiply-connected polygonal regions;

Voronoi diagram

13. AVAILABILITY

r~Al Unlimited

I I For Official Distribution. Do Not Release to NTIS

r~Xl Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402. Q^ Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

33

15. Price

USCOMM-DC 6043-P80


U.S. Department of Commerce

National Institute of Standards and Technology

(formerly National Bureau of Standards)

Gaithersburg, MD 20899

Official Business

Penalty for Private Use $300

'"^^O S^«^^

Stimulating America s Progress 1913-1988