Here 4 cases are shown with number of sites of 100, 2000, 5964, 525383,  respectively. The third row of the table shows the CPU time used, and also shown in the figure.  The X axis of the figure is the number of sites, the Y axis is the CPU time in seconds.  The computational time basically follows the O(n*log(n)) with 7e-10 as the front coefficient.

Download and running

mkdir somewhere/quad_edge and download everything in directory quad_edge.
1. source files: Delaunay_java.tar.gz
2. class files: Delaunay_class.tar.gz
3. library: jdsl.tar.gz
If   you do not bother to complile, download 2 and 3 and set CLASSPATH as follow on Unix.
set quad = "somewhere/quad_edge"
setenv CLASSPATH $quad/../:.:$quad/lib:$quad/cs252/asgn4:$quad/cs252/:$quad/:lib/jdsl.jar
you can modify setPath which I set my path by  "source setPath".

To compile and run the code

javac DelaunayDC.java
java DelaunayDC -d -r100
java DelaunayDC -d rt1.dat
java DelaunayDC -b -r
 
-d is Divide-and-Conquer method. -b is Brutal  Force method. -r100 means generat 100 random
sites in the code without input file. rt1.dat means your provide the site file rt1.dat

To run with Embedded Planar Graph (EPG) output

cd EPG
javac Voroinoi.java
java Asgn4

than click on the pop up window. Since the supporting code Asgn4 of drawing is constructed with integer, the double random number may result in a duplicate integer, and throw out Duplicate points Exception. But they are really not duplicate and I only have the complile Asgn4.class, thus nothing can be done.

There is also this stand alone adda.java self contained everything in one file that
 portrait1 can be reproduced in less than 4 minutes following these steps.

step 1. download adda.java and rt1.dat and gunzip rt1.dat.gz on unix.
step 2. javac adda.java
step 3. java adda -d rt1.dat   (if path not set on unix, type "setenv CLASSPATH ." ) It generate the triangle data.
step 4. drawrt.m read in the triangle data in step 3 and draw it. download drawrt.m  in the same director, open a Matlab window, type drawrt, portrait1 will appear on the window.

javadoc

Also here is CS252Project.tar.gz that tar everything with subdirectory
structure including this page.

Algorithm

Embedded planar graph
There maybe many ways to construct an embedded planar graph (EPG). Since the graph must be connected all the time during the construction and when inserting
an edge the  next(or previous) edge of both vertex ring must be part of the
EPG already, the edges could not be added simply one by one. I figure out two ways.

• Brutal force method:

step 1:

Save all the QuadEdge in an array with a method which can  remove any element of the array.

step 2:

Add edge from the left most convex hull as the first edge, since this first edge is supposed to connect to an infinite point in the EPG and the DivideAndConquer algorithm give the handle for free (automatically). Remove this QuadEdge from the array of step 1 and mark in QuadEdge that it has been added. Sine QuadEdge keeps all the reference of its neighborhood, we proceed to add the next edge in the ring until we could not proceed (next edge is marked added).

step 3:

while (there are still elements in the array of step 1)
{
  Try each edge to see if it could be inserted to attached.
  If it could , then add the edge to the graph(EPG) and remove
  it from the aforementioned array in step 1.
  if i++> n^2, break;
}

step 3 could be a dead loop if something is wrong, for example the the graph is tangled and not planar. Thus I put the *break* there.

• Recursive Traverse:

I create a method (addToEPG) in QuadEdge class to add  *this.QuadEdge*  (itself) to the embedded planar graph. Imagine a tree of vessel, from the root, the fluid filling up each branch and leaves. Here is the pseudo code.

addToEPG()
{
  if(this edge is not added)
  {
    add this edge with the embedding information;
    mark it as added
    Onext.addToEPG();
    lnext.addToEPG();
  }
  return
}

The Onext and lnext exhaust all the possibility of the other edges that have a path to this edge, thus it can traverse all the edges recursively, since we are supposed to deal with a connected graph. I do not need to check if Onext or lnext is null, since they are never null even on a boundary of any sense.

I also put a proxy of Vertex in the Geometric data to decide  to attachVertex  or to insertEdge. This is no problem of attaching an edge, but for inserting an edge,  when the other end does not have an Order.AFTER or Order.BEFORE information, it needs special care. The method is O(n*log(n)) with a much smaller coefficient in front, thus is a match algorithm with Divide and Conquer. As you see the pseudo code is very simple and real code is just half a page. To call it

leftMostConvexHull.addToEPG();

leftMostConvexHull is the QuadEdge as its name implys and is a free handle as explain in step 2 of my first method.

The quad_edge data structure is a natrual computer implementation of hte corresponding edge algebra.
The main advantage of it is that the construction and modification of arbitrary diagrams can be effected
by as few as two basic topological operators, makeEdge and splice. The quad-edge data structure contains no separate records for vertice or faces; a vertex is implicitly defined as a ring of edges, and the standard way to refer to it is to specify one of its out going edges. All the operation on the edge  can be
reduced to two operation, Onext and rot. Thus there are only two private variable reference in
each QuadEdge class, Onext and rot. All the other neighbors can be find by these two gateway *Onext*,
*rot*.  The topological relation is as shown in the following code:

  public QuadEdge onext() { return Onext; }
  public QuadEdge rot() { return rot; }
  public QuadEdge sym() { return rot.rot(); }
  public QuadEdge rotSym() {return rot.sym(); }
  public QuadEdge oprev() { return rot.onext().rot(); }
  public QuadEdge lnext() { return rotSym().onext().rot(); }
  public QuadEdge lprev() { return onext().sym(); }
  public QuadEdge rnext() { return rot.onext().rotSym(); }
  public QuadEdge rprev() { return sym().onext(); }

QuadEdge also keeps a geometric data field:

  public GeomSites orig() { return data; }
  public GeomSites dest() { return sym().orig(); }
  public GeomSites right(){ return rot.orig(); }
  public GeomSites left() { return rot.sym().orig(); }

makeEdge() and splice() is very simple. makeEdge() create four instance of QuadEdge
and set Onext and rot as

    q0.rot = q1;
    q1.rot = q2;
    q2.rot = q3;
    q3.rot = q0;

    q0.Onext = q0;
    q1.Onext = q3; //on the sphere, left=right
    q2.Onext = q2;
    q3.Onext = q1;

return q0 as handle.

By exchanging two pair of next reference, splice() can perform the  splitting and joining
of graphs. The idea are really simple and smart.

splice()
{
    a.Onext = bOnext
    b.Onext = a.Onext   (the value before a.Onext is overide by bOnext)

    alpha.Onext = beta.Onext;
     beta.Onext = alpha.Onext;
}
Refer figure 10 of [1] for meaning of alpha and beta.

Connect(a,b) is not an independent function, it is a combination of makeEdge() and splice().
It add a new edge q connecting the detination of a to the origin of b, in such a way
that a.left()=e.left()=b.left().

First sort sites in increasing x-coordinate. When there are ties we resolve them by sorting in incrasing y-coordinate. and throwing away duplicate points. This make future splittings constant time operations. After spliting in the middle and recursively triangulating L and R untiil L and R end up with 2 or 3 sites and be handled separately. When the recursive function coming back (withdraw), we must consider the merge step. There are two lemma which solve the difficulty[1].

1. Let L and R be two sets of points. Any edge of the Delaunay diagram of L U R whose  endpoints
are both in L is in the Delaunay diagram of L.
In other words, the addition of new points does not introduce new edges between the old points.

2. Let T_L and T_R be Delaunay triangulations with vertex sets L and R. Then we can always construct a Delanauy triangulation T for the set L U  R such that every edge of T that is not in T_L or in T_R has one endpoint in L and one in R.

Thus when we perform merge, no new L-L or R-R edges will be added. In other words, merge only add cross edge. All these edges must cross a line parallel to the y-axis and  placed at the splitting x value. This established a linear ordering of the cross edges! Cross edges will produce incrementally in ascending y direction. A rising bubble with the inCircle test will decide which edge to delete. The rising bubble only has a single degree of freedom.
The algorithm runs in time O(n log n) and uses linear storage.

References

1.  Guibas, L. and Stolfi, J., "Primitives for the Manipulation of General Subdivisions and the Computation of Voronoi Diagrams" ACM transactions on Graphics, 4:74-23, April 1985.

2. "Graph Drawing - Algorithms for the Visualization of Graphs", Giuseppe D.B, Peter E. , Roberto Tamassia, Ioannis G.T.
 

3. "Data Structures and Algorithms in Java", Michael T.G. and Roberto Tamassia .