This thesis for the Master of Science
degree by
Qibin Tao
has been approved
Christopher Smith
/ 2~] /00
Date
Tao, Qibin (M.S., Computer Science)
Delaunay Triangulation Implemented in Java
Thesis directed by associate Professor McConnell, Ross
ABSTRACT
In this thesis, I research the Delaunay Triangulation problem of Computational
Geometry, and mainly discuss the data structures and algorithms used to compute it.
The Randomized Incremental Algorithm is implemented in Java, and the expected
running time is bounded in 0(n Inn ).
This abstract accurately represents the content of the candidates thesis. I recommend
its publication.
Signed
McConnell, Ross
in
ACKNOWLEDGMENT
My thanks to my advisor, Ross McConnell, for his patience and help during the pass
two years.
CONTENTS
Chapter
1. Introduction........................................................1
2. The Properties of Delaunay Triangulation...........................3
2.1 Triangulations of Planar Point Set................................3
2.2 Delaunay Triangulation............................................4
2.3 The Relationship of Voronoi Diagram and Delaunay Triangulation....7
3. General Geometry Data Structure....................................11
3.1 DoublyConnected Edge List.........................................11
3.2 Objects Represented in Java.......................................14
4. Algorithms.........................................................16
4.1 Randomized Incremental Algorithm..................................16
4.2 Randomly Generating Points........................................18
4.3 Initial Incremental Triangle......................................20
4.4 Point Location Query..............................................22
v
4.5 Flipping.........................................................22
4.6 The Correctness of the Algorithm.................................29
5. Analysis...........................................................36
6. Application........................................................40
6.1 Representation of Maps............................................40
6.2 Finite Element Analysis..........;...............................41
6.3 The Nearest Neighbor.............................................41
6.4 Minimum Spanning Tree............................................42
Appendix
A. Delaunay Triangulaion.............................................45
B. Graph.............................................................67
References............................................................71
vi
DELAUNAY TRIANGULATION IMPLEMENTED IN JAVA
B.S., Nanjing University of Aeronautics and Astronautics, 1985
M.S., Nanjing University of Aeronautics and Astronautics, 1988
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
by
Qibin Tao
2000
1. Introduction
In computational geometry, the triangulation of planar point set is a very important
topic like the convex hull. What is the triangulation of planar point set? For a set of
points P, we add edges to connect two points (vertices) until no edge can be added
without crossing other edges. It is obvious that there exist many triangulations for a
set of P, but some triangulations cant provide any information about the relations
among these points such as the closest neighbors and the shortest path passing
through P. So it is an interesting problem that what kind triangulation can provide
enough information about the relations of sites (points).
Voronoi [2] first discussed a problem as the following description: if there is a set of
points (sites) P in the plane, for an arbitrary point q, can we find a point peP that is
closest to q? Voronoi assigned very point to the nearest site in a region, so that every
site was bounded in a certain region. This configuration is called Voronoi Diagram,
and it provides all kind information about trading areas of sites, and their relations.
Delaunay in 1935 [2] found that when the dual graph is drawn with straight lines, it
produces a planar triangulation of sites P, now called the Delaunay Triangulation.
1
Delaunay Triangulation from another viewpoint researches the problem and more
directly provides the information of point relations. In Delaunay Triangulation, we
use fat triangles to triangulate a set of points, so that no triangle with more sharp and
long shape can be found in the triangulation. This triangulation keeps the points on a
triangle have a more strong relationship, so that it provides more effective
applications.
During the development of computational geometry, computer scientists try to find
more effective algorithms for this problem. There are several algorithms to compute
Delaunay Triangulation such as Intersection of Halfplanes, Incremental
Construction, Divide and Conquer. Most of the algorithms need to compute Voronoi
Diagram first and then translate it into Delaunay Triangulation, so these algorithms
make their computing ineffective or make their implementing complicated.
The algorithm used in this thesis is a simple and clever way, and it flips some edges
of a triangulation so that there is no more sharp and long triangle existing. In order to
make the running time not worse than the algorithms mentioned above, we must
build a suitable data structure. In the algorithm, we introduce Incremental
Configuration to build a data structure, and bound the expected running time
2
in O (n Inn).
This algorithm is implemented in Java version in my thesis, and the software can be
developed to possess real applications.
3
2. The properties of Delaunay Triangulation
2.1 Triangulation of Planar Point Sets
We first see a problem of planar point triangulation. Let P = {Pi, P2,,,,,, Pn }be a set
of points in a plane, can we get a subdivision from this set of points so that in this
subdivision no edge connecting two vertices can be added without intersecting other
edges. We call this subdivision as a maximal planar subdivision. For a set of points
P, if there is a maximal planar subdivision whose vertex set is P, this subdivision is
called a triangulation of P.
It is obvious that there exist triangulations for any point set. It can be proved that a
triangulation of P consists of triangles. For a set of points, a convex hull can be
computed that bounds all points, and the edges constructing the boundary of convex
hull keep the same in any triangulation. It is obvious that the segment connecting two
consecutive points on the convex hull does not cross any other edges of
triangulations, so the edges on the boundary of convex hull must appear in any
triangulation. The points in the convex hull can be triangulated into several polygons,
and a polygon must consist of triangles [1]. Theorem 2.1 indicates that the number of
triangles created in any triangulation keep the same and is bounded in 0(n). This
conclusion will be used in the analysis of running time afterwards.
4
Theorem 2.1 Let P be a set of n points in the plane, not all collinear, and let k denote
the number of points in P that lie on the boundary of the convex hull of P. Then any
triangulation of P has 2n2k triangles and 3n3k edges.
The proof is easy, and Euler's formula tells that:
nne+nf=2.
n, ne, and nf respectively present the number of points, edges and faces. Let m denote
the number of triangles of a triangulation, and then the number of faces of the
triangulation is m+1 (including unbounded faces). Each triangle has three edges, and
adjacent faces use each incident edge as a boundary twice. So the total number of
edges should be (3m + k) / 2. Taking these values into the formula, we get the
conclusion of Theorem 2.1.
2.2 The Delaunay Triangulation
There are many ways to triangulate a set of points according to different rules, but
some triangulations produce some triangles that have long and sharp shapes. These
triangulations are not suitable to some applications. So scientists try to find an
optimal triangulation that meets a special goal. Delaunay Triangulation is a good
triangulation according to the angle optimal goal to avoid appearing of these
5
triangles. Let us introduce what is an illegal or legal edge, and how to flip an illegal
edge in the following contents.
Suppose there is a triangulation, and it produces m triangles. We arrange the 3m
angles of these triangles by increasing values. Let A = { ai, a2,,, an } be the
sequence of angles; hence ai < aj, if i
another triangulation that produce an anglevector A = { ar,a r,,, an }, and there
exists an index i with 1 < i ^ 3m such that
aj = aj for all j ai.
We say that A is larger than A lexicographically. If A is larger than the anglevector
of any other triangulation, we call this triangulation with angleoptimal. Figure 2.1
illustrates how to obtain an angle optimal triangulation from an arbitrary
triangulation.
6
Pi
Pi
Figure 2.1 An Illegal Edge
In Figure 2.1 (a), the edge is PiPj incident to two triangles PiPjPk and PiPjPi, and the
two triangles have six angles {ai, a2, ,a6}. In Figure 2.1 (b), the edge PkPi is
incident to two triangles PiPkPi and PkPiPj, and the two triangles have six angles {ar,
or, ,a6}. We call the edgePiPj an illegal edge if:
min ai < min ai, 1< i < 6
ai presents the angles of the left quadrilateral, and ai presents the angles of the right
one. Before giving a simple and clear way to decide whether an edge is illegal,
we induce a theorem called Thales Theorem [1].
7
Theorem 2.2 Let C be a circle. Suppose that p, q, a, and b lie on C, that r lies inside
C, and that s lies outside C. Then
Z arb > Z apq = Z aqb > Zasb.
s
Figure 2.2 Thales Theorem
So we can get Lemma 2.3 by combining the definition of illegal edge and Theorem
2.2, and will prove its correctness in Chapter 2.3.
Lemma 2.3 Let edge PiPj be incident to triangles PiPjPk and PiPjPi, and C be the
circle through Pi, Pj, and Pk (see Figure 2.3). The edge PiPj is illegal if only if the
point Pi lies in the interior of C.
8
It needs to be demonstrated that Lemma is symmetric, and the conclusion is the same
if a circle is through Pi, Pj, and Pi. Note that the two triangles conform quadrilateral,
and otherwise Pi is certainly out of C.
2.3 The relationship of Voronoi Diagram
and Delaunay Triangulation
In order to introduce some properties of Delaunay Triangulation, we first discuss the
relationship of Voronoi Diagram and Delaunay Triangulation. The Voronoi diagram
is the subdivision of the plane divided into n regions for n points, one for each site,
such that the region of a site contains all points for which this site is the closest site.
9
o
o
o
Figure 2.4 Voronoi Diagram
In Figure 2.4, there are 4 sites, and 4 regions bounded by some edges are generated.
These edges are called the edges of Voronoi Diagram. If we link two sites between
which there exists an edge to bisect them in Voronoi Diagram, then we can produce
& Delaunay Triangulation, showing as Figure 2.5. So Voronoi Diagram and
Delaunay Triangulation are dual graphs, and they can translate into each other.
Figure 2.5 Delaunay Triangulation
10
From the relationship of the Voronoi Diagram and Delaunay Triangulation, we can
get a lot of properties of Delaunay Triangulation. Lemma 2.3 gives a verification for
illegal edges, but in the implementation we must use math formula to represent it.
Usually complicated matrix operations are needed to compute this verification, but
from the relationship of voronoi Diagram and Delaunay triangulation we can give a
simple and direct verification for illegal edges.
Pi
Figure 2.6 Illegal Edge Verification
Seeing Figure 2.6, we draw a perpendicular bisector Li of the edge PiPk, and L2 of the
PkPj. We can get an intersection O of Li and L2, and the point O is the center of circle
that is through the points Pi, Pj, and Pk. Then we draw a bisector L3 of PkPi. If the
11
point Pi is outside of the circle C, the points O and Pk must be on the same side of L3.
L3 is certainly not an edge of the Voronoi Diagram. So there do not exist an edge
between points Pi and Pk, and the edge PiPj does not need to be flipped. If the point Pi
is in the interior of the circle C, the point O and Pk must be on the different side of
L3. Then L3 must be an edge of Voronoi Diagram bounding the site Pk. So there must
be an edge between points Pk and Pi, and we must flip the edge PiPj.
Through demonstrating this verification, we also give a proof for Lemma 2.3.
12
3. General Geometry Data Structure
3.1 Doubly Connecting Edge List
When dealing with the geometrical problems, we must build a suitable data structure.
For example, a map is a planar subdivision composed of a set of segments. In order
to decide the location of a point, we need some operations to support point location
query. In general case, when we start from any one segment, we can walk to any one
segment of this subdivision. So a suitable data structure must be built to ensure these
operations.
In a planar subdivision, we use three elements to describe its configuration, vertices,
edges and faces. The embedding of a node of graph is called a vertex, the embedding
of an arc is called an edge, and a maximal connected subset of the plane that doesn't
contain a point on an edge or a vertex is called a face. Usually we use the doubly
linked list to support the operations mentioned above. A doublyconnected edge list
contains three classes to record the information.
A vertex object includes the coordinates of a point, and in a planar subdivision, these
coordinates are the values of x and y. An edge object includes the direction of an
edge. In subdivision an edge is called an oriented edge that always appears as a
13
boundary of a face. An edge usually bounds two faces so that we can regard the
different sides of an edge as two distinct halfedges, and each of these halfedges is
responsible to one face. Suppose the counterclockwise traversal as our walking
direction. When an observer walks along an edge, the face which lies to its left of the
observer is respect to this halfedge and the direction on which the observer walks is
the direction of this halfedge. So each halfedge is an oriented edge and has an
original vertex and destination vertex.
For each face, although it can be bounded or unbounded by some edges, it must have
at least one edge as its boundary. For this reason, a face can be represented by the
halfedges, and we know one boundary edge of a face, then we can get all of the
boundary edges. In order to achieve this goal, we build a doublyconnected edge list
on which each edge has two pointers to its next edge and previous edge of the same
face boundary. So we can walk along the boundary of a face from any its
representing edge, and keep the face on our left always.
The method mentioned above only resolves the problem of walking around one face.
How to walk from one face to another face is another important problem, and it can
ensure the searching operation going through the whole map. An edge can be
14
represented as two half edges, and each of them bounds one of two separated faces.
These two faces are adjacent so that we call these two half edges the twins. They
actually belong to one edge but appear as different direction. So each edge has a
pointer to its twin. The content above can be illustrated in Figure 3.1.
Vs
Figure 3.1 Doublyconnected Edge List
In the Figure 3.1, it can be seen that how a planar subdivision is represented in
vertices, edges and faces.
Vertices: Vi, V2, V3, V4, Vs.
Edges: ei2, e23, e3i, ei3, e34, e4s, esi.
Faces: fi, fz.
15
The face fi can be represented in anyone of edges ei2, e23, and e3i, and the face h can
be represented by anyone of edges ei3, es4, e45, esi. The edge e3i is the twin of edge
ei3. The original vertex of edge e3i is V3, and the destination vertex is Vi. Its twin
edge ei3's original vertex is V3 the destination vertex is Vi. The edge e3i's next edge
is ei2, and its previous edge is e23.
3.2 Objects Represented in Java
In Java language, we can use three classes to record the information discussed above:
Class vertex {
intx;
inty;
}
Class edge{
vertex original;
vertex destination;
edge next;
16
edge previous;
edge twin;
}
Class triangle{
edge currentEdge;
}
Delaunay Triangulation, each face has a fixed configuration: triangle, so its
boundary is only composed of three edges. In this view, we can simplify Class
triangle as the following:
Class triangle {
edge triangleEdgel;
edge triangleEdge2;
edge triangleEdge3;
}
When a triangle is created in Delaunay Triangulation, we only update its three edges.
17
4. Algorithms
In the computing Delauny Triangulation, there are two main ways used frequently:
computing Voronoi Diagram and translating it into Delaunay Triangulation, or
computing it directly by Randomized Incremental Algorithm. Though there are other
algorithms, their running time is not effective or their implementations are too
complicated. In my paper, the randomized incremental algorithm is chosen for
implementing.
4.1 Randomized Incremental Algorithm
In Computational Geometry, there are two algorithms used frequently: the sweeping
line and randomized incremental algorithms. The sweeping line strategy is that: a
line is moved to the map from a far place, and the corresponding data structure is
updated when this line meets the events (the points or segments of map). The
randomized incremental algorithm uses another kind strategy to build its data
structure. Its strategy is as the following description: it starts from the simplest
configuration, and the parts constructing the whole configuration are added randomly
step by step. These parts can be points, segments or configurations of low dimension.
On each step, the current configuration need to be adjusted in order to meet some
constrains. Because the parts are added randomly, in some cases the running time for
18
adjustment is in the worst case, in some cases in the best case, or between the both of
them. Usually we use the expected running time to estimate them. The randomized
incremental algorithm has its advantages because it keeps the information of
constructing configuration. So it is often used to resolve the problems involved in
point location query.
For Delauny Triangulation, it starts from a triangle that is built far from the sample
points, and then we add all of sample points step by step. We define the sample point
set as the following: P{Pi, P2,.Pn), and Delaunary Triangulation created for the
ith step as Dl
Algorithm main( )
input: A sample point array P[1,2,....n].
output Vertices of each triangle in Delaunay Triangulation.
1. GetSartTriangleO generates a start triangle.
2. For i from 1 to n, n is the number of points.
3. RandomPop() pops a point randomly.
4. Query the location of the point, and find the triangle that contains this
point.
19
5. Split the triangle that contains the point.
6 Flip all illegal edges.
7. Discard the edges incident to the start triangle.
The following content will discuss the details of each step.
4.2 Randomly Generating Points
When we get a map, we usually collect sample points from its surface into an array
according to a certain rule. For example we maybe measure its heights along a row
by row, so the arrangement of these points may have a relationship. In order to
destroy this relationship and make the arrangement of points seem random, we need
to arrange the order of these points again. To ensure the running time of arrangement
not to exceed the total running time, we can use an algorithm as the following
description.
For n sample points, we put these points into an array Vertex [n]. Then we use a
random number generator to produce an index, and pop the element with this index
from this array.
20
Algorithm RandomPop( AQ, BQ, n)
1. ForifromOton1.
2. A[i] = random( ), produces a random real number between 0 and 1.
3. QuickSort(float[] A, int[] B, 0, n1) sorts A[] according to its increasing
order, and B[] records the index of the array after sorting.
This algorithm pops a point from the point array and the index of this point is
generated randomly. This algorithm is showed as the following Table 4.1.
Before Sorting
Index 1 2 3 4
Random Number 0.985 0.523 0.784 0.236
After Sorting
Index 4 2 3 1
Random Number 0.236 0.523 0.784 0.985
Table 4.1 Random Number Generator
21
For the input of point array, we generate a random real number corresponding to
each index of the array, and then sort these numbers according to their increasing
order. So we get a new arrangement of die index, and this index has a random order.
In order to keep the running time on (9 (n Inn), we use the randomized Quick
Sorting, and its expected running time is O (n Inn).
4.3 Initial Incremental Triangle
In the incremental algorithms, we generate an initial configuration at first from which
a final configuration will be developed. Because in Delaunay Triangulation, an
incremental algorithm is involved in the operation of point location query, we must
give a boundary in which all of points locate. We choose a triangle as the start
configuration that bounds all of points.
This initial triangle is added extra, and it must be discarded when Delaunary
Triangulation is finished. So we should choose this triangle correctly. In a theoretical
view, it is better that the top points of the triangle are farther from the points which
will be inserted inside it, but it causes huge coordinates comparing with the
coordinates of the points. We explain how to choose an initial triangle with suitable
coordinates from the following theorem.
22
Theorem 4.1 Let P be a set of n points in the plane, not collinear, then these points
must lie on the boundary of the convex hull of P [1].
Pi
Figure 4.2 Initial Triangle
It seems that the introduction of the initial triangle does not affect Delaunay
Triangulation of P. It can be seen that all points are bounded in a convex hull from
the Figure 4.2. If we set a convenience start triangle to bound this Delaunay
Triangulation, its vertices only link the vertices on the convex hull boundary to
construct a new Delaunay Triangulation. The flipping of edges does not occur on the
23
configuration composed of P, so we can discard those edges outside Delaunay
Triangulation of P after finishing the construction.
4.4 Point Location Query
Point Location Query is a typical and difficult topic in computational geometry, and
many problems involve in it. Given a point in a planar subdivision, we want to know
we need to look for a
in a face is illustrated
in Figure 4.3.
where it is. It means that if we know the coordinate of a point,
face that contains this point. How to verify wheather a point is
Figure 4.3 Point Location Query
24
When we walk along the boundary edges of face f in counterclockwise, the point is
always on our left. We can use a math formula to present this verification. For an
edge e, its direction is from point Pi to P2, and the coordinates of Pi and P2 are Xi,
Xi, X2, Y2. A point P's coordinates are X and Y. Then we walk along the edge e in its
direction to keep the point P on the left, and the coordinates of P, Pi and P2 must
meet the condition as the following.
( XXi) / (YYi) > ( X2 Xi) / (Y2 Yi).
This formula includes the case of point on the edge when the two sides of equation
are equal each other.
In Delaunay Triangulation, when we insert a point in the current configuration, we
first look for a triangle that contains this point. Because the incremental algorithm
constructs a configuration step by step, we have a chance to store the information of
configuration for each step. We show how to construct a data structure for querying
by the following graphs in Figure 4.4 and Figure 4.5.
25
Figure 4.4 Data Structure for Querying
Du is an directed acyclic graph when Pii is inserted. Its bottom nodes (leaves)
represent all of triangles of the current Delauray Triangulation. Then we inserted
point Pi, and query which triangle contains this point. For example it lies on the
triangle di, and di will be split into three triangles. In this case di will become an
internal node of the graph, and three new leaves d3, d4, and ds will be created which
correspond to di. Then we check each edge of di to decide wheather they are illegal.
If an edge is illegal, we need to flip it. (The next section will discuss how to flip a
edge).
26
For example the edge PiPj of dl is illegal (see Figure 4.5), so we need flip the edge
PiPj. We link point Pr to Pk to create a new edge PrPk instead of the edge PiPj. This
flipping involves two triangles ds and d2, and creates two new triangles d6 and d7. So
the two triangles make the leaves ds and 62 become internal nodes and themselves
become two new leaves. This operation continues until all of edges are legal. At last
this flipping create a new Delaunay Triangulation represented by a data structure Di.
Figure 4.5 Building a Data Structure
27
Pi
di
ch
Figure 4.5 Building a Data Structure
This data structure is a directed acyclic graph, and its leaves represent all of triangles
of the current Delaunary triangulation. Its internal nodes represent those triangles
created in the early stages and destroyed now. Because the current configuration is
built from the initial triangle, this data structure makes point location query easy.
When we insert a point into the configuration, we search from the root of data
structure (the initial triangle) to its leaves to find a leaf that contains this point. The
searching process is very easy in this graph, and on each internal node we decide
which child contains this point until we achieve a leaf.
28
In Java, we can make a class triangle like the following:
class triangle {
edge triEdgel;
edge triEdge2;
edge triEdge3;
triangle child 1;
triangle child2;
triangle child3;
}
In this class we define three children for the triangles created by it, because the
splitting operation creates three triangles and the flipping operation creates two
triangles. When we update the class triangle, we start from its child 1 to child3 no
matter what there is the splitting operation or flipping operation.
This data structure might search a node repeatedly because one node may have
several parents to point to it. Avoiding this case, we add a flag in the class triangle
29
present the information of visiting. flag=0 means this node is not visited until now,
and flag=l mean it has been visited.
class triangle {
int flag;
}
4.5 Flipping
We have discussed in the chapter 2.2, a simplest way of Delaunay Triangulation is
that we flip all illegal edges until no illegal edge existing. In the incremental
algorithm, a point is inserted inside the Delaunay Triangulation created on the early
stage, so that it may causes some edges illegal and destroys the property of Delaunay
Triangulation. Let us analyze the procedure of this way. When a new point locates in
a triangle of the configuration, we know that there are at least three new edges to link
this point (three degree). We will prove in the next chapter that the edges created
during the flipping operations must be the edges of the current Delaunay
Triangulation. Showing as Figure 4.6.
30
Pi
Pi
Pk
Pi
Pj
Pj
Figure 4.6 Filpping
When the point Pr is inserted in the triangle PiPjPk, it creates the three new edges PrPi,
PrPj, and PrPk. We need to check wheather the edges PiPj, PjPk and PkPi are illegal.
For example the edge PiPj is illegal, and a new edge PrPi will be created through
flipping it Then we need to check the edges PiPi and PjPi. The iteration continues
until the respected edges are legal.
Algorithms:LegalizedEdge(Pr, PiPj, T)
1. Pr is inserted into T, and PiPj is the edge of T that may need to be flipped.
2. if PiPj is illegal.
31
3. then Let PiPjPi be the triangle adjacent to PiPjPk along the edge PiPj.
4. Replace PiPj with PrPi.
5. LegalizedEdge(Pr,PiPi,T).
6. LegalizedEdge(Pr,PjPi,T).
4.6 The Correctness of the Algorithm
Concluding the discussions above, we can give a precise description of the main
algorithm.
Pi
(a)
Pi
Figure 4.6 Splitting
Algorithm main()
1. GetStartTriangle() gives a suitable initial triangle that contains all of sample
points, and it can be discarded at the last stage.
32
2. for i from 0 to n1.
3. RandomPop(), it pops an element Pr from the array of points.
4. Insert Pr into the current configuration T.
5. Query a triangle of T that contains this point' Pr.
6. if Pr lies in the interior of triangle PiPjPk(see Figure 4.7 (a)),
7. split the triangle PiPjPk into three triangles.
8. LegalizedEdge (Pr, PiPj, T).
9. LegalizedEdge (Pr, PjPk, T ).
10. LegalizedEdge (Pr, PkPi, T ).
11. else Pr lies on a edge of the triangle (PiPj), create two new edges PrPk and
PrPi, so that the two triangles incident to the edge PiPj are split into
four triangles(see Figure 47 (b)).
12. LegalizedEdge (Pr, PiPi, T).
13. LegalizedEdge (Pr, PiPj, T).
14. LegalizedEdge ( Pr, PjPk,T).
15. LegalizedEdge ( Pr, PkPi,T).
16. Discard the initial triangle.
33
This recursive algorithm seems so beautiful, but we need to prove its correctness.
The following content shows the proofs.
1. When the point Pr is inserted in the current Delaunay Triangulation, what we need
to do is embedding some edges on this point and ensuring the properties of Delaunay
Triangulation. In LegalizedEdge( Pr, PiPj, T ), it creates two new edges PrPi and PrPj
if there is not any point in the interior of the circle C through Pi, Pj and Pr.
Otherwise, it will create three new edges PrPi, PrPj, and PrPi (dash lines in Figure 4.8)
and delete the edge PiPj.
Pk
Pj
Figure 4.8
34
For edges PiPi and PjPi, we recursively call LegalizedEdge(Pr, PiPi, T) and
LegalizedEdge(Pr, PjPi, T), and it will produce a series of new edges until it gets a
legal edge. This algorithm verifies if the point Pi incident to the edge PiPj lies in the
interior of C. If Pi is inside C, the two triangles must form a convex quadrilateral. So
convex quadrilaterals make the new edges created by flipping not across other edges,
and make these new edges be the edges of Delaunay Triangulation created newly. In
the following, we will prove that the flipping operation will stop until it meets a legal
edge.
2. We prove that if a edge in LegalizedEdge( , ) is legal edge, this recursion
certainly stops. Seeing Figure 4.9, if PiPj is legal, Pi (the point Pi is a point of the
triangle pointed by the edge PiPjs twin) should be outside of the circumcircle C.
We can prove that there is not any other point in C. We assume that there is a point
Pg existing in C, also know that Pi is outside of C. According Thales' Theorem, the
angle PiPgPj > PiPiPj. If we shrink circumcircle C to C' which is through points Pi, Pj,
and Pi, we find that point Pg is still in the interior of C' because of angle PiPgPj >
PiPiPj. Otherwise the triangle PiPjPi is a triangle of Delaunay Triangulation, and the
circumcircle C' should not exist any other point So the assumption above accuses
contradictable to the properties of Delaunay Triangulation.
35
Pi
Pi
Pj
Figure 4.9
3. The last thing we need to prove is that the iteration LegalizedEdge(,) is not
infinite loop. It is clear because Delaunay Triangulation give O (n) number triangles
from Theorem 2.1. So the recursion function searches at most O (n) triangles.
Concluding the above proofs, we can regard this algorithms is correct.
Implementing the algorithm by Java, the difficulty is how to look for Pi according to
the current edge PiPj. Actually PiPi is the next edge of PiPj's twin and PiPj is the
previous edge of PiPjs twin.
36
So in the iteration of LegalizedEdge(,), we must update each edge data structure to
make it store new information of Delaunay Triangulation. It makes this algorithm
complicated.
37
5. Analysis
What is the running time of this randomized version of this incremental algorithm?
When we add a point into Delaunay Triangulation created on early stage, the
procedure includes several operations: querying, splitting, and flipping. The running
time of each step is depended on the order of point permutation. In some cases, when
a point is input, it only splits a triangle that contains it. We call the cases the best
case, because it only deals with one triangle and adds one for the depth of data
structure D. In the worstcase, it deals with all current triangles of Tr ( Delaunay
Triangulation created on the r stage).
The other cases are between these two extremes. So we must use the expected
running time to bound this algorithm running time that is a average running time
over all possible permutations.
Before we analyses the running time bound Delaunayay Triangulation, we introduce
an important theorem.
Theorem 5.1 The expected number of triangles created in Delaunay Triangulation is
at most 9n+l [1].
38
In order to prove this Theorem and next conclusions, we introduce a method
backwards analysis. Assume that the r1 step has been finished, we want to get some
results on the r step from the r1 step.
When we analyze the running time of Delaunay Triangulation, we can find it
composes of the running time of several operations: querying, splitting and flipping.
Because the expected number of the triangles of Delaunay Triangulation is at most
9n+l, for every insertion of point we at most flip illegal edge 9n+l times. The
expected running time of operation flipping is bounded by 0(n).
Now we mainly pay attention to the expected running time of location query. We
discuss this problem from data structure D. When a point is inserted into Delaunay
Triangulation, we will search from the root of D to its leaves to decide which
triangle contains this point. So the running time is depended on the depth of D, and
we can use backwards analysis to obtain the expected running time.
Assume that the expected depth of D ri is H ri, we try to get the expected depth of
Dr. See Figure 5.1
39
Figure 5.1
From Dri to D r, the expected number of triangles is created is 9(r)+l from
Theorem 5.1, and we simply present it as r. We have known that one triangle of
Delaunay Triangulation can produce two or three triangles by splitting or flipping
operation. The triangle in which the inserted point locates will only be split, or its
edges will be flipped continually. So in the best case it will add one to the depth of D
ri, and in the worst case it will flip all triangles of Delaunay Triangulation to add r.
We use the average value of depth that is increased. From Theorem 5.1, in order to
produce r (the expected number of triangles) triangles from one triangle, it must do
flipping the expected number of In r times. So it add the depth of Dri In r. Now the
40
depth of Dr becomes H ri + In r. We can get the following formula of the expected
depth of D n.
Hn = In 1+ In 2+,,, + In n1 +ln n.
Hn
Hn > In n/2 + In n/2 + ,,,+ In n/2 + In n/2 = n/2 In n/2.
So we can bound the value of H n by 0(n Inn), and get the expected running time of
location query is 0( n In n).
From the discussion above, we can conclude that the expected running time of this
incremental algorithm is 0( n In n).
41
6. Applications
6.1 Representation of Maps
In graphics, we need to represent a surface of the earth in three dimensions image.
Usually we only measure the heights of some sample points for a certain region, and
we must use these sample points to make our graph more natural. These sample
points can be linked with edges so that the surface is separated many discrete planar
elements that may be triangles, rectangles or other shapes. Usually we use triangles
as our element shapes.
The important problem is how to triangulate these sample points, because for some
points there are many ways to triangulate them. But some triangulations will cause
long and sharp triangles and make the representation of surface not natural.
Delaunay Triangulation is a effective method to solve this problem. We can project
the sample points to a plane, and then triangulate corresponding points on the plane
by Delaunay Triangulation. In this way we can get the edges that do not intersect
each other, and make every point have at least two edges to link it. We still use these
edges to link the sample points on the surface so that the surface will be separate
many triangles. So we can use these triangles to represent the surface.
42
6.2 Finite Element Analysis
Analyzing the structural property of complex shapes is often accomplished by a
technique called finite element analysis. This is used, for example, by aerodynamics
to model airplane wings. The domain to be studied is partitioned into a grid of finite
element, and then the relevant differential equations modeling dynamics are solved
by discrediting over the partition. The stability of the numerical procedures used
depends on the quality of the partition. Delaunay Triangulation is a good partition
for this problem.
6.3 The Nearest Neighbor Problem
We have know that Voronoi Diagram can be used to decide the nearest neighbor
problem. We prove that for a set of n points Voronoi Diagram can be computed from
Delaunay Triangulation by 0(n) running time.
Proof: For very edge of Delaunay Triangulation, we draw a perpetual bisector, and
then for each point we can a get a region in which all of points are closest this point.
Each region is bounded in several bisectors. Because the number of edges in
Delaunay Triangulation is 0(n), the running time of doing perpetual bisectors of the
edges is 0(n).
43
6.4 Minimum Spanning Tree
A minimum Spanning Tree of a set of points is a minimum length tree that spans all
the Points. When the length of tree is measured by Euclidean length, the tree is often
called the Euclidean minimum spanning tree (EMST). EMST has many applications,
for example, the computer network uses it to minimizes total wire length.. Kruskals
algorithm (1956) is a simple way to solve this problem: a greedy strategy finds
EMST based on that a shortest tree must be built up by incrementally by adding the
shortest edges, not explored, also maintaining ( acyclicity ). Because for a set of n
points, it can generates nxn edges between very two points, the sorting these edges
according to their length needs 0(nxn Inn) running time.
We can prove that EMST is a subset of Delaunay Triangulations edges. So we only
find MSTs edges from Delaunay Triangulations edge. The sorting of these edges
lengths need 0(n Inn) running time, so the total running time for this algorithm is
Q(nlnn ).
44
Appendix A. Delaunay Triangulation
import java.io.*;
import java.util.*;
import java.awt.*;
import java.applet. *;
import java.util.Random;
// This class includes the coordinates of a point in the plane.
class vertex {
public int x;
public int y;
public vertexO{;}
public vertex(int x, int y) {
this.x=x;
this.y=y;
}
}
// This class includes the information of an edge,
class edge{
public vertex origin; //The start point of an edge,
public vertex destin; //The end point of an edge,
public edge next; //The next edge of the current edge,
public edge prev; //The previous edge of the current edge,
public edge twin; //The twin edge of the current edge,
public edge(){;}
}
// The triangle created in Delaunay Triangulation,
class triangle {
45
public int flag; //The status of a triangle visited,
public int flagl;
public triangleO{flag=0;flagl=0;}// Initialize the value of flag.
public edge triEdgel;//The edges of a triangle,
public edge triEdge2;
public edge triEdge3;
public triangle child 1;// The triangles created by this triangle,
public triangle child2;
public triangle child3;
}
// This class includes the methods to build the data structure,
// and all of operations for Delaunay triangulation.
class tree {
static PrintWriter screen=new PrintWriter(System.out,true);
static BufferedReader keyboard
=new BufferedReader(new InputStreamReader(System.in));
static FileWriter file2;
static PrintWriter outputFile;
static int triangleNumber=0;//Count the number of the triangles
//created in Delaunay Triangulation.
static int triangleNumberl=0;
static int vertexNumber; //The number of the points.
static vertex[] node; //The array of the points.
static triangle root; //The root of the tree.
static edge e;
//The method provides an output file that stores the edges of Delaunay
// Triangulation.
46
static void initl() throws IOException{
file2=new FileWriter("resultsl .txt");
outputFile=new PrintWriter(file2);
}
// The method reads the coordinates of the points and puts them
//in the array of vertices.
static void init() throws IOException{
FileReader file=new FileReader("t.txt");//Open a input file named "t.tet".
BufferedReader inputFile=new BufferedReader(file);
//Read the number of the points.
vertexNumber=newInteger(inputFile.readLineO).intValueO;
node=new vertex[vertexNumber];
inti;
for(i=0;i
node[i]=new vertexO;
node[i].x=newInteger(inputFile.readLine()).intValueO;
node[i].y=newInteger(inputFile.readLineO).intValueO;
}
}
// Create a initial triangle that includes all of points, and let its top points far from
// the points triangulated.
static void getStart(){
int max_x,max_y,m;
47
int i,j;
max_x=Math.abs(node[0] .x);
max_y=Math.abs(node[0] .y);
// Look for the maximal absolute value of the coordinates of the points.
for(i=l ;i
if(max_x
if(max_y
}
m=max_x;
if(max_y>max_x)m=max_y;
vertex[] bound;
bound=new vertex[3];
bound[0]=new vertex(0,3*m); // Set the coordinates of the three
// triangles top points.
bound[l]=new vertex(3*m,0);
bound[2]=new vertex(3*m,3*m);
edge el=new edgeO; //Create the three edges of the initial triangle,
edge e2=new edge();
edge e3=new edgeO;
e 1 .origin=bound[0];
el ,destin=bound[l];
e2.origin=bound[l ];
e2.destin=bound[2];
e3 .origin=bound[2];
e3 .destin=bound[0];
//Give each edge the next, previous and twin edges.
el .next=e2;
el.prev=e3;
48
el.twin=null; //The initial triangle's twin is set to null.
e2.next=e3;
e2.prev=el;
e2.twin=null;
e3.next=el;
e3.prev=e2;
e3.twin=null;
triangle temp=new triangle();
temp.triEdge 1 =e 1;
temp.triEdge2=e2;
temp.triEdge3 =e3;
temp.childl=null; //Its children are null before triangulation.
temp.child2=null;
temp.child3=null;
root=temp;
}
//The method prints all leaves of the tree, and these leaves represent the triangles of
//Delaunay Triangulation.
static void printLeaf(triangle myTriangle) throws IOException{
triangle temp=new triangleO;
myTriangle.flag=l;//Mean that this triangle has been visited.
temp=myTriangle;
if(temp.child 1 =null) {
triangleNumber=triangleNumber+l;
screen.println(temp.triEdgel .origin.x+","+temp.triEdge 1 .origin.y+">"+
temp.triEdgel .destin.x+","+temp.triEdgel .destin.y+">"+
49
temp.triEdgel .next.destin.x+",''+temp.triEdgel .next.destin.y);
outputFile.println(temp.triEdgel .origin.x);// Put the coordinates of the
outputFile.println(temp.triEdgel.origin.y);//triangles into the output file,
outputFile.println(temp.triEdgel.destin.x);//and arrange them in the order
outputFile.println(temp.triEdgel.destin.y);//by one by one.
outputFile.println(temp.triEdge 1 .next.destin.x);
outputFile.println(temp.triEdge 1 .next.destin.y);
return;
}
if(temp.childl !=null&&temp.childl.flag=0) printLeaf(temp.childl);
if(temp.child2 !=null&&temp.child2.flag==0) printLeaf(temp.child2);
if(temp.child3 !=null&&temp.child3.flag=0) printLeaf(temp.child3);
}
// The method counts the number of the triangles created, and it searches
//from the graph.
static void getNumberOfTriangle(triangle myTriangle){
triangle temp=new triangleO;
myTriangle.flagl=l;
temp=myTriangle;
if(temp.child 1 =null) {
triangleNumber 1 =triangleNumber 1+1;
return;
}
if(temp.childl!=null&&temp.childl.flagl=0)getNumberOfTriangle(temp.childl);
50
ifl[temp.child2!=null&&temp.child2.flagl=0)getNumberOfrriangle(temp.child2);
if(temp.child3!=null&&temp.child3.flagl=0)getNumberOfTriangle(temp.child3);
}
// The method includes the operations of splitting and flipping.
static void splitAndFlip(vertex myNode,triangle myTriangle){
inti,j;
intdl,d2,d3;
triangle temp =new triangle();
triangle templ=new triangleQ;
triangle temp2=new triangle();
triangle temp3=new triangle();
triangle temp4=new triangle();
// Give the value of verification to decide if a point lies on an edge,
dl =criter(myTriangle.triEdge 1 .origin.x,
myTriangle.triEdge 1 .origin.y,
myTriangle.triEdgel .destin.x,
myTriangle.triEdgel .destin.y,
myNode.x,myNode.y);
d2=criter(myTriangle.triEdge2.origin.x,
myTriangle.triEdge2.origin.y,
myTriangle.triEdge2.destin.x,
myTriangle.triEdge2.destin.y,
myNode.x,myNode.y);
d3=criter(myTriangle.triEdge3.origin.x,
myTriangle.triEdge3 .origin.y,
myTriangle.triEdge3 .destin.x,
myTriangle.triEdge3 .destin.y,
51
myNode.x,myNode.y);
//If a point does not lie on any edge of a triangle, it must lie inside the triangle.
//Then we split this triangle into the three triangles.
if(dl !=0&&d2!=0&&d3!=0){
//Give all information of these three triangles, and create
//the three new triangles.
//Put these three triangles as the children of the triangle split.
edge el=new edge();
edge twinl=new edgeO;
edge e2=new edgeO;
edge twin2=new edge();
edge e3=new edgeO;
edge twin3=new edgeO;
e 1 .origin=myNode;
el .destin=myTriangle.triEdgel .origin;
twinl.origin=el.destin;
twinl .destin=el .origin;
el.twin=twinl;
twinl .twin=el;
e2.origin=myNode;
e2.destin=myTriangle.triEdgel .destin;
twin2.origin=e2.destin;
twin2.destin=e2.origin;
e2.twin=twin2;
twin2.twin=e2;
e3 .origin=myNode;
e3.destin=myTriangle.triEdge2.destin;
52
twin3 .origin=e3 .destin;
twin3 .destin=e3 .origin;
e3.twin=twin3;
twin3.twin=e3;
e 1 .next=myTriangle.triEdge 1;
el.prev=twin2;
myTriangle.triEdge 1 .next=twin2;
myTriangle.triEdge 1 .prev=el;
twin2.next=el;
twin2.prev=myTriangle.triEdgel;
temp 1 .triEdge 1 =e 1;
tempi .triEdge2=myTriangle.triEdgel;
temp 1 .triEdge3=twin2;
updateEdge(temp 1);
e2.next=myTriangle.triEdge2;
e2.prev=twin3;
myTriangle.triEdge2.next=twin3;
myTriangle.triEdge2.prev=e2;
twin3.next=e2;
twin3 .prev=myTriangle.triEdge2;
temp2.triEdgel=e2;
temp2.triEdge2=myTriangle.triEdge2;
temp2.triEdge3=twin3;
updateEdge(temp2);
e3 .next=myTriangle.triEdge3;
e3.prev=twinl;
myTriangle.triEdge3 .next=twinl;
myTriangle.triEdge3 ,prev=e3;
twinl.next=e3;
twinl .prev=myTriangle.triEdge3;
53
temp3 .triEdge 1 =e3;
temp3 .triEdge2=myT dangle.triEdge3;
temp3.triEdge3=twinl;
updateEdge(temp3);
tempi .child 1 =temp 1 .child2=templ .child3=null;
temp2.childl=temp2.child2=temp2.child3=null;
temp3 .child 1 =temp3 .child2=temp3 ,child3=null;
myTriangle.childl=templ;
myTriangle.child2=temp2;
myTriangle.child3=temp3;
//Flip each edge of the triangle until we get a legal edge.
legalizeEdge(myNode,my Triangle. triEdge 1);
legalizeEdge(myNode,myTriangle.triEdge2);
legalizeEdge(myNode,myTriangle.triEdge3);
}
// If the point lies on the one edge of triangle, it will split the triangle into
// the four triangles
else{
edge tempEdge=new edge();
edge el =new edgeO;
edge twinl=new edgeO;
edge e2 =new edgeO;
edge twin2=new edgeO;
edge e3 =new edgeO;
edge twin3=new edgeO;
edge e4 =new edgeO;
edge twin4=new edgeO;
if(di=0)tempEdge=myTriangle.triEdgel;// Look for the edge
if(d2=0)tempEdge=myTriangle.triEdge2;//containing
if(d3=0)tempEdge=myTriangle.triEdge3 ;//the point.
54
temp=fmdTriangle(root,tempEdge.twin);//Look for a triangle
// corresponding to
//the edge's twin.
// Create four triangles, and give them all information.
//Then put them as the children of the current triangle.
el .origin=tempEdge.origin;
el .destin=myNode;
twinl .origin=el .destin;
twinl .destin=el .origin;
el.twin=twinl;
twinl .twin=el;
e2.origin=myNode;
e2.destin=tempEdge.destin;
twin2.origin=e2.destin;
twin2.destin=e2.origin;
e2.twin=twin2;
twin2.twin=e2;
e3 .origin=myNode;
e3 .destin=tempEdge.next. destin;
twin3 .origin=e3 .destin;
twin3 .destin=e3 .origin;
e3.twin=twin3;
twin3.twin=e3;
e4.origin=myNode;
e4.destin=tempEdge.twin.next.destin;
twin4.origin=e4.destin;
twin4.destin=e4.origin;
e4.twin=twin4;
twin4.twin=e4;
el.next=e3;
e 1 .prev=tempEdge.prev;
55
e3 .next=tempEdge.prev;
e3.prev=el;
tempEdge.prev.next=el;
tempEdge.prev.prev=e3;
temp 1 .triEdge 1 =e 1;
temp 1 .triEdge2=e3;
temp 1 .triEdge3=tempEdge.prev;
updateEdge(temp 1);
e2.next=tempEdge.next;
e2.prev=twin3;
tempEdge.next.next=twin3;
tempEdge.next.prev=e2;
twin3.next=e2;
twin3 .prev=tempEdge.next;
temp2.triEdge 1 =e2;
temp2.triEdge2=tempEdge.next;
temp2.triEdge3=twin3;
updateEdge(temp2);
twinl .next=tempEdge.twin.next;
twinl .prev=twin4;
tempEdge.twin.next.next==twin4;
tempEdge.twin.next.prev=twinl;
twin4.next=twin 1;
twin4.prev=tempEdge.twin.next;
temp3 .triEdge 1 =twin 1;
temp3 .triEdge2=tempEdge.twin.next;
temp3.triEdge3=twin4;
updateEdge(temp3);
56
twin2.next=e4;
twin2.prev=tempEdge.twin.prev;
e4.next=tempEdge.twin.prev;
e4.prev=twin2;
tempEdge.twin.prev.next=twin2;
tempEdge.twin.prev.prev=e4;
temp4.triEdge 1 =twin2;
temp4.triEdge2=e4;
temp4.triEdge3=tempEdge.twin.prev;
updateEdge(temp4);
temp 1 .child 1 =temp 1 .child2=temp 1 .child3=null;
temp2.childl=temp2.child2=temp2.child3=null;
myTriangle.child 1 =temp 1;
myTriangle.child2=temp2;
myTriangle.child3=null;
temp 1 .child 1 =temp 1 .child2=temp 1 .child3=null;
temp2.childl=temp2.child2=temp2.child3=null;
temp3 .child 1 =temp3 .child2=temp3 .child3=null;
temp4.childl=temp4.child2=temp4.child3=null;
temp.childl=temp3;
temp.child2=temp4;
temp.child3=null;
//Flip the edges of quadrilateral that bound these four triangles.
legalizeEdge(myNode,tempEdge.prev);
legalizeEdge(myNode,tempEdge.next);
legalizeEdge(myNode,tempEdge.twin.next);
legalizeEdge(myNode,tempEdge.twin.prev);
}
57
}
// Give a point and find a triangle from the current Delaunay Triangulation,
//that contains this point.
static triangle query(vertex myNode){
hit ij;
intdl,d2,d3;
triangle temp =new triangleO;
triangle temp4=new triangleO;
// Search the root of graph until the leaves, and find this triangle
//from the leaves.
temp=root;
while(temp.childl !=null){
temp4=temp;
for(i=l;i<=3;i++){
if(i==l&&temp4.childl !=null)temp=temp4.childl;
if(i=2&&temp4.child2!=null)temp=temp4.child2;
if(i=3&&temp4.child3 !=null)temp=temp4.child3;
// If this point lies on the left of very edge of a triangle, this triangle must
// contains this point.
d 1 =criter(temp.triEdge 1 .origin.x,temp.triEdge 1 .origin.y,
temp.triEdgel .destin.Xjtemp.triEdgel .destin.y,
myNode.x,myNode.y);
d2=criter(temp.triEdge2.origin.x,temp.triEdge2.origin.y,
temp.triEdge2.destin.x,temp.triEdge2.destin.y,
myNode.x,myNode.y);
d3=criter(temp.triEdge3 .origin.x,temp.triEdge3 .origin.y,
58
temp.triEdge3 .destin.x,temp.triEdge3 .destin.y,
myNode.x,myNode.y);
if(d 1 <=0&&d2<=0&&d3<=0)break;
}
}
return temp;
}
// The recursion function flips a edge until getting a legal edge,
static void legalizeEdge(vertex p,edge e){
int cir;
if(e.twin==null)retum;// It reaches the initial triangle's edges.
else{
// Get the verification value to decide if a edge is illegal.
cir=verifyCircle(e.origin,e.destin,p,e.twin.next.destin);
screen.println(cir);
if(cir>0)retum; // The edge is legal, the recursion is finished.
// Otherwise, the recursion will continue.
else{
triangle triangle l=new triangle();
triangle triangle2=new triangleO;
triangle templ=new triangle();
triangle temp2=new triangle();
edge newEdge=new edgeO;
edge newTwin=new edge();
59
edge tempEdgel=new edgeO;
edge tempEdge2=new edge();
templ=findTriangle(root,e);// Look for a triangle incident to the edge e.
temp2=findTriangle(root,e.twin);//Look for a triangle incident to the edge
// e's twin.
//Update the triangles created due to flipping.
newEdge.origin=p;
newEdge.destin=e.twin.next.destin;
newTwin.origin=e.twin.next.destin;
newTwin.destin=p;
newEdge.next=e.twin.prev;
newEdge.prev=e.next;
e.twin.prev.next=e.next;
e.twin.prev.prev=newEdge;
e.next.next==newEdge;
e.next.prev=e.twin.prev;
newTwin.next=e.prev;
newTwin.prev=e.twin.next;
e.prev.next=e.twin.next;
e.prev.prev=newTwin;
e.twin.next.next=newTwin;
e.twin.next.prev=e.prev;
newEdge.twin=newTwin;
newT win.twin=newEdge;
triangle 1 .triEdge 1 =e.prev;
triangle 1 .triEdge2=e.twin.next;
triangle 1 .triEdge3=newTwin;
tempEdgel=trianglel.triEdge2;
triangle2.triEdgel=newEdge;
60
triangle2.triEdge2=e.twin.prev;
triangle2.triEdge3=e.next;
tempEdge2=triangle2.triEdge2;
/*triangle 1 .triEdge3 .twin=triangle2.triEdge 1;
triangle2.triEdge 1 .twin=triangle 1 .triEdge3; */
// Put these new triangles as the children of the triangles involved in flipping.
trianglel .childl=trianglel .child2=trianglel .child3=null;
triangle2.childl=triangle2.child2=triangle2.child3=null;
tempi .childl=trianglel;
temp 1 .child2=triangle2;
temp 1 .child3=null;
temp2.childl=trianglel;
temp2.child2=triangle2;
temp2.child3=null;
legalizeEdge(p,tempEdge 1);
legalizeEdge(p,tempEdge2);
}
}
}
// The method is used to decide if a point lies on the left of an edge.
static int criter(int xl ,int yl ,int x2,int y2,int x3,int y3){
intd;
d=x2*y3+xl *y2+x3 *y 1 x2*y 1 x3 *y2xl *y3;
return d;
61
// The method is used to update the relationship of the edges of a triangle.
static void updateEdge(triangle myTriangle){
myTriangle.triEdgel.next=myTriangle.triEdge2;
myTriangle.triEdge 1 .prev=myTriangle.triEdge3;
myTriangle.triEdge2.next=myTriangle.triEdge3;
myTriangle.triEdge2.prev=myTriangle.triEdgel;
myTriangle.triEdge3.next=myTriangle.triEdgel;
myTriangle.triEdge3.prev=myTriangle.triEdge2;
}
// The method is used to calculate the verification value, and decide
//if an edge is illegal.
static int verifyCircle(vertex x,vertex y,vertex z, vertex t){
int xl ,x2,y 1 ,y2,zl ,z2,tl ,t2,cir;
xl=x.x; x2=x.y;
yl=y.x; y2=y.y;
zl=z.x; z2=z.y;
tl=t.x; t2=t.y;
cir=( xl *(y2*(zl *zl+z2*z2)+t2*(y 1 *y l+y2*y2)+z2*(tl *tl+t2*t2)
t2 (z 1 z 1 +z2 z2)y2 (t 1 *t 1 +t2*t2)z2 (y 1 *y 1 +y2 *y2)) )
( x2*(y 1 *(zl *zl+z2*z2)+tl *(y 1 *y l+y2*y2)+zl *(tl *tl+t2*t2)
tl *(zl *zl+z2*z2)y 1 *(tl *tl+t2*t2)zl *(y 1 *y l+y2*y2)))
+
((xl*xl+x2*x2)*(yl*z2+y2*tl+t2*zltl*z2t2*yly2*zl) )
( y 1 *z2*(tl *tl+t2*t2)+y2*tl *(zl *zl+z2*z2)+zl *t2*(y 1 *y l+y2*y2)
(yl*yl+y2*y2)*z2*tl(zl*zl+z2*z2)*t2*yl(tl*tl+t2*t2)*y2*zl);
62
}
return cir;
// The method is used to look for a triangle with the specified edge from the current
// delaunay triangulation.
static triangle findTriangle(triangle myTriangle,edge my Edge) {
triangle temp=new triangle();
triangle templ=new triangleO;
temp=myTriangle;
if(temp==null)retum temp;
if((temp.triEdge 1 =myEdge temp.triEdge2=myEdge temp.triEdge3=myEdge)
&&(temp.childl=null))retum temp;
else{
temp 1 =findTriangle(temp.childl ,myEdge);
if(templ==null){
temp 1 =findTriangle(temp.child2,myEdge);
if(temp 1 =null) {
templ=findTriangle(temp.child3,myEdge);
}
}
}
return tempi;
// The method is used to sort a array according to the increasing order.
// The algorithm is Quick Sorting.
static void QuickSort(float[] a,int[] b,int m,int n){
if(m
63
int p=partition(a,b,m,n);
QuickSort(a,b,m,p1);
QuickSort(a,b,p+l ,n);
}
}
static int partition(fIoat[] a,int[] b,int i,int j){
float pivot,temp;
int k,middle,p,ch;
middle=(i+j)/2;
pivot=a[middle]; a[middle]=a[i] ;a[i]=pivot;
p=i;
for(k=i+l ;k<=j ;k++) {
if(a[k]
temp=a[++p] ;a[p]=a[k] ;a[k]=temp;
ch=b[p];b[p]=b[k];b[k]=ch;
}
}
temp=a[i] ;a[i]=a[p] ;a[p]=temp;
retump;
static public void main(String[] args) throws IOException{
int k,indexAfterRandom;
triangle temp=new triangleO;
int i,pointNumber;
int maxNumberOfPoint=l 00;
initlQ;// Open an output file to store the vertices of Delaunay Triangulation.
64
init(); // Read the coordinates of all points that will be triangulated.
getStartO;// Get the initial triangle.
screen.println("Please input the number of points(ln):");
pointNumber=new Integer(keyboard.readLine()) .intV alueO;
float[] randomNumber=new float[maxNumberOf Point];
int[] randomIndex=new int[maxNumberOf Point];
Random value=new RandomO;
// Create n (the number of points) random number between 0 and 1.
for(i=0;i
randomNumber[i]=value.nextFloatO;
randomIndex[i]=i;
}
QuickSort(randomNumber,randomIndex,0,pointNumber1);
for(i=0;i
screen.println(randomIndex[i]);
}
for(i=0;i
mdexAfterRandom=randomIndex[i];// The index is pop randomly.
temp=query(node[indexAfterRandom]);// Query the location of this point.
splitAndFlip(node[indexAfterRandom],temp);// Split and Flip the triangle
// that contains this point.
65
}
getNumberOfTriangle(root);
outputFile.priritln(triangleNumberl);
printLeaf(root);
outputFile.closeO;
/* temp=findTriangle(root,temp.triEdgel .twin);
screen.println(temp.triEdgel.origin.x+" "+temp.triEdgel.origin.y+""
+temp.triEdgel.destin.x+" "+temp.triEdgel.destin.y+" "
+temp.triEdge2.destin.x+" "+temp.triEdge2.destin.y);
i=verifyCircle(temp.triEdge 1 .origin,temp.triEdgel .destin,
temp.triEdgel .next.destin,temp.triEdgel .twin.next.destin);
screen.println(i);*/
for(i=0;i
screen.printIn(randomIndex[i]);
}
}
66
t
Appendix B. Graph
import java.io.*;
import java.awt.*;
import java.awt.event. *;
class Gui extends Frame implements WindowListener,MouseListener{
public int myNumber;
public int myMax;
static PrintWriter screen= new PrintWriter(System.out,true);
static int[J myX;
static int[] myY;
public Gui(String s,intO x,int[] y, int number, int max){
super(s);
myNumber=number;
myMax=max;
myX=new int[number];
myY=new int [number];
for(int i=0;i
myX[i]=x[i];
myY[i]=y[i];
}
screen.println(myX[0]);
setBackground(Color. white);
addMouseListener(this);
addWindowListener(this);
}
public void mouseClicked(MouseEvent event) {}
public void mouseEntered(MouseEvent event) {}
public void mouseExited(MouseEvent event)]}
public void mouseReleased(MouseEvent event)]}
67
public void mousePressed(MouseEvent event) {
inti;
Graphics g=getGraphics();
for(i=0;i
g.drawLine((myX[i]+myMax)* 10+100,(my Y[i]+myMax)* 10+100,
(myX[i+l ]+myMax)* 10+100,(myY[i+l ]+myMax)* 10+100);
g.drawLine((myX[i+l]+myMax)* 10+100,(my Y[i+1 ]+myMax)* 10+100,
(myX[i+2]+myMax)* 10+100,(myY[i+2]+myMax)* 10+100);
g.drawLine((myX[i+2]+myMax)* 10+100,(myY[i+2]+myMax)* 10+100,
(myX[i]+myMax)* 10+100,(my Y[i]+myMax)* 10+100);
}
screen.println(myX[0]);
}
public void windowClosed(WindowEvent event) {}
public void windowDeiconified(WindowEvent event) {}
public void windowIconified(WindowEvent event) {}
public void windowActivated(WindowEvent event) {}
public void windowDeactivated(WindowEvent event) {}
public void windowOpened(WindowEvent event) {}
public void windowClosing(WindowEvent event) {
System.exit(0);
}
68
}
class Exl {
static PrintWriter screen=new PrintWriter(System.out,true);
public static void main(String[J args) throws IOException {
int numberOfTriangle;
intQ x;
intQy;
int maxX,maxY,max;
FileReader file=new FileReader("resultsl.txt");
BufferedReader inputFile=new BufferedReader(file);
numberOfTriangle=newInteger(inputFile.readLineO)intValue();
numberOfT riangle=numberOfTriangle* 3;
x=new intfnumberOfTriangle];
y=new int[numberOfTriangle];
inti;
for(i=0;i
x[i]=newInteger(inputFile.readLine()).intValue();
y[i]=newInteger(inputFile.readLineO).intValue();
}
maxX=Math.abs(x[0]); maxY=Math.abs(y[0]);
for(i=0; i
69
if(Math.abs(x[i])>maxX)maxX=Math.abs(x[i]);
if(Math.abs(y[i])>maxY)maxY=Math.abs(y[i]);
}
max=maxX;
if(maxY>maxX)max=maxY;
max=3*max;
screen.println(x[0]);
Gui screen=new Gui("Delaunay
Triangulation" ,x,y,numberOfTriangle,max);
screen.setSize(400,400);
screen.setVisible(true);
}
}
70
References
[1] M. de. berg, M. Van. Kreveld, M. Overmars, O.Schwarzkopf. Computational
Geometry Algorithms and Applications, Addision Wesley Longman, Inc,1998.
[2] Joseph. ORourke. Computational Geometry In C, second edition, Cambridge
University Press, 1998.
[3] Roger T. Stevens. Graphics Programming in C, M&T Publishing, Inc, Redwood
City, California, 1998.
[4] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest. Introduction to
Algorithms, the MIT Press, Cambridge, Massachusetts, London, England.
[5] Thomas A. Standish. Data Structure in Java, Addision Wesley Longman, Inc,
1998.
[6] Barry Holmes. Programming with Java, Jones and Bartlett Publisher, Inc, 1998.
[7] Vartan Piroumian. Java GUI Development, A Division of Macmillan Computer,
1999.
71
