NBSIR 87-3679

Expected 0(N) and 0(N^3

) Algorithms for Constructing

Voronoi Diagrams in Two and

Three Dimensions

Javier Berna!

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Applied Mathematics

Mathematical Analysis Division

Gaithersburg, MD 20899

November 1987

U.S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS

NBSIR 87-3679

EXPECTED 0(N) AND 0(N4/3

)

ALGORITHMS FOR CONSTRUCTING

VORONOI DIAGRAMS IN TWO AND

THREE DIMENSIONS

Javier Bernal

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Applied Mathematics

Mathematical Analysis Division

Gaithersburg, MD 20899

November 1987

U.S. DEPARTMENT OF COMMERCE, C. William Verity, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Expected 0(N) and 0(N4 / 3

)

algorithms for constructing Voronoi

diagrams in two and three

dimensions

Javier Bernal

National Bureau of Standards, Gaithersburg, MD 20899.

Abstract. Bentley, Weide and Yao have proposed an expected 0(N) cell technique for computing Voronoi diagrams in two dimensions that does not generalize

readily to three. In this paper their work is further developed and generalized

to produce expected 0(N ) and 0(N4 / 3

) algorithms for constructing Voronoi diagrams in two and three dimensions, respectively. Computational experience is presented for the algorithm in two dimensions.

Key words, algorithm, computational geometry, computational complexity,

Voronoi diagram, expected time analysis

AMS(MOS) subject classifications. 68U05

1. Introduction

Consider a finite set 5 of points in Euclidean space, and let P be any point

in 5. The Voronoi polyhedron V(P,S) of P relative to S is the set of all

points in the space such that P is as close to any point in this set as is

any other point in S. V(P, S) is the result of a bisection process for P

with respect to S; i. e., V[P,S) is the intersection of the closed half-spaces

that contain P and that are determined by the perpendicular bisectors of

the line segments connecting P with the other points in 5. It follows that

1
V(P, S) is a closed convex polyhedron, possibly unbounded, which contains

P in its interior. The Voronoi polyhedra of the points in S (relative to 5 if

not stated otherwise) fill the space without common interior points, and the

union of their boundaries forms a diagram, the Voronoi diagram for 5, which

partitions the space into the Voronoi polyhedra. A point in the space that

is a vertex of some Voronoi polyhedron is said to be a vertex of the Voronoi

diagram for 5. In ^-dimensional Euclidean space a vertex of the Voronoi

diagram for a set S is always the common vertex of at least k + 1 Voronoi

polyhedra of points in 5. It is called a degenerate vertex whenever it is a vertex of more than k + 1 Voronoi polyhedra.

Bentley, Weide and Yao [l] have shown how the Voronoi diagram for a

set 5 of N points, chosen independently from a uniform distribution on a

square in the plane, can be constructed in linear expected time. With M

defined as the largest integer less than or equal to V1 / 2

, i. e. the floor of TV 1 / 2

, Bentley, et al. divide the square into M2 equal-sized square cells. With L(N) defined as the floor of log N (here log denotes the natural logarithm), cells not

contained in the outermost L(N) layers of cells of the square are called ‘inner

cells’, the rest ‘outer cells’. Essentially, the algorithm takes three steps. In

the first step, each point in S is assigned (in constant time) to a cell in which

it is contained. In the second step, the Voronoi polygons of points assigned

to inner cells are constructed. Given a point P assigned to an inner cell, a search is conducted through each of the layers of cells surrounding P for

points assigned to cells in these layers. This search procedure, called a spiral

search around P, starts with the cell assigned to P, and then proceeds in

outward direction to each of the layers of cells surrounding this cell. V

r (P, 5)

is progressively built through a bisection process for P with respect to the

set of points in 5 found through the search. A geometric test is available

for deciding whether V(P, 5) has been obtained. In most cases V

r (P, 5) is obtained after examining only a small number of cells and points so that the

expected time for the second step is O(N). Finally, in the third step, the

Voronoi polygons of those points assigned to outer cells are constructed by

applying to them an 0(K log K) worst case algorithm, e. g. Shamos’ [3].

The expected time for this last step is also O(JV).

In this paper, we present a version of this algorithm in which the third

step is replaced by steps that also require expected time O(N) but that have

the advantage of generalizing to three dimensions. Again, with M defined as

2
the floor of TV 1 / 2

, the square is divided into M2 equal-sized cells. We define

inner cells and outer cells differently from [ 1] . Namely, cells contained in the

outermost two layers of cells of the square are now called outer cells, the

rest inner cells. As in [l], each point in 5 is assigned to a cell in which it is

contained. Voronoi polygons of points assigned to inner cells are constructed

by applying the Bentley, et al. cell-based method together with a generalization of its geometrical test. Finally, a modified version of Bowyer’s insertion

process [2] is used to handle points assigned to outer cells.

An insertion process is a method for updating the Voronoi diagram for a

set of k — 1 points for an additional k

th point in order to obtain the Voronoi

diagram for all k points. The first step in Bowyers algorithm consists of

identifying a vertex of the Voronoi diagram for the k — 1 points that will not

be a vertex of the Voronoi diagram for the k points. To do this, Bowyer uses

a ‘walk’ that starts near the centroid of the k — 1 points and that ends at a

point in whose Voronoi polygon the k

th point lies. The point thus found is

clearly a point among the k — 1 points that is closest to the k

th point. We

have modified this step by taking advantage of the cell structure. In place of a walk we use a spiral search around the k

th point and, as described in [l],

find a point among the k — 1 points that is closest to it. The remainder of

Bowyer’s algorithm, in which vertices of the Voronoi diagram for the k — 1 points which lie in the Voronoi polygon of the k

th point are deleted, and the

vertices of the Voronoi polygon of the k

th point are added, is left essentially

unaltered.

In Section 2, the two-dimensional algorithm is outlined and justified. In

Section 3, we prove the algorithm is of expected O(N) complexity. In Section 4, the algorithm is generalized to three dimensions and shown to be of

expected 0(TV4//3

) complexity. Finally, in Section 5, computational experience with the implementation of the two-dimensional algorithm is presented.

2. The two-dimensional algorithm

In the following, cells are divided into four classes of cells (Figure l). Inner

cells not contained in any of the outermost L(N) layers of cells of the square

are called class 1 cells. Inner cells within T(TV) layers of cells from exactly

one side of the square are called class 2 cells. Inner cells within T(TV) layers

3
2 ceils

L(N)

ceils

co

CO JS

O

Class 4

Figure 1: Regions of square that contain the four classes of cells. Class 1, 2,

and 3 cells are inner cells, and class 4 cells are outer cells.

4
of cells from two sides of the square are called class 3 cells. Finally, all outer

cells are called class 4 cells. We assume that N is large enough so that none

of the four classes is empty. We define 54 C 5 as the set of points assigned

to class 1 cells. 52, 53 , 54 are analogously defined with respect to class 2,

class 3, class 4 cells, respectively.

Throughout the outline of the algorithm that follows, we say a point in

5 has been processed if its Voronoi polygon has been constructed. We say

a cell is empty if no points of 5 are assigned to it. We say an inner cell

has been activated if it has been found to be empty or if each point of 5

assigned to the cell has been processed. We select a centermost class 1 cell.

The order in which inner cells are activated by our algorithm is determined

by proceeding through the layers of cells in the square in a spiral-like fashion

around the centermost cell. We say an inner cell is the currently active cell if

the points assigned to it are currently being processed. We say a point in the

currently active cell is the currently active point if the point is currently being

processed. We say a cell has been searched during a spiral search around a

point in 5 if all points assigned to the cell have already been found through

the search.

We say that points P and P' in 5 are Voronoi neighbors in S if V(P, 5)

and V(P',S) are contiguous in the Voronoi diagram for 5. For each P in

5 we let V

r (P, 5) represent the set of Voronoi neighbors of P in 5, and we

assume V(P, 5) can be readily extracted from V(P, 5). During the execution

of the algorithm, for each P in 5 that has not been processed we define Z

V (P)

as the set of points assigned to inner cells that are known Voronoi neighbors

of P in 5; i. e., the set of points in 5\S4 that have been processed and found

to be Voronoi neighbors of P in 5. Every point assigned to an inner cell is

processed by our algorithm with a spiral search around the point. Given

P, one such point, and assuming it is the currently active point, we let dt represent the radius of the largest circle centered at P whose interior only

intersects cells that have been searched up to the current time. We also let Z*

be the set of points in 5 that have been found up to the current time during

the search. Then, as in [1], we define the tentative Voronoi polygon V1 of P

at the current time as the Voronoi polygon of P relative to {P} U Z

l UZ*(P)

(the plane if Z

l and ZV (P) are empty). Accordingly, when the search starts

and Z

l is still empty, we define the initial tentative Voronoi polygon V

1 of P

as the tentative Voronoi polygon of P at that time. Also as in [l], each V

l is

5
constructed by updating a previous one in such a way that each V

1 is exactly

the result of a bisection process for P with respect to Z

l U Zv(P)-

Given points 0 and O' in the plane, we let QQ' and dist(Q ,0') represent, respectively, the closed linear segment that connects 0 and Q', and the

distance from 0 to O'. In the outline of the algorithm, given Q, a point in

the plane, and W, a finite or infinite subset of the plane, we let dmax(Q , W)

represent the maximum value of dist(Q ,Q') for Q' in W. Given rays rfi and

f2 of common origin, we let 771

( 7

* 1 , 7* 2 ) represent the size of the angle produced

by a clockwise rotation from rfi to f2

. Given points Q and O’ in the plane, we

let QQ

'

represent the ray through O' of origin 0. Given a point P in 52 , we

let both 1\{P) and 1 2 {P ) represent the ray of origin P that is perpendicular

to the side of the square closest to P. Given a point P in S3 , we let li(P)

and l 2 (P) represent the two rays of origin P that are perpendicular to the

two sides of the square closest to P with m(/1 (P), l 2 (P)) equal to 90°. Given

P, a point in S2 or S3 , and assuming it is the currently active point, we

let mt i and mt2 represent at the current time the smallest positive values

of m(PP', l\(P)) and m(/2 (P), PP'), respectively, for P' in Z

l U ZV {P). We

say points Pi and P2 in Z

l U ZV (P) determine mtl and mt2 , respectively, if

mn equals m(P

P

1 .l 1 (P)) and m<2 equals m(/2 (P), P

P

2 ). Assuming points

Pi and P2 determine rri tl and m(2 , respectively, we let C

l represent the interior of that region of the plane obtained by a clockwise rotation from PPi to

PP2. We are now ready to formulate the algorithm.

Start of algorithm.

Step 1 . Assign points to cells and select first class 1 cell to be activated.

Let M be the floor of N1 ^ 2. Partition the square into M2 equal-sized square cells.

Determine inner and outer cells.

Assign each point in 5 to a cell.

For each cell, list the points assigned to it. Determine the centermost cell.

If the centermost cell is empty then go to Step 2.

Else designate this cell as the currently active cell.

Go to Step 3.

Step 2. Select next inner (class 1, 2, or 3) cell.

6
If all inner cells have been activated then go to Step 8 . Else choose the next inner cell to be activated.

If this cell is empty then go to Step 2 . Else designate this cell as the currently active cell.

Determine class for currently active cell.

If class 1 then go to Step 3.

Else go to Step 4.

Step 3. Construct Voronoi polygon of a point in Si. Let P be a point assigned to the currently active cell that has not been

processed.

Designate P as the currently active point.

Start spiral search around P and construct V

1

. Update V1 and each subsequent V1 as appropriate.

For each \n compute Dt = dmax(P,Vt ) and dt

. Terminate search when one of the following criteria is met.

1 . 2 Dt < dt .

2. All cells in the square have been searched.

Upon termination go to Step 7.

Step 4. Begin construction of Voronoi polygon of a point in S2 or S3 . Let P be a point assigned to the currently active cell that has not been

processed.

Designate P as the currently active point.

Determine l\(P) and Z2 (P).

Start spiral search around P and construct V1

. Update Ul and each subsequent Vf as appropriate.

For each U

f compute Dt = dmax{P.Vt ). dtl mfl , and 77i f2 , and determine CL

Terminate search when one of the following criteria is met.

1 . 2 Dt < dt .

2. All cells in the square have been searched.

3. All cells that intersect C

l have been searched.

7
Upon termination, if neither criterion 1 nor criterion 2 has been met

then go to Step 5.

Else go to Step 7.

Step 5. Continue construction of Voronoi polygon of a point in S2 or S3.

Let C = CV

Let Pi and P2 be points in S that determine mtl and mt 2 , respectively.

Compute d! = dmax(P, {Pi, i^})-

Resume spiral search around P.

Update each V

rt as appropriate.

For each V1 compute Dt = dmax(P,Vt ) and dt. Terminate search when one of the following criteria is met.

1 . 2 Dt < dt .

2. All cells in the square have been searched.

3. d! < dt Upon termination, if neither criterion 1 nor criterion 2 has been met

then go to Step 6.

Else go to Step 7.

Step 6. Complete construction of Voronoi polygon of a point in S2 or S3.

Resume spiral search around P

.

Update each V1 as appropriate.

For each \n compute Dt = dmax(P, V1 \ C) and dt

. Terminate search when one of the following criteria is met.

1. 2 Dt < dt.

2. All cells in the square have been searched.

Upon termination go to Step 7.

Step 7. Save Voronoi polygon of a point assigned to an inner cell. Identify U(P, S) with V

1

. Mark P as processed and save V(P, S).

For each P' in V(P, S) that has not been processed let

Z

V (P') = ZV (P') U {P}.

8
Determine whether currently active cell has been activated.

If activated then go to Step 2 . Else if P is in Si then go to Step 3.

Else go to Step 4.

Step 8. Construct and save Voronoi polygons of points in S4 . Determine S4

.

If S4 is empty then stop.

Else perform insertion process on S4 . Perform for each P in S4 a bisection process with respect to

V(P, S4 ) U Z

V (P) and identify V(P, S) with result of process.

For each P in S4 mark P as processed and save V(P, 5).

Stop.

End of algorithm.

Justification of algorithm. As established in [l], the Voronoi polygons of

points in Si can be constructed with Step 3 of the algorithm. Let P, V

4

, Z)t , dt be as defined in Step 3. Let P' be a point in 5 with 2 Dt < dist(P, P'). It follows that P' can not affect V

1 since the perpendicular bisector of P

P' does

not intersect V

1

. Thus, during the spiral search around P, we may conclude

that V(P, S) is equal to V

1 as soon as 2 Dt < dt

. We justify, with the aid of Figure 2 , that we can produce the Voronoi

polygons of points in S2 or S3 with Steps 4, 5, 6 of the algorithm. Let P be

as defined in Step 4, and let C, Pi, P2 , d! be as defined in Step 5. Let V

be the portion of V(P, {P, Pi, P2 }) that is contained in C (shaded region in

Figure 2 ). Let P' be a point in S that does not lie in C or in the interior of

the circles with centers (P + Pi)/2, (P + P2 )/ 2 , and diameters dist(P1 Px ), dist(P, P2 ), respectively, as shown in Figure 2 . It follows that P' does not

affect V as indicated in Figure 2 by the perpendicular bisector b of PPL

We note that points in S contained in C are found through the spiral search

around P with Step 4, and that points in S contained in the circle with center

P and radius d! are found with Step 5 (this circle contains the two circles

mentioned above and is easier to search). Thus, with Step 6 , in which V is not considered during the geometrical test and through which only points

that do not affect V are found, we can produce V(P, 5).

9
Figure 2: P\ a point in 5 that does not affect V (shaded region).

10
Finally, we can construct the Voronoi polygons of points in S4 with Step 8

of the algorithm since two points in S4 that are Voronoi neighbors in 5 must

be Voronoi neighbors in S4 .

3. Proof of complexity

In this section, we assume that the points in 5 have been chosen from a uniform distribution on the square, and prove that the algorithm presented

in Section 2 has expected O(N) execution time. First, we state the following

theorem by Bentley, Weide and Yao [1]. Here, the time involved is defined

as the number of cells plus the number of points examined.

Theorem Let P be a point in the square, and let P' be a point in S closest to

P . Then the expected time required to find P

1 through a spiral search around

P is constant.

By modifying the proof of the theorem, Bentley, et al. also establish the

following two observations crucial to our proof of optimality.

1 . Let P be a point in S such that at least one point in S is contained in each

of the octants around P shown in Figure 3. Then the expected time required

to find at least one point in each octant through a spiral search around P is

constant.

2. Under similar assumptions, let Pl} i = 1 , . .

. , 8, be the first eight points

in S obtained through a spiral search around P such that they are contained in the octants around P, one per octant. Let d and D be the values

of dmax(P,{Pi, ..., P»}) and dmax(P. V(P, {P, . .

.

,

P8 })) ; respectively.

V(P, S) can be constructed by searching only those cells intersecting the interior of the circle with center P and radius 2 D (Figure 3). Since 2 D <y/2d,

it follows from Observation 1 that the expected time required to search all of

these cells through a spiral search around P is constant.

Assigning the points in 5 to the appropriate cells can be accomplished in

0(N) time [l]. Thus, it will suffice to show that the expected time involved in

constructing wdth the algorithm all Voronoi polygons of points in S;, for each

i, i = 1, . .

.

,

4, is bounded above by O(N). Let P be a point in 5 \ S4 so that

V(P, S) is constructed with the algorithm through a spiral search around P.

In what follows, we let w denote the time involved in constructing V(P, 5),

11
Figure 3: Points Pt , i = contained in octants around P. P' is outside the circle of radius 2 D so it does not affect V[P, {P, Px ,

. .

.

,

P&})

as indicated by the perpendicular bisector b of PP'.

12
and use the fact that w is bounded above by

o(j +

i=i

where j

is the number of cells examined with the search, k is the number

of points in those j

cells, and Vi is the number of vertices in the tentative

Voronoi polygon of P when the i th point is found through the search. Finally,

we let E(w) denote the expected value of w, i. e. the expected time involved

in constructing V

r (P, S).

Proof for S\. Let P be a point in Si. As in [1], we say P is closed if at

least one point in 5 is contained in each of the octants around P as shown

in Figure 3.

Let pi be the probability that P is closed, and t\ the expected number of

points examined while constructing V(P, S) with the algorithm when P is

closed. P2 and t 2 are similarly defined, respectively, for P not closed.

If t is the expected number of points examined while constructing V(P, 5)

with the algorithm, then t = pi • t\ + p2

* ^ 2

- Of course pi < 1, and from Observation 2 above, ti = 0(1). Since at most

all points are examined when P is not closed, it follows that t 2 < 0(N).

Next, we find an upper bound for p2 using an argument of [1]. If no points

are found in a given octant, then at least 0(L(N)2

) cells must be empty. The

probability of the octant being empty is then bounded above by \

It follows that p2 < 8

Therefore,

t = 1 • 0(1) + 8e~0(t(Ar)2)

• O(N) = 0(1).

A similar argument can be used to show that the expected number of cells

examined while constructing V(P, 5) with the algorithm is constant. Finally,

since the number of vertices in any tentative Voronoi polygon of P is at most

the number of points examined while constructing V(P, 5), the expected

number of vertices in any tentative Voronoi polygon of P is also constant.

It follows, then, that

E{w) < 0(1) + 0(1) • 0(1) = 0(1)

for each point in Si.

Since at most N points are contained in 51 then the expected time required

13
for Si is

N • 0(1) = O(N).

Proof for S2

. Let P be a point in S2 and, without any loss of generality,

assume P is within L(N) layers of cells from the right-hand side of the

square. Let l represent the ray that both h(P) and Z 2 (P) represent. As

shown in Figure 4, we say P is closed if within the first L{N) layers of cells

that surround P at least one point in S is contained in each of six octants

around P, octants I through V/, and at least one point in 5 (which may be

one of the points in octants I or VI) is found in each of the upper and lower

portions of the outermost layer of the square.

If P is closed let Pi and P2 be points in 5 within the first L(N) layers of cells

that surround P in the upper and lower portions of the outermost layer of

—* — the square, respectively, with the smallest positive values of m(PPi,l) and

m(/, PP2 ). Let C be the interior of the region of the plane obtained by a

clockwise rotation from PPi to PP2

. Clearly, since P is further than two cells from all sides of the square, PPi and

PP2 intersect the boundary of the square at points within t'he first 2 L(N)

layers of cells that surround P. Therefore, by examining the first 2 L(N)

layers of cells that surround P, all cells intersecting C are examined.

Let d' be the value of dmax(P, {Pi, P2 }). Then the circle of radius d' with

center P is also contained in the first 2 L(N) layers of cells that surround P.

Finally, let U be the Voronoi polygon of P relative to those points within

the first L(N) layers of cells that surround P. Note that points outside C

and the circle of radius d' with center P do not affect the part of U that is

contained in C. Accordingly, let D' be the value of dmax(P,U \ C). Then,

since P is closed, the circle of radius 2 D' with center P is also contained in

the first 2 L(N) layers of cells that surround P.

It follows from these observations that the Voronoi polygon of a closed point

can be constructed with the algorithm by examining no more than the first

2 L(N) layers of cells that surround the point.

Let pi be the probability that P is closed, and ti the expected number of

points examined while constructing V(P, 5) with the algorithm when P is

closed. p2 and t 2 are similarly defined, respectively, for P not closed.

If t is the expected number of points examined while constructing V(P,S) with the algorithm, then t = pi • ti -f p2

• t 2.

14
Outermost

2 layers

L(N) cells

/" pr

\

\ Octant II

J2

o

o

\

\\

#

\

\

\

Octant III \

\

\

\

/

Octant IV /

/

/

/

/

/

/ Octant V

/

/

Id

eta nt I /

/

/

/

/

/

Clipper portion

i

\

\

\

\

\

> Lower portion

\

\

eta

VI

nt

\

Figure 4 : A closed point in S2. One point is contained in each of the octants

I through VI. Pi and P2 are contained in the upper and lower portions of

the outermost layer of the square, respectively.

15
Of course pi < 1 , t 2 < O(N), and from the above discussion 0(L(N)2

) is an

upper bound for t j.

In order to find an upper bound for p2 , we argue as follows. If no points are found in one of the octants II through V, then at least 0(L(N)2

) cells must

be empty. If no points are found in one of the octants I and VI, then at

least 0(L(N)) cells must be empty. Finally, if either the upper or the lower

portion of the outermost layer of the square is empty, then 0(L(N)) cells

must be empty. It follows that p2 1 4 e

-o(l'(-^)2

) 4. 4 e -°{ Li N)) m Thus,

t < 1 • 0{L(Nf) + (4e- 0{L<- N),) + 4e~°WN))

) O(N) = 0(L(N)2

) + 0(1) = 0{L(Nf).

An argument similar to that used for points in S\ can be used to show that

the expected number of cells examined while constructing V(P, S) with the

algorithm is 0{L{N)2

), and that the expected number of vertices in any

tentative Voronoi polygon of P is also 0(L(N)2

). Thus,

E(w) < 0(L(N)2

) + 0(L(N)2

) 0(L(N)2

) = 0(L(N)4)

for each point in S2

. The expected number of points in S2 is 0{Nl ^ 2 L[ N)). Hence, the expected

time required for S2 is

0(N1/2 l(N)) 0(L(N

)

4

) = O(N).

Proof for S3 . Let P be a point in S3 and, without any loss of generality,

assume P is within L(N) layers of cells from the right-hand and bottom

sides of the square. Let /j and / 2 represent the rays that l\(P) and / 2 (P)

represent, respectively. As shown in Figure 5, we say P is closed if within

the first L(N) layers of cells that surround P at least one point in S is

contained in each of four octants around P, octants / through IV, and the

right-hand and bottom portions of the outermost layer of the square.

That O(N) is an upper bound for the expected time required for S3 now

follows by an argument similar to the one used for S2. Prooffor S4. The expected time required by the insertion process of Step 8 of

the algorithm is at most proportional to the product of the expected number

16
Outermost

2 layers

L(N) cells

r

rFT

<n

“55o

\ Octant II \

\

Octant

\ \

\

\ I \ i /

Octant III \ 1 / \ I / 4-1 v

/

/

• /

Octant IV /

1 Pi

• . /

/

/

t /

/

/

\ Right-hand

' portion

/

/

/

/

/

/

jL

Bottom portion ^2

Figure 5: A closed point in 53

. One point is contained in each of the octants /

through I\ . Pi and Pn are contained in the right-hand and bottom portions

of the outermost layer of the square, respectively.

17
of points in 54 and the expected maximum number of vertices in the Voronoi

diagram for any subset of S4 . Since the expected number of points in S4 is OfiV1 /

2 ), it follows from the Euler-Poincare formula that this time is at most

0(N1/2

) 0(N1/2

) = O(N).

Finally, it suffices to show that the expected time required by the bisection

process of Step 8 is also at most O(N). Let r be the number of points in S4

. Let i = 1, ... ,r, be the points in S4

. For each i, i = 1 , . .

.

,r, define as the number of points in V’(Pn S4 ), Ui as the final number of points in Z

V (P;),

and Vi as the maximum number of vertices of any polygon obtained during

the bisection process for Pt.

It follows that the time required by the bisection process is at most proportional to

r r r £K- + Ui) • Vi < £(u>* + Ui)

2 = £(w? + 2wiUi + u2

).

i = l i = 1 i=l

Again, r has expected value 0(N1 ' 2

), so that by the Euler-Poincare formula,

J2i=i wi ^ as expected value 0(iV 1//2 ). In order to calculate an upper bound for the expected value of each uj,

i = 1, . .

.

,r, we proceed as follows. Given z, 1 < i< r, let P' be a point in

5 not contained in 54 such that P' is outside the first 2 L(N) layers of cells

that surround P;. As previously proven, P' is a Voronoi neighbor of Pl in

S with probability at most proportional to so that the expected

value for Ui is bounded above by

0(L(N)2

)

• 1 + N- e - 0[L{N)) = 0(L(N)2

) + 0(1) = 0(L(N)2

).

It follows now that the expected value of Ya=i{w{ + 2 WiUi -f u\) is at most

0{N) + 0{N1/2

)

• 0(L(N)2

) + 0{

N

1/2

)

• 0{L{N)

4

) = 0(N);

••

i. e., the expected time required by the bisection process is at most 0(N).

18
4. The three-dimensional algorithm

We now present an algorithm for constructing the Voronoi diagram for a

set 5 of N points contained in a cube in three-dimensional Euclidean space.

First, with M defined as the floor of N

1//3 we divide the cube into M3 equalsized cubic cells. Cells contained in the outermost two layers of cells of the

cube we call outer cells, the rest inner cells. Next, as in the two-dimensional

case, each point in 5 is assigned to a cell in which it is contained. Finally,

the Voronoi polyhedron of each point in S is constructed according to its cell

assignment by generalizing the two-dimensional algorithm of Section 2.

In order to outline the algorithm, cells are further divided into five classes

of cells. With L(N) defined again as the floor of log TV, inner cells not con- tained in any of the outermost L(N) layers of cells of the cube are called

class 1 cells. Inner cells within L(N) layers of cells from exactly one face of

the cube are called class 2 cells. Inner cells within L(N) layers of cells from

exactly two faces of the cube ared called class 3 cells. Inner cells within L(N) layers of cells from three faces of the cube are called class 4 cells. Finally, all

outer cells are called class 5 cells. We assume N is large enough so that none

of the five classes is empty. We define Si C 5 as the set of points assigned to

class 1 cells. 52 , S3 , S4 , S5 are analogously defined with respect to class 2,

class 3, class 4, class 5, respectively.

Throughout the following, definitions and meaning of terminology, such as Z

V (P ), V(P, 5), are as in the two-dimensional case with the words polyhedron and space replacing the words polygon and plane, respectively, when

necessary. However, points in S2 , S3 , S4 require some additional definitions

and terminology which we present separately for the purpose of clarity. Most

importantly, we define symbols C"(P) and d\ for each point P in S2 US3 US4 , and describe what it means to say that P ‘has been closed.’

Definitions and terminology for S2 . Let P be a point in S2. Let F be the face

of the cube closest to P. Let l be the ray with origin P that is perpendicular

to F. Let P' be the point at which / intersects F. Let m be a line through

P' that is perpendicular to an edge of the cube in F

.

Let H" be the plane

parallel to F that contains the point [P + P')/ 2. We define C"{P) as that

closed half-space determined by H" that contains P

.

Assume P is the currently active point, and P[, i = 1 ,..., 8 , are points

in F contained in the octants around P' as shown in Figure 6. We say that

19
at the current time P has been closed and that {P/,z = closes P

if there exist points Pi, i = 1, . .

.

,

8, in Z

l U ZV (P) such that the rays PP;,

i = 1 ,. .

. , 8, intersect F at the points P/, i — 1 ,

. .

. , 8, respectively. Assuming

P has been closed we define d[ at the current time as the smallest value of

dmax(P,{Pl,i = 1,...,8}) for {P/,z = 1,...,8} in the family of sets that

close P.

Definitions and terminology for S3

. Let P be a point in S3 . Let Pj, j = 1,2,

be the two faces of the cube closest to P. For each j, j = 1,2, let lj be the

ray with origin P that is perpendicular to Fj. For each j, j = 1,2, let Pj be

the point at which intersects Fj. Let Pjo, j = 1,2, be the vertices of the

cube common to P4 and P2 in the order shown in Figure 7. Let m be the

line that contains the edge of the cube common to Pj and P2

. For each j,

j — 1,2, let m3 be the line through Pj perpendicular to m. Let ra0 represent

the same line that m2 represents. For each j , j = 1,2, let Pj_i,i and Pj0 be

the closed half-planes determined by mJ _ 1 and mj, respectively, that contain

Pjq. For each y, j = 1,2, let H'f be the plane parallel to Fj that contains

(P + Pj)/2. We define C"(P) as the intersection of the closed half-spaces

determined by FT" and FFf that contain P.

Assume P is the currently active point, and Pji? j = 1,2, z = 0, . .

.

,

7, are points such that with Pj7 = Pj7 , for each j, j = 1,2, PjJ? z = 1,...,6, are points in Pj contained in the six octants around Pj as shown in Figure 7, and

Pj ., 7 and Pj0 are points in Pj_i,i and Pj0 , respectively. We say that at the

current time P has been closed and that {Pj;,j = 1,2 , z = 0, . .

.

,

7} closes P

if there exist points Pj;, j = 1,2, i = 0, . .

.

,

6, in Z

l U ZV (P ) such that for

each j, j = 1,2, the rays PPj;, i — 1, . .

.

,

6, intersect F3 at the points Pji5 i =

1, . .

.

,

6, respectively, and the ray PPjo intersects Pj_i,i and Ejo at the points

Pj j 7 and Pj0 , respectively. Assuming P has been closed we define dj at the

current time as the smallest value of dmax(P, {P-^ j — 1,2, i = 0,...,7})

for (Pjj, j = 1, 2, i = 0, . .

.

,

7} in the family of sets that close P.

Definitions and terminology for S4 . Let P be a point in 54

. Let Fj, j — 1, 2, 3,

be the three faces of the cube closest to P in the order shown in Figure 8.

Let P0 represent the same face that P3 represents. For each j, j = 1,2,3, let

lj be the ray with origin P perpendicular to Pj. For each j, j = 1,2,3, let

Pj be the point at which lj intersects Pj. Let Pq represent the same point

that Pj represents. Let Pj' be the vertex of the cube common to P4 , P2 , and

20
Figure 6: View of the face closest to a point in S2 that has been closed.

21
m

Figure 7 : View of the two faces closest to a point in 53 that has been closed.

22
/

/

Figure 8: View of the three faces closest to a point in S4 that has been closed.

23
F3. For each j, j = 1,2,3, let rrij be the line that contains the edge of the

cube common to Fj_ l and Fj. For each j, j — 1,2,3, let rrij _ 1 , 1 and mj0 be

the lines through Pj_ 1 and Pj, respectively, perpendicular to mr For each j,

j — 1,2,3, let Ej_i' 1 and Ej0 be the closed half-planes determined by 77ij_i,i

and mjo, respectively, that do not contain P^ , and that are contained in the

planes that contain FJ _ 1 and Fj, respectively. For each j , j = 1,2,3, let H"

be the plane parallel to Fj that contains (P + Pj)/ 2. We define C"{P ) as the

intersection of the closed half-spaces determined by H", H", and H" that

contain P.

Assume P is the currently active point, and Pfi, j = 1,2,3, i = 0,...,5,

are points such that with Pg5 = P^, for each j , j = 1,2,3, P-i?

i = 1,...,4, are points in Fj contained in the four octants around P- as shown in Figure 8, and Pj_ 1

5

and Pj0 are points in Pj- 1,1 and Ej0 , respectively. We

say that at the current time P has been closed and that {Pjiij — 1>2,3,

i = 0, . .

. , 5} closes P if there exist points Pji, j = 1 , 2, 3, i = 0, . .

. , 4, in

U Z

V {P

)

such that for each j, j = 1,2,3, the rays PPji , i = 1,...,4,

intersect Fj at the points Pj-, i = 1,...,4, respectively, and the ray PPjo

intersects Pj_i,i and Pj0 at the points Pj j 5 and P^, respectively. Assuming P has been closed we define d[ at the current time as the smallest value

of dmax(P, {P'i,j = 1,2, 3, i = 0, ... ,5}) for {Pj;, j = 1, 2, 3, i = 0, . .

.

,

5} in

the family of sets that close P.

A modified version of Bowyer’s three-dimensional insertion process [2] is

used in what follows. It is the obvious generalization to three dimensions of

the modified version of Bowyer’s two-dimensional insertion process.

Start of algorithm.

Step 1. Assign points to cells and select first class 1 cell to be activated.

Let M be the floor of N1 ^ 3

. Partition the cube into M3 equal-sized cubic cells.

Determine inner and outer cells.

Assign each point in 5 to a cell.

For each cell, list the points assigned to it. Determine the centermost cell.

If the centermost cell is empty then go to Step 2.

Else designate this cell as the currently active cell.

24
Go to Step 3.

Step 2. Select next inner (class 1 , 2, 3, or 4) cell.

If all inner cells have been activated then go to Step 8.

Else choose the next inner cell to be activated.

If this cell is empty then go to Step 2.

Else designate this cell as the currently active cell.

Determine class for currently active cell.

If class 1 then go to Step 3.

Else go to Step 4.

Step 3. Construct Voronoi polyhedron of a point in S\.

Let P be a point assigned to the currently active cell that has not been

processed.

Designate P as the currently active point.

Start spiral search around P and construct V1

. Update V1 and each subsequent V

1 as appropriate.

For each \n compute Dt = dmax(P.Vt ) and dt

. Terminate search when one of the following criteria is met.

1. 2 Dt < dt.

2. All cells in the cube have been searched.

Upon termination go to Step 7.

Step 4. Begin construction of Voronoi polyhedron of a point in Sn, S3 , or

£4.

Let P be a point assigned to the currently active cell that has not been

processed.

Designate P as the currently active point.

Start spiral search around P and construct V1

. Update V1 and each subsequent Vl as appropriate.

For each \

rt compute Dt = dmax(P,Vt ) and dt

. Terminate search when one of the following criteria is met.

1. 2 Dt < dt.

2. All cells in the cube have been searched.

25
3. P has been closed.

Upon termination, if neither criterion 1 nor criterion 2 has been met

then go to Step 5.

Else go to Step 7.

Step 5. Continue construction of Voronoi polyhedron of a point in S2 , S3 , or S4

. Resume spiral search around P.

Update each V

1 as appropriate.

For each V

1 compute Dt = dmax(P,\n

), dt , and d[.

Terminate search when one of the following criteria is met.

1 . 2 Dt < dt .

2. All cells in the cube have been searched.

3. y/2 d[ < dt

.

Upon termination, if neither criterion 1 nor criterion 2 has been met

then go to Step 6 . Else go to Step 7.

Step 6 . Complete construction of Voronoi polyhedron of a point in So, S3 , or S4

. Determine C"(P).

Resume spiral search around P.

Update each V1 as appropriate.

For each \n compute Dt = dmax(P.Vt D C"(P)) and dt

. Terminate search when one of the following criteria is met.

1 . 2 Dt < dt

.

2. All cells in the cube have been searched.

Upon termination go to Step 7.

Step 7. Save Voronoi polyhedron of a point assigned to an inner cell.

Identify U(F, S) with Vt. Mark P as processed and save V(P, S).

For each P' in U(P, S) that has not been processed let

26
Z

V {P') = Z

v {P')u{P}.

Determine whether currently active cell has been activated.

If activated then go to Step 2.

Else if P is in Si then go to Step 3.

Else go to Step 4.

Step 8. Construct and save Voronoi polyhedra of points in S5. Determine S5

.

If S5 is empty then stop.

Else let Z5 = UPesB Z

v (P).

Perform insertion process on S5 U Z5 . For each P in S5 identify V(P,S) with V

r (P, S5 U Z5 ). For each P in S5 mark P as processed and save V(P, 5).

Stop.

End of algorithm.

Justification of algorithm. Because of similarities with the two dimensional

method, we need not justify that the above algorithm constructs the Voronoi

polyhedra of points in Si or S5. In addition, because of similarities among

52, S3 , and 54 , we only justify that it constructs the Voronoi polyhedra of

points in 52. Let F, l, P'

,

m, H", C”{P ) be as defined above for a point P in 52. Assume

that P is the currently active point and that it has been closed. In addition,

assume that at the current time y/2 d[ < dt and that a set {P/,z = 1, . . .

,

8}

of points in F closes P with d[ = dmax(P, {P/, i = 1 ,

. .

. , 8}). Let V be that

part of V1 that is not contained in C"(P). Assume S \ Z

l U Z

V [P ) U {P} is not empty and O' is a point in this set. Let H' be the plane that perpendicularly bisects O'

P

, and let C' be the open half-space determined by H'

that contains P. We show that C' contains V, so that O' does not affect V, and thus the termination criteria may change from those in Step 5 to those

in Step 6.

Let Pi, i = 1 , . . .

,

8, be points in Z

1 U Z

V (P) such that the rays PP,,

i = 1,...,8, intersect P at the points P-, i = 1,...,8, respectively. We

assume P/ = Pt for each i,i — 1, . .

.

,

8, since, at worst, V is just a larger set. We assume P/ ^ P' , i = 1,. .

.

,8, since, otherwise, V equals the empty set.

Finally, we assume P is not a convex combination of P' and O'.

27
We first argue that we may assume that O' is contained in F

.

For suppose

O' is not in F

.

Let Q'

0 be the point in the plane that contains F such that

O' and Q'0 are equidistant from P, and Q'q lies in the half-plane that contains

Q' and which is determined by a line through P perpendicular to F. Such a O'0 exists because the perpendicular distance from P to F is less than d[ and

dist(Q'

, P) > y/2 d'

t

. Let H'Q be the plane that is the perpendicular bisector

of Q'0 P. Let C'q be the open half-space determined by H'Q that contains P

.

It follows, as shown in Figure 9, that if C'Q contains V, which lies above H'\

then so does C'. Thus, without any loss of generality, we assume that Q'

is contained in F. Further, we assume that P' P[ / P' P2 and that O' is

contained in the region in F obtained by a clockwise rotation along F from

P' P2 to P' P[ (refer back to Figure 6).

Let H[ and H2 be the planes that are the perpendicular bisectors of P[P

and P2 P, respectively. Let R be the region which is the intersection of the

closed half-space determined by H" that does not contain P and the closed

half-spaces determined by H[ and H2 that contain P. Clearly R contains V.

We will show that C' contains V by showing that R is the convex hull of a region (A"") and a ray (

r

"), both of which lie in C'

.

Since C’ is convex, C'

then contains R and the result follows.

The following definitions will be useful. Let H"’ be the plane that contains P

and is parallel to H"

.

Let Q"\ P"\ and P2" be the perpendicular projections

onto H"’ of O', Pj, and P2', respectively. Let d

'0 , d'1? d2 , d'^'

, d"\ and d'" be the

values of dist(Q'

, P), dist(P(, P), dist(P2 , P), dist(Q'", P), dist(P"', P), and

dist(P"'

, P), respectively. Let /q, and l" be the lines that are the intersections of H" with if', H[, and id2 , respectively. Let and V2 be the

lines in H'" that perpendicularly bisect 0'"P, P{"

P

, and P"'P<, respectively.

Let P" be the perpendicular projection of P onto H" . The region K" is defined to be the wedge obtained by intersecting the half-planes in H" determined by l" and l" that contain P" (Figure 10). Another region, K"\ is similarly defined with respect to iP", /2", and P (Figure 11).

We begin the task of showing that C' contains the region K" by arguing

that dq > y/2 dmax(P, {P"'> P”'})-

(O' + P)/2, (P[ + P)/2, (P' + P)/2, (O'" + P)/2, (PI" + P)/2, (P'" + P)/2

are contained in 1'0

', l ", / 2 , /q", V"

,

respectively, so that from the definitions

of 0'", P"', and P2", it follows by similar triangles that and l2 are the

perpendicular projections onto PL'" of /J, /", and /", respectively. If h' is the

28
Figure 9: Perpendicular view of unique plane that contains P, Q'

0 , and Q

'

.

29
/

/

/

Figure 10: Perpendicular view of H"

.

The shaded area is the wedge K .

30
p'" /

/

/

Figure 11: Perpendicular view of H"'. The shaded area is the wedge K"'

.

31
value of dist{P\ P) then, since d'0 > \/2d't , we have

«T ) 2 + (*')

2 = K)2 > 2 «)2 > 2 K)2 = 2(Ki")

2 + W2

) > 2 (d'?y + (h')

2

.

Thus, d"' > y/2 d'"

.

Similarly dg" > \/2d2", which shows

d"' > \/2 dmax(P: {P”' > P"'})-

Since Pf and P2 are in contiguous octants, it follows that I'f does not intersect

K

1 ". In addition, since K’" is the perpendicular projection onto H

1" of if",

it also follows that /J does not intersect K" (Figures 10 and 11). In other

words, by the definition of C' contains K".

Next, we define r" and show it is contained in C’

.

Let H* be the unique plane that contains P, Pj, and P2. Let C* be the

closed half-space determined by H* that contains PL Since, in particular,

dj, > dmax(P, {P[, P2 })? if follows that some convex combination, say T', of

P[ and P2 is also a convex combination of P' and O' with T' ^ O' (Figure 12).

Thus, since P' is an interior point of C*, and T' is contained in if*, it follows

that O' does not belong to C*.

Let /

be the line that is the intersection of the planes H[ and if2

. Let H**

be the plane that contains P and O' and that is perpendicular to H*

. Let

/** be the perpendicular projection onto H** of l". Since l" is perpendicular

to H* (Figure 13), so is /**. Thus, as shown in Figure 14, /** contains a ray

r** that lies wholly in C' D C*

.

Therefore, by the definition of /’*, it follows

that l" contains a ray r" that is also contained in C' D C*

.

Since both K" and r" lie in C', their convex hull R lies in the convex region

C', and our assertion is proved.

Proof of complexity. We now assume that the points in 5 have been chosen from a uniform distribution in the cube. Under this assumption, we

claim that the expected time involved with the above algorithm is O(N)

and 0(N4 ^ 3

) for obtaining all Voronoi polyhedra of points in 5 \ 5s and S5 , respectively. Again, because of similarities with the two-dimensional algorithm, we only present the proof for S5 .

32
Figure 12: Perpendicular view of face of cube closest to P. P{P^ and P'Q'

intersect at the point T

1

.

33
p

;

Figure 13: Perpendicular view of plane H*

. /" is perpendicular to Hthrough

the point T*

.

34
Figure 14 : Perpendicular view of plane H**

. The shaded area is C' H C*.

The ray r*m is the solid portion of the line /**.

35
Let Z5 be as defined in Step 8 of the algorithm. We first show that the

expected number of points in 55 U Z5 is 0(N2 ! 3

). Clearly, the expected number of points in 55 is 0(X2 / 3

). Thus, it suffices to

prove that the expected number of points in Z5 is 0(N2//3

). In what follows, given G, a finite nonempty set in Euclidean space, and

Gi, G2 , nonempty subsets of G, we define X(G, Gi,G2 ) as the number of

points in G2 that are Voronoi neighbors in G of points in Gi. Accordingly,

E(N(G, Gi, G2 )) is defined as the expected value of X(G, Gi, G2 ). The following observation is crucial to our proof:

Let X, Y be finite nonempty sets in Euclidean space. Let X' be a nonempty

subset of X . Then

X(X,X',X) < N(X U V,X',X) + N(X U Y,Y,X).

We prove a stronger result. Let P and P' be points in X and X', respectively.

We show that if P' is a Voronoi neighbor in A

r of P then either P' is a Voronoi

neighbor in A" U Y of P or there exists at least one point in Y which is a

Voronoi neighbor in A" U Y of P. If not, all Voronoi neighbors in X U Y of P

lie in X, none of which is P\ But this implies that the Voronoi polyhedron

of P relative to X coincides with the Voronoi polyhedron of P relative to

X U V, so that P' cannot be a Voronoi neighbor in A" of P, a contradiction.

Let B be the cube from which the set S has been chosen. Let B' be the cube

obtained by surrounding B with L(N) + 1 additional layers of cells. Define

cells in B' \ B as exterior cells. Let Z be the set whose members are the

centroids of the exterior cells.

It follows from the observation above that

N{S, 5S , S) < N(S U Z, S5 , S) + N(S U Z, Z, 5).

Thus, since each term above is a random variable that depends solely on the

choice of 5, we must have

E(N(S, S5 , 5)) < E{N(S U Z, S5 , S)) + E{N(S U Z, Z, 5)).

Let Z' be the set of points in Z contained in the first layer of exterior cells

that surrounds B. From the definition of Z, points in Z' are the only points

in Z that can be Voronoi neighbors in 5 U Z of points in 5. Clearly, points

36
in S5 U Z

are not contained in any of the outermost L(N) layers of cells of

B\ so that the expected number of Voronoi neighbors in 5 U Z of a point in

S5 U Z

is 0(1) [l]. Since the expected number of points in 55 U Z' is 0(iV2//3 ), this implies

E(N{S u Z, S5 , 5)) = E(N(S U Z,Z,S)) = 0{

N

2 ' 3

),

so that

E(N(S, Ss , 5)) < 0(N2/1

) + 0(N2/3) = 0(N2/3

).

Since E(N(S, S5 , 5)) is by definition an upper bound on the expected number

of points in Z5 , our assertion follows.

Finally, we show that the expected time required to perform the insertion

process of Step 8 is bounded above by 0(iV4//3

). Assume, inductively, that the Voronoi diagram for k — 1 points in S5 U has

been constructed. Since the expected number of points in S5 UZ5 is 0(V2 / 3

),

it follows that 0(N2 / 3

) is an upper bound on the expected time required to

find that one of the k — 1 points closest to an additional k

th point, and that

the expected number of vertices of any polyhedron in the Voronoi diagram for

the k — 1 points is at most 0( V

2y/3

). Thus, the expected time required to find

a vertex of the Voronoi diagram for the k — 1 points that does not belong

to the Voronoi diagram for all k points is at most 0(vV2 / 3

). This implies

that the expected time required by the first step of the insertion process, as described in the introduction, is 0(N2//3

) for each point in S5 U Z5 , and

0(iV2/3 ).0(iV2/3

) = 0(iV4/3

)

for all of 55 U Z5. Next, we show that the expected time required by the remainder of the

insertion process is also 0(N4 ^ 3

) for all of S5 U Z5

. As described in the

introduction, and again assuming that the Voronoi diagram for k — 1 points

in S 5 UZ5 has been constructed, the remainder of the insertion process consists

of deleting those vertices in the Voronoi diagram for the k — 1 points that lie

in the Voronoi polyhedron of the k

th point, and adding to the diagram the

vertices of the Voronoi polyhedron of the k

th point. Let r be the number of

points in S5 U and for each fc, 1 < k < r, let d(k ) and a(k) be the number

of vertices that are deleted and added, respectively, when updating for the

k

th point in S5 U Z5 the Voronoi diagram for the first k — 1 points. Clearly,

37
for each k, 1 < k < r, the expected value of a(k) is at most 0(iV2 ^ 3

), so that the expected value of Y7k=i a W mosf 0(N4)/3

). Also, since deleted

vertices come from the set of previously added vertices, we must have

Y.d(k) < Y,a(k),

k — 1 k —

1

and the desired result follows.

5. Computational experience

The two-dimensional algorithm was implemented in standard FORTRAN

on a Control Data Cyber 2051 at the National Bureau of Standards. The

storage of the data structure and the treatment of degenerate vertices were

based on considerations from [2]. Table 1 shows the computing times per

N Proposed algorithm Bowyer’s algorithm

8,100 1.27 X 1(T3 1.07 X 10"3

14,400 1.37 X 10~ 3 1.29 X 10'3

52,900 1.30 X 10~ 3 1.97 X 10" 3

102,400 1.28 X 10"3 2.59 X 10~ 3

122,500 1.28 X 10~ 3 2.80 X 10" 3

Table 1: Computing time per Voronoi polygon.

Voronoi polygon for the algorithm when applied to several regular nontrivial

configurations. The times obtained when applying the modified version of

Bowyer’s algorithm to the same sets are also shown. The unit of computing

time is given in CPU seconds.

Acknowledgments

The author wishes to express his deep gratitude to Dr. Sally E. Howe who

strongly influenced this work and who provided invaluable assistance in re- viewing several versions of the manuscript. He also thanks Dr. Christoph

x Any references to commercial products do not imply recommendation or endorsement

by the National Bureau of Standards of these products.

38
Witzgall for reading the manuscript and offering his criticism and helpful

suggestions. Finally, his gratitude also goes to Dr. Francis Sullivan for his

encouragement, and to Mrs. Rosemary Massengill and Mrs. Ellen Sullivan

for their help in typing the first version of the manuscript.

39
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, 4

(1980), pp. 563-580.

[2] A. Bowyer, Computing Dirichlet tessellations, Comput. J., 24, 2 (1981),

pp. 162- 166 .

[3] M. I. Shamos, Computational geometry, Ph.D. dissertation, Yale Univ.,

New Haven, Conn., 1978.

40
BS.114A (REV. 2-80

-

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR 87-3679

2. Performing Organ. Report No. 3. Publ ication Date

NOVEMBER 1987

4. TITLE AND SUBTITLE

4/3

Expected 0 ( N ) and 0(N ) algorithms for constructing Voronoi diagrams in

two and three dimensions

5. AUTHOR(S)

Javier Bernal

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions) 7. Contract/Grant No.

NATIONAL bureau of standards

DEPARTMENT OF COMMERCE 8. Type of Report & Period Covered

WASHINGTON, D.C. 20234

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

10. SUPPLEMENTARY NOTES

J J Document describes a computer program; SF-185, FIPS Software Summary, is attached.

11. ABSTRACT (A 200~word or less factual summary of most significant information. If document includes a significant

bi bl iography or literature survey, mention it here)

Bentley, Weide and Yao have proposed an expected 0(N) cell technique for computing

Voronoi diagrams in two dimensions that does not generalize readily to three. In

this paper their work is further developed and generalized to produce expected

0 ( N ) and 0 (n 4/3) algorithms for constructing Voronoi diagrams in two and three

dimensions, respecti vely, Computational experience is presented for the

algorithm in two dimensions.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) algorithm; computational geometry; computational complexity; Voronoi diagram;

expected time analysis; Voronoi polyhedra

13. AVAILABILITY

[X] Uni imited

1 j For Official Distribution. Do Not Release to NTIS

i Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

20402.

fXl Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

43

15. Price

$11.95
-

ROMpI