OPTIMIZATION ALGORITHMS FOR COMMUNICATION AND CONTROL
NETWORKS
by
Asaad Jabri
B.S., Metropolitan State College of Denver, 1991
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Computer Science
1999
This thesis for the Master of Science
degree by
Asaad Charles Jabri
has been approved
by
Tom Altman
Mike Radenkovic
Date
ori&fttilman
Jabri, Asaad Charles (M.S., Computer Science)
Optimization Algorithms for Communication and Control Networks
Thesis directed by Professor Tom Altman
ABSTRACT
In this thesis we discuss various optimization algorithms for computer /control
networks whose architectures can be represented via directed or undirected
graphs.
We present one algorithm that assigns users (whose physical coordinates are
determined by a G.S.P. to the closest network node, i.e., a server). We then derive
two important (and original) NPcompleteness results dealing with the design of
network spanning trees. In particular, we show that the problems of finding a
spanning tree of a weighted graph G = (V, E) of weight equal to integer K, as well
as of weight bounded by integers C and D, are both NP complete.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Tom Altman
IV
ACKNOWLEDGEMENT
I would like to thank my thesis supervisor, Professor Tom Altman for his help and
advice to make the thesis possible. I would like to thank the thesis committee,
Professor Mike Radenkovic and Professor Boris Stilman for their efforts to
review my thesis.
CONTENTS
Chapter
1. Introduction.......................................................1
2. Literature Survey..................................................3
2.1 Voronoi Diagrams Literature Survey.................................3
2.1.1 Convex Hull Of a Point Set.........................................4
2.1.2 Finding Nearest Neighbors..........................................8
2.1.3 Planar Subdivision Searching......................................13
2.1.4 Range searching...................................................16
2.1.5 Intersection problems.............................................19
2.2 Spanning Trees literature survey..................................23
3. MinimumCost Spanning Tree Algorithms.............................25
3.1 Kruskals algorithm ..............................................25
3.2 The set Merging Problem...........................................27
3.2.1 Naive SetMerging Algorithm ......................................28
3.2.2 Linked List SetMerging Algorithm .............................. 28
3.2.3 Tree structured set Merging Algorithm Weighted Union
with Path compression ...........................................29
4. NPCompleteness results ..........................................30
4.1 Subset Sum........................................................30
4.2 First Theorem.....................................................30
4.2.1 Proof.............................................................30
4.3 Second Theorem................................................... 31
4.3.1 Proof.............................................................31
5. Sweep line algorithm..............................................32
5.1 Voronoi Algorithm using sweep line technique......................32
6. Conclusion........................................................34
Appendix
A. Source code for Neighborhood Search..... ............................35
References..............................................................45
1. Introduction
During the eighties and nineties, the number of transitions from mainframe
computing to clientserver computing has changed dramatically. This has caused a
change in the architecture of networks, which we now call distributed computing.
During the past few years, the popularity of the web has increased dramatically. The
Internet traffic is dominated by the usage of the web. This explosive Web traffic has
wide ranges of reasons behind it.
First, graphical user interfaces such as browser for navigating the Web are more
abundant than ever, they are very easy to install and use. Internet Explorer,
Netscape Navigator, and NCSA Mosaic are examples of such browsers. E
commerce is growing in popularity, also.
Second, languages and protocols used for information exchange are machine
independent. In cyberspace, there is a big increase in Web pages. On the other hand,
educational institutions and people in the fields of research use the Internet
extensively. Internet traffic has grown significantly, and sparked much research
activities for improvements. The Internet usage has revolutionized the way Earth
does business. Things have been accelerated and compressed. Promotion of
business or educational transactions on the Internet can be achieved by improving
many components that make up the Internet (hardware and software), and by
reducing bottleneck situations, if we want to achieve faster client server response.
First, we can start by improving the performance of the Internet servers. Second, we
can reduce unnecessary Internet traffic by reducing the amount of time a data
packet has to spend in the physical layer of the Internetworks, and the distance it
has to travel (only one data packet can be on a cable at any one time). Networking
devices are able to control amount of traffic and speed up the flow of data. Such
devices could be repeaters, hubs, bridges, and routers. Routers play very important
roles, because they use Layer three of the O.S.I model to determine the best path for
delivery of data, and they help control broadcasts. What is meant by best path is the
shortest path that is available. In order to achieve this goal, routers keep a table
internally that has all the available paths and uses spanning trees to view all the
available routes.
A third way to improve the performance of the client server response time on the
Web is to assign users represented by their computers or workstations to the closest
1
server, in their area. Thus, we decrease the amount of unnecessary traffic, increase
productivity, and data packets avoid costly collisions.
The purpose of this thesis is to first present a discussion of KruskaTs algorithm and
its implementation. An important part of the implementation of Kruksals algorithm
is handling of the Set Merging Problem. Three different techniques for handling the
Set Merging Problem will be presented and analyzed.
Second. we present two computational complexity results addressing problems in
network design and construction. Specifically, we will show that finding spanning
trees, whose cost is restricted via an equality restraint or via upper and lower
bounds, is NPcomplete. This result has significant implications in that it justifies
the use of heuristics or approximation algorithms to address the abovementioned
problems.
Third, we propose an algorithm for solving the closest network node problem (i,e,.
assigns users represented by their computers or workstations to the closest server, in
their area). The algorithm uses the Voronoi diagram, the sweep line approach, and
the run time is 0(n log n).
2
2. Literature survey:
Before we dive into the sea of algorithms methodologies that we will present and
explain, we explore the development of solutions to seven key problems. The
problems considered were chosen on the basis of their longevity, their signifi
cance in Communication and Control Networks, and because of their relationship
to the main thesis topics. In each case, the problem has been actively considered
by both practitioners and theoreticians. The set of problems considered gives an
accurate overview of the problems and methods of computational geometry (in
the plane). In most cases, the solutions are practical and we will discuss some im
plementation issues. Where appropriate, a discussion to higher dimension is in
cluded. It is important to mention that the first five algorithms are related to the
first topic which is the Voronoi diagrams, and the last two are related to the
minimum spanning tree problem.
2.1 Voronoi diagrams
The problems we consider are the computation of convex hulls, the searching of
planar subdivisions, the finding of nearest neighbors, range searching and inter
section detection and computation. We choose these problems for three reasons.
First, they give an accurate representation of the stateoftheart in computational
geometry. Second, the development of solutions to each over the past few years
provides evidence of the success of computational geometry and also provides a
good survey of the techniques used throughout the field. Finally, the problems
are sufficiently simple to be easily studied.
As we describe in greater detail below, much of what is done in computational
geometry arises from the study of methods of sorting and searching in dimensions
higher than one. Sorting and searching are well understood for totally ordered
sets in 1dimension [Kn73]. However, beyond this simple case it is often difficult
to define the exact sorting/searching problems under consideration. When it is
useful in what follows, we will characterize problems as sorting and searching
problems and identify the ordering information that is available to us. This is use
ful in establishing connections among geometric problems that superficially ap
pear dissimilar. These connections are described as a means to provide the reader
with techniques that might be useful in solving his/her next geometric problem.
The planar convex hull problem described below is a good example of a higher
dimensional sorting problem. Next, the problem of planar subdivision searching is
considered. Appropriate planar subdivisions can be considered as convex hulls of
3
3dimensional polyhedra. Here, we are organizing the points on the convex hull
so that the faces of the hull, corresponding to regions in the planar subdivision,
can be efficiently searched. We will see this representation arise again when we
discuss hierarchical representations of convex polygons and polyhedra that can be
used for efficient intersection calculation and searching.
The problems of finding nearest neighbors and planar subdivision searching are
considered next. Nothing in the problem statements suggests that the problems
are related. However, the algorithms proposed to solve the problems are quite
similar. This is a situation where the algorithms proposed unite the problems
rather than vice versa, as is more common. Ultimately, finding nearest neighbors
is best done by building a Voronoi diagram. This Voronoi diagram is a planar
subdivision. Processing of this planar subdivision is done by the general methods
introduced for processing general planar subdivisions. This involves the same
techniques used for processing convex polyhedra in 3 dimensions.
The connections mentioned above are examples of the ideas we intend to develop
below. Where appropriate, implementation details are discussed. Connections
between the problems are described to document how computational geometry is
becoming unified.
We now describe the five problems that form the basis of literature survey for this
thesis. Each demonstrates central ideas of the field in both its statement and solu
tion. The reader is encouraged to pause after reading each problem to ponder rea
sonable methods that might lead to a solution.
2.1.1 Convex Hull Of a Point Set
The input is a set of points P,,P2 P in ^dimensional Euclidean space. The
desired output is the smallest convex set containing all of the points. In general
this' set will be a polyhedron in the space. Specializing to the plane, in the non
degenerate case the solution will be a polygon P with vertices from the original
set. Each vertex of P will be external in some direction. That is, if Pj = (x., y.) is
a vertex of the convex hull, then there exist real numbers a and b such that for all
M
ax, + by, > aXj + byf
This suggests a relationship with linear programming problems where the feasible
region is specified as the intersection of a set of half planes (in the plane or half
spaces in higher dimension). Indeed, as we will see below, this relationship has
4
been exploited in establishing relationships between the convex hull problem and
the half space intersection problem.
Another characterization of the convex hull problem is given for an analog
computing model. We might consider a board corresponding to the plane with
nails hammered in at points of the point set. Of a large rubber band has a circum
ference that includes of all the points; the convex hull is the set of points (or the
polygon) that impede the motion of the rubber band when it is released [Sh74].
Computing convex hulls is the oldest and most popular problem of computational
geometry. Many results on the problem date back to the late 1960s and seem to
have their roots in statistical applications. Connections between convex hulls and
mathematical programming seem to have led to the strategies behind the initial
algorithms. The two earliest published algorithms are due to Graham [Gr71 ] in
dimension 2 and to Chand and Kapur [CK70] in higher dimensions.
At the time, Chand and Kapur were working at Lockheed and the convex hull
problem arose in their work. They describe an algorithm based upon the principle
or giftwrapping for computing convex hulls of spatial point sets. The key idea
of their algorithm is the observation that in any dimension, an edge (of dimension
1) of the convex hull belongs to exactly two facets of the convex hull. The di
mension of the facets will be one less than the dimension of the resulting hull of a
set of 1000 points in 6 dimensions.
A 2 dimensional rendering of their algorithm might go as follows:
Find a point, p, of the convex hull and make it the current point
While you havent returned to p, do
Begin
find the point that is most counterclockwise from p
add that point to the convex hull and make it the current point.
End
The initial point could be the point of minimal y coordinate. The most counter
clockwise point is found by scanning the points and doing tests to see if triples are
counterclockwise. To see whether B is more counterclockwise than C with re
spect to A, we ask whether the triangle ABC has a counterclockwise orientation.
The Graham algorithm was discovered independently and arose in response to a
request from scientists at Bell Labs who were doing the computation for statistical
applications. His algorithm consists of two steps. First, the points are sorted with
respect to an interior point. Next, a Graham scan is performed. During this
step, points are traced about the perimeter of the polygon determining whether a
left turn or right turn was made at each vertex. These correspond to tests to see
5
whether we are moving clockwise or counterclockwise. If the motion is counter
clockwise, the point is tentatively accepted and the scan moves forward. At a
clockwise point, the point is rejected and the scan moves back. Since the scan can
have no more than n3 rejections, the process halts in a linear number of steps.
The behavior of this algorithm on a set of points is shown in Figure 3. As a 2 di
mensional algorithm, it possesses a simplicity that is not present in the gift
wrapping algorithm; however, it does not naturally extend to higher dimensions.
In the years following these algorithms, numerous improvements were suggested
often proposing a special model where ones algorithm was the best. Of the prob
abilistic algorithms and analyses, those of Eddy and Floyd are probably the most
significant ([Ed77, FI76]). In most cases, linear behavior will result and per
formance will be superior to that of a sorting based method. Preparata and Hong
proposed an algorithm for 3 dimension convex hull computations in 1977 [PH77].
This algorithm is based upon divideandconquer method. They compute left and
right hulls and then determine common tangents (in 2 dimensions) and a wrapping
(in 3 dimensions) that can be used to merge the hulls.
Already by 1977, it had been determined that sorting must play a central role in
the derivation of a convex hull algorithm. It was seen that orientation would be
the significant test in the inner loop. Optimal algorithms were known for 2 and 3
dimensions. Open questions remained, however, of these central question asked
for extensions to higher dimensions and for methods that took advantage of the
information they had on input and/or had running times proportional to their input
size. For example, suppose we were given a starshaped polygon, could we com
pute its convex hull in linear time? If we determined in our computation that few
points were on the convex hull, could we compute the hull in faster asymptotic
time? How hard would it be to just identify points of the convex hull (in higher
dimensions) without building all of its structure?
By 1977, most of the important questions about convex hulls had been asked and
many had already been answered. The connections between computing convex
hulls and sorting were already well established. What remained were questions
about optimal convex hull algorithms in cases where the output size was small.
For example, the giftwrapping algorithms in cases where the output size was
small, or the giftwrapping algorithm on random examples. Its running time is
O(nh) where n is the number of points and h is the number of points on the con
vex hull. Typically h is O(log n) and so both the giftwrapping and Graham algo
rithms have cleaner implementation. However, for a set of n points all on the
hull, it gives quadratic performance. The situation is worse in higher dimensions.
Here, there is a lower bound of O(ri^2) (resp. 0(n(d'l)/2) was applied in even (resp.
6
odd) dimensions. This bound comes from worst cases where it could take this
amount of effort to report the output. However, in practice the output size is sig
nificantly smaller. It is known that computing components of the convex hull is
polynomially reducible to linear programming [DR80]. Three key questions
arose. First, if we know something about the input can we use this information to
compute the convex hull. Second, in two dimensions, is it possible to have an al
gorithm with running time optimal in the number of points on the hull even if this
number is not known in advance.
Finally, can lineartime linear programming algorithms [Me83,Me84] and other
techniques be used to determine all components of the convex hull in time linear
in their number. Progress has been made on all of these questions in the past dec
ade. Perhaps the most interesting instance of the first question is the problem of
determining the convex hull of a set of points given as the vertices of a simple
polygon. Here, the issue is that the points have already been sorted in some order.
That is, by virtue of forming a polygon, they have and order. We would be able
to do a Graham scan in linear time if we knew an ordering about an interior point.
However, the order we are given on the boundary may not correspond to such an
ordering. Hence, the Graham scan must be modified to determine when there has
been and inappropriate backup. Here, the algorithms proposed typically compute
an upper hull and a lower hull, keeping track of turns to assure that vertices are
not incorrectly deleted [GY83,Le83]. It remains open to characterize other situa
tions where the input contains sorting information that can be used to simplify the
hull computation.
Next we turn to the situation where not all points are on the convex hull. Here,
the convex hull could potentially be computed in o(n log n) operations if o(n)
points were on the hull. To do so requires not sorting the points since this re
quires 0(n log n) operations. Kirkpatrick and Seidel [KS86] were able to do so.
Their algorithm finds a method of avoiding the recursive technique of finding 2
halves and merging. That technique has a running time which satisfies the recur
rence relation
T(N) = 2 T(N/2) + 0(N).
This method can be slow because it will compute many hull edges that are re
moved in later mergings. Instead, they find a merging edge in linear time (via a
variant of the linear time linear programming method [Dy84,Me83]) and reduce
to the problem of solving two sub problems. This leads to a time bound deter
mined by the recurrence relation
7
T(N,h) = T(N/2,h) + T(N/2,h2) +0(N)
where h1 and h2 represent the number of points on the two smaller hulls and h' +
h2 h. This yields an algorithm of running time 0(n log h) that is aptly titled the
ultimate convex hull algorithm.
In higher dimensions, two methods have evolved to compute convex hulls. One
is an extension of the giftwrapping technique of Chand and Kapur [Sw85]. The
other is a method named beneathbeyond [Se81], Using either method, near op
timal algorithms can be found in the case where the convex hull is as bad as pos
sible. However, the running time is exponential in the dimension even if the out
put is small. In his thesis, Siedel has improved the situation [Se86], By using the
lineartime is linear programming algorithms, he computes convex hulls in di
mension d. If the input consisted of m points and the output had F facets, he can
enumerate all facets in time 0(m2f(d) + d 4 Flogm). Unfortunately, f(d) is the
constant term that arises in the Megiddo linear programming algorithm and
grows exponentially with d. This algorithm is not practical, though it yields in
sights that should lead to practical implementations during the next decade of
computational geometry.
Thus, we see that significant progress has been made on the convex hull problem
since the time of the ChandKapur and Graham results. It is now the case that
theoretically optimal (or nearly optimal) algorithms are known in all cases. Also,
practical algorithms are understood for many situations. There remain a few im
plementation issues that become relevant for applications where convex hull of
large point sets were to be computed or where only part of the structure of the
convex hull is needed. Many of the techniques given in this and the previous
convex hull section find use in real problems of statistics, patter recognition, and
classification algorithms( see e.g., [Ba84]) in describing and classifying point sets.
The convex hull is and especially desirable indicator to be used in recognizing
patterns because of its rotational invariance.
2.1.2 Finding Nearest Neighbors
The input is a set of points PI,P2,...,Pn in ^/dimensional Euclidean space. The
desired output is the nearest neighbor of each point or possibly the pair of points
that are closest. The answer will consist of a directed (possibly disconnected)
graph with an edge from / to j if Pj is the nearest neighbor of P. Note that the
graph is directed since nearest neighbor is not a symmetric relation. Or, the an
swer will consist of the pair of points that realize the relationship
8
min d(Pf P.).
l
A closely related problem asks for furthest neighbors. Here again the output is a
graph. The problem where we request the pair of points that realize the maximum
interpoint distance in the set is known as the diameter problem since this distance
is called the diameter of the set. We show below that the diameter of a set is de
termined by two points that lie on its convex hull. This establishes a potential re
lationship between this problem and the convex hull problem.
In our study of the nearest and furthest neighbors, we study a structure commonly
known as the closest point Voronoi diagram. This is a subdivision of space as
signing to each point P. a region Rt such that points in R. such that points in R. are
closer to P. than to any of the other points. Each of the regions Ri is a convex
polytope. Furthermore, the only unbounded regions correspond to points on the
convex hull of the original set. We will also consider a furthest point Voronoi
diagram. Weve mentioned these structures here because they provide interesting
examples of planar subdivisions for consideration in our next problem.
The next problems to be considered on point sets were computations involving
distance properties. Convex hull algorithms found points that were the extrema of
a set. It was natural to next focus attention on finding pairs of points that realized
extremal distances. The first problems that arose were to find the closest and fur
thest pair of a set of n points.
First, the 1 dimensional versions of the problem were considered. Here, connec
tions could be made with results from Kenneth [Kn73] on sorting and selection.
It was observed that the closest pair of points must be neighbors and that finding
all pairs of neighbors would require information similar to a sort [Sh75]. The fur
thest pair of points would be the maximum and minimum of the set. Identifying
them would require an algorithm for finding the max and min of a set of n num
bers. Thus were derived informal arguments that n log n was asymptotically cor
rect for the closest pair of points and linear time was necessary and sufficient for
the furthest pair.
Next, attention was focused on higher dimensions. In 2 dimensions, more geome
try was necessary in order to establish which points were neighbors and had to be
considered as potential closest pairs in the set. Here, two methods ultimately
emerged for solving the problem.
9
The first method was based on building a structure that would contain all potential
nearest neighbors. This structure assigned to each point a region containing points
closer to it than to any of the others. As such, it could be constructed from the
perpendicular bisectors generated from each pair of points. Under further study, it
became clear that in the plane the structure was actually a planar subdivision. This
followed since the region corresponding to each point was convex and hence a lin
ear number of vertices and edges. Furthermore, it was discovered that this struc
ture could be built by a divideandconquer method in 0(n log n) operations. Fi
nally, it was learned that this structure had been previously defined in other con
texts. Among its other names were Voronoi diagram and Dirichlet tessellation.
From the structure arose the first subquadratic algorithms for finding nearest
neighbors in 2 dimensions. The algorithms consisted of finding the Voronoi dia
gram as a planar graph. Then, all potential nearest neighbors were given by con
sidering edges of the dual of this graph. This dual was known as the Delaunay tri
angulation and had been previously discovered in the 1930s. Since there were
only linearly many edges in the Delaunay triangulation, the minimum distance
could be found in linear time. This method was discovered independently be
ShamodHoey and Reiss at Yale in 1975. The details appear in [SH75]. Reiss im
plemented his algorithm as part of a BLISS program to perform other functions.
However, there was no sense at the time of reasonable data structures for doing
such a task and his implementation was quite complex.
This algorithm seemed to be an overkill to solve so simple a problem and alter
natives were sought. Strong at IBM and Bentley [Be76] discovered a simpler di
vide and conquer method to solve the problem. Their algorithms divided the
points into a left set L and a right set R. For each set, the nearest neighbors were
found along with their distance d(L) and d( R). Next, they observed that a point in
R could only be within distance d(L) of a constant number of points of L. Hence, a
merge of two sets with shortest distances known could be cone by considering only
a linear number of distances and hence linear time would suffice.
The Voronoi diagram method also solves numerous other problems on planar
point sets (see [PS85]). The StrongBentley method extends more naturally to
higher dimensions while the Voronoi diagram on n points can be of size 0(nd2) in
dimension d. However, by 1977 neither method was sufficiently well understood
to lead to a clean implementation. Even for planar point sets, it was not under
stood what data structure ought to be used to implement the Voronoi diagram.
Determining the furthest distance was equivalent to determining the diameter of
the set. The diameter of the set was shown to be the diameter of the convex hull.
This followed from the observation that if one or both of the furthest points did
not lie on the hull, a larger distance was possible. Next, it was shown that poten
10
tial endpoints of a diameter must be antipodal. That is, each edge could be con
sidered to have a range of angles that tangents to it (or one of its endpoints) might
make. When the ranges for 2 vertices involve values differing by n they are con
sidered to be antipodal. Connections with work of Redoes[Er46] established that
there were only linearly many pairs of antipodal points. This yielded the linear
time bound on the scan step of the diameter algorithm.
In 3 dimensions, it was again shown that the diameter would be determined by a
pair of points on the convex hull. It is known that the diameter can be realized at
most a linear number of times[Gr56]. This corresponds to having only a linear
number of antipodal pairs in the plane. However, after a few false starts, the
problem of finding the diameter in fewer than a quadratic number of operations
remained open.
Thus, nearest neighbors required n log n operations by the element uniqueness
lower bounds. N log n algorithms were known in all dimensions by the divide
and conquer method. There was also a method in 2 dimensions based upon the
Voronoi diagram that had numerous other applications. However, no clean
method of implementing the Voronoi diagram or handling degeneracies of point
sets (e.g., 4 cocicular points) was known and there was no accepted data structure
for the problem. Curiously, at about this time, Baumgart was writing a thesis in
computer vision [Ba74] that would provide such a structure. With regards to fur
thest distances, linear time was known to be necessary and sufficient in 1 dimen
sion, n log n was known to be necessary and sufficient in 2 and the problem was
open in higher dimensions.
The rediscovery of the Voronoi diagram provided a framework for discussing
nearest neighbor problems. Furthermore, by 1977, there were known algorithms
requiring 0(n log n) operations for building Voronoi diagrams. While on such
algorithm had been implemented [Re77], implementation was a nontrivial task.
In addition to the problem of finding a data structure upon which to build the im
plementation, there were questions of numerical accuracy and degeneracy. The
divideand conquer nature of all known algorithms was cause for concern in terms
of resolution deficiency and numerical stability of algorithms. What was needed
was a simple yet complete explanation of the Voronoi diagram. There also re
mained the question higher dimensional Voronoi diagrams. However, Brown es
tablished a relationship between convex hulls and Voronoi diagrams [Br79a] that
related convex hull computation in dimension d+1 to Voronoi diagram construc
tion in dimension d. As a result, the results of Seidel [ Se86] given in the previ
ous section can be applied to higher dimensional cases.
11
Significant work was done of the development of good Voronoi diagram algo
rithms during the decade 19771987. Authors sought clean formulations of algo
rithms that would be implementable, extensions to the Voronoi diagram that
would also handle line segments and curves and data structures to use in doing the
computation. This work is presented in two papers [Fo86,GS85] that resolve
many of the outstanding issues and represent the current state of the art.
Guibas and Stolfi[GS85] present a generalization of the wingededge data struc
ture of Boumgart [Ba74] to a quadedge structure that can represent planar subdi
visions and their duals simultaneously. Their work can even represent subdivi
sions on nonorientable manifolds. Their development is complete justifying the
choice of the quadedge as the data structure of choice for the representation of
subdivisions of 2manifolds (ie planar subdivisions and 3dimensional polyhedra).
They are able to give a onepage implementation of the Voronoi diagram ( and its
dual, the Delaunay triangulation) based on their data structure. Since the manipu
lation of their data structure is straightforward, the result is a simple implementa
tion of the Voronoi diagram that is reasonably robust. They can identify a set of
(more than 3) cocircular points efficiently. This identification would be desirable
in constructing Delaunay triangulations since the resulting regular polygon should
not be triangulated.
Fortune [F086] presents the first optimal algorithm for computing the Voronoi
diagram that is not based upon divide and conquer. He uses a sweepline tech
nique. In sweeping across the plane, the Voronoi region of a point will be en
countered before the point is reached. Any sweepline algorithm must have a
method of identifying the regions beginning as it occurs. An algorithm that does
otherwise is destined to have quadratic running time. Fortune introduces a trans
formation of the plane that causes a points Voronoi region to be first encountered
at the point when doing a planar sweep. The result is that he is able to insert the
region at that point, in a manner similar to the plane sweeps of [SH76, B078] de
scribed above. On insertion, he can identify when adjacent regions will change
discontinuously in linear number of events and hence the algorithm has 0(n log n)
running time. His transformation is a conceptual device for understanding the
algorithm. By undoing the transformation, Voronoi regions are identified as they
appear. He gives further transformations for cases where the Voronoi diagrams of
line segments or weighted point sets are to be computed. *Note that this differs
from the sweepline for finding all segment intersections shape or disappear add
ing an event to the event queue. There can only be a where there can be a quad
ratic number of events.
12
The Voronoi diagram has been transformed from a theoretical device used to
solve problems into a tool that is used in practical situations. This progress obvi
ously arises from the development of algorithms that can be used in the practical
situations. Surprisingly, it also arises from the dissemination of some of the early
papers on the subject to engineers in related fields who have problems that need
such a structure. Thus, it seems likely that these applications will derive the next
decade of work on the Voronoi diagram and Delaunay triangulation. A step in
this direction is the thesis work of Mike Laszlo [DL87,La87]. This work provides
an extension of the work of Guibas and Stolfi to spatial subdivisions. There ap
pears to be the potential for the application of this work to numerous problems
involving computational geometry as well as applications to finite applications to
finite element calculations (see e.g.[JB87]).
2.1.3 Planar Subdivision Searching
The input is a subdivision of the plane into regions that are convex polygons that
only overlap in edges and vertices. For example, these polygons may have arisen
from a closest point Voronoi diagram as described above. For convenience, let us
assume that the polygons have all been triangulated so that the input is set of tri
angles Tj, T2,Tn. We will assume further that two triangles overlap in an edge if
the complete edge belongs to both of them.
The problem here is to efficiently search the planar subdivision. Searching a pla
nar subdivision assumes that the subdivision has been sorted. This sorting is done
during a preprocessing phase. The result of the sorting is to create a data structure
that can be searched. The result of sorting is to create a data structure that can be
searched. It is assumed that the sorting is done once and supports numerous
searches. Hence, the cost of the sorting can be amortized over the searches. Al
gorithms for planar subdivision searching have their complexities specified by a
triple of numbers. This triple involves the preprocessing time for the sort, the
space needed to store the searching data structure and the time for each search.
Planar subdivision searching can be viewed as a multidimensional extension of
binary search. In the latter problem, we are given a set of numbers that can be
viewed as a set of points on a line (say the xaxis with the numbers being x coor
dinates). Then, the problem becomes one of organizing the points, which we do
by sorting x coordinates, and searching by standard binary search (see e.g.,
[Kn73]). This operation requires n log n time for the sort, linear space to store the
search data structure and log n time for each search. Just as planar subdivision
searching is an extension into higher dimensions. We will mention below exten
13
sions of the problem into d dimensions and consider the general problem of spa
tial subdivision searching. Our next problem is a different version of the multi
dimensional searching problem.
Planar subdivision by polygons in dimension 2 is the generalization of the subdi
vision of a line by a set of points in dimension 1. There are numerous applica
tions of an algorithm for searching planar subdivisions. For example, the Voronoi
diagram presents applications where we might want to locate the region of a point
in a planar subdivision. Similar problems arise in using finite element analyses to
solve partial differential equations. Also, the problem arises in the postoffice
problem as stated by Knuth. Finally, in higher dimensions, the knapsack prob
lem can be posed as a problem of searching spatial subdivisions. Planar subdivi
sions for applications mentioned here might arise in either of two similar contexts.
Either they arise from the regions formed from the arrangement of a set of n infi
nite lines. Or, they might arise from the Euclidean embedding of a planar graph.
In either case, the structure to be searched is the same. And, we can assume that
the structure can be built once and will remain static for a large number of
searches. Hence, it is reasonable to think of preprocessing to build a usable
search structure that can be rapidly searched.
The central problem here is that of finding an ordering on the regions if the planar
subdivision. The earliest methods [DL74,DL76] did this by reducing the problem
dimension by one using slabbing methods. Here the observation was that the
planar subdivision could be divided into a set of slabs that could be easily
searched. This was done by creating slabs in which the regions were monotoni
cally ordered. Furthermore, the slabs were created in monotone order as well.
Thus, the unknown problem of how to do binary search in two dimensions was
reduced to the problem of doing 2 binary searches. The first was a binary search
through an ordered set of slabs to find the correct one and the second was a search
through an ordered set of regions in the slab to find the correct one. However, the
storage required can be such that a segment has to be stored within every slab.
The slab algorithm resulted in an algorithm requiring 0(n2log n) preprocessing
time and O(n2) preprocessing space on the worst case. However, the retrieval
time was 2 log n accounting for two binary searches of n items each.
Thus, a solution existed. The challenge now became to improve upon the running
times. For such algorithms, complexity is measured via triples of numbers
(PT(n),PS(n),ST(n)) where PT(n) (resp. PS(n)) represents the preprocessing time
(resp. space) and ST(n) represents the search time. It became clear that PT (n) and
PS(n) could be made linear (e.g., by doing nothing), but that if ST(n) were to be
14
sublinear then PT(n) would have to be at least n log n. Then the question arose as
to whether, it was also possible to achieve a solution that was (n log n, n, log ri).
Lee and Preparata [LP79] reduced the bestknown algorithm to in log n, n, log*n).
Finally, in an elegant fashion, Lipton and Tarjan [LT77a,LT77b] achieved a solu
tion of (rt log n, n, log ri). Their methodology was to show that divide and con
quer could be applied to the problem of splitting a planar graph. In linear time
they were able to find a separator that divided the graph into two parts each con
taining a positive fraction of the original vertices. Furthermore, the separator was
small. While their method was not easily implementable, it did provide the first
proof that divide and conquer could be made to work in 2 dimensions as it had in
1. This supported the belief that sorting could be applied to higher dimensional
geometries.
The situation regarding planar subdivision searching in 1977 was similar to that of
Voronoi diagrams. The LiptonTarjan result had solved the theoretical problem
but was not of practical significance. It remained for the practitioners to propose
techniques that would be simpler and hence likely to apply in practice. Two
themes emerged in the development of practical techniques.
The first is the notion of hierarchical search. Kirkpatrick [Ki83] proposed this as
a technique for representing a planar subdivision. The basic idea is that when
searching a planar subdivision, more and more detail are needed about less and
less of the planar subdivision as the search progresses. This suggests a technique
whereby the planar subdivision consisting of a constant number vertices and re
gions. This can be searched in constant time. When the appropriate region has
been identified, only this region of the subdivision is refined by adding vertices
and so further subdividing the region. This is shown in Figure 7. Note that we
only visit one of the subdivided regions.
This process continues for long n steps. Then, in O(log ri) time, the search has
been completed. The hierarchical technique appears to be a basic technique a
computational geometry [DS87] having also found application to polyhedral in
tersection detection [DK85].
The other approach follows [EGS87] the topology of the arrangement and applies
a divide and conquer technique. It is based upon ideas that led to some of the ini
tial methods for solving this problem. The idea is simple. A chain of edges is
found that divides the subdivision in half. The chain also has the property that it
is monotone and hence can be searched quickly. Furthermore, future chains build
upon the initial chain so that all searches can be completed in time O{log ri). A
similar technique is given in [ST86].
15
The progress on planar subdivision searching provides the insights necessary to
lift searching from a onedimensional primitive to a well understood operation in
2 dimensions. The techniques described above have made the problem suffi
ciently practical to allow its use in application systems. Indeed, as graphical edi
tors and window systems become more complex, naive algorithms are demon
strating behaviors that are insufficient and these techniques will be needed. There
remain open the problems of extending searching to spatial subdivisions. The re
sults in [DL74,DL76] stood unchallenged for almost a decade. Recent work [
Ch85a, Cl85] has improved upon the space bounds and gave probabilistic algo
rithms that vastly improve the older results. Extensions beyond dimension 3 re
main open.
The subsequent development of rage searching techniques after quad and kd trees
is interesting within the context of this thesis. As opposed to the other problems
we consider, this is an example where the first algorithm found was the practical
one and led work of theoretical interest. In all other cases, the theoretical work
typically happened before 1977 and the work of the last decade has involved re
fining theoretical insights and looking for simplified algorithms that are imple
mentable. For range searching, the initial algorithms were driven by practical con
siderations and are quite implementable. The algorithms that followed were bom
out of theoretical considerations. Most of these theoretical considerations sought
to fix the bad worst case possibilities of range searching and extending rage
searching to work for nonrectangular search domains.
2.1.4 Range searching
Next we consider a problem that is similar in spirit to the planar/spatial subdivi
sion search. The input is a set of points P,,P2,...,Pn in Jdimensional Euclidean
space. We again want to preprocess the points to allow for efficient searches.
Here the searches involve ranges in the ^dimensional Euclidean space and ask
which of the points lie in the range. Here a range can be viewed as a box, as a
simples or as a halfspace. For example, a range might be described as a set of
ranges ml,m2,...,m(tMIM2,Md where all points P( = {xjI,xj2,...,xj^ with mk
ik
positive halfspace determined by the hyperplane. A hyperplane search asks for
subset of the points that satisfy the hyperplane. A simples search asks for the sub
set of the points that simultaneously satisfy all of a set of hyperplanes.
Two versions of this problem must be distinguished, these are the counting prob
lem and the reporting problem. In the former, the output is a count of all points
that satisfy the conditions (i,e. lie in the range). The latter requires not only for
16
the number of points but also for the names of the points in the range. An algo
rithm for the reporting problem will obviously solve the counting problem. How
ever, there may be a more efficient solution if only a count is desired, this distinc
tion is important in what follows.
In considering this problem in conjunction with the previous problem, it appears
as though a dual transformation should be possible (see e.g., [Br79b]). That is, it
should be possible to transform this problem into one of point location in a spatial
subdivision. In practice, this appears to not be the case as the algorithms devel
oped for the problems have significantly different forms.
The range searching problem has its roots in the implementation of databaselike
queries for determining ranges of values. We might be given a set of attributes
(e.g. salary, age, height and weight) and represent employees of a company as
points in R4 corresponding to their 4tuple of these attributes. Then, we might
probe in order to determine all employees satisfying a set of ranges simultane
ously. This might consist of finding all employees within a certain age range,
having salary between specified limits, ... Usually, the database is reasonably
static so that preprocessing of the points can be done off line in support of faster
queries. Knuth posed the problem of designing a data structure [Kn73] for multi
dimensional searching problems. The search queries to be supported are similar
to those considered in the previous section, however subtle differences lead to al
gorithms of a very different form.
Bentley was the first to propose algorithms of doing rectangular range queries
[Be75,BS75,FB74]. He proposed a structure called the quad tree for doing toe
search in 2 dimensions. Knuth later named the generalization of this search struc
ture the kD tree. The quad tree is relatively simple. The problem of searching
in 2 dimensions is reduced to two methods of searching in 1 dimension just as we
saw in the slab methods of the previous section. However, here the reduction oc
curs in an iterative fashion that is symmetric. The resulting algorithm is one that
has good average case behavior, though its worst case can be bad [LW77].
The quad tree is a 2 dimensional binary search tree that alternates search direc
tions. So, given a set of points in the plane, a divider is found parallel to the y
axis. Next, each of the two regions is divided by a divider parallel to the xaxis.
The alternating process continues until regions contain only 1 point.
Searching of this structure now can be cone by following a path down the binary
search tree corresponding to the subdivision given. If we wish to determine how
many points lie within a rectangle, we determine which of the two regions at a
node the rectangle intersects and use this information to proceed. As opposed to
17
other binary trees, it could be the case that both sub trees need to be followed.
Quad tree based algorithms were already popular within computer graphics and
pattern recognition by 1977[Hu78,HS79]. Indeed, the method of traversing a
quad tree uses techniques similar to those proposed by Wamock [Wa69] for use
in hidden surface removal. However, the discoveries appear to have been inde
pendent.
As opposed to the other problems we discuss in this section, range searching is a
situation where the practical has driven the theoretical. The quad tree was pro
posed as an interesting search structure and remains the data structure of choice in
numerous applications [Sa84]. However, the problem of optimal range searching
has attracted significant theoretical interest. Here, the search problem is the half
space searching problem. So, points are to be preprocessed so that those lying in
a halfspace can be determined efficiently. The search structure is to require only
linear space and the search time is to be sublinear. It has been shown that this
searching model embodies all polyhedral range searches [DE84]. It has also been
shown that with more storage, all searches can be achieved in 0(1og n) time
[EMK82]. Finally, a lower bound of O(nl/3) is known on the problem [Fr81].
The first observation that sublinear halfplane searches were possible is due to
Willard [Wi82] who produced an algorithm of search time O(n774). The method
of Willard consisted of creating 6 regions in the plane each containing an equal
share of the points and such that a line would intersect no more than 4 of them.
Then, a search could be done by reducing a problem to 4 problems of 1/6 of its
original size in constant time. This yields the recurrence relation
T(n) = T(n/2) + T(n/4) +c
with solution 0(ri695).
In higher dimensions, there was a flurry of activity resulting in and upper bound
of O(n899) [DEY84[ in 3 dimensions and an observation that the commonly used
technique would not extend past 4 dimensions {Av85]. Next, Cole [Co85] pro
posed an alternative technique using portions of multiple hyperplanes for the sub
division rather than simply building quadrants or octants as had been cone in di
mensions 2 and 3. He was able to achieve an upper bound of 0(n977) in dimen
sion 4. Next, there followed a construction built upon Coles that a query time of
0(na) [YY85]. Finally, Haussler and Welzl [HW86] showed that a probabilistic
technique could be used to achieve a query time of 0(na). Their construction is
related to the VapnikChervconenkis dimension of probability theory. Curiously,
their construction is nonconstructive. They show that an algorithm exists but
18
cannot demonstrate one. Their enet construction appears likely to have signifi
cant application to related problems of region counting in arrangements [GS87].
A similar construction was discovered independently by Clarkson [Cl 86].
As mentioned above, range searching has proceeded in a different direction from
that taken in the other four problems we discuss. There has been some develop
ment of techniques for using range searching and new applications. However, the
significant advances during the past decade have been theoretical. Elegant
mathematics has resulted from the analysis of extensions to what is basically an
intuitive method. It is likely that these extensions will generate new results of
practical interest during the next decade.
2.1.5 Intersection problems
Our final problem is the problem of polygon intersection computation. We con
sider here a version of the problem for which the input is a set of line segments
s,.s2, ...,sn. There are two versions of the problem. In one, the output is a listing
of all pairs of segments that intersect. For the second, a YES/NO answer is de
sired tilling whether any pair of segments intersect. The analog to the report
ing/counting problems above is obvious Clearly, determining whether any pair of
segments intersect is no harder than enumerating all pairs that do.
Intersection is a basic problem of geometry. This problem forms a basis for many
others where the line segments are replaced by axisoriented rectangles or convex
polygons. The distinction between detection and computation is significant since
for many applications, a fast detector is sufficient.
Initially, intersection problems seem different from the other problems we con
sider here. In particular, there is no obvious sorting or searching in the problem.
But, we shall see otherwise below. Convex hull computations are the analogs of
sorting. In the same way, intersection of convex bodies is the analog of merging
two sorted lists. This will become clear in the next sections. The other lower
bound we consider will be on the element uniqueness problem.
In this problem the input is the set xrx^...,xn of real numbers. The output is a
YES/NO answer corresponding to whether there exist l
We assume a model of computation where queries are linear forms in the inputs
19
(i.e., tests comparing lax. to c for real numbers a /=1, ..n and c). In this case,
it is known [DL79,Be83] that O(n log n) queries are necessary and sufficient.
No study of geometry would be complete without a discussion of intersection
problems. Much of the motivation for intersection comes from a desire to dis
play geometries on a graphics screen. In 2 dimensions, we are often concerned
with a screen full of overlapping windows and our goal is to determine which
windows need to repainted or a paint order for windows. In 3 dimensions, inter
section detection and computation forms the basis of most hidden surface algo
rithms. Both of these problems involve determining intersections among line
segments and identifying which (if any) pairs from a set of line segments inter
sect.
This problem was first considered by Shamos and Hoey [SH76] who developed a
method based upon the notion of plane sweep to trace across a set of segments
and determine whether any pair intersected. Their technique builds on ideas simi
lar to the slabs methods for planar subdivision searching. Its implementation gen
eralizes the implementation of scanline based hidden surface routines as devel
oped by Watkins [Wa70].
The plane sweep uses slabs defined by endpoints of the segments. We consider
these endpoints to be projected on the xaxis and slabs to extend vertically. Seg
ments at any value of x can be ordered by their y coordinates. And, within a slab,
this order can change only if an intersection has occurred. On the boundary be
tween slabs, a new line segment is inserted or deleted. For n line segments, there
a 2n slabs. A very naive algorithm is to determine the slabs and sweep across
slabs determining for each slab if it contains an intersection. This can be done in
time linear in the number of segments within the slab since the order of the seg
ments at the start of the slab is known. This leads to a quadratic algorithm.
Shamos and Hoey improve upon this algorithm by observing that if two segments
intersect, they must have been adjacent in the ordering of segments at some point
(not necessarily at a slab boundary).
1. When the second segment was inserted, they were adjacent in the ordering of
segments
2. When the second segment was inserted, they were not adjacent in the ordering
of segments, but all segments occurring between them either intersect one (or
both) of the segments or were deleted before they intersected.
20
3. When the first segment was deleted, they were adjacent in the ordering of seg
ments.
4. When the first segment was deleted, they were not adjacent in the ordering of
segments, but all segments occurring between them either intersect one (or
both) of the segments or were inserted after they intersected.
As they do their sweep, Shamos and Hoey are careful to always maintain the order
of the line segments. This is done by using a balanced tree scheme so that inser
tions and deletions can be done at a cost of 0(n log ri) operations. Whenever a
segment is about to be inserted, it is checked against its. neighbors for intersection.
And, whenever a segment is bout to be deleted, its neighbors are checked to see if
they intersect each other. Thus, any two segments that will ever be adjacent in
some ordering are checked for intersection. By appealing to the list of situations
given above, it can the be shown that if there is an intersection it will be detected.
Furthermore, the running time of the algorithm is O(n log ri) because the x coordi
nates of the endpoints must be sorted and because 0(h) insertions and deletions
must be done. Finally, this is optimal by appealing to the element uniqueness
problem and showing that a sort must be done.
The only failing of the ShamosHoey algorithm is that it will not necessarily find
all intersections. Bentley and Ottman [B078] proposed an extension that does
compute all intersections. For this algorithm to work, more slabs are necessary.
Every time an intersection is found, an old slab must be subdivided to create a two
new ones. With appropriate modification, it can be shown that the ShamosHoey
algorithm always detects intersections before they occur and therefore, we will
always know in advance that slabs will have to be created. When we reach an in
tersection point corresponding segments that intersected. This requires O(log ri)
operations. Putting this together yields an algorithm with running time 0(h log n
+ I log ri) for reporting all intersections where / is the number of intersections re
ported. In case where I is small, this represents a significant improvement over
previous methods. However, in cases where / is large, this algorithm will actually
be slower than the naive algorithm.
Since the work of Bentley and Ottman, there has been significant progress on in
tersection problem. There have been both theoretical and practical improvements
to their algorithm. Also, numerous other directions in intersection have been pur
sued. Furthermore, work on significant practical problems, most notably hidden
surface removal has begun to occur at both the theoretical and practical levels.
On the basic problem of finding all intersections, improving the upper bound of
0(h log n + I log ri) has been a nagging open problem. Mairson and Stolfi [MS82]
21
showed that if there were two sets of line segments, n painted red and n painted
green such that all intersecting pairs involved a red segment and a green segment,
then O(n log n + 1) and was efficient in practice. Finally, Chazelle [Ch84] im
proved the upper bound to O(n log2n/loglogn + I). Kalin and Ramsey [KR84] did
significant empirical studies of algorithms for the problem and concluded that
Myers was the fastest in practice and that further improvements of his algorithm
were possible. Their conclusion was that for randomly generated line segments, a
technique that wisely combines enumeration and bucketing will give the best be
havior. While, most pictures sent to a hidden surface remover are not random, it
is reasonable to believe that their observations will extend.
Most of the intersection work is geared towards extending methods towards ob
jects of increasing complexity. The holy grail of intersection problems is an algo
rithm that will operate on a database of arbitrary objects and allow for insertions
and deletions of objects as well as determining all objects that a new object inter
sects. The applications of such a technique are manifold. Towards this end, there
has been work on understanding intersections of convex polygons [Sh75], convex
polyhedra [MP78,DK83,DK85,CD87] and curved polygons [DSV86]. Also,
searching techniques have been developed to strike at the problem of a dynamic
database [Ch86]. In all cases, the central problem is the same. One is always
seeking a method of sorting objects that do not have a total order. The goal of
such sorting is to be able to later search a database of such objects to determine
overlap. An interesting technique that gets at the essence of this problem is the
work of Chazelle [Ch85b] on embedding intersection problems in higher dimen
sional spaces and solving them by subdivision searching techniques.
The next step is to relate this work to hidden surface removal problems in a mean
ingful fashion. This has begun with theoretical work [De86] giving lower bounds
for worst case hidden surface removal. Of more immediate potential practical im
portance are results that use theoretical techniques to study key problems. Nota
ble in this category is the work of Hubschman and Zucker on frametoframe co
herence [HZ82] and the forthcoming thesis of Panduranga [Pa87]. The former
paper is a first attempt at using object space techniques to represent an object for
multiple hidden surface removals. The latter considers the problem of ray tracing
in object space by tracing back reflections onto primary objects. As weve seen
above, intersection problems represent an area where theory has applied in prac
tice and the current progress suggests that this indeed continue.
22
2.2 Spanning Trees literature survey:
Input description: A graph G = (V, E) with weighted edges.
Problem description: The subset of E c E of minimum weight forming a tree on
V. The minimum spanning tree of a graph defines the cheapest subset of edges
that keeps the graph in one connected component. Telephone companies and
networking companies are particularly interested in minimum spanning trees,
because the minimum spanning tree of a set of sites defines the wiring scheme
that connects the sites using as little wire as possible. It is the mother of all
network problems.
Minimum spanning trees prove important for several reasons:
They can be computed quickly and easily, and they create a sparse subgraph
that reflects a lot about the original graph.
They provide a way to identify clusters in sets of points. Deleting the long
edges from a minimum spanning tree leaves connected components that
define natural clusters in the data set.
They can be used to give approximate solutions to hard problems such as
Steiner tree and traveling salesman.
As an educational tool, minimum spanning tree algorithms provided graphic
evidence that greedy algorithms can yield provably optimal solutions.
Two classical algorithms efficiently construct minimum spanning trees, namely
Prims and Kruskals. We refer the reader to Chapter five of this thesis for a more
detailed implementation and discussion.
Good expositions on Prims [Pri57] and Kruskals [Kru56] algorithms will appear
in any textbook on algorithms. The fastest implementations of Prims and
Kruskals algorithms use Fibonacci heaps [FT 87].
A recent break through on the minimum spanning tree problem is the lineartime
randomized algorithm of Karger, Klein, and Tarjan [KKT95].
The history of the minimum spanning tree problem dates back at least to Boruvka,
in 1926 and is presented in [GH85].
23
3. MinimumCost Spanning Tree Algorithms
In this chapter, we present a discussion of KruskaFs algorithm and its implemen
tation. An important part of the implementation of KruksaFs algorithm is han
dling of the Set Merging Problem. Three different techniques for handling the Set
Merging Problem, will be presented and analyzed.
3.1 Kruskals algorithm
A spanning tree for a graph G is a subgraph of G that contains every vertex of G
and is a tree. But every connected graph has a spanning tree, and any two span
ning trees for a graph have the same number of edges.
A weighted graph is a graph for which each edge has an associated real number
weight. The sum of the weights of all the edges is the total weight of the graph A
minimal spanning tree for a weighted graph is a spanning tree that has the least
possible total weight compared to all other spanning trees for the graph. If G is a
weighted graph and e is an edge of G then w(e) denotes the weight of e and w(G)
denotes the total weight of G.
The problem of finding a minimal spanning tree for a graph is certainly solvable.
One solution is to list all spanning trees for the graph, compute the total weight of
each, and choose one for which this total is minimal. (Note that the well ordering
principle guarantees the existence of such a minimal total.). This solution, how
ever, is inefficient in its use of computing time because the number of distinct
spanning trees is so large. For instance, a computer graph with n vertices has n(n'2)
spanning trees. For graphs with n vertices and m edges, Kruskals algorithm can
be implemented so as to have worst case orders of mlog1m and n2 respectively.
In other words, KruskaFs algorithm is devised to find a Minimum Cost Spanning
Tree of the Graph, where a spanning tree is an unconnected tree that connects all
of the vertices of the graph. The cost of the spanning tree is the sum of the costs
of the edges of the spanning tree, and a minimum cost spanning tree is simply
(one of) the lowest cost such spanning trees for a given graph.
Therefore, in Kruskals Algorithm, the edges of a weighed graph are examined
one by one in order of increasing weight. At each stage the edge being examined
is added to what will become the minimal spanning tree provided that this addi
tion does not create a circuit. After n1 edges have been added (where n is the
25
number of vertices of the graph), these edges together with the vertices of the
graph form a minimal spanning tree for the graph.
The following is a short description of the algorithm:
KruskaPs Algorithm for obtaining the minimum tree T for graph G(v,e) is as fol
lows:
initialize T to empty
repeat the following for each edge e in order of nondecreasing
weight:
if T + e is acyclic then add e to T
until  T\ = n\
To better understand the algorithm, and to be able to show how KruskaPs works
for a graph, the following is a more detailed version:
Input: G is a weighted graph with n vertices
[build a subgraph T of G and to consist of all the vertices of G with edges
added in
order of increasing weight. At each stage, let m be the number of edges of
TJ
1. Initialize T to have all the vertices of G and no edges.
2. Let E be the set of all edges of G and let m :=
3. [Precondition: G is connected]
while (m
3.a. Find and edge e in E of least weight
3.b. Delete e from E
3.c. if addition of e to the edge set of T does not produce a circuit
then add e to the edged set of T and set m := (m+1)
end while [Postcondition: T is a minimal spanning tree for G],
KruskaPs Algorithm is an example of a greedy algorithm, in that each choice that
is made optimizes the current situation without regard for future choices. It is a
remarkably simple algorithm, in that it minimizes an arithmetic sum, and yet no
addition takes place.
Kruksals is based on the following theorem: For any set of vertices S, any small
est edge incident to S is in some minimum Cost Spanning Tree. This theorem can
be proved by a typical swapping argument.
26
An important step in the implementation of Kruksals is the determination of
whether and edge being evaluated for addition to T would create a cycle. This is
an example of the problem of Set Merging, which is discussed at length below. In
analyzing Kruksals Algorithm, two activities have to be examined. One is the
priority queue that offers the edges in order of nondecreasing weight. The other
is the algorithm used to handle set merging. Below, it will be shown that set
merging can be handled in almost linear time, so it turns out that the dominating
activity is creating and maintaining the priority queue. This requires at most 0(e
log e) time, where e equals the number of edges in the graph, so that is the worst
case timing of Kruskals Algorithm. If the edges were to be given in sorted order,
then the set merging activity would dominate, and Kruskals would be essentially
a linear time algorithm.
Another algorithm for finding a Minimum Cost Spanning Tree is Prims algo
rithm. Prims works by starting at a specific node of the graph and building a tree
out from that node. At each step, the edge of lowest cost incident to the subtree
under construction is added to the subtree. After all vertices are connected, an
minimum spanning tree is obtained.
The major differences between Kruskals and Prims is this: In Kruskals, many
different, unconnected, subtrees are built and united until a single tree spans the
entire graph. In Prims, only a single subtree is created and expanded until tree
spans the graph.
3.2 The set Merging Problem
The Set Merging Problem is as follows: Set S is divided into subsets such that no
subsets overlap, and the union of all of the subsets is equal to S For any element
in S, we need to be able to find out which subset it is in, and for any two subsets
of S, we need to be able to merge them into a single subset.
The process that finds out and returns which subset an element e is in will be
called FIND(e). The process of merging two subsets will be called UNION(a,6);
this is called a destructive merge because the history of the subsets in which an
element resides is lost. A sequence of FINDs and UNIONs must be processed on
line, that is, the FINDs and UNIONs will be executed sequentially, and the re
sults of each FIND and UNION must be determined before the next process is
executed.
For any set S on n elements, there can be at most h1 UNIONs, for at that point
the set will be contained in a single subset. To completely describe the status of
set S would require n FINDs between each union, so in any nontrivial applica
27
tion of Set Merging, there wold be at most 0(n2) Finds in the sequence. In gen
eral, for m FINDs and n UNIONs, n
This Set Merging Process applies to Kruksals Algorithm as follows: As edges are
added to T, subtrees of G are created. The vertices in these subtrees are subset
of the entire set V of vertices. When edge(v, w) is being evaluated to see whether it
will cause a cycle in T, it must be determined in which subsets v and w reside.
That is, FIND(v) and FIND(w) must be executed. If they are in the same subset,
then a cycle would be created, so edge(v, w) is discarded. If they are in different
subsets, then edge(v,w) will be added to T, and, in addition, the two sets contain
ing v and w must be merged into a new set by the process UNION(FIND(v),
FIND(w). The different algorithms for handling the Set Merging Problem will be
presented.
3.2.1 Naive SetMerging Algorithm
The naive setmerging algorithm consists of an array of size n (the number of ver
tices) called VertexSet. The array entry i contains the name of the subset in which
vertex i resides. Thus the process FIND(i) can be done in constant time. The UN
ION^, 6) operation is complete by going through the entire array and changing all
as to b's. Thus a single UNION is 0(n) steps, and because there can be at most
n1 UNIONs, the total time required for all of the UNIONs is 0(n2). So for the
sequence of m FINDs and n1 UNIONs, the total time is 0(m+ n2). If m is close
to n2 then this algorithm is optimal according to asymptotic analysis.
3.2.2 Linked List SetMerging Algorithm
A simple improvement to the naive algorithm is to maintain the subsets as linked
lists. An array VertexSet is still maintained which contains the name of the set for
each vertex, so FIND is still done in constant time. The linked lists for the subset
members can be implemented by using two additional arrays of size n called
FIRST and NEXT. The FIRST array contains the first element of a subset.
The NEXT array contains the element following element / in the subset. When
element UNION instruction does not need to traverse an entire array of size n. To
complete UNION(a, b), then i is the last member of a subset, the entry NEXT(/)
is zero. Now the e linked list for only one of the subsets a or b needs to be trav
ersed and its elements in the array VertexSet are changed. Then the two linked
lists are concatenated to form a single linked list for the new merged subset. An
other significant improvement can be made if the UNION process merges the
smaller set into the larger; This is called weighted union.
28
For example, if the size of set a is one thousand, and the size of set b is two, much
less work is done if b is merged into a and not the other way around. In order to
implement this another array is maintained that keeps track of the sizes of the
various subsets.The analysis of weighted union is as follows: assign the amount of
work of the UNION process to each element whenever it is merged into another
set. Now, by merging a smaller subset into a larger, the resultant subset is at least
twice the size of the smaller subset. Obviously a subset can double in size only
log2 times before it is as big as the entire set, so each element i can be merged at
most log n times. Therefore the entire sequence of UNION operations using
weighted UNION is O(n log n). Thus the entire sequence of m FINDs and n1
UNIONs is 0(m + n log ri), and if m = 0(n log n) this algorithm is optimal.
3.2.3 Tree structured set Merging Algorithm Weighted Union with Path
compression
In this algorithm, the various subsets are represented by a rooted, undirected tree.
The element at the root is the name of the set, and all other members of the set are
children of the root. A UNION(a, b) operation can now be done in constant time
by making the root of one of the subsets a child of the root of the other subset.
This tree structure is implemented by using an array in which the entry for each
element / is the parent of i Thus the FIND( i) instruction simply follows the
path from element i up to the root of the tree containing i. The time required is
dependent on the length of this path, which of course is no longer a constant as it
was in the first two algorithms. In worst case, this path could have length n
1 .However, this path length can be kept "very small by using weighted union and
path compression.
As before, weighted union simply means making the root of the smaller subset a
child of the larger subset. Path compression is accomplished during a FIND op
eration as follows: each vertex on the path from element i to the root of the subset
containing i is made a child of the root. The combined effect of weighted union
and path compression is to reduce the time required for a series of FIND opera
tions to almost linear. The degree to which it is not linear is the inverse of the
Ackerman function, which for all practical purposes is a constant less than five.
Therefore (simplifying somewhat by treating the effects of the inverse Ackerman
function as a constant) the treestructured setmerging algorithm with weighted
union and path compression can process m FINDs and n UNIONs in 0(m + n),
which is essentially linear.
29
4. NPCompleteness results
In this chapter we present two computational complexity results addressing prob
lems in network design and construction. Specifically, we will show that finding
spanning trees, whose cost is restricted via an equality restraint or via upper and
lower bounds, is NPcomplete. This result has significant implications in that it
justifies the use of heuristics or approximation algorithms to address the above
mentioned problems.
First, let us formally define the Subset Sum problem. This problem is NP
complete and it is closely related to one of the original six basic NPcomplete
problems identified by Karp in [Karp72], the Set Partition problem.
4.1 Subset Sum
INSTANCE: Set .4 = {a/, a 2, a3l ..... a} of positive integers and a positive integer
B.
QUESTION: Is there a subset of A, say A', such that the sum of the element val
ues in A is exactly B?
Let G = (V,E) be a weighed graph, with each edge having a nonenegative integer
weight, and K be a positive integer.
4.2 First Theorem
Given G and K, the problem of determining if G has a spanning tree of weight
equal to K is NPcomplete.
4.2.1 Proof
We can nondeterministically guess a spanning tree T, and verify in deterministic
linear time if its weight is equal to K. Hence the problem is clearly in NP. Let us
take an instance of the subset sum problem A= { aj, as, ..... a}, B. We shall
construct a graph G that will have a spanning tree of weight B if and only if the
subset sum problem has A c A of weight B. First, construct one node and label
it as vo For each at in A construct a pair of nodes v,/ and v,2 and connect them
with an edge of weight 0. Next, connect each node v,7 to v<) with an edge of
weight 0. Finally, connect each node vi2 to vo with an edge of weight a,. Now, if
A has a subset A' of weight B then G will have a spanning tree T of the same
weight, by simply choosing those nonzero weight edges (vq, vj that correspond
30
to each a, in A The remaining edges of T would, of course, be the zero weight
edges (v0, vji) for each ay not in A and the nzero weight edges from v,/ to vi2 .
Conversely, if G has a spanning tree of weight B then A of the same weight obvi
ously exists.
4.3 Second Theorem
Given a weighted graph G and two positive integers C and D, the problem of de
termining if G has a spanning tree of weight bounded by C and D is NPcomplete.
4.3.1 Proof
By restriction, follows directly from the first Theorem with C = D = K.
31
5. The Sweep Line technique
When space sweep is used, it reduces an wdimensional static problem to a (n1)
dimensional dynamic problem. Let us say that we have one or more objects lying
in some Cartesian space, and space being completely traversed by a moving hy
perplane, (the sweeping plane). For any time T, the sweeping plane intersects
some subset of the elements, the active ones. These intersections evolve continu
ously in time, except at certain times where discrete events occur: new objects are
added to the active set, old objects get deleted, or the configuration of the inter
sections on the sweeping plane undergoes a major change. The sweep algorithm is
a simulation of this process, is discrete, and uses two data structures: a queue of
future events, and a representation of the active set and its cross section by the
current sweep plane. Each iteration removes the next event from the queue, and
updates the cross section data structure to mirror the effect of the sweeping plane
advancing past that point. Depending on the algorithm, each iteration may also
reveal some future events that were not known at the beginning of the sweep;
these must be inserted into the event queue, in the appropriate order, before pro
ceeding to the next iteration. This technique has been used in algorithms for com
putation of Voronoi diagrams.
5.1 Voronoi Algorithm using sweep line technique
In this section we will describe an optimal sweep algorithm for computing the
Voronoi diagrams of n sites in the plane. For simplicity, we assume the sites are
in general position, so that there are no two sites with same x or y coordinate, and
no four sites are cocircular:
Input: A setP: {pi, .... p} of point sites in the plane
Output: The voronoi diagram Vor(P) given inside a bounding box in a doubly
connected edge list structure
1. Initialize the event Q with all site events
2. While Q is not empty
3. do Consider the event with largest ^coordinate in Q
4. if the event is a site event, occurring at sit /?,
5. then HANDLESITEEVENT(p,)
6. else HANDLECIRCLEEVENT(pl), where pL is the
lowest point of the circle causing the event
32
7. Remove the event from Q
8. The internal nodes still present in T correspond to the half infinite edges of the
voronoi diagram. Compute a bounding box that contains all vertices of the Voro
noi diagram in its interior, and attach the half infinite edges to the bounding box
by updating the doubly connected edge list appropriately
9. Traverse the half edges of the doubly connected edge list to add the cell records
and the pointers to and from them .
HANDLE SET EVENT (p,)
1 .Search in T for the arc A vertically above p,, and delete all circles events involv
ing A from Q.
2. Replace the leaf of T that represents A with a subtree having three leaves: the
middle leaf stores the new site p, and the other two leaves store the site pj that was
originally stored with ^4. Store the tuples [p, ,pj\ and [p7 pi] representing the new
breakpoints at the two new internal nodes. Perform rebalancing operations on T if
necessary.
3. Create new records in the Voronoi diagram structure for the two half edges
separating V(p,) and V(pj) .which will be traced out by the two new breakpoints.
4. Check the triples of consecutive arcs involving one of the three new arcs. In
sert the corresponding circle event only if the circle intersects the sweep line and
the circle event isnt present yet in Q.
HANDLECIRCLEEVENTfc,)
1. Search in T for the arc A vertically above pt that is about to disappear, and de
lete all circle events that involve A from Q.
2. Delete the leaf that represents A fromT. Update the tuples representing the
breakpoints at the internal nodes. Perform rebalancing operations on T if neces
sary.
3. Add the center of the circle causing the event as a vertex record in the Voronoi
diagram structure and create two halfedge records corresponding to the new
breakpoint of the voronoi diagram. Set the pointers between them appropriately.
4. Check the new triples of consecutive arcs that arise because of the disappear
ance of A. Insert the corresponding circle event into Q only if the circle intersects
the sweep line and the circle event is not present yet in Q.
33
6. Conclusion
The literature is growing at such a rate that it was only possible to scratch the sur
face of developments on the seven problems. Beyond these problems, there has
been significant work in a number of other directions that went unmentioned here.
This decade has seen computational geometry pulled in two opposing directions.
On the one hand, the quest for improved theoretical results has caused the field to
make contact with a mathematical community that may provide the tools for these
improvements. In the other direction, the quest for relevance has continued con
tact with researchers in continuous computational geometry, computer graphics,
and solid modeling robotics. The combination of these pulls has led to an
interesting and applicable intellectual discipline. What is most significant about
the improvements is the level of technology transfer. As the results mentioned
here were developed, they became known beyond the theoretical computer
science community from which discrete computational geometry first emerged.
Practitioners were quick to recognize the potential relevance of this field and
crossfertilization was a key theme in the second decade. The momentum of the
field is such that this trend will clearly continue. As the field grows, the
community of scholars involved reaches a size and excitement that neither Robin
Forrest nor Mike Shamos could have realized when they named it. Obviously,
more researchers will be attracted to the field, old problems will be solved, new
problems will emerge, the process of turning theoretical results into practical
algorithms will continue,... as is the nature of intellectual disciplines. Beyond the
obvious, there appears to be one need that is just being recognized and deserves
mention. This is the concern about numerical issues when translating efficient
algorithms such as those mentioned above into practice. There is a need to not
only apply existing techniques but also to identify the limitations of such
techniques and find methods of extending them.
34
Appendix
A. Source code for Neighborhood Search
#include
#include
#include
#define IS Full(ptr) (!(ptr)>
typedef struct point *LIST;
struct point
{
double px;
double py;
LIST next;
};
double x = 0, y = 0;
LIST pointA = NULL, pointB = NULL, shourtest list = NULL;
U* H'***********!'!tll!3l:*>t:=l'*!ft!l'Jl'*** X**************************//
$ sfc $ $& $$$$$$ ae if: $$ liesie if: ^sicjoicitcy^
void MakeList(LIST get)
{
LIST pNewCell;
pNewCell = (LIST) malloc(sizeof(struct point));
pNewCell>px = get>px;
pNewCell>py = get>py;
pNewCell>next = shortest list;
shortest list = pNewCell;
}
35
^***************************************************************^
^^**************************** ******************************yy
U*t * ***** ***** ** ******** * * ****** * * * * ** * *** 3je * * ** * ff.jl
LIST merge(LIST listl, LIST list2)
{
if (listl = = NULL)
return list 2;
else if (list 2 = = NULL)
return listl;
else if (listl >px <= list2>px)
{
listl>next = merge(listl>next, list2);
return listl;
}
else
{
list2~>next = merge(listl, list2>next);
return list2;
}
jj sfe * * 9fe * * * * * $ * $ * * $ afc % * $ $ * * * $ * * * )fc * $ * * * * *
//* *************************** sp[l(* *******************************
j j***************************************************************
LIST split(LIST list)
{
LIST pSecondPoint;
if (list = = NULL)
return NULL;
else if (list^next = = NULL)
return NULL;
else
{
pSecondPoint = list>next;
list> next = pSecondPoint>next;
pSecondPoint>next = split(pSecondPoint>next);
return pSecondPoint;
}
}
36
^****************************************************************
/I****************************MergeSort****************************
//*************************************>***************************
LIST MergeSort(LIST list)
{
LIST SecondList;
If (list = = NULL)
return NULL;
else if (list^next = = NULL)
return list;
else
{
SecondList = split(list);
return merge(MergeSort(list), MergeSort(SecondList));
}
37
jjif. 9(Ã‚Â£9)e>)c9f:9)c3fc3c9^9)ca)c9fc9ff34cdef:9)c9fc3)c3f:9c3f;9)cae9f:>)c9}c9f:f:%sJ<}c9f:9fc9jc9c>fcafÃ‚Â£a(Ã‚Â£3f:)f:34c94c9)c9ic9)c>}ca(;a(cd}:3^9fc9(:df:%94c3tc9f;9fe9tc9f:df:9{c3};
^****3E***?ic$***#$*9ic$$$)e)Ã‚Â£)e}c3t)e$*3t:3c3c9lc$)cj<3c)c:t:3<3'}c:cjc3cjt:jcj:j<9i:3tc3<3f3fe3';c;c;cjc3cj<):)t:fc3f'
void print list(LIST target)
{
printf (\nThe List Contains: );
while (target != Null)
{
printf (\n);
printf ((%g, %g,) \n, target>px, target*>py);
target = target>next;
}
printf (\n\n);
}
U* ****************************iti****>t**>t***:t:!c::ilc*****>:i:***:t:*******
I I* s************************* getpointB *****************************
//****************************************************************
int getpointB()
{
LIST pNewPoint;
int i =0, c = 0 ;
do
{
pNewPoint = (LIST) malloc(sizeof(struct point));
if (IS FULL(pNewPoint))
{
fprintf(stderr, The memory is full! !\n);
exit (1);
}
printf(\nPlease enter a value for x:);
scanf(%lf,&x);
printf(\nPlease enter a value for y:);
scanf(%lf,&y);
pNewPoint>px = x;
pNewPoint>py = y;
pNewPoint>next = pointB;
pointB = pNewPoint;
printf(\nDo you want to enter next value?);
printf(\nPlease enter 0 for no or 1 for yes: );
38
scanf(%d, &c);
i ++;
}while (c != 0);
return i;
!/**** *********************** gefpointA *****************************
int getpointA()
{
LIST pNewPoint;
int i = 0, c = 0;
do
{
pNewPoint = (LIST) malloc(sizeof(struct point));
if (IS FULL(pNewPoint))
{
fprintf (stderr, The memory is full!! \n);
exit(l);
}
printf(\nPlease enter a value for x:);
scanf(%lf,&x);
printf( \ nPlease enter a value for y:):
scanf( %lf &y);
pNewPoint^px = x;
pNewPoint>py = y;
pNewPoint^ next = pointA;
pointA = pNewPoint;
printf( \ nDo you want to enter next value? );
printf( \ nPlease inter 0 for no or 1 for yes : );
scanf(%d &c);
i + +;
}
while (c = 0);
return i;
39
II* **************** *********>(: J^XtfClCh BCk *************** * * ******** *
LIST Extra Check(LIST targetl, LIST current, double d, double m)
{
int hit = 0;
double d of hit = 0;
LIST best = NULL;
while (current! = NULL)
{
if ((targetl>pxm) <=current>px && current>px <= (targetl>px+m)
&&
(targetl >py m) < =current>py && current>py <= (targetl>py +
m))
{
d of hit = sqrt ((current>px targetl>px) (current>px targetl>px
) + (
current>py targetl >py )*( current>py targetl >py));
if (d of hit < d)
{
d = d of hit;
best = current;
hit = 1 ;
}
current = currents next;
}
else
{
if (hit)
return (best);
else
return NULL;
}
}
return NULL ;
40
^y****************************************************************
//* ******************* Neighborhood Search ************************
^y* ***************************************************************
void Neighborhood Search(LIST target)
{
int i = 1, hit = 0 ;
double d of hit = 0 d of shortest = 1000 ;
LIST ptemp, shortest node, update ;
Ptemp = pointA;
While (ptemp != NULL)
{
if ((target>px i)<=ptemp>px && ptemp>px<= (target>px+i) &&
(target>pyi)<=ptemp>py && ptemp>py<= (target>py+i))
{
d of hit = sqrt ((ptemp>px target>px)* (ptemp>px target>px) +
(ptemp>py target ^ py) (ptemp> py target>py));
if (d of hit < d of shortest)
{
d of shortest = d of hit;
printf (= = % g = = \n, d of shortest);
shortest node = ptemp;
hit = 1;
}
ptemp = ptemp> next;
}
else
{
if (! hit)
i++;
else
{
if ((update = Extra Check(target, ptemp, d of shortest, i*sqrt (2))) =
Null)
shortest node = update;
break;
}
}
}
MakeList(shortest node);
}
41
y/y^ ^ 5^ 5(c >}c 5}e 5js s)c 9(c s}c sft sfs sjs s(e sc >fc sjc j}c s^c 4c s4c a{c y}^f*^"^ Jl4*lf*i*s^*l{5fss^HCJi{Sis*l{*l!5l{,l*5i{5lejls9{Hs*ls^t:*i!Sie*ls*i{*ff,lS5ls,i*>i{,KH**ll!,i<,ie}i{JiS3fc
II* Jje **** ac *** ** 3(c ** aje * * * * **** ************ * *************** *
LIST invert(LIST lead)
{
LIST middle, trail;
middle = NULL;
while (lead = NULL)
{
trail = middle;
middle = lead;
lead = lead^next;
middle>next = trail;
}
return middle;
}
^****************************************************************
*
!I$t ************5fc*5fc***:t: ^ ^5j4^tHe3(E^t^SS((^<)j!5jC>C!)t)<5jÃ‚Â£4')()(c3jcHe^^(%4'^*^5l*^*^*
*
^y****************************************************************
*
void Release Memory(LIST list)
{
LIST dump;
while ((dump = list) != NULL)
{
list = list^next;
free(dump);
}
}
42
jI*Is sfc sfe sjc 3jc sjc aj< sfs 9fe sc sc 9}tf 9}c sfc sjc 9^s 3c ?)c s{e sjc 9fc sjc *1* sjc s]c sc sfc dfc sfc sfe sfe ?fe s)e 9)c sf< 3fc 9fs a]c )fe ac afc ){c sj< sjc ae s{c sfc s{c d)c 9}e
jjif. ^J^*************************************************************
void main ()
{
int numA = 0, numB = 0, index = 0, j = 0;
LIST current;
/* input and print how many points in the set A */
numA = getpointA ();
printf ( \nThere are % d points in the A set. \n, numA);
/* input and print how many points in the set B */
numB = getpointB();
printf (\nThere are % d points in the B set. \n, numB);
/* Sort the points in set A & set B */
pointA = MergeSort(pointA);
pointB = MergeSort(pointB);
/* Search the nearest neighbors from set A of all points of set B */
current = pointB;
while (current != NULL)
{
Neighborhood Search(current);
current = current^ next;
}
shortest list = invert(shortest list);
/ Output of each points in set B & its nearest neighbor /
while (pointB != NULL && shortest list != NULL)
{
printf (B SET (%g, %g)...,pointB>px, pointB>py);
pointB = pointB>next;
printf (NEAREST (%g, %g,)...\n, shortest list>px, shortest
list>py);
shortest list = shortest list>next;
}
/ Garbage Collection /
Release Memory (pointA);
Release Memory (pointB);
43
Release Memory (shortest list);
}
44
References
[Av85]
[Ba84]
[Ba84
[Be83]
[Be75]
[B078]
[BS75]
[Be76]
Avis, K., On the partitionability of point sets in space, ACM
symposium on computational geometry, 1985, 116120
Baird, H. Modelbased image matching using location, Ph.D.
Thesis, Princeton University, May, 1984.
Baumagart, B. G., Geometric modeling for computer vision,
Stanford U., Computer Science Dept. Technical Report, STAN
CS74463.
BenOr, M,. Lower bounds for algebraic computation trees,
Proceedings of the 15th ACM Symposium on the Theory of
computing, 8086, 1983.
Bentley, J. L., Multidimensional binary search trees used for
associative searching, Communications of the ACM, 18, 509517,
1975.
Bentley, J.L. and Ottman, Th., Algorithms for reporting and
counting geometric intersections, CarnegieMellon University,
August 1978.
Bentley, J.L. and Stanat, D.F., Analysis of range searches in quad
trees, Information processing Letters, 3, 170173, 1975.
Bentley, J.L., Divide and conquer algorithms for closest point
45
[Br79a]
[Br79b]
[CK70]
[Ch84]
[Ch85a]
[Ch86b]
[Ch86]
[CD87]
[C185]
problems in multidimensional spaces, PhD thesis, University of
North Carolina, 1976.
Brown, K. Q., Voronoi diagrams from convex hulls, Information
Processing Letters,.9, 223228,1979.
Brown, K. Q., Geometric transformations for fast geometric
algorithms, PhD thesis, Carnegie Mellon University, December,
1979.
Chand, D.R. and Kapur, S.S., An algorithm for convex polytopes,
JACM, 17, 7886, 1970.
Chazelle, B., Intersecting is easier than sorting, Proceedings of the
16th ACM Symposium on the Theory of computing, 125134, 1984
also JCSS, .
Chazelle, B., How to search in history, Information control, 64,
7799, 1985.
Chazelle, B., Fast searching in a real algebraic manifold with
applications to geometric complexity, Proc. CAAP85, LNCS,
SpringerVerlag, 1985.
Chazelle, B., Filtering search: a new approach to query
answering, SIAMJ. on Computing, 15, 703722, 1986.
Chazelle, B. and Dobkin, D., Intersection of convex objects, JACM,
34, 127,1987.
Clarkson, K., A probabilistic algorithm for the post office problem,
46
[Cl 86]
[Co85]
[De86]
[DE84]
[DEY84]
[DK83]
[DK85]
[DL87]
Proceedings of the 17th ACM Symposium on the Theory of
computing, 175185, 1985.
Clarkson, K., Further applications of random sampling to
computational geometry, Proceedings of the 18th ACM
Symposium on the Theory of Computing, 414423, 1986.
Cole, R., Partitioning point sets in 4 dimensions, ICALP 85,
Springer Verlag Lecture Notes in Computer Science, 194, 111
119,1985.
Devai, F., Quadratic bounds for hidden line elimination, ACM
symposium on computational geometry, 269275, 1986.
Dobkin, D.P. and Edelsbrunner. H.E., Space searching for
intersecting objects, IEEE Symposium on foundations of computer
science, 387392, 1984, J. of Algorithms, .
Dobkin, D.P., Edelsbrunner, H.E., and Yao, F.F., A 3space
partition and its applications, unpublished manuscript.
Dobkin, D.P. and Kirkpatrick, D.G., Fast detection of polyhedral
intersections, Theoretical computer science, 27, 241352, 1983.
Dobkin, D.P. and Kirkpatrick, D.G., Alinear slgorithm for
determining the separation of convex polyhedra, J. of Algorithms,
6,381392, 1985.
Dobkin, D.P. and Laszlo, M.J., Primitives for the manipulation of
threedimensional subdivisions, ACM Symposium on the Theory
47
[DL74]
[DL76]
[DL79]
[DR80]
[DS87]
[DSV86]
[Dy84]
[Ed77]
[EGS87]
[EMK82]
of Computing, geometry, 1987, .
Dobkin, D.P. and Lipton, R.J., On some generalizations of binary
search, Proceedings of the 6th ACM Symposium on Theory of
Computing, 310316, 1974.
Dobkin, D.P. and Lipton, R.J., Multidimensional searching, SIAM
Journal on computing, 5, 181186, 1976.
Dobkin, D.P. and Lipton, R.J., On the complexity of computations
under varying sets of primitives, JCSS, 18, 1(1979),8691.
Dobkin, D.P. and Reiss, S.P., On the complexity of linear
programming, Theoretical Computer Science, 11,3, 118, 1980.
Dobkin, D.P. and Souvaine, D., Computational geometry a
users guide, in Advances in Robotics, Vol. 1. Edited by J. T.
Schwartz and C.K. Yap, Lawrence Erlbaum Associates, 1987.
Dobkin, D.p., Souvaine, D., and Van Wyk, C., Decomposition and
intersection f simple splinegons, submitted for publication.
Dyer, M., Linear time algorithms for two and three variable linear
programs, SIAM Journal on Computing, 13, 3145, 1984.
Eddy, W. Anew convex hull algorithm for planar sets, ACM
Transactions on Mathematical Software, 3,4, 398403, 1977.
Edelsbrunner, H.E., Guibas, L.J. and Stolfi, J., Optimal point
location in a monotone subdivision, SIAM J on Computing, .
Edelsbrunner, H., Maurer, H. and Kirkpatrick, D., Polygonal
48
intersection searching, Information Processing Letters, 14, 7479,
1982.
[EW82] Edelsbrunner, H. and Welzl, E., Halfplane range search in linear space and O(n0.695) query time, unpublished manuscript.
[Er46] Erdos, P., On sets of distances of n points, American Mathematic! Monthly, 53, 248250, 1946.
[FP81] Faux, I.E. and Pratt, M.J., Computational geometry for design and manufacture, Ellis Horwood Ltd., Chichester, 1981.
[FB74] Finkel, R.A. and Bentley, J.L., Quad trees: a data structure for retrieval on composite keys, Acta informatica, 4, 19,1974.
[FI 76] Floyd, R., private communication.
[Fo68] Forrest, A.R., Curves and surfaces for computeraided design, PhD thesis, Cambridge University, July, 1968.
[Fo86] Fortune, S., A sweepline algorithm for Voronoi diagrams, ACM symposium on computational geometry, 313322, 1986.
[Fr81] Fredman, M.L., Lower bounds on the complexity of some optimal data structures, SIAM Jon computing, 10,110,1981.
[Gr71] Graham, R.L., An efficient algorithm for determining the convex hull of a finite planar set, Information Processing Letters, 1, 132133, 1972.
[Gr83] Graham, R.L. and Yao, F.F., Finding the convex hull of a simple polygon, J. of Alogrithms, 4, 324331, 1983.
49
[Gr56]
[GS87]
[GS87]
[HW86]
[HZ82]
[Hu78]
[HS79]
[JB87]
[KR84]
Grunbaum, B., Aproofof Vaszonyis conjecture, Bulletin Res,
Council of Israel, 6, 7778, 1956.
Guibas, L.J. and Stolfi, J., Primitives for the manipulation of
general subdivisions and the computation of Voronoi diagrams,
ACM Trans. On Graphics, 2, 75123, 1985.
Guibas, L.J. and Sharir, M., private communication.
Haussler, K. and Welzl, E., Epsilonnets and simplex range
queries, ACM symposium on computational geometry,
1986, 6171, Discrete and Computational Geometry, .
Hubschman, H. and Zucker, S., Frametoframe coherence
and the hidden surface computation: constraints for a convex
world, ACM Transactions on Graphics, 1,2, April, 1982, 129
162.
Hunter, G.M., Efficient computation and data structures for
graphics, PhD thesis, Princeton University, 1978.
Hunter, G.M. and Steiglitz, K., Operations on images using
quad trees, IEEE Transactions on Pattern Analysis and
Machine Intelligence, 1, 145153, 1979.
Jameson, A. and Baker, T., Improvements to the aircraft Euler
method, AIAA 25th Aerospace sciences meeting, paper AIAA
870452, 1987.
Kalin, R. and Ramsey, N., A new analysis and comparison
50
of algorithms for the planar segment intersection problem,
unpublished manuscript.
[Ka72]
[Ki83]
[KS86]
[Kn73]
[La87]
[Ke83]
[LP79]
Karp, R.M. [1972], Reducibility among combinatorial problems,
in R.E. Miller and J.W. Thatcher (eds.), Complexity of
Computer Computations, Plenum Press, New York, 85103.
(1.5; 3.1; 5.2; 7.1; 7.4; Al.l; A1.2; A1.3; A2.1; A2.2; A3.1;
A3.2; A5.1; A6; A10.1)
Kirkpatrick, D.G., Optimal search in planar subdivisions,
SIAMJon Computing, 12,2835, 1983.
Kirkpatrick, D.G. and Seidel, R., The ultimate convex hull
algorithm?, S1AMJon computing, 15,287299, 1986.
Knuth, D.E., The art of computer programming, Vol. 3 
Sorting and Searching, AddisonWesley, Reading, MA, 1973.
Laszlo, M.J., Manipulating threedimensional subdivisions, PhD
thesis, Princeton University, 1987, to appear.
Lee, D.T., On finding the convex hull of a simple polygon,
International Journal of Computing and Information Sciences,
12,8798, 1983.
Lee, D.T. and Preparata, F.P., Location of a point in a planar
subdivision and its applications, SIAM Journal on Computing,
6,3,594606, 1979.
51
[LW77]
[LT77a]
[LT77b]
[MS82]
[Me83]
[Me84]
[MP78]
[My85]
Lee, D.T. and Wong, C.K., Worstcase analysis for region and
partial region searches in multidimensional binary trees and
balanced quad trees, Acta Informatica, 9, 2329, 1977.
Lipton, R.J. and Tarjan, R.E., A separator theorem for
planar graphs, Conference on theoretical computer
science, 110,1977.
Lipton, R.J. and Tarjan, R.E., Applications of a planar separator
theorem, IEEE Symposium on foundations of computer science
162170, 1977.
Mairson, H.G. and Stolfi, J., unpublished manuscript.
Megiddo, N., Linear time algorithm for linear programming in
R3 and related problems, SIAM J on Computing, 12,759776,
1983.
Megiddo, N., Linear programming in linear time when the
dimension is fixed, JACM, 31, 114127, 1984.
Muller, D.E. and Preparata, F,P., Finding the intersection of
two conves polyhedra, Theoretical Computer Science, 7,217
236, 1978.
Myers, W.W., An 0(E log E + I) expected time algorithm for
the planar segment intersection problem, SIAM J on,
Computing, 14 625637, 1985.
52
[Pa87]
Panduranga, E.S., Reflections in curved surfaces, PhD Thesis,
[PH77]
[PS85]
[Re77]
[Ro85]
[Sa84]
[ST86]
[Se81]
[Se86]
Princeton U.,.
Preparata, F.P. and Hong, S.J., Convex hulls of finite sets of
points in two and three dimensions, Communications of the
ACM, 20, 2, 8793, 1977.
Preparata, F.P. and Shamos, M.I., Computational Geometry
SpringerVerlag, 1985.
Reiss, S.P., private communication.
Rosenstiehl, P., private communication.
Samet, H., The quadtree and related hierarvhical data structures,
ACM Computing Surveys, 16, 187260, 1984.
Samak, N. and Tarjan, R.E., Planar point location using
persistent search trees, Communications of the ACM, 29, 669
679, 1986.
Seidel, R., A convex hull algorithm optimal for point sets in even
dimensions, Technical Report, University of British Columbia,
1981.
Seidel, R., Outputsize sensitive algorithms for constructive
problems in computational geometry, PhD thesis, Cornell
University, 1986.
53
[Sh74]
[Sh75]
[SH75]
[Sw85]
[Wa69]
[Wa70]
[Wi82]
[YY85]
Shamos, M.I., private communication.
Shamos, M.I., Geometric complexity, Proceedings of the 7th
ACM Symposium on the Theory of Computing, 224233, 1975.
Shamos, M.I. and Hoey, D., Geometric intersection problems,
IEEE Symposium on foundations of computer science, 208
215,1976.
Swart, G., Finding the convex hull facet by facet, J. of ,
Algorithm 6, 1748, 1985.
Wamock, J.E., A hiddensurface algorthim for computer
generated halftone pictures, University of Utah Computer
Science Department, TR 415, 1969.
Watloms. G.S., A realtime visible surface algorithm,
University of Utah Computer Science Department, UTECCSc
70101, June, 1970.
(
Wilard, D., Polygon retrieval, SIAM Journal on Computing,
11, 149165,1982.
Yao, A. and Yao, F., A general approach to ddimensional
geometric queries, Proceedings of the 17th ACM Symposium
on the Theory of Computing, 163169, 1985.
54

PAGE 1
OPTIMIZATION ALGORITHMS FOR COMMUNICATION AND CONTROL NETWORKS by Asaad Jabri B.S., Metropolitan State College ofDenver, 1991 A thesis submitted to the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Computer Science 1999
PAGE 2
This thesis for the Master of Science degree by Asaad Charles Jabri has been approved by TomAttman Mike Radenkovic /2 /1o( Cf 'f Date
PAGE 3
Jabri, Asaad Charles (M.S., Computer Science) Optimization Algorithms for Communication and Control Networks Thesis directed by Professor Tom Altman ABSTRACT In this thesis we discuss various optimization algorithms for computer /control networks whose architectures can be represented via directed or undirected graphs. We present one algorithm that assigns users (whose physical coordinates are determined by a G.S.P. to the closest network node, i.e., a server). We then derive two important (and original) NPcompleteness results dealing with the design of network spanning trees. In particular, we show that the problems of finding a spanning tree of a weighted graph G = ( V. E) of weight equal to integer K, as well as of weight bounded by integers C and D, are both NPcomplete. This abstract accurately represents the content of the candidate's thesis. I recommend its publication. Signed Tom Altman lV
PAGE 4
ACKNOWLEDGEMENT I would like to thank my thesis supervisor Professor Tom Altman for his help and advice to make the thesis possible. I would like to thank the thesis committee, Professor Mike Radenkovic and Professor Boris Stilman for their efforts to review my thesis.
PAGE 5
CONTENTS Chapter I. Introduction ........................................................................ I 2. Literature Survey .................................................................. 3 2.1 Voronoi Diagrams Literature Survey ........................................... 3 2.1.1 Convex Hull Of a Point Set ...................................................... 4 2.1.2 Finding Nearest Neighbors ....................................................... 8 2.1.3 Planar Subdivision Searching ................................................... 13 2.1.4 Range searching .................................................................. 16 2.1.5 Intersection problems ............................................................ 19 2.2 Spanning Trees literature survey ................................................ 23 3. MinimumCost Spanning Tree Algorithms ................................... 25 3.1 Kruskal's algorithm ............................................................. 25 3.2 The set Merging Problem ....................................................... 27 3.2.1 Naive SetMerging Algorithm ................................................. 28 3.2.2 Linked List SetMerging Algorithm ..... : .................................... 28 3.2.3 Tree structured set Merging AlgorithmWeighted Union with Path compression .......................................................... 29 4. NPCompleteness results ....................................................... 30 4.1 Subset Sum ........................................................................ 30 4.2 First Theorem ..................................................................... 30 4.2.1 Proof ............................................................................... 30 4.3 Second Theorem .................................................................. 31 4.3.1 Proof ............................................................................... 31 5. Sweep line algorithm ............................................................ 32 5.1 Voronoi Algorithm using sweep line technique ............................. .32 6. Conclusion ........................................................................ 34 Appendix A. S9urce code for Neighborhood Search ........ ...................................... 35 References ................................................................................... 45 v
PAGE 6
1. Introduction During the eighties and nineties, the number of transitions from mainframe computing to clientserver computing has changed dramatically. This has caused a change in the architecture of networks, which we now call distributed computing. During the past few years, the popularity of the web has increased dramatically. The Internet traffic is dominated by the usage of the web. This explosive Web traffic has wide ranges of reasons behind it. First, graphical user interfaces such as browser for navigating the Web are more abundant than ever, they are very easy to install and use. Internet Explorer, Netscape Navigator, and NCSA Mosaic are examples of such browsers. commerce is growing in popularity, also. Second, languages and protocols used for information exchange are machine independent. In cyberspace, there is a big increase in Web pages. On the other hand, educational institutions and people in the fields of research use the Internet extensively. Internet traffic has grown significantly, and sparked much research activities for improvements. The Internet usage has revolutionized the way Earth does business. Things have been accelerated and compressed. Promotion of business or educational transactions on the Internet can be achieved by improving many components that make up the Internet (hardware and software), and by reducing bottleneck situations, if we want to achieve faster client server response. First, we can start by improving the performance of the Internet servers. Second, we can reduce unnecessary Internet traffic by reducing the amount of time a data packet has to spend in the physical layer of the Internetworks, and the distance it has to travel (only one data packet can be on a cable at any one time). Networking devices are able to control amount of traffic and speed up the flow of data. Such devices could be repeaters, hubs, bridges, and routers. Routers play very important roles, because they use Layer three of the O.S.I model to determine the best path for delivery of data, and they help control broadcasts. What is meant by best path is the shortest path that is available. In order to achieve this goal, routers keep a table internally that has all the available paths and uses spanning trees to view all the available routes. A third way to improve the performance of the client server response time on the Web is to assign users represented by their computers or workstations to the closest
PAGE 7
server, in their area. Thus, we decrease the amount of unnecessary traffic, increase productivity, and data packets avoid costly collisions. The purpose of this thesis is to first present a discussion ofK.ruskal s algorithm and its implementation. An important part of the implementation ofKruksal's algorithm is handling of the Set Merging Problem. Three different techniques for handling the Set Merging Problem will be presented and analyzed. Second, we present two computational complexity results addressing problems in network design and construction. Specifically, we will show that finding spanning trees, whose cost is restricted via an equality restraint or via upper and lower bounds, is NPcomplete. This result has significant implications in that it justifies the use of heuristics or approximation algorithms to address the abovementioned problems. Third, we propose an algorithm for solving the closest network node problem (i,e,. assigns users represented by their computers or workstations to the closest server, in their area). The algorithm uses the Voronoi diagram, the sweepline approach, and the run time is O(n log n). 2
PAGE 8
2. Literature survey: Before we dive into the sea of algorithms methodologies that we will present and explain, we explore the development of solutions to seven key problems. The problems considered were chosen on the basis of their longevity, their signifi cance in Communication and Control Networks, and because of their retationship to the main thesis topics. In each case, the problem has been actively considered by both practitioners and theoreticians. The set of problems considered gives an accurate overview of the problems and methods of computational geometry (in the plane). In most cases, the solutions are practical and we will discuss some im plementation issues. Where appropriate, a discussion to higher dimension is in cluded. It is important to mention that the first five algorithms are related to the first topic which is the Voronoi diagrams, and the last two are related to the minimum spanning tree problem. 2.1 Voronoi diagrams The problems we consider are the computation of convex hulls, the searching of planar subdivisions, the finding of nearest neighbors, range searching and inter section detection and computation. We choose these problems for three reasons. First, they give an accurate representation of the stateoftheart in computational geometry. Second, the development of solutions to each over the past few years provides evidence of the success of computational geometry and also provides a good survey of the techniques used throughout the field. Finally, the problems are sufficiently simple to be easily studied. As we describe in greater detail below, much of what is done in computational geometry arises from the study of methods of sorting and searching in dimensions nigher than one. Sorting and searching are well understood for totally ordered sets in !dimension [Kn73]. However, beyond this simple case it is often difficult to define the exact sorting/searching problems under consideration. When it is useful in what follows, we will characterize problems as sorting and searching problems and identify the ordering information that is available to us. This is useful in establishing connections among geometric problems that superficially ap pear dissimilar. These connections are described as a means to provide the reader with techniques that might be useful in solving his/her next geometric problem. The planar convex hull problem described below is a good example of a higher dimensional sorting problem. Next, the problem of planar subdivision searching is considered. Appropriate planar subdivisions can be considered as convex hulls of 3
PAGE 9
3dimensional polyhedra. Here, we are organizing the points on the convex hull so that the faces of the hull, corresponding to regions in the planar subdivision, can be efficiently searched. We will see this representation arise again when we discuss hierarchical representations of convex polygons and polyhedra that can be used for efficient intersection calculation and searching. The problems of finding nearest neighbors and planar subdivision searching are considered next. Nothing in the problem statements suggests that the problems are related. However, the algorithms proposed to solve the problems are quite similar. This is a situation where the algorithms proposed unite the problems rather than vice versa, as is more common. Ultimately, finding nearest neighbors is best done by building a Voronoi diagram. This Voronoi diagram is a planar subdivision. Processing of this planar subdivision is done by the general methods introduced for processing general planar subdivisions. This involves the same techniques used for processing convex polyhedra in 3 dimensions. The connections mentioned above are examples of the ideas we intend to develop below. Where appropriate, implementation details are discussed. Connections between the problems are described to document how computational geometry is becoming unified. We now describe the five problems that form the basis of literature survey for this thesis. Each demonstrates central ideas of the field in both its statement and solu tion. The reader is encouraged to pause after reading each problem to ponder rea sonable methods that might lead to a solution. 2.1.1 Convex Hull Of a Point Set The input is a set of points P 1 P 2 .. P n in ddimensional Euclidean space. The desired output is the smallest convex set containing all of the points. In general this' set will be a polyhedron in the space. Specializing to the plane, in the non degenerate case the solution will be a polygon P with vertices from the original set. Each vertex of P will be external in some direction. That is, if P; = (x;, Y;) is a vertex of the convex hull, then there exist real numbers a and b such that for all ax.+ b1J. >ax.+ bu. I '1 J 'J This suggests a relationship with linear programming problems where the feasible region is specified as the intersection of a set of half planes (in the plane or half spaces in higher dimension). Indeed, as we will see below, this relationship has 4
PAGE 10
been exploited in establishing relationships between the convex hull problem and the half space intersection problem. Another characterization of the convex hull problem is given for an analog computing model. We might consider a board corresponding to the plane with nails hammered in at points of the point set. Of a large rubber band has a circum ference that includes of all the points; the convex hull is the set of points (or the polygon) that impede the motion of the rubber band when it is released [Sh74]. Computing convex hulls is the oldest and most popular problem of computational geometry. Many results on the problem date back to the late 1960's and seem to have their roots in statistical applications. Connections between convex hulls and mathematical programming seem to have led to the strategies behind the initial algorithms. The two earliest published algorithms are due to Graham [Gr71] in dimension 2 and to Chand and Kapur [CK70] in higher dimensions. At the time, Chand and Kapur were working at Lockheed and the convex hull problem arose in their work. They describe an algorithm based upon the principle or "giftwrapping" for computing convex hulls of spatial point sets. The key idea of their algorithm is the observation that in any dimension, an edge (of dimension I) of the convex hull belongs to exactly two facets of the convex hull. The di mension of the facets will be one less than the dimension of the resulting hull of a set of 1 000 points in 6 dimensions. A 2 dimensional rendering of their algorithm might go as follows: Find a point, p, of the convex hull and make it the current point While you haven't returned top, do Begin find the point that is most counterclockwise from p add that point to the convex hull and make it the current point. End The initial point could be the point of minimal y coordinate. The most counter clockwise point is found by scanning the points and doing tests to see if triples are counterclockwise. To see whether B is more counterclockwise than C with re spect to A, we ask whether the triangle ABC has a counterclockwise orientation. The Graham algorithm was discovered independently and arose in response to a request from scientists at Bell Labs who were doing the computation for statistical applications. His algorithm consists of two steps. First, the points are sorted with respect to an interior point. Next, a "Graham scan" is performed. During this step, points are traced about the perimeter of the polygon determining whether a left turn or right turn was made at each vertex. These correspond to tests to see 5
PAGE 11
whether we are moving clockwise or counterclockwise. If the motion is counter clockwise, the point is tentatively accepted and the scan moves forward. At a clockwise point, the point is rejected and the scan moves back. Since the scan can have no more than n3 rejections, the process halts in a linear number of steps. The behavior of this algorithm on a set of points is shown in Figure 3. As a 2 di mensional algorithm, it possesses a simplicity that is not present in the gift wrapping algorithm; however, it does not naturally extend to higher dimensions. In the years following these algorithms, numerous improvements were suggested often proposing a special model where one's algorithm was the best. Ofthe prob abilistic algorithms and analyses, those of Eddy and Floyd are probably the most significant ([Ed77, F176]). In most cases, linear behavior will result and per formance will be superior to that of a sorting based method. Preparata and Hong proposed an algorithm for 3 dimension convex hull computations in 1977 [PH77]. This algorithm is based upon divideandconquer method. They compute left and right hulls and then determine common tangents (in 2 dimensions) and a wrapping (in 3 dimensions) that can be used to merge the hulls. Already by 1977, it had been determined that sorting must play a central role in the derivation of a convex hull algorithm. It was seen that orientation would be the significant test in the inner loop. Optimal algorithms were known for 2 and 3 dimensions. Open questions remained, however, of these central question asked for extensions to higher dimensions and for methods that took advantage of the information they had on input and/or had running times proportional to their input size. For example, suppose we were given a starshaped polygon, could we com pute its convex hull in linear time? If we determined in our computation that few points were on the convex hull, could we compute the hull in faster asymptotic time? How hard would it be to just identify points of the convex hull (in higher dimensions) without building all of its structure? By 1977, most of the important questions about convex hulls had been asked and many had already been answered. The connections between computing convex hulls and sorting were already well established. What remained were questions about optimal convex hull algorithms in cases where the output size was small. For example, the giftwrapping algorithms in cases where the output size was small, or the giftwrapping algorithm on random examples. Its running time is O(nh) where n is the number of points and h is the number of points on the con vex hull. Typically his 0(/og n) and so both the giftwrapping and Graham algo rithms have cleaner implementation. However, for a set of n points all on the hull, it gives quadratic performance. The situation is worse in higher dimensions. Here, there is a lower bound of O(nd/2 ) (resp. O(n(dt)/2 ) was applied in even (resp. 6
PAGE 12
odd) dimensions. This bound comes from worst cases where it could take this amount of effort to report the output. However, in practice the output size is sig nificantly smaller. It is known that computing components of the convex hull is polynomially reducible to linear programming [DR80]. Three key questions arose. First, if we know something about the input can we use this information to compute the convex hull. Second, in two dimensions, is it possible to have an al gorithm with running time optimal in the number of points on the hull even if this number is not known in advance. Finally, can lineartime linear programming algorithms [Me83,Me84] and other techniques be used to determine all components of the convex hull in time linear in their number. Progress has been made on all of these questions in the past dec ade. Perhaps the most interesting instance of the first question is the problem of determining the convex hull of a set of points given as the vertices of a simple polygon. Here, the issue is that the points have already been sorted in some order. That is, by virtue of forming a polygon, they have and order. We would be able to do a Graham scan in linear time if we knew an ordering about an interior point. However, the order we are given on the boundary may not correspond to such an ordering. Hence, the Graham scan must be modified to determine when there has been and inappropriate backup. Here, the algorithms proposed typically compute an upper hull and a lower hull, keeping track of turns to assure that vertices are not incorrectly deleted [GY83,Le83]. It remains open to characterize other situa tions where the input contains sorting information that can be used to simplify the hull computation. Next we tum to the situation where not all points are on the convex hull. Here, the convex hull could potentially be computed in o(n log n) operations if o(n) points were on the hull. To do so requires not sorting the points since this re quires O(n log n) operations. Kirkpatrick and Seidel [KS86] were able to do so. Their algorithm finds a method of avoiding the recursive technique of finding 2 halves and merging. That technique has a running time which satisfies the recur rence relation T(N) = 2 T(N/2) + O(N). This method can be slow because it will compute many hull edges that are removed in later mergings. Instead, they find a merg.ing edge in linear time (via a variant of the linear time linear programming method [Dy84,Me83]) and reduce to the problem of solving two sub problems. This leads to a time bound deter mined by the recurrence relation 7
PAGE 13
T(N,h) = T(N/2,h1 ) + T(N/2,h2 ) +O(N) where h1 and h2 represent the number of points on the two smaller hulls and h1 + h2 = h. This yields an algorithm of running time O(n log h) that is aptly titled the ultimate convex hull algorithm. In higher dimensions, two methods have evolved to compute convex hulls. One is an extension of the giftwrapping technique of Chand and Kapur (Sw85]. The other is a method named beneathbeyond [Se81]. Using either method, near op timal algorithms can be found in the case where the convex hull is as bad as pos sible. However, the running time is exponential in the dimension even iftheout put is small. In his thesis, Siedel has improved the situation [Se86]. By using the lineartime is linear programming algorithms, he computes convex hulls in di mension d. If the input consisted of m points and the output had F facets, he can enumerate all facets in time O(mJ"(d) + d 4 Flogm). Unfortunately, f(d) is the "constant" term that arises in the Megiddo linear programming algorithm and grows exponentially with d, This algorithm is not practical, though it yields in sights that should lead to practical implementations during the next decade of computational geometry. Thus, we see that significant progress has been made on the convex hull problem since the time of the ChandKapur and Graham results. It is now the case that theoretically optimal (or nearly optimal) algorithms are known in all cases. Also, practical algorithms are understood for many situations. There remain a few im plementation issues that become relevant for applications where convex hull of large point sets were to be computed or where only part of the structure of the convex hull is needed. Many of the techniques given in this and the previous convex hull section find use in real problems of statistics, patter recognition, and classification algorithms( see e.g., [Ba84]) in describing and classifying point sets. The convex hull is and especially desirable indicator to be used in recognizing patterns because of its rotational invariance. 2.1.2 Finding Nearest Neighbors The input is a set of points P1,Prpn in ddimensional Euclidean space. The desired output is the nearest neighbor of each point or possibly the pair of points that are closest. The answer will consist of a directed (possibly disconnected) graph with an edge from ito j if P; is the nearest neighbor of Pi" Note that the graph is directed since nearest neighbor is not a symmetric relation. Or, the an swer will consist of the pair of points that realize the relationship 8
PAGE 14
min d( P,, P1 ). A closely related problem asks for furthest neighbors. Here again the output is a graph. The problem where we request the pair of points that realize the maximum interpoint distance in the set is known as the diameter problem since this distance is called the diameter of the set. We show below that the diameter of a set is de termined by two points that lie on its convex hull. This establishes a potential re lationship between this problem and the convex hull problem. In our study of the nearest and furthest neighbors, we study a structure commonly known as the closest point Voronoi diagram. This is a subdivision of space as signing to each point Pi a region Ri such that points in Ri such that points in Ri are closer to Pi than to any of the other points. Each of the regions Ri is a convex polytope. Furthermore, the only unbounded regions correspond to points on the convex hull of the original set. We will also consider a furthest point Voronoi diagram. We've mentioned these structures here because they provide interesting examples of planar subdivisions for consideration in our next problem. The next problems to be considered on point sets were computations involving distance properties. Convex hull algorithms found points that were the extrema of a set. It was natural to next focus attention on finding pairs of points that realized extremal distances. The first problems that arose were to find the closest and fur thest pair of a set of n points. First, the 1 dimensional versions of the problem were considered. Here, connec tions could be made with results from Kenneth [Kn73] on sorting and selection. It was observed that the closest pair of points must be neighbors and that finding all pairs of neighbors would require information similar to a sort [Sh75]. The fur thest pair of points would be the maximum and minimum of the set. Identifying them would require an algorithm for finding the max and min of a set of n num bers. Thus were derived informal arguments that n log n was asymptotically cor rect for the closest pair of points and linear time was necessary and sufficient for the furthest pair. Next, attention was focused on higher dimensions. In 2 dimensions, more geome try was necessary in order to establish which points were neighbors and had to be considered as potential closest pairs in the set. Here, two methods ultimately emerged for solving the problem. 9
PAGE 15
The first method was based on building a structure that would contain all potential nearest neighbors. This structure assigned to each point a region containing points closer to it than to any of the others. As such, it could be constructed from the perpendicular bisectors generated from each pair of points. Under further study, it became clear that in the plane the structure was actually a planar subdivision. This followed since the region corresponding to each point was convex and hence a lin ear number of vertices and edges. Furthermore, it was discovered that this struc ture could be built by a divideandconquer method in O(n log n) operations. Fi nally, it was learned that this structure had been previously defined in other con texts. Among its other names were Voronoi diagram and Dirichlet tessellation. From the structure arose the first subquadratic algorithms for finding nearest neighbors in 2 dimensions. The algorithms consisted of finding the Voronoi dia gram as a planar graph. Then, all potential nearest neighbors were given by con sidering edges of the dual of this graph. This dual was known as the Delaunay tri angulation and had been previously discovered in the 1930's. Since there were only linearly many edges in the Delaunay triangulation, the minimum distance could be found in linear time. This method was discovered independently be ShamodHoey and Reiss at Yale in 1975. The details appear in [SH75]. Reiss im plemented his algorithm as part of a BLISS program to perform other functions. However, there was no sense at the time of reasonable data structures for doing. such a task and his implementation was quite complex. This algorithm seemed to be an overkill to solve "so simple a problem" and alter natives were sought. Strong at IBM and Bentley [Be76] discovered a simpler di vide and conquer method to solve the problem. Their algorithms divided the points into a left set L and a right set R. For each set, the nearest neighbors were found along with their distance d(L) and d( R). Next, they observed that a point in R could only be within distance d(L) of a constant number of points of L. Hence, a merge of two sets with shortest distances known could be cone by considering only a linear number of distances and hence linear time would suffice. The Voronoi diagram method also solves numerous other problems on planar point sets (see [PS85]). The StrongBentley method extends more naturally to higher dimensions while the Voronoi diagram on n points can be of size O(nd 2 ) in dimension d. However, by 1977 neither method was sufficiently well understood to lead to a clean implementation. Even for planar point sets, it was not under stood what data structure ought to be used to implement the Voronoi diagram. Determining the furthest distance was equivalent to determining the diameter of the set. The diameter of the set was shown to be the diameter of the convex hull. This followed from the observation that if one or both of the furthest points did not lie on the hull, a larger distance was possible. Next, it was shown that poten10
PAGE 16
tial endpoints of a diameter must be antipodal. That is, each edge could be con sidered to have a range of angles that tangents to it (or one of its endpoints) might make. When the ranges for 2 vertices involve values differing by n they are con sidered to be antipodal. Connections with work of Redoes[Er46] established that there were only linearly many pairs of antipodal points. This yielded the linear time bound on the scan step of the diameter algorithm. In 3 dimensions, it was again shown that the diameter would be determined by a pair of points on the convex hull. It is known that the diameter can be realized at most a linear number of times[Gr56]. This corresponds to having only a linear number of antipodal pairs in the plane. However, after a few false starts, the problem of finding the diameter in fewer than a quadratic number of operations remained open. Thus, nearest neighbors required n log n operations by the element uniqueness lower bounds. N log n algorithms were known in all dimensions by the divide and conquer method. There was also a method in 2 dimensions based the Voronoi diagram that had numerous other applications. However, no clean method of implementing the Voronoi diagram or handling degeneracies of point sets (e.g., 4 cocicular points) was known and there was no accepted data structure for the problem. Curiously, at about this time, Baumgart was writing a thesis in computer vision [Ba74] that would provide such a structure. With regards to fur thest distances, linear time was known to be necessary and sufficient in 1 dimen sion, n log n was known to be necessary and sufficient in 2 and the problem was open in higher dimensions. The rediscovery of the Voronoi diagram provided a framework for discussing nearest neighbor problems. Furthermore, by 1977, there were known algorithms requiring O(n log n) operations for building Voronoi diagrams. While on such algorithm had been implemented [Re77], implementation was a nontrivial task. In addition to the problem of finding a data structure upon which to build the im plementation, there were questions of numerical accuracy and degeneracy. The divideand conquer nature of all known algorithms was cause for concern in terms of resolution deficiency and numerical stability of algorithms. What was needed was a simple yet complete explanation of the Voronoi diagram. There also re mained the question higher dimensional Voronoi diagrams. However, Brown es tablished a relationship between convex hulls and Voronoi diagrams [Br79a] that related convex hull computation in dimension d+ I to Voronoi diagram construetion in dimension d. As a result, the results of Seidel [ Se86] given in the previ ous section can be applied to higher dimensional cases. 11
PAGE 17
Significant work was done of the development of good Yoronoi diagram algo rithms during the decade 19771987. Authors sought clean formulations of algo rithms that would be implementable, extensions to the Yoronoi diagram that would also handle line segments and curves and data structures to use in doing the computation. This work is presented in two papers [Fo86,GS85] that resolve many of the outstanding issues and represent the current state of the art. Guibas and Stolfi[GS85] present a generalization of the wingededge data struc ture of Baumgart [Ba74] to a quadedge structure that can represent planar subdi visions and their duals simultaneously. Their work can even represent subdivi sions on nonorientable manifolds. Their development is complete justifying the choice of the quadedge as the data structure of choice for the representation of subdivisions of2manifolds (ie planar subdivisions and 3dimensional polyhedra). They are able to give a onepage implementation of the Yoronoi diagram ( and its dual, the Delaunay triangulation) based on their data structure. Since the manipu lation of their data structure is straightforward, the result is a simple implementa tion of the Vorohoi diagram that is reasonably robust. They can identify a set of (more than 3) cocircular points efficiently. This identification would be desirable in constructing Delaunay triangulations since the resulting regular polygon should not be triangulated. Fortune [Fo86] presents the first optimal algorithm for computing the Voronoi diagram that is not based upon divide and conquer. He uses a sweepline tech nique. In sweeping across the plane, the Voronoi region of a point will be en countered before the point is reached. Any sweepline algorithm must have a method of identifying the region's beginning as it occurs. An algorithm that does otherwise is destined to have quadratic running time. Fortune introduces a trans formation of the plane that causes a point's Voronoi region to be first encountered at the point when doing a planar sweep. The result is that he is able to insert the region at that point, in a manner similar to the plane sweeps of [SH76, 8078] de scribed above. On insertion, he can identify when adjacent regions will change discontinuously in linear number of events and hence the algorithm has O(n log n) running time. His transformation is a conceptual device for understanding the algorithm. By undoing the transformation, Voronoi regions are identified as they appear. He gives further transformations for cases where the Voronoi diagrams of line segments or weighted point sets are to be computed. *Note that this differs from the sweepline for finding all segment intersections shape or disappear add ing an event to the event queue. There can only be a where there can be a quad ratic number of events. 12
PAGE 18
The Voronoi diagram has been transformed from a theoretical device used to solve problems into a tool that is used in practical situations. This progress obvi ously arises from the development of algorithms that can be used in the practical situations. Surprisingly, it also arises from the dissemination of some of the early papers on the subject to engineers in related fields who have problems that need such a structure. Thus, it seems likely that these applications will derive the next decade of work on the Voronoi diagram and Delaunay triangulation. A step in this direction is the thesis work of Mike Laszlo [DL87,La87]. This work provides an extension of the work of Guibas and Stolfi to spatial subdivisions. There ap pears to be the potential for the application of this work to numerous problems involving computational geometry as well as applications to finite applications to finite element calculations (see e.g.[JB87]). 2.1.3 Planar Subdivision Searching The input is a subdivision of the plane into regions that are convex polygons that only overlap in edges and vertices. For example, these polygons may have arisen from a closest point Voronoi diagram as described above. For convenience, let us assume that the polygons have all been triangulated so that the input is set of tri angles T1 T2 ... Tn. We will assume further that two triangles overlap in an edge if the complete edge belongs to both of them. The problem here is to efficiently search the planar subdivision. Searching a pla nar subdivision assumes that the subdivision has been sorted. This sorting is done during a preprocessing phase. The result of the sorting is to create a data structure that can be searched. The result of sorting is to create a data structure that can be searched. It is assumed that the sorting is done once and supports numerous searches. Hence, the cost of the sorting can be amortized over the searches. Al gorithms for planar subdivision searching have their complexities specified by a triple of numbers. This triple involves the preprocessing time for the sort, the space needed to store the searching data structure and the time for each search. Planar subdivision searching can be viewed as a multidimensional extension of binary search. In the latter problem, we are given a set of numbers that can be viewed as a set of points on a line (say the xaxis with the numbers being x coor dinates). Then, the problem becomes one of organizing the points, which we do by sorting x coordinates, and searching by standard binary search (see e.g., [Kn73]). This operation requires n log n time for the sort, linear space to store the search data structure and log n time for each search. Just as planar subdivision searching is an extension into higher dimensions. We will mention below exten13
PAGE 19
sions of the problem into d dimensions and consider the general problem of spa tial subdivision searching. Our next problem is a different version of the multi dimensional searching problem. Planar subdivision by polygons in dimension 2 is the generalization of the subdi vision of a line by a set of points in dimension 1. There are numerous applica tions of an algo;ritlun for searching planar subdivisions. For example, the Voronoi diagram presents applications where we might want to locate the region of a point in a planar subdivision. Similar problems arise in using finite element analyses to solve partial differential equations. Also, the problem arises in the "postoffice problem" as stated by Knuth. Finally, in higher dimensions, the knapsack prob lem can be posed as a problem of searching spatial subdivisions. Planar subdivi sions for applications mentioned here might arise in either of two similar contexts. Either they arise from the regions formed from the arrangement of a set of n infi nite lines. Or, they might arise from the Euclidean embedding of a planar graph. In either case, the structure to be searched is the same. And, we can assume that the structure can be built once and will remain static for a large number of searches. Hence, it is reasonable to think of preprocessing to build a usable search structure that can be rapidly searched. The central problem here is that of finding an ordering on the regions if the planar subdivision. The earliest methods [DL 74,DL 76] did this by reducing the problem dimension by one using "slabbing methods". Here the observation was that the planar subdivision could be divided into a set of slabs that could be easily searched. This was done by creating slabs in which the regions were monotoni. cally ordered. Furthermore, the slabs were created in monotone order as well. Thus, the unknown problem of how to do binary search in two dimensions was reduced to the problem of doing 2 binary searches. The first was a binary search through an ordered set of slabs to find the correct one and the second was a search through an ordered set of regions in the slab to find the correct one. However, the storage required can be such that a segment has to be stored within every slab. The slab algoritlun resulted in an algoritlun requiring O(n2/og n) preprocessing time and O(n2 ) preprocessing space on the worst case. However, the retrieval time was 2 log n accounting for two binary searches of n items each. Thus, a solution existed. The challenge now became to improve upon the running times. For such algorithms, complexity is measured via triples of numbers (PT(n).PS(n),ST(n)) where PT(n) (resp. PS(n)) represents the preprocessing time (resp. space) and ST(n) represents the search time. It became clear that PT (n) and PS(n) could be made linear (e.g., by doing nothing), but that if ST(n) were to be 14
PAGE 20
sublinear then PT(n) would have to be at least n log n. Then the question arose as to whether, it was also possible to achieve a solution that was (n log n, n, log n). Lee and Preparata [LP79] reduced the bestknown algorithm to (n log n, n, lofln). Finally, in an elegant fashion, Lipton and Tarjan [LT77a,L T77b] achieved a solu tion of (n log n, n, log n). Their methodology was to show that divide and con quer could be applied to the problem of splitting a planar graph. In linear time they were able to find a separator that divided the graph into two parts each con taining a positive fraction of the original vertices. Furthermore, the separator was small. While their method was not easily implementable, it did provide the first proof that divide and conquer could be made to work in 2 dimensions as it had in 1. This supported the belief that sorting could be applied to higher dimensional geometries. The situation regarding planar subdivision searching in 1977 was similar to that of Voronoi diagrams. The LiptonTarjan result had solved the theoretical problem but was not of practical significance. It remained for the practitioners to propose techniques that would be simpler and hence likely to apply in practice. Two themes emerged in the development of practical techniques. The first is the notion of hierarchical search. Kirkpatrick [Ki83] proposed this as a technique for representing a planar subdivision. The basic idea is that when searching a planar subdivision, more and more detail are needed about less and less of the planar subdivision as the search progresses. This suggests a technique whereby the planar subdivision consisting of a constant number vertices and re gions. This can be searched in constant time. When the appropriate region has been identified, only this region of the subdivision is refined by adding vertices and so further subdividing the region. This is shown in Figure 7. Note that we only visit one of the subdivided regions. This process continues for long n steps. Then, in 0(/og n) time, the search has been completed. The hierarchical technique appears to be a basic technique a computational geometry [DS87] having also found application to polyhedral in tersection detection [DK85]. The other approach follows (EGS87] the topology of the arrangement and applies a divide and conquer technique. It is based upon ideas that led to some of the ini tial methods for solving this problem. The idea is simple. A chain of edges is found that divides .the subdivision in half. The chain also has the property that it is monotone and hence can be searched quickly. Furthermore, future chains build upon the initial chain so that all searches can be completed in time O(log n). A similar technique is given in [ST86]. 15
PAGE 21
The progress on planar subdivision searching provides the insights necessary to lift searching from a onedimensional primitive to a well understood operation in 2 dimensions. The techniques described above have made the problem suffi ciently practical to allow its use in application systems. Indeed, as graphical edi tors and window systems become more complex, naive algorithms are demon strating behaviors that are insufficient and these techniques will be needed. There remain open the problems of extending searching to spatial subdivisions. The re sults in [DL74,DL76] stood unchallenged for almost a decade. Recent work [ Ch85a, C185] has improved upon the space bounds and gave probabilistic algo rithms that vastly improve the older results. Extensions beyond dimension 3 re mam open. The subsequent development of rage searching techniques after quad and kd trees is interesting within the context of this thesis. As opposed to the other problems we consider, this is an example where the first algorithm found was the practical one and led work of theoretical interest. In all other cases, the theoretical work typically happened before 1977 and the work of the last decade has involved re fining theoretical insights and looking for simplified algorithms that are imple mentable. For range searching; the initial algorithms were driven by practical con siderations and are quite implementable. The algorithms that followed were born out of theoretical considerations. Most of these theoretical considerations sought to fix the bad worst case possibilities of range searching and extending rage searching to work for nonrectangular search domains. 2.1.4 Range searching Next we consider a problem that is similar in spirit to the planar/spatial subdivi sion search. The input is a set of points P1,P2 ... ,Pn in ddimensional Euclidean space. We again want to preprocess the points to allow for efficient searches. Here the searches involve ranges in the ddimensional Euclidean space and ask which of the points lie in the range. Here a range can be viewed as a box, as a simples or as a halfspace. For example, a range might be described as a set of ranges m1,m2 .. ,mciMI'M2 ... Md where all points P; = (x;/'x;2 .. ,x;d) with mk
PAGE 22
the number of points but also for the names of the points in the range.. An algo rithm for the reporting problem will obviously solve the counting problem. How ever, there may be a more efficient solution if only a count is desired, this distinc tion is important in what follows. In considering this problem in conjunction with the previous problem, it appears as though a dual transformation should be possible (see e.g., [Br79b]). That is, it should be possible to transform this problem into one of point location in a spatial subdivision. In practice, this appears to not be the case as the algorithms devel oped for the problems have significantly different forms. The range searching problem has its roots in the implementation of databaselike queries for determining ranges of values. We might be given a set of attributes (e.g. salary, age, height and weight) and represent employees of a company as points in R4 corresponding to their 4tuple of these attributes. Then, we might probe in order to determine all employees satisfying a set of ranges simultane ously. This might consist of finding all employees within a certain age range, having salary between specified limits, . Usually, the database is reasonably static so that preprocessing of the points can be done offline in support of faster queries. Knuth posed the problem of designing a data structure [Kn73] for multi dimensional searching problems. The search queries to be supported are similar to those considered in the previous section, however subtle differences lead to al gorithms of a very form. Bentley was the first to propose algorithms of doing rectangular range queries [Be75,BS75,FB74]. He proposed a structure called the "quad tree" for doing toe search in 2 dimensions. Knuth later named the generalization of this search struc ture the "kD tree". The quad tree is relatively simple. The problem of searching in 2 dimensions is reduced to' two methods of searching in 1 dimension just as we saw in the slab methods of the previous section. However, here the reduction oc curs in an iterative fashion that is symmetric. The resulting algorithm is one that has good average case behavior, though its worst case can be bad [LW77]. The quad tree is a 2 dimensional binary search tree that alternates search direc tions. So, given a set of points in the plane, a divider is found parallel to the yaxis. Next, each of the two regions is divided by a divider parallel to the xaxis. The alternating process continues until regions contain only 1 point. Searching of this structure now can be cone by following a path down the binary search tree corresponding to the subdivision given. If we wish to determine how many points lie within a rectangle, we determine which of the two regions at a node the rectangle intersects and use this information to proceed. As opposed to 17
PAGE 23
other binary trees, it could be the case that both sub trees need to be followed. Quad tree based algorithms were already popular within computer graphics and pattern recognition by 1977[Hu78,HS79]. Indeed, the method of traversing. a quad tree uses techniques similar to those proposed by Warnock [Wa69] for use in hidden surface removal. However, the discoveries appear to have been inde pendent. As opposed to the other problems we discuss in this section, range searching is a situation where the practical has driven the theoretical. The quad tree was pro posed as an interesting search structure and remains the data structure of choice in numerous applications [Sa84]. However, the problem of optimal range searching has attracted significant theoretical interest. Here, the search problem is the half space searching problem. So, points are to be preprocessed so that those lying in a halfspace can be determined efficiently. The search structure is to require only linear space and the search time is to be sublinear. It has been shown that this searching model embodies all polyhedral range searches [DE84]. It has also been shown that with more storage, all searches can be achieved in O(log n) time [EMK82]. Finally, a lower bound ofO(n113) is known on the problem [Fr81]. The first observation that sublinear halfplane searches were possible is due to Willard [Wi82] who produced an algorithm of search time O(n774). The method of Willard consisted of creating 6 regions in the plane each containing an equal share of the points and such that a line would intersect no more than 4 of them. Then, a search could be done by reducing a problem to 4 problems of 116 of its original size in constant time. This yields the recurrence relation T(n) = T(n/2) + T(n/4) +c with solution O(n'695). In higher dimensions, there was a flurry of activity resulting in and upper bound of [DEY84[ in 3 dimensions and an observation that the commonly used technique would not extend past 4 dimensions {Av85]. Next, Cole [Co85] pro posed an alternative technique using portions of multiple hyperplanes for the sub division rather than simply building quadrants or octants as had been cone in di mensions 2 and 3. He was able to achieve an upper bound of O(n977 ) in dimen sion 4. Next, there followed a construction built upon Cole's that a query time of O(n) [YY85]. Finally, Haussler and Welzl [HW86] showed that a probabilistic technique could be used to achieve a query time of O(n). Their construction is related to the VapnikChervconen.kis dimension of probability theory. Curiously, their construction is nonconstructive. They show that an algorithm exists but 18
PAGE 24
carmot demonstrate one. Their enet construction appears likely to have signifi cant application to related problems of region counting in arrangements [GS87]. A similar construction was discovered independently by Clarkson [C 186]. As mentioned above, range searching has proceeded in a different direction from that taken in the other four problems we discuss. There has been some develop ment of techniques for using range searching and new applications. However, the significant advances during the past decade have been theoretical. Elegant mathematics has resulted from the analysis of extensions to what is basically an intuitive method. It is likely that these extensions will generate new results of practical interest during the next decade. 2.1.5 Intersection problems Our final problem is the problem of polygon intersection computation. We con sider here a version of the problem for which the input is a set of line segments s ts2'" .. ,sn. There are two versions of the problem. In one, the output is a listing of all pairs of segments that intersect. For the second, a YES/NO answer is de sired tilling whether any pair of segments intersect. The analog to the report ing/counting problems above is obvious Clearly, determining whether any pair of segments intersect is no harder than enumerating all pairs that do. Intersection is a basic problem of geometry. This problem forms a basis for many others where the line segments are replaced by axisoriented rectangles or convex polygons. The distinction between detection and computation is significant since for many applications, a fast detector is sufficient. Initially, intersection problems seem different from the other problems we con sider here. In particular, there is no obvious sorting or searching in the problem. But, we shall see otherwise below. Convex hull computations are the analogs of sorting. In the same way, intersection of convex bodies is the analog of merging two sorted lists. This will become clear in the next sections. The other lower bound we consider will be on the element uniqueness problem. In this problem the input is the set x1,x2' ... ,xn of real numbers. The output is a YES/NO answer corresponding to whether there exist 1
PAGE 25
(i.e., tests comparing .Ea .x. to c for real numbers a., i=l, ... nand c). In this case, , it is known [DL79,Be83] that O(n log n) queries are necessary and sufficient. No study of geometry would be complete without a discussion of intersection problems. Much of the motivation for intersection comes from a desire to dis play geometries on a graphics screen. In 2 dimensions, we are often concerned with a screen full of overlapping windows and our goal is to determine which windows need to repainted or a paint order for windows. In 3 dimensions, inter section detection and computation forms the basis of most hidden surface algo rithms. Both of these problems involve determining intersections among line segments and identifying which (if any) pairs from a set of line segments inter sect. This problem was first considered by Shamos and Hoey [SH76] who developed a method based upon the notion of plane sweep to trace across a set of segments and determine whether any pair intersected. Their technique builds on ideas simi lar to the slabs methods for planar subdivision searching. Its implementation gen eralizes the implementation of scanline based hidden surface routines as devel oped by Watkins [Wa70]. The plane sweep uses slabs defined by endpoints of the segments. We consider these endpoints to be projected on the xaxis and slabs to extend vertically. Seg ments at any value of x can be ordered by their y coordinates. And, within a slab, this order can change only if an intersection has occurred. On the boundary be tween slabs, a new line segment is inserted or deleted. For n line segments, there a 2n slabs. A very naive algorithm is to determine the slabs and sweep across slabs determining for each slab if it contains an intersection. This can be done in timelinear in the number of segments within the slab since the order of the seg ments at the start of the slab is known. This leads to a quadratic algorithm. Shamos and Hoey improve upon this algorithm by observing that if two segments intersect, they must have been adjacent in the ordering of segments at some point (not necessarily at a slab boundary). 1. When the second segment was inserted, they were adjacent in the ordering of segments 2. When the second segment was inserted, they were not adjacent in the ordering of segments, but all segments occurring between them either intersect one (or both) of the segments or were deletedbefore they intersected. 20
PAGE 26
3. When the first segment was deleted, they were adjacent in the ordering of seg ments. 4. When the first segment was deleted, they were not adjacent in the ordering of segments, but all segments occurring between them either intersect one (or both) of the segments or were inserted after they intersected. As they do their sweep, Shamos and Hoey are careful to always maintain the order of the line segments. This is done by using a balanced tree scheme so that inser tions and deletions can be done at a cost of O(n log n) operations. Whenever a segment is about to be inserted, it is checked against its. neighbors for intersection. And, whenever a segment is bout to be deleted, its neighbors are checked to see if they intersect each other. Thus, any two segments that will ever be adjacent in some ordering are checked for intersection. By appealing to the list of situations given above, it can the be shown that if there is an intersection it will be detected. Furthermore, the running time of the algorithm is O(n log n) because the x coordi nates of the endpoints must be sorted and because O(n) insertions and deletions must be done. Finally, this is optimal by appealing to the element uniqueness problem and showing that a sort must be done. The only failing of the ShamosHoey algorithm is that it will not necessarily find all intersections. Bentley and Ottman [B078] proposed an extension that does compute all intersections. For this algorithm to work, more slabs are necessary. Every time an intersection is found, an old slab must be subdivided to create a two new ones. With appropriate modification, it can be shown that the SharnosHoey algorithm always detects intersections before they occur and therefore, we will always know in advance that slabs will have to be created. When we reach an in tersection point corresponding segments that intersected. This requires O(log n) operations. Putting this together yields an algorithm with running time O(n log n + I log n) for reporting all intersections where I is the number of intersections re ported. In case where I is small, this represents a significant improvement over previous methods. However, in cases where I is large, this algorithm will actually be slower than the naive algorithm. Since the work of Bentley and Ottman, there has been significant progress on in tersection problem. There have been both theoretical and practical improvements to their algorithm. Also, numerous other directions in intersection have been pur sued. Furthermore, work on significant practical problems, most notably hidden surface removal has begun to occur at both the theoretical and practical levels. On the basic problem of finding all intersections, improving the upper bound of O(n log n + I log n) has been a nagging open problem. Mairson and Stolfi [MS82] 21
PAGE 27
showed that if there were two sets of line segments, n painted red and n painted green such that all intersecting pairs involved a red segment and a green segment, then O(n log n + 1) and was efficient in practice. Finally, Chazelle [Ch84] im proved the upper bound to O(n lofln!loglogn + 1). Kalin and Ramsey [KR84] did significant empirical studies of algorithms for the problem and concluded that Myer's was the fastest in practice and that further improvements of his algorithm were possible. Their conclusion was that for randomly generated line segments, a technique that wisely combines enumeration and bucketing will give the best be havior. While, most pictures sent to a hidden surface remover are not random, it is reasonable to believe that their observations will extend. Most of the intersection work is geared towards extending methods towards ob jects of increasing complexity. The holy grail of intersection problems is an algo rithm that will operate on a database of arbitrary objects and allow for insertions and deletions of objects as well as determining all objects that a new object inter sects. The applications of such a technique are manifold. Towards this end, there has been work on understanding intersections of convex polygons [Sh75], convex polyhedra [MP78,DK83,DK85,CD87] and. curved polygons [DSV86]. Also, searching techniques have been developed to strike at the problem of a dynamic database (Ch86]. In all cases, the central problem is the same. One is always seeking a method of sorting objects that do not have a total order. The goal of such sorting is to be able to later search a database of such objects to determine overlap. An interesting technique that gets at the essence of this problem is the work of Chazelle [Ch85b] on embedding intersection problems in higher dimen sional spaces and solving them by subdivision searching techniques. The next step is to relate this work to hidden surface removal problems in a mean ingful fashion. This has begun with theoretical work [De86] giving lower bounds for worst case hidden surface removal. Of more immediate potential practical im portance are results that use theoretical techniques to study key problems. Nota ble in this category is the work of Hubschman and Zucker on frametoframe co herence [HZ82] and the forthcoming thesis of Panduranga [Pa87]. The former paper is a first attempt at using object space techniques to represent an object for multiple hidden surface removals. The latter considers the problem of ray tracing in object space by tracing back reflections onto primary objects. As we've seen above, intersection problems represent an area where theory has applied in prac tice and the current progress suggests that this indeed continue. 22
PAGE 28
2.2 Spanning Trees literature survey: Input description: A graph G = (V, E) with weighted edges. Problem description: The subset of E' c E of minimum weight forming a tree on V. The minimum spanning tree of a graph defines the cheapest subset of edges that keeps the graph in one connected component. Telephone companies and networking companies are particularly interested in minimum spanning trees, because the minimum spanning tree of a set of sites defines the wiring scheme that connects the sites using as little wire as possible. It is the mother of all network problems. Minimum spanning trees prove important for several reasons: They can be computed quickly and easily, and they create a sparse subgraph that reflects a lot about the original graph. They provide a way to identify clusters in sets of points. Deleting the long edges from a minimum spanning tree leaves connected components that define natural clusters in the data set. They can be used to give approximate solutions to hard problems such as Steiner tree and traveling salesman. As an educational tool, minimum spanning tree algorithms provided graphic evidence that greedy algorithms can yield provably optimal solutions. Two classical algorithms efficiently construct minimum spanning trees, namely Prim's and Kruskal's. We refer the reader to Chapter five of this thesis for a more detailed implementation and discussion. Good expositions on Prim's [Pri57] and Kruskal's [Kru56] algorithms will appear in any textbook on algorithms. The fastest implementations of Prim's and Kruskal's algorithms use Fibonacci heaps [FT 87]. A recent break through on the minimum spanning tree problem is the lineartime randomized algorithm of Karger, Klein, and Tarjan [KKT95]. The history of the minimum spanning tree problem dates back at least to Boruvka, in 1926 and is presented in [GH85]. 23
PAGE 29
3. MinimumCost Spanning Tree Algorithms In this chapter, we present a discussion of Kruskal's algorithm and its implemen tation. An important part of the implementation of Kruksal's algorithm is han dling of the Set Merging Problem. Three different techniques for handling the Set Merging Problem, will be presented and analyzed. 3.1 Kruskal's algorithm A spanning tree for a graph G is a subgraph of G that contains every vertex of G and is a tree. But every connected graph has a spanning tree, and any two span ning trees for a graph have the same number of edges. A weighted graph is a graph for which each edge has an associated real number weight. The sum of the weights of all the edges is the total weight of the graph A minimal spanning tree for a weighted graph is a spanning tree that has the least possible total weight compared to all other spanning trees for the graph. If G is a weighted graph and e is an edge of G then w(e) denotes the weight of e and w(G) denotes the total weight of G. The problem of finding a minimal spanning tree for a graph is certainly solvable. One solution is to list all spanning trees for the graph, compute the total weight of each, and choose one for which this total is minimal. (Note that the well ordering principle guarantees the existence of such a minimal total.). This solution, how ever, is inefficient in its use of computing time because the number of distinct spanning trees is so large. For instance, a computer graph with n vertices has n(n2 ) spanning For graphs with n vertices and m edges, Kruskal's algorithm can be implemented so as to have worst case orders of mlog2m and n2 respectively. In other words, Kruskal's algorithm is devised to find a Minimum Cost Spanning Tree of the Graph, where a spanning tree is an unconnected tree that connects all of the vertices of the graph. The cost of the spanning tree is the sum of the costs of the edges of the spanning tree, and a minimum cost spanning tree is simply (one of) the lowest cost such spanning trees for a given graph. Therefore, in Kruskal's Algorithm, the edges of a weighed graph are examined one by one in order of increasing weight. At each stage the edge being examined is added to what will become the minimal spanning tree provided that this addi tion does not create a circuit. After n1 edges have been added (where n is the 25
PAGE 30
. number of vertices of the graph), these edges together with the vertices of the graph form a minimal spanning tree for the graph. The following is a short description of the algorithm: Kruskal's Algorithm for obtaining the minimum tree T for graph G(v,e) is as fol lows: initialize T to empty repeat the following for each edge e in order of nondecreasing weight: if T + e is acyclic then add e to T until I Tl = n1 To better understand the algorithm, and to be able to show how Kruskal's works for a graph, the following is a more detailed version: Input: G is a weighted graph with n vertices [build a subgraph T of G and to consist of all the vertices of G with edges added in T} order of increasing weight. At each stage, let m be the number of edges of 1. Initialize T to have all the vertices of G and no edges. 2. Let E be the set of all edges of G and let m := 3. [Precondition: G is connected] while (m
PAGE 31
An important step in the implementation of Kruksal' s is the determination of whether and edge being for addition to T would create a cycle. This is an example of the problem of Set Merging, which is discussed at length below. In analyzing K.ruksal's Algorithm, two activities have to be examined. One is the priority queue that offers the edges in order of nondecreasing weight. The other is the algorithm used to handle set merging. Below, it will be shown that set merging can be handled in almost linear time, so it turns out that the dominating activity is creating and maintaining the priority queue. This requires at most O(e log e) time, where e equals the number of edges in the graph, so that is the worst case timing of K.ruskal' s Algorithm. If the edges were to be given in sorted order, then the set merging activity would dominate, and K.ruskal's would be essentially a linear time algorithm. Another algorithm for finding a Minimum Cost Spanning Tree is Prim's algo rithm. Prim's works by starting at a specific node of the graph and building a tree out from that node. At each step, the edge of lowest cost incident to the subtree under construction is added to the subtree. After all vertices are cmmected, an minimum spanning tree is obtained. The major differences between K.ruskal's and Prim's is this: In K.ruskal's, many different, unconnected, subtrees are built and united until a single tree spans the entire graph. In Prim's, only a single subtree is created and expanded until tree spans the graph. 3.2 The set Merging Problem The Set Merging Problem is as follows: Set S is divided into subsets such that no subsets overlap, and the union of all of the subsets is equal to S. For any element in S we need to be able to find out which subset it is in, and for any two subsets of S, we need to be able to merge them into a single subset. The process that finds out and returns which subset an element e is in will be called FIND(e). The process of merging two subsets will be called UNION(a,b); this is called a destructive merge because the history of the subsets in which an element resides is lost. A sequence of FINDs and UNIONs must be processed on line, that is, the FINDs and UNIONs will be executed sequentially, and the re sults of each FIND and UNION must be determined before the next process is executed. For any setS on n elements, there can be at most n1 UNIONs, for at that point the set will be contained in a single subset. To completely describe the status of set S would require n FINDs between each union, so in any nontrivial applica27
PAGE 32
tion of Set Merging, there wold be at most O(n2 ) Finds in the sequence. In gen eral, form FINDs and n UNIONs, n sm sn2 This Set Merging Process applies to Kruksal's Algorithm as follows: As edges are added to T, subtrees of G are created. The vertices in these subtrees are subset of the entire set Vofvertices. When edge(v, w) is being evaluated to see whether it will cause a cycle in T, it must be determined in which subsets v and w reside. That is, FIND(v) and FIND(w) must be executed. If they are in the same subset, then a cycle would be created, so edge(v, w) is discarded. If they are in different subsets, then edge(v, w) will be added to T, and, in addition, the two sets contain ing v and w must be merged into a new set by the process UNION(FIND(v), FIND(w). The different algorithms for handling the Set Merging Problem will be presented. 3.2.1 Nai"ve SetMerging Algorithm The naive setmerging algorithm consists of an array of size n (the number ofvertices) called VertexSet. The array entry i contains the name of the subset in which vertex i resides. Thus the process FIND(i) can be done in constant time. The UNION(a,b) operation is complete by going through the entire array and changing all a's to b's. Thus a single UNION is O(n) steps, and because there can be at most nI UNIONs, the total time required for all of the UNIONs is O(n2). So for the sequence of m FINDs and n1 UNIONs, the total time is O(m+ n2). If m is close to n2 then this algorithm is optimal according to asymptotic analysis. Linked List SetMerging Algorithm A simple improvement to the naive algorithm is to maintain the subsets as linked lists. An array Vertex Set is still maintained which contains the name of the set for each vertex, so FIND is still done in constant time. The linked lists for the subset members can be implemented by using two additional arrays of size n called FIRST and NEXT. The FIRST array contains the "first" element of a subset. The NEXT array contains the element following element i in the subset. When element UNION instruction does not need to traverse an entire array of size n. To complete UNION(a, b), then i is the last member of a subset, the entry NEXT(i) is zero. Now the e linked list for only one of the subsets a or b needs to be trav ersed and its elements in the array V ertexSet are changed. Then the two linked lists are concatenated to form a single linked list for the new merged subset. An other significant improvement can be made if the UNION process merges the smaller set into the larger; This is called weighted union, 28
PAGE 33
For example, if the size of set a is one thousand, and the size of set b is two, much less work is done if b is merged into a and not the other way around. In order to implement this another array is maintained that keeps track of the sizes of the various subsets.The analysis of weighted union is as follows: assign the amount of work of the UNION process to each element whenever it is merged into another set. Now, by merging a smaller subset into a larger, the resultant subset is at least twice the size of the smaller subset. Obviously a subset can double in size only log2n times before it is as big as the entire set, so each element i can be merged at most log n times. Therefore the entire sequence of UNION operations using weighted UNION is O(n log n). Thus the entire sequence of m FINDs and n1 UNIONs is O(m + n log n), and if m = O(n log n) this algorithm is optimal. 3.2.3 Tree structured set Merging AlgorithmWeighted Union with Path compression In this algorithm, the various subsets are represented by a rooted, undirected tree. The element at the root is the name of the set, and all other members of the set are children of the root. A UNION(a, b) operation can now be done in constant time by making the root of one of the subsets a child of the root of the other subset. This tree structure is implemented by using an array in which the entry for each element i is the parent of i Thus the FIND( i) instruction simply follows the path from element i up to the root of the tree containing i. The time required is dependent on the length of this path, which of course is no longer a constant as it was in the first two algorithms. In worst case, this path could have length nl.However, this path length can be kept very small by using weighted union and path compression. As before, weighted union simply means making the root of the smaller subset a child of the larger subset. Path compression is accomplished during a FIND op eration as follows: each vertex on the path from element ito the root ofthe subset containing i is made a child of the root. The combined effect of weighted union and path compression is to reduce the time required for a series of FIND opera tions to almost linear. The degree to which it is not linear is the "inverse" of the Ackerman function, which for all practical purposes is a constant less than five. Therefore (simplifying somewhat by treating the effects of the inverse Ackerman function as a constant) the treestructured setmerging algorithm with weighted union and path compression can process m FINDs and n UNIONs in O(m + n), which is essentially linear. 29
PAGE 34
4. NPCompleteness results In this chapter we present two computational complexity results addressing prob lems in network design and construction. Specifically, we will show that finding spanning trees, whose cost is restricted via an equality restraint or via upper and lower bounds, is NPcomplete. This result has significant implications in that it justifies the use of heuristics or approximation algorithms to address the above mentioned problems. First, let us formally define the Subset Sum problem. This problem is NP complete and it is closely related to one of the original six basic NPcomplete problems identified by Karp in [Karp72], the Set Partition problem. 4.1 Subset Sum INSTANCE: Set A = (a,, a2, a3, .... an} of positive integers and a positive integer B. QUESTION: Is there a subset of A, say A', such that the sum of the element val ues in A is exactly B? Let G = (V, E) be a weighed graph, with each edge having a nonenegative integer weight, and K be a positive integer. 4.2 First Theorem Given G and K, the problem of determining if G has a spanning tree of weight equal to K is NPcomplete. 4.2.1 Proof We can nondeterministically guess a spanning tree T, and verify in deterministic linear time if its weight is equal to K. Hence the problem is clearly in NP. Let us take an instance of the subset sum problem A= {a,, a2, a3, .... ,an}, B. We shall construct a graph G that will have a spanning tree of weight B if and only if the subset sum problem has A c A of weight B. First, construct one node and label it as v 0 For each a; in A construct a pair of nodes V;J and v;2 and connect them with an edge of weight 0. Next, connect each node v;1 to vo with an edge of weight 0. Finally, connect each node v;2to vo with an edge of weight a;. Now, if A has a subset A of weight B then G will have a spanning tree T of the same weight, by simply choosing those nonzero weight edges (vo, vJ that correspond 30
PAGE 35
to each a; in A'. The remaining edges of T would, of course, be the zero weight edges (v0 v11) for each a1 not in A' and the nzero weight edges from ViJ to v;2 Conversely, ifG has a spanning tree ofweight B then A' ofthe same weight obvi ously exists. 4.3 Second Theorem Given a weighted graph G and two positive integers C and D, the problem of de termining if G has a spanning tree of weight bounded by C and D is NPcomplete. 4.3.1 Proof By restriction, follows directly from the first Theorem with C = D = K. 31
PAGE 36
5. The Sweep Line technique When space sweep is used, it reduces anndimensional static problem to a (n1) dimensional dynamic problem. Let us say that we have one or more objects lying in some Cartesian space, and space being completely traversed by a moving hyperplane, (the sweeping plane). For any time T, the sweeping plane intersects some subset of the elements, the active ones. These intersections evolve continu ously in time, except at certain times where discrete events occur: new objects are added to the active set, old objects get deleted, or the configuration of the inter sections on the sweeping plane undergoes a major change. The sweep algorithm is a simulation of this process, is discrete, and uses two data structures: a queue of future events, and a representation of the active set and its cross section by the current sweep plane. Each iteration removes the next event from the queue, and updates the cross section data structure to mirror the effect of the sweeping plane advancing past that point. Depending on the algorithm, each iteration may also reveal some future events that were not known at the beginning of the sweep; these must be inserted into the event queue, in the appropriate order, before pro ceeding to the next iteration. This technique has been used in algorithms for com putation of Voronoi diagrams. 5.1 Voronoi Algorithm using sweep line technique In this section we will describe an optimal sweep algorithm for computing the Voronoi diagrams of n sites in the plane. For simplicity, we assume the sites are in general position, so that there are no two sites with same x or y coordinate, and no four sites are cocircular: Input: A set P: {p1 .. Pn} of point sites in the plane Output: The voronoi diagram Vor(P) given inside a bounding box in a doubly connected edge list structure 1. Initialize the event Q with all site events 2. While Q is not empty 3. do Consider the event with largest ycoordinate in Q 4. if the event is a site event, occurring at sit p; 5. then HANDLESITEEVENT(p;) 6. else HANDLECIRCLEEVENT(pL), where PL is the lowest point of the circle causing the event 32
PAGE 37
7. Remove the event from Q 8. The internal nodes still present in T correspond to the half infinite edges of the voronoi diagram. Compute a bounding box that contains all vertices of the Voro noi diagram in its interior, and attach the half infinite edges to the bounding box by updating the doubly connected edge list appropriately 9. Traverse the half edges of the doubly connected edge list to add the cell records and the pointers to and from them HANDLESETEVENT (pi) !.Search in T for the arc A vertically above Pi, and delete all circles events involv ingA from Q. 2. Replace the leaf of T that represents A with a subtree having three leaves: the middle leaf stores the new site Pi and the other two leaves store the site pj that was originally stored with A. Store the tuples [pi pj] and [pj. Pi] representing the new breakpoints at the two new internal nodes. Perform rebalancing operations on T if necessary. 3. Create new records in the Voronoi diagram structure for the two half edges separating V(pi) and V(p), which will be traced out by the two new breakpoints. 4. Check the triples of consecutive arcs involving one of the three new arcs. In sert the corresponding circle event only if the circle intersects the sweep line and the circle event isn't present yet in Q. HANDLECIRCLEEVENT(p1 ) 1. Search in T for the arc A vertically above p1 that is about to disappear, and de lete all circle events that involve A from Q. 2. Delete the leaf that represents A fromT. Update the tuples representing the breakpoints at the internal nodes. Perform rebalancing operations on T if neces sary. 3. Add the center of the circle causing the event as a vertex record in the Voronoi diagram structure and create two halfedge records corresponding to the new breakpoint of the voronoi diagram. Set the pointers between them appropriately. 4. Check the new triples of consecutive arcs that arise because of the disappear ance of A. Insert the corresponding .circle event into Q only if the circle intersects the sweep line and the circle event is not present yet in Q. 33
PAGE 38
6. Conclusion The literature is growing at such a rate that it was only possible to scratch the sur face of developments on the seven problems. Beyond these problems, there has been significant work in a number of other directions that went unmentioned here. This decade has seen computational geometry pulled in two opposing directions. On the one hand, the quest for improved theoretical results has caused the field to make contact with a mathematical community that may provide the tools for these improvements. In the other direction, the quest for relevance has continued con tact with researchers in continuous computational geometry, computer graphics, and solid modeling robotics. The combination of these pulls has led to an interesting and applicable intellectual discipline. What is most significant about the improvements is the level of "technology transfer". As the results mentioned here were developed, they became known beyond the theoretical computer science community from which discrete computational geometry first emerged. Practitioners were quick to recognize the potential relevance of this field and crossfertilization was a key theme in the second decade. The momentum of the field is such that this trend will clearly continue. As the field grows, the community of scholars involved reaches a size and excitement that neither Robin Forrest nor Mike Shamos could have realized when they named it. Obviously, more researchers will be attracted to the field, old problems will be solved, new problems will emerge, the process of turning theoretical results into practical algorithms will continue, ... as is the nature of intellectual disciplines. Beyond the obvious, there appears to be one need that is just being recognized and deserves mention. This is the concern about numerical issues when translating efficient algorithms such as those mentioned above into practice. There is a need to not only apply existing techniques but also to identify the limitations of such techniques and find methods of extending them. 34
PAGE 39
Appendix A. Source code for Neighborhood Search #include #include #include #define IS Full(ptr) (!(ptr)) typedef struct point *LIST; struct point { double px; double py; LIST next; } ; double x = 0, y = 0; LIST pointA = NULL, pointB = NULL, shourtest list = NULL; V****************************************************************a I I**************************** MakeList* * * * * * * * * * * *******I I V***************************************************************a void MakeList(LIST get) { LIST pNewCell; pNewCell =(LIST) malloc(sizeof(struct point)); pNewCell?px = get?px; pNewCell?py = get?py; pNewCell?next =shortest list; shortest list = pNewCell ; } 35
PAGE 40
H***************************************************************H merge****************************** II H***************************************************************H LIST merge(LIST list 1, LIST list2) { } if (list I ==NULL) return list 2; else if (list 2 ==NULL) return list 1 ; else if (list I <= { } list I merge(listl list2) ; return list!; else { } = merge(listl, return list2; H**************************************************************** I I**************************** split******************************** H*************************************************************** LIST split(LIST list) { } LIST pSecondPoint; if (list = = NULL) return NULL; else if next== NULL) return NULL; else { pSecondPoint = next = = return pSecondPoint; } 36
PAGE 41
1/****************************MergeSort**************************** LIST MergeSort(LIST list) { LIST SecondList; If (list = = NULL) return NULL; else if (list? next== NULL) return list; else { SecondList = split(list); return merge(MergeSort(list), MergeSort(SecondList)); } } 37
PAGE 42
fi**************************************************************** II************************* printlist* * * * * * * * * * * * * * * * * fi**************************************************************** void print list(LIST target) { } printf ("\nThe List Contains: "); while (target !::::: Null) { } printf ("\n"); printf("(%g, %g,) \n", target7px, target7py); target ;:::: target7next; printf ("\n\n") ; fi**************************************************************** II*************************** getpointB *** * * * * * * * * * * * * * fi**************************************************************** int getpointB( ) { LIST pNewPoint ; int i :::::0, c ;:::: 0 ; do { pNewPoint =(LIST) malloc(sizeof(struct point)); if (IS FULL(pNewPoint)) { fprintf(stderr, "The memory is full! !\n") ; exit (1) ; } printf("\nPlease enter a value for x:"); scanf("%1 f',&x); printf("\nPlease enter a value for y:"); scanf("% 1 f',&y); pNewPoint7px ;:::: x; pNewPoint7py = y; pNewPoint7next = pointB; pointB ;:::: pNewPoint; printf("\nDo you want to enter next value?"); printf("\nPlease enter 0 for no or 1 for yes: "); 38
PAGE 43
} scanf("%d", &c); i ++; }while (c != 0); return i; fl**************************************************************** 1/***************************getpointA ***************************** fl**************************************************************** int getpointA( ) { } LIST pNewPoint; int i = 0, c = 0; do { pNewPoint =(LIST) malloc(sizeof(struct point)); if (IS FULL(pNewPoint) ) { } fprintf (stderr, "The memory is full ! \n"); exit(l); printf("\nPlease enter a value for x:" ); scanf("% If' ,&x) ; printf( "\ nPlease enter a value for y:") : scanf( "%If", &y); pNewPoint7px = x ; pNewPoint7py = y ; pNewPoint7next = pointA; pointA = pNewPoint; printf(" \ nDo you want to enter next value?") ; printf(" \ nPle ase inter 0 for no or I for yes : ") ; scanf("%d ", &c) ; i + +; } while(c!=O); return i; 39
PAGE 44
a**************************************************************** a*************************** Extracheck******* * * * * * ** * * ** ** //**************************************************************** LIST Extra Check(LIST target], LIST current, doubled, double m) { int hit= 0; double d of hit = 0; LIST best = NULL; while (current!= NULL) { if ((targetl7pxm) <=current7px && current7px <= (targetl7px+m) && m)) )+( } { } (targetl7pym) < =current7py && current7py <= (targetl7py + d of hit= sqrt ((current7pxtargetl7px) (current7pxtargetl7px current7pytargetl7py )*( current7pytargetl7py)); if (d of hit< d) { } d = d of hit; best = current; hit= 1 ; current= current7next; else { } if (hit) return (best) ; else return NULL; return NULL; } 40
PAGE 45
V**************************************************************** //******* ***** *********Neighborhood Search************************ V**************************************************************** void Neighborhood Search(LIST target) { int i = 1, hit = 0 ; double d of hit = 0 d of shortest = 1 000 ; LIST ptemp, shortest node update ; Ptemp = pointA; While (ptemp = NULL) { } else { if ((target7pxi)<=ptemp7px && ptemp7px<= (target7px+i) && (target7pyi)<=ptemp7py && ptemp?py<= (target7py+i)) { d of hit = sqrt (( ptemp7px target7px)*(ptemp7pxtarget?px) + (ptemp7pytarget7py)*(ptemp7pytarget7py)); if ( d of hit < d of shortest) { } d of shortest = d of hit ; printf ("= = % g = = \n", d of shortest); shortest node = ptemp; hit= 1; ptemp = ptemp?next; if (!hit) i++; else { if((update =Extra Check(target, ptemp, d of shortest, i*sqrt (2)))! = Null) } } } } shortest node = update; break; MakeList( shortest node); 41
PAGE 46
a**************************************************************** II********************* invert*************************************** a**************************************************************** LIST invert(LIST lead) { } LIST middle, trail; middle = NULL; while (lead =NULL) { } trail = middle; middle =.lead; lead= lead?next; middle? next= trail; return middle; a**************************************************************** II********************* ReleaseMemory* * * ** * * * * ** * * * * * ** a**************************************************************** void Release Memory(LIST list) { } LIST dump; while ((dump= list)!= NULL) { } list = list? next; free( dump); 42
PAGE 47
II**************************************************************** ll*****************************tnain******************************* V**************************************************************** void main () { int numA = 0, numB = 0, index = 0, j = 0; LIST current; I* input and print how many points in the set A *I numA = getpointA ( ) ; printf(" \nThere are% d points in the A set. \n", numA); I* input and print how many points in the set B *I numB = getpointB( ) ; printf("\nThere are% d points in the B set. \n", numB); I* Sort the points in set A & set B *I pointA = MergeSort(pointA); pointB = MergeSort(pointB); I* Search the nearest neighbors from set A of all points of set B *I current = pointB; while (current!= NULL) { } Neighborhood Search( current); current = current?next ; shortest list= invert(shortest list); I Output of each points in set B & its nearest neighbor I while (pointB !=NULL && shortest list!= NULL) { printf ("B SET (%g, o/og ) ... ",pointB?px, pointB?py) ; pointB = pointB?next ; printf ("NEAREST ( %g, %g,) ... \n", shortest list?px, shortest list?py); shortest list= shortest list7next; } I Garbage Collection I Release Memory (pointA) ; Release Memory (pointB) ; 43
PAGE 48
Release Memory (shortest list); } 44
PAGE 49
References [Av85] [Ba84] [Ba84 [Be83] [Be75] [B078] [BS75] [Be76] Avis, K., On the partitionability of point sets in space, ACM symposium on computational geometry, 1985, 116120 Baird, H. Modelbased image matching using location, Ph.D. Thesis, Princeton University, May, 1984. Baumagart, B. G., Geometric modeling for computer vision, Stanford U., Computer Science Dept. Technical Report, STAN CS74463. BenOr, M,. Lower bounds for algebraic computation trees, Proceedings ofthe 151h ACM Symposium on the Theory of computing, 8086, 1983. Bentley, J. L., Multidimensional binary search trees used for associative searching, Communications of the ACM, 18, 509517, 1975. Bentley, J.L. and Ottman, Th., Algorithms for reporting and counting geometric intersections, CarnegieMellon University, August 1978. Bentley, J.L. and Stanat, D.F., Analysis of range searches in quad trees, Information processing Letters, 3, 170173, 1975. Bentley, J.L., Divide and conquer algorithms for closest point 45
PAGE 50
[Br79a] [Br79b] [CK70] [Ch84] [Ch85a] (Ch86b] [Ch86] [CD87] [Cl85] problems in multidimensional spaces, PhD thesis, University of North Carolina, 1976. Brown, K. Q., Voronoi diagrams from convex hulls, Information Processing Letters,_9, 223228,1979. Brown, K. Q., Geometric transformations for fast geometric algorithms, PhD thesis, Carnegie Mellon University, December, 1979. Chand, D.R. and Kapur, S.S., An algorithm for convex polytopes, JACM, 17,7886, 1970. Chazelle, B., Intersecting is easier than sorting, Proceedings of the 161h ACM Symposium on the Theory of computing, 125134, 1984 also JCSS,. Chazelle, B., How to search in history, Information control, 64, 7799, 1985. Chazelle, B., Fast searching in a real algebraic manifold with applications to geometric complexity, Proc. CAAP'85, LNCS, SpringerVerlag, 1985. Chazelle, B., Filtering search: a new approach to query answering, SIAM J on Computing, 15, 703722, 1986. Chazelle, B. and Dobkin, D., Intersection of convex objects, JACM, 34, 127,1987. Clarkson, K., A probabilistic algorithm for the post office problem, 46
PAGE 51
[C186] [Co85] [De86] [DE84] [DEY84] [DK83] [DK85] [DL87] Proceedings of the 17th ACM Symposium on the Theory of computing, 175185, 1985. Clarkson, K., Further applications of random sampling to computational geometry, Proceedings ofthe 18th ACM Symposium on the Theory of Computing, 414423, 1986. Cole, R., Partitioning point sets in 4 dimensions, ICALP '85, SpringerVerlag Lecture Notes in Computer Science, 194, 111119,1985. Devai, F., Quadratic bounds for hidden line elimination, ACM symposium on computational geometry, 269275, 1986. Dobkin, D.P. and Edelsbrunner. H.E., Space searching for intersecting objects, IEEE Symposium on foundations of computer science, 387392, 1984, J of Algorithms, Dobkin, D.P., Ede1sbrunner, H.E., and Yao, F.F., A 3space partition and its applications, unpublished manuscript. Dobkin, D.P. and Kirkpatrick, D.G., Fast detection ofpolyhedral intersections, Theoretical computer science, 27, 241352, 1983. Dobkin, D.P. and Kirkpatrick, D.G., A1inear slgorithm for determining the separation of convex polyhedra, J of Algorithms. 6, 381392, 1985. Dobkin, D.P. and Laszlo, M.J., Primitives for the manipulation of threedimensional subdivisions, ACM Symposium on the Theory 47
PAGE 52
[DL74] [DL76] [DL79) [DR80] [DS87] [DSV86] [Dy84] [Ed77] [EGS87] [EMK82] of Computing, geometry, 1987, Dobkin, D.P. and Lipton, R.J., On some generalizations of binary search, Proceedings ofthe 61h ACM Symposium on Theory of Computing, 31 0316, 197 4. Dobkin, D.P. and Lipton, R.J., Multidimensional searching, SIAM Journal on computing, 5, 181186, 1976. Dobkin, D.P. and Lipton, R.J., On the complexity of computations under varying sets of primitives, JCSS, 18, 1 ( 1979),8691. Dobkin, D.P. and Reiss, S.P On the complexity of linear programming, Theoretical Computer Science, 11 ,3, 118, 1980. Dobkin, D.P. and Souvaine, D., Computational geometrya user's guide, in Advances in Robotics, Vol. 1. Edited by J. T. Schwartz and C.K. Yap, Lawrence Erlbaum Associates, 1987. Dobkin, D.p., Souvaine, D., and VanWyk, C., Decomposition and intersection f simple splinegons, submitted for publication. Dyer, M., Linear time algoritluns for twoand threevariable linear programs, SIAM Journal on Computing, 13, 3145, 1984. Eddy, W. Anew convex hull algoritlun for planar sets, ACM Transactions on Mathematical Software, 3,4, 398403, 1977. Edelsbrunner, H.E., Guibas, L.J. and Stolfi, J. Optimal point location in a monotone subdivision, SIAM Jon Computing, Edelsbrunner, H., Maurer, H. and Kirkpatrick, D., Polygonal 48
PAGE 53
[EW82] [Er46] [FP81] [FB74] [F176] [Fo68] [Fo86] [Fr81] [Gr71] [Gr83] intersection searching, Information Processing Letters, 14, 7479, 1982. Edelsbrunner, H. and Welzl, E., Halfplane range search in linear space and O(n0.695) query time, unpublished manuscript. Erdos, P., On sets of distances of n points, American Mathematic! Monthly, 53, 248250, 1946. Faux, I.E. and Pratt, M.J ., Computational geometry for design and manufacture, Ellis Horwood Ltd., Chichester, 1981. Finkel, R.A. and Bentley, J.L., Quad trees: a data structure for retrieval on .composite keys, Acta informatica, 4, 19,1974. Floyd, R., private communication. Forrest, A.R., Curves and surfaces for computeraided design, PhD thesis, Cambridge University, July, 1968. Fortune, S., A sweepline algorithm for Voronoi diagrams, ACM symposium on tomputational geometry, 313322, 1986. Fredman, M.L., Lower bounds on the complexity of some optimal data structures, SIAM Jon computing, I 0, 110,1981. Graham, R.L., An efficient algorithm for determining the convex hull of a finite planar set, Information Processing Letters, 1, 132133, 1972. Graham, R.L. and Yao, F.F., Finding the convex hull of a simple polygon, J. of Alogrithms, 4, 324331, 1983. 49
PAGE 54
[Gr56] [GS87] [GS87] [HW86] [HZ82] [Hu78] [HS79] (JB87] [KR84] Grunbaum, B., Aproof ofVaszonyi's conjecture, Bulletin Res, Council of Israel, 6, 7778, 1956. Guibas, L.J. and Stolfi, J., Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams, ACM Trans. On Graphics, 2, 75123, 1985. Guibas, L.J. and Sharir, M., private communication. Haussler, K. and Welzl, E., Epsilonnets and simplex range queries, ACM symposium on computational geometry, 1986, 6171, Discrete and Computational Geometry, Hubschman, H. and Zucker, S., Frametoframe coherence and the hidden surface computation: constraints for a convex world, ACM Transactions on Graphics, 1,2, April, 1982, 129162. Hunter, G.M., Efficient computation and data structures for graphics, PhD thesis, Princeton University, 1978. Hunter, G.M. and Steiglitz, K., Operations on images using quad trees, IEEE Transactions on Pattern Analysis and Machine Intelligence, I, 145153, 1979. Jameson, A. and Baker, T., Improvements to the aircraft Euler method, AIAA 251h Aerospace sciences meeting, paper AIAA870452, 1987. Kalin, R. and Ramsey, N., A new analysis and comparison 50
PAGE 55
[Ka72] [Ki83] [KS86] [Kn73] [La87] [Ke83] [LP79] of algorithms for the planar segment intersection problem, unpublished manuscript. Karp, R.M. [1972], "Reducibility among combinatorial problems," in R.E. Miller and J.W. Thatcher (eds.), Complexity of Computer Computations, Plenum Press, New York, 85103. (1.5; 3.1; 5.2; 7.1; 7.4; Al.1; Al.2; Al.3; A2.1; A2.2; A3.1; A3.2; A5.1; A6; Al0.1) Kirkpatrick, D.G., Optimal search in planar subdivisions, SIAM Jon Computing, 12,2835, 1983. Kirkpatrick, D.G. and Seidel, R., The ultimate convex hull algorithm?, SIAM Jon computing, 15, 287299, 1986. Knuth, D.E., The art of computer programming, Vol. 3 Sorting and Searching, AddisonWesley, Reading, MA, 1973. Laszlo, M.J., Manipulating threedimensional subdivisions, PhD thesis, Princeton University, 1987, to appear. Lee, D.T., On finding the convex hull of a simple polygon, International Journal of Computing and Information Sciences, 12, 8798, 1983. Lee, D.T. and Preparata, F.P., Location of a point in a planar subdivision and its applications, SIAM Journal on Computing, 6,3, 594606, 1979. 51
PAGE 56
[LW77] [LT77a] [LT77b] [MS82] [Me83] [Me84] [MP78] [My85] Lee, D.T. and Wong, C.K., Worstcase analysis for region and partial region searches in multidimensional binary trees and balanced quad trees, Acta Informatica, 9, 2329, 1977. Lipton, R.J. and Tarjan, R.E., A separator theorem for planar graphs, Conference on theoretical computer science, 11 0, 1977. Lipton, R.J. and Tarjan, R.E., Applications of a planar separator the'orem, IEEE Symposium on foundations of computer science 162170, 1977. Mairson, H.G. and Stolfi, J., unpublished manuscript. Megiddo, N., Linear time algorithm for linear programming in R3 and related problems, SIAM Jon Computing, 12,759776, 1983. Megiddo, N ., Linear programming in linear time when the dimensionisfixed,JACM, 31,114127,1984. Muller, D.E. and Preparata, F,P., Finding the intersection of two conves polyhedra, Theoretical Computer Science, 7, 217236, 1978. Myers, W.W., An O(E logE+ I) expected time algorithm for the planar segment intersection problem, SIAM Jon, Computing, 14 625637, 1985. 52
PAGE 57
[Pa87] [PH77] (PS85] [Re77] (Ro85] [Sa84] [ST86] (Se81] [Se86] Panduranga, E.S., Reflections in curved surfaces, PhD Thesis, Princeton U.,. Preparata, F.P. and Hong, S.J., Convex hulls of finite sets of points in two and three dimensions, Communications of the ACM, 20, 2, 8793, 1977. Preparata, F.P. and Shamos, M.l., Computational Geometry SpringerVerlag, 1985. Reiss, S.P., private communication. Rosenstiehl, P., private communication. Samet, H., The quadtree and related hierarvhical data structures, ACMComputingSurveys, 16, 187260, 1984. Samak, N. and Tarjan, R.E., Planar point location using persistent search trees, Communications of the ACM, 29, 669679, 1986. Seidel, R., A convex hull algorithm optimal for point sets in even dimensions, Technical Report, University of British Columbia, 1981. Seidel, R., Outputsize sensitive algorithms for constructive problems in computational geometry, PhD thesis, Cornell University, 1986. 53
PAGE 58
[Sh74] [Sh75] [SH75] [Sw85] [Wa69] [Wa70] [Wi82] [YY85] Shamos, M.I., private communication. Shamos, M.I., Geometric complexity, Proceedings of the ih ACM Symposium on the Theory ofComputing, 224233, 1975. Shamos, M.l. and Hoey, D., Geometric intersection problems, IEEE Symposium on foundations of computer science, 208215,1976. Swart, G., Finding the convex hull facet by facet, J of Algorithm 6, 1748, 1985. Warnock, J .E., A hiddensurface algorthim for computer generated halftone pictures, University ofUtah Computer Science Department, TR 415, 1969. Watloms. G.S., A realtime visible surface algorithm, University ofUtah Computer Science Department, UTECCSc70101, June, 1970. Wilard, D., Polygon retrieval, SIAM Journal on Computing, 11, 149165, 1982. Yao, A. and Yao, F., A general approach to ddimensional geometric queries, Proceedings of the 1 ih ACM Symposium on the Theory of Computing, 163169, 1985. 54
