/*
A D implementation of the Voronoi Olden benchmark. Voronoi
generates a random set of points and computes a Voronoi diagram for
the points.
L. Guibas and J. Stolfi. "General Subdivisions and Voronoi Diagrams"
ACM Trans. on Graphics 4(2):74-123, 1985.
The Java version of voronoi (slightly) differs from the C version
in several ways. The C version allocates an array of 4 edges and
uses pointer addition to implement quick rotate operations. The
Java version does not use pointer addition to implement these
operations.
Run it with:
time voronoi1_d -n 100000 -v
*/
version (Tango) {
import tango.stdc.stdio: printf, fprintf, sprintf, stderr;
import tango.stdc.stdlib: exit;
import tango.math.Math: sqrt, exp, log;
import tango.stdc.time: CLOCKS_PER_SEC, clock;
import Integer = tango.text.convert.Integer;
alias Integer.parse toInt;
} else {
import std.c.stdio: printf, fprintf, sprintf, stderr;
import std.c.stdlib: exit;
import std.math: sqrt, exp, log;
import std.c.time: CLOCKS_PER_SEC, clock;
import std.conv: toInt;
}
double myclock() {
return clock() / cast(double)CLOCKS_PER_SEC;
}
class Stack(T) {
T[] data;
int length() { return data.length; }
bool empty() { return data.length == 0; }
void push(T item) { data ~= item; }
T pop() {
assert(data.length);
T last = data[$-1];
data.length = data.length - 1;
return last;
}
}
/**
* A class that represents a wrapper around a double value so
* that we can use it as an 'out' parameter. The java.lang.Double
* class is immutable.
**/
final class MyDouble {
public double value;
this(double d) {
value = d;
}
public char* toCString() {
auto repr = new char[64];
sprintf(repr.ptr, "%f", value);
return repr.ptr;
}
}
/**
* Vector Routines from CMU vision library.
* They are used only for the Voronoi Diagram, not the Delaunay Triagulation.
* They are slow because of large call-by-value parameters.
**/
class Vec2 {
double x,y;
double norm;
this() {}
this(double xx, double yy) {
x = xx;
y = yy;
norm = x * x + y * y;
}
final public double X() {
return x;
}
final public double Y() {
return y;
}
final public double Norm() {
return norm;
}
final public void setNorm(double d) {
norm = d;
}
final public char* toCString() {
auto repr = new char[256];
sprintf(repr.ptr, "%f %f", x, y);
return repr.ptr;
}
final Vec2 circle_center(Vec2 b, Vec2 c) {
Vec2 vv1 = b.sub(c);
double d1 = vv1.magn();
vv1 = sum(b);
Vec2 vv2 = vv1.times(0.5);
if (d1 < 0.0)
// there is no intersection point, the bisectors coincide.
return vv2;
else {
Vec2 vv3 = b.sub(this);
Vec2 vv4 = c.sub(this);
double d3 = vv3.cprod(vv4) ;
double d2 = -2.0 * d3 ;
Vec2 vv5 = c.sub(b);
double d4 = vv5.dot(vv4);
Vec2 vv6 = vv3.cross();
Vec2 vv7 = vv6.times(d4 / d2);
return vv2.sum(vv7);
}
}
/**
* cprod: forms triple scalar product of [u,v,k], where k = u cross v
* (returns the magnitude of u cross v in space)
**/
final double cprod(Vec2 v) {
return x * v.y - y * v.x;
}
/* V2_dot: vector dot product */
final double dot(Vec2 v) {
return x * v.x + y * v.y;
}
/* V2_times: multiply a vector by a scalar */
final Vec2 times(double c) {
return new Vec2(c * x, c * y);
}
/* V2_sum, V2_sub: Vector addition and subtraction */
final Vec2 sum(Vec2 v) {
return new Vec2(x + v.x, y + v.y);
}
final Vec2 sub(Vec2 v) {
return new Vec2(x - v.x,y - v.y);
}
/* V2_magn: magnitude of vector */
final double magn() {
return sqrt(x * x + y * y);
}
/* returns k X v (cross product). this is a vector perpendicular to v */
final Vec2 cross() {
return new Vec2(y, -x);
}
} // Vec2 ends
/**
* A class that represents a voronoi diagram. The diagram is represnted
* as a binary tree of points.
**/
final class Vertex : Vec2 {
/**
* The left and right child of the tree that represents the voronoi diagram.
**/
Vertex left, right;
/**
* Seed value used during tree creation.
**/
static int seed;
this() { }
this(double x, double y) {
super(x, y);
left = null;
right = null;
}
public void setLeft(Vertex l) {
left = l;
}
public void setRight(Vertex r) {
right = r;
}
public Vertex getLeft() {
return left;
}
public Vertex getRight() {
return right;
}
/**
* Generate a voronoi diagram
**/
static Vertex createPoints(int n, MyDouble curmax, int i) {
if (n < 1 )
return null;
Vertex cur = new Vertex();
Vertex right = Vertex.createPoints(n / 2, curmax, i);
i -= n/2;
cur.x = curmax.value * exp(log(Vertex.drand()) / i);
cur.y = Vertex.drand();
cur.norm = cur.x * cur.x + cur.y * cur.y;
cur.right = right;
curmax.value = cur.X();
Vertex left = Vertex.createPoints(n / 2, curmax, i - 1);
cur.left = left;
return cur;
}
/**
* Builds delaunay triangulation.
**/
Edge buildDelaunayTriangulation(Vertex extra) {
EdgePair retVal = buildDelaunay(extra);
return retVal.getLeft();
}
/**
* Recursive delaunay triangulation procedure
* Contains modifications for axis-switching division.
**/
EdgePair buildDelaunay(Vertex extra) {
EdgePair retval = null;
if (getRight() !is null && getLeft() !is null) {
// more than three elements; do recursion
Vertex minx = getLow();
Vertex maxx = extra;
EdgePair delright = getRight().buildDelaunay(extra);
EdgePair delleft = getLeft().buildDelaunay(this);
retval = Edge.doMerge(delleft.getLeft(), delleft.getRight(),
delright.getLeft(), delright.getRight());
Edge ldo = retval.getLeft();
while (ldo.orig() != minx) {
ldo = ldo.rPrev();
}
Edge rdo = retval.getRight();
while (rdo.orig() != maxx) {
rdo = rdo.lPrev();
}
retval = new EdgePair(ldo, rdo);
} else if (getLeft() is null) {
// two points
Edge a = Edge.makeEdge(this, extra);
retval = new EdgePair(a, a.symmetric());
} else {
// left, !right three points
// 3 cases: triangles of 2 orientations, and 3 points on a line. */
Vertex s1 = getLeft();
Vertex s2 = this;
Vertex s3 = extra;
Edge a = Edge.makeEdge(s1, s2);
Edge b = Edge.makeEdge(s2, s3);
a.symmetric().splice(b);
Edge c = b.connectLeft(a);
if (s1.ccw(s3, s2)) {
retval = new EdgePair(c.symmetric(), c);
} else {
retval = new EdgePair(a, b.symmetric());
if (s1.ccw(s2, s3))
c.deleteEdge();
}
}
return retval;
}
/**
* Print the tree
**/
void print() {
Vertex tleft, tright;
printf("X=%f Y=%f\n", X(), Y());
if (left is null)
printf("NULL\n");
else
left.print();
if (right is null)
printf("NULL\n");
else
right.print();
}
/**
* Traverse down the left child to the end
**/
Vertex getLow() {
Vertex temp;
Vertex tree = this;
while ((temp=tree.getLeft()) !is null)
tree = temp;
return tree;
}
/****************************************************************/
/* Geometric primitives
****************************************************************/
bool incircle(Vertex b, Vertex c, Vertex d) {
// incircle, as in the Guibas-Stolfi paper
double adx, ady, bdx, bdy, cdx, cdy, dx, dy, anorm, bnorm, cnorm, dnorm;
double dret;
Vertex loc_a,loc_b,loc_c,loc_d;
int donedx,donedy,donednorm;
loc_d = d;
dx = loc_d.X(); dy = loc_d.Y(); dnorm = loc_d.Norm();
loc_a = this;
adx = loc_a.X() - dx; ady = loc_a.Y() - dy; anorm = loc_a.Norm();
loc_b = b;
bdx = loc_b.X() - dx; bdy = loc_b.Y() - dy; bnorm = loc_b.Norm();
loc_c = c;
cdx = loc_c.X() - dx; cdy = loc_c.Y() - dy; cnorm = loc_c.Norm();
dret = (anorm - dnorm) * (bdx * cdy - bdy * cdx);
dret += (bnorm - dnorm) * (cdx * ady - cdy * adx);
dret += (cnorm - dnorm) * (adx * bdy - ady * bdx);
return (0.0 < dret) ? true : false;
}
bool ccw(Vertex b, Vertex c) {
// TRUE iff this, B, C form a counterclockwise oriented triangle
double dret;
double xa,ya, xb, yb, xc, yc;
Vertex loc_a, loc_b, loc_c;
int donexa, doneya, donexb, doneyb, donexc, doneyc;
loc_a = this;
xa = loc_a.X();
ya = loc_a.Y();
loc_b = b;
xb = loc_b.X();
yb = loc_b.Y();
loc_c = c;
xc = loc_c.X();
yc = loc_c.Y();
dret = (xa-xc) * (yb-yc) - (xb-xc) * (ya-yc);
return (dret > 0.0)? true : false;
}
/**
* A routine used by the random number generator
**/
static int mult(int p,int q) {
int p1, p0, q1, q0;
int CONST_m1 = 10000;
p1 = p / CONST_m1; p0 = p % CONST_m1;
q1 = q / CONST_m1; q0 = q % CONST_m1;
return ((p0 * q1 + p1 * q0) % CONST_m1) * CONST_m1 + p0 * q0;
}
/**
* Generate the nth random number
**/
static int skiprand(int seed, int n) {
for (; n != 0; n--)
seed = random(seed);
return seed;
}
static int random(int seed) {
int CONST_b = 31415821;
seed = mult(seed, CONST_b) + 1;
return seed;
}
static double drand() {
double retval = (cast(double)(Vertex.seed = Vertex.random(Vertex.seed))) /
cast(double)2147483647;
return retval;
}
} // Vertex ends
/**
* A class that represents the quad edge data structure which implements
* the edge algebra as described in the algorithm.
* <p>
* Each edge contains 4 parts, or other edges. Any edge in the group may
* be accessed using a series of rotate and flip operations. The 0th
* edge is the canonical representative of the group.
* <p>
* The quad edge class does not contain separate information for vertice
* or faces; a vertex is implicitly defined as a ring of edges (the ring
* is created using the next field).
**/
final class Edge {
/**
* Group of edges that describe the quad edge
**/
Edge quadList[];
/**
* The position of this edge within the quad list
**/
int listPos;
/**
* The vertex that this edge represents
**/
Vertex vertex;
/**
* Contains a reference to a connected quad edge
**/
Edge next;
/**
* Create a new edge which.
**/
this(Vertex v, Edge[] ql, int pos) {
vertex = v;
quadList = ql;
listPos = pos;
}
/**
* Create a new edge which.
**/
this(Edge[] ql, int pos) {
this(null, ql, pos);
}
/**
* Create a string representation of the edge
**/
public char* toCString() {
if (vertex !is null)
return vertex.toCString();
else
return "None";
}
public static Edge makeEdge(Vertex o, Vertex d) {
Edge ql[] = new Edge[4];
ql[0] = new Edge(ql, 0);
ql[1] = new Edge(ql, 1);
ql[2] = new Edge(ql, 2);
ql[3] = new Edge(ql, 3);
ql[0].next = ql[0];
ql[1].next = ql[3];
ql[2].next = ql[2];
ql[3].next = ql[1];
Edge base = ql[0];
base.setOrig(o);
base.setDest(d);
return base;
}
public void setNext(Edge n) {
next = n;
}
/**
* Initialize the data (vertex) for the edge's origin
**/
public void setOrig(Vertex o) {
vertex = o;
}
/**
* Initialize the data (vertex) for the edge's destination
**/
public void setDest(Vertex d) {
symmetric().setOrig(d);
}
Edge oNext() {
return next;
}
Edge oPrev() {
return this.rotate().oNext().rotate();
}
Edge lNext() {
return this.rotateInv().oNext().rotate();
}
Edge lPrev() {
return this.oNext().symmetric();
}
Edge rNext() {
return this.rotate().oNext().rotateInv();
}
Edge rPrev() {
return this.symmetric().oNext();
}
Edge dNext() {
return this.symmetric().oNext().symmetric();
}
Edge dPrev() {
return this.rotateInv().oNext().rotateInv();
}
Vertex orig() {
return vertex;
}
Vertex dest() {
return symmetric().orig();
}
/**
* Return the symmetric of the edge. The symmetric is the same edge
* with the opposite direction.
* @return the symmetric of the edge
**/
Edge symmetric() {
return quadList[(listPos + 2) % 4];
}
/**
* Return the rotated version of the edge. The rotated version is a
* 90 degree counterclockwise turn.
* @return the rotated version of the edge
**/
Edge rotate() {
return quadList[(listPos + 1) % 4];
}
/**
* Return the inverse rotated version of the edge. The inverse
* is a 90 degree clockwise turn.
* @return the inverse rotated edge.
**/
Edge rotateInv() {
return quadList[(listPos + 3) % 4];
}
Edge nextQuadEdge() {
return quadList[(listPos + 1) % 4];
}
Edge connectLeft(Edge b) {
Vertex t1,t2;
Edge ans, lnexta;
t1 = dest();
lnexta = lNext();
t2 = b.orig();
ans = Edge.makeEdge(t1, t2);
ans.splice(lnexta);
ans.symmetric().splice(b);
return ans;
}
Edge connectRight(Edge b) {
Vertex t1, t2;
Edge ans, oprevb, q1;
t1 = dest();
t2 = b.orig();
oprevb = b.oPrev();
ans = Edge.makeEdge(t1, t2);
ans.splice(symmetric());
ans.symmetric().splice(oprevb);
return ans;
}
/****************************************************************/
/* Quad-edge manipulation primitives
****************************************************************/
void swapedge() {
Edge a = oPrev();
Edge syme = symmetric();
Edge b = syme.oPrev();
splice(a);
syme.splice(b);
Edge lnexttmp = a.lNext();
splice(lnexttmp);
lnexttmp = b.lNext();
syme.splice(lnexttmp);
Vertex a1 = a.dest();
Vertex b1 = b.dest();
setOrig(a1);
setDest(b1);
}
void splice(Edge b) {
Edge alpha = oNext().rotate();
Edge beta = b.oNext().rotate();
Edge t1 = beta.oNext();
Edge temp = alpha.oNext();
alpha.setNext(t1);
beta.setNext(temp);
temp = oNext();
t1 = b.oNext();
b.setNext(temp);
setNext(t1);
}
bool valid(Edge basel) {
Vertex t1 = basel.orig();
Vertex t3 = basel.dest();
Vertex t2 = dest();
return t1.ccw(t2, t3);
}
void deleteEdge() {
Edge f = oPrev();
splice(f);
f = symmetric().oPrev();
symmetric().splice(f);
}
static EdgePair doMerge(Edge ldo, Edge ldi, Edge rdi, Edge rdo) {
while (true) {
Vertex t3 = rdi.orig();
Vertex t1 = ldi.orig();
Vertex t2 = ldi.dest();
while (t1.ccw(t2, t3)) {
ldi = ldi.lNext();
t1=ldi.orig();
t2=ldi.dest();
}
t2 = rdi.dest();
if (t2.ccw(t3, t1))
rdi = rdi.rPrev();
else
break;
}
Edge basel = rdi.symmetric().connectLeft(ldi);
Edge lcand = basel.rPrev();
Edge rcand = basel.oPrev();
Vertex t1 = basel.orig();
Vertex t2 = basel.dest();
if (t1 == rdo.orig())
rdo = basel;
if (t2 == ldo.orig())
ldo = basel.symmetric();
while (true) {
Edge t = lcand.oNext();
if (t.valid(basel)) {
Vertex v4 = basel.orig();
Vertex v1 = lcand.dest();
Vertex v3 = lcand.orig();
Vertex v2 = t.dest();
while (v1.incircle(v2,v3,v4)){
lcand.deleteEdge();
lcand = t;
t = lcand.oNext();
v1 = lcand.dest();
v3 = lcand.orig();
v2 = t.dest();
}
}
t = rcand.oPrev();
if (t.valid(basel)) {
Vertex v4 = basel.dest();
Vertex v1 = t.dest();
Vertex v2 = rcand.dest();
Vertex v3 = rcand.orig();
while (v1.incircle(v2, v3, v4)) {
rcand.deleteEdge();
rcand = t;
t = rcand.oPrev();
v2 = rcand.dest();
v3 = rcand.orig();
v1 = t.dest();
}
}
bool lvalid = lcand.valid(basel);
bool rvalid = rcand.valid(basel);
if ((!lvalid) && (!rvalid))
return new EdgePair(ldo, rdo);
Vertex v1 = lcand.dest();
Vertex v2 = lcand.orig();
Vertex v3 = rcand.orig();
Vertex v4 = rcand.dest();
if (!lvalid || (rvalid && v1.incircle(v2, v3, v4))) {
basel = rcand.connectLeft(basel.symmetric());
rcand = basel.symmetric().lNext();
} else {
basel = lcand.connectRight(basel).symmetric();
lcand = basel.rPrev();
}
}
assert(0);
}
/**
* Print the voronoi diagram and its dual, the delaunay triangle for the
* diagram.
**/
void outputVoronoiDiagram() {
Edge nex = this;
// Plot voronoi diagram edges with one endpoint at infinity.
do {
Vec2 v21 = nex.dest();
Vec2 v22 = nex.orig();
Edge tmp = nex.oNext();
Vec2 v23 = tmp.dest();
Vec2 cvxvec = v21.sub(v22);
Vec2 center = v22.circle_center(v21, v23);
Vec2 vv1 = v22.sum(v22);
Vec2 vv2 = vv1.times(0.5);
Vec2 vv3 = center.sub(vv2);
double ln = 1.0 + vv3.magn();
double d1 = ln / cvxvec.magn();
vv1 = cvxvec.cross();
vv2 = vv1.times(d1) ;
vv3 = center.sum(vv2);
printf("Vedge %s %s\n", center.toCString(), vv3.toCString());
nex = nex.rNext();
} while (nex != this);
// plot delaunay triangle edges and finite VD edges.
Stack!(Edge) edges = new Stack!(Edge);
bool[Edge] seen;
pushRing(edges, seen);
printf("no. of edges = %d\n", edges.length);
while (!edges.empty()) {
Edge edge = edges.pop();
if (edge in seen && seen[edge]) {
Edge prev = edge;
nex = edge.oNext();
do {
Vertex v1 = prev.orig();
double d1 = v1.X();
Vertex v2 = prev.dest();
double d2 = v2.X();
if (d1 >= d2) {
printf("Dedge %s %s\n", v1.toCString(), v2.toCString());
Edge sprev = prev.symmetric();
Edge snex = sprev.oNext();
v1 = prev.orig();
v2 = prev.dest();
Vertex v3 = nex.dest();
Vertex v4 = snex.dest();
if (v1.ccw(v2, v3) != v1.ccw(v2, v4)) {
Vec2 v21 = prev.orig();
Vec2 v22 = prev.dest();
Vec2 v23 = nex.dest();
Vec2 vv1 = v21.circle_center(v22, v23);
v21 = sprev.orig();
v22 = sprev.dest();
v23 = snex.dest();
Vec2 vv2 = v21.circle_center(v22, v23);
printf("Vedge %s %s\n", vv1.toCString(), vv2.toCString());
}
}
seen[prev] = false;
prev = nex;
nex = nex.oNext();
} while (prev != edge);
}
edge.symmetric().pushRing(edges, seen);
}
}
void pushRing(Stack!(Edge) stack, ref bool[Edge] seen) {
Edge nex = oNext();
while (nex != this) {
if (!(nex in seen)) {
seen[nex] = true;
stack.push(nex);
}
nex = nex.oNext();
}
}
void pushNonezeroRing(Stack!(Edge) stack, ref bool[Edge] seen) {
Edge nex = oNext();
while (nex != this) {
if (nex in seen) {
seen.remove(nex);
stack.push(nex);
}
nex = nex.oNext();
}
}
} // Edge ends
/**
* A class that represents an edge pair
**/
final class EdgePair {
Edge left;
Edge right;
this(Edge l, Edge r) {
left = l;
right = r;
}
public Edge getLeft() {
return left;
}
public Edge getRight() {
return right;
}
} // EdgePair ends
struct Voronoi {
/**
* The number of points in the diagram
**/
private static int points = 0;
/**
* Set to true to print informative messages
**/
private static bool printMsgs = false;
/**
* Set to true to print the voronoi diagram and its dual,
* the delaunay diagram
**/
private static bool printResults = false;
/**
* The main routine which creates the points and then performs
* the delaunay triagulation.
* @param args the command line parameters
**/
public static void main(char[][] args) {
parseCmdLine(args);
if (printMsgs)
printf("Getting %d points\n", points);
auto start0 = myclock();
Vertex.seed = 1023;
Vertex extra = Vertex.createPoints(1, new MyDouble(1.0), points);
Vertex point = Vertex.createPoints(points-1, new MyDouble(extra.X()),
points-1);
auto end0 = myclock();
if (printMsgs)
printf("Doing voronoi on %d nodes\n", points);
auto start1 = myclock();
Edge edge = point.buildDelaunayTriangulation(extra);
auto end1 = myclock();
if (printResults)
edge.outputVoronoiDiagram();
if (printMsgs) {
printf("Build time %f\n", end0 - start0);
printf("Compute time %f\n", end1 - start1);
printf("Total time %f\n", end1 - start0);
}
printf("Done!\n");
}
/**
* Parse the command line options.
* @param args the command line options.
**/
private static final void parseCmdLine(char[][] args) {
int i = 1;
while (i < args.length && args[i][0] == '-') {
char[] arg = args[i++];
if (arg == "-n") {
if (i < args.length)
points = toInt(args[i++]);
else
throw new Exception("-n requires the number of points");
} else if (arg == "-p") {
printResults = true;
} else if (arg == "-v") {
printMsgs = true;
} else if (arg == "-h") {
usage();
}
}
if (points == 0)
usage();
}
/**
* The usage routine which describes the program options.
**/
private static final void usage() {
fprintf(stderr, "usage: voronoi_d -n <points> [-p] [-m] [-h]\n");
fprintf(stderr, " -n the number of points in the diagram\n");
fprintf(stderr, " -p (print detailed results/messages - the voronoi diagram>)\n");
fprintf(stderr, " -v (print informative message)\n");
fprintf(stderr, " -h (this message)\n");
exit(0);
}
} // Voronoi ends
void main(char[][] args) {
Voronoi.main(args);
}