source/blog/posts/2010/08/14-memory-allocation/voronoi.d

   1 /*
   2 Translated by Leonardo Fantascienza, downloaded from http://codepad.org/xGDCS3KO
   3 Modified by Leandro Lucarella to be really quiet when -v is not used.
   4
   5 A D implementation of the Voronoi Olden benchmark. Voronoi
   6 generates a random set of points and computes a Voronoi diagram for
   7 the points.
   8
   9 L. Guibas and J. Stolfi. "General Subdivisions and Voronoi Diagrams"
  10 ACM Trans. on Graphics 4(2):74-123, 1985.
  11
  12 The Java version of voronoi (slightly) differs from the C version
  13 in several ways.  The C version allocates an array of 4 edges and
  14 uses pointer addition to implement quick rotate operations.  The
  15 Java version does not use pointer addition to implement these
  16 operations.
  17
  18 Run it with:
  19 time voronoi1_d -n 100000 -v
  20 */
  21
  22 version (Tango) {
  23     import tango.stdc.stdio: printf, fprintf, sprintf, stderr;
  24     import tango.stdc.stdlib: exit;
  25     import tango.math.Math: sqrt, exp, log;
  26     import tango.stdc.time: CLOCKS_PER_SEC, clock;
  27     import Integer = tango.text.convert.Integer;
  28     alias Integer.parse toInt;
  29 } else {
  30     import std.c.stdio: printf, fprintf, sprintf, stderr;
  31     import std.c.stdlib: exit;
  32     import std.math: sqrt, exp, log;
  33     import std.c.time: CLOCKS_PER_SEC, clock;
  34     import std.conv: toInt;
  35 }
  36
  37
  38 double myclock() {
  39     return clock() / cast(double)CLOCKS_PER_SEC;
  40 }
  41
  42
  43 class Stack(T) {
  44     T[] data;
  45     int length() { return data.length; }
  46     bool empty() { return data.length == 0; }
  47     void push(T item) { data ~= item; }
  48     T pop() {
  49         assert(data.length);
  50         T last = data[$-1];
  51         data.length = data.length - 1;
  52         return last;
  53     }
  54 }
  55
  56
  57 /**
  58 * A class that represents a wrapper around a double value so
  59 * that we can use it as an 'out' parameter.  The java.lang.Double
  60 * class is immutable.
  61 **/
  62 final class MyDouble {
  63     public double value;
  64
  65     this(double d) {
  66         value = d;
  67     }
  68
  69     public char* toCString() {
  70         auto repr = new char[64];
  71         sprintf(repr.ptr, "%f", value);
  72         return repr.ptr;
  73     }
  74 }
  75
  76
  77 /**
  78 * Vector Routines from CMU vision library.
  79 * They are used only for the Voronoi Diagram, not the Delaunay Triagulation.
  80 * They are slow because of large call-by-value parameters.
  81 **/
  82 class Vec2 {
  83     double x,y;
  84     double norm;
  85
  86     this() {}
  87
  88     this(double xx, double yy) {
  89         x = xx;
  90         y = yy;
  91         norm =  x * x + y * y;
  92     }
  93
  94     final public double X() {
  95         return x;
  96     }
  97
  98     final public double Y() {
  99         return y;
 100     }
 101
 102     final public double Norm() {
 103         return norm;
 104     }
 105
 106     final public void setNorm(double d) {
 107         norm = d;
 108     }
 109
 110     final public char* toCString() {
 111         auto repr = new char[256];
 112         sprintf(repr.ptr, "%f %f", x, y);
 113         return repr.ptr;
 114     }
 115
 116     final Vec2 circle_center(Vec2 b, Vec2 c) {
 117         Vec2 vv1 = b.sub(c);
 118         double d1 = vv1.magn();
 119         vv1 = sum(b);
 120         Vec2 vv2 = vv1.times(0.5);
 121         if (d1 < 0.0)
 122           // there is no intersection point, the bisectors coincide.
 123             return vv2;
 124         else {
 125             Vec2 vv3 = b.sub(this);
 126             Vec2 vv4 = c.sub(this);
 127             double d3 = vv3.cprod(vv4) ;
 128             double d2 = -2.0 * d3 ;
 129             Vec2 vv5 = c.sub(b);
 130             double d4 = vv5.dot(vv4);
 131             Vec2 vv6 = vv3.cross();
 132             Vec2 vv7 = vv6.times(d4 / d2);
 133             return vv2.sum(vv7);
 134         }
 135     }
 136
 137     /**
 138     * cprod: forms triple scalar product of [u,v,k], where k = u cross v
 139     * (returns the magnitude of u cross v in space)
 140     **/
 141     final double cprod(Vec2 v) {
 142         return x * v.y - y * v.x;
 143     }
 144
 145     /* V2_dot: vector dot product */
 146     final double dot(Vec2 v) {
 147         return x * v.x + y * v.y;
 148     }
 149
 150     /* V2_times: multiply a vector by a scalar */
 151     final Vec2 times(double c) {
 152         return new Vec2(c * x, c * y);
 153     }
 154
 155     /* V2_sum, V2_sub: Vector addition and subtraction */
 156     final Vec2 sum(Vec2 v) {
 157         return new Vec2(x + v.x, y + v.y);
 158     }
 159
 160     final Vec2 sub(Vec2 v) {
 161         return new Vec2(x - v.x,y - v.y);
 162     }
 163
 164     /* V2_magn: magnitude of vector */
 165     final double magn() {
 166         return sqrt(x * x + y * y);
 167     }
 168
 169     /* returns k X v (cross product).  this is a vector perpendicular to v */
 170     final Vec2 cross() {
 171         return new Vec2(y, -x);
 172     }
 173 } // Vec2 ends
 174
 175
 176 /**
 177 * A class that represents a voronoi diagram.  The diagram is represnted
 178 * as a binary tree of points.
 179 **/
 180 final class Vertex : Vec2 {
 181     /**
 182     * The left and right child of the tree that represents the voronoi diagram.
 183     **/
 184     Vertex left, right;
 185
 186     /**
 187     * Seed value used during tree creation.
 188     **/
 189     static int seed;
 190
 191     this() { }
 192
 193     this(double x, double y) {
 194         super(x, y);
 195         left = null;
 196         right = null;
 197     }
 198
 199     public void setLeft(Vertex l) {
 200         left = l;
 201     }
 202
 203     public void setRight(Vertex r) {
 204         right = r;
 205     }
 206
 207     public Vertex getLeft() {
 208         return left;
 209     }
 210
 211     public Vertex getRight() {
 212         return right;
 213     }
 214
 215     /**
 216     * Generate a voronoi diagram
 217     **/
 218     static Vertex createPoints(int n, MyDouble curmax, int i) {
 219         if (n < 1 )
 220             return null;
 221
 222         Vertex cur = new Vertex();
 223         Vertex right = Vertex.createPoints(n / 2, curmax, i);
 224         i -= n/2;
 225         cur.x = curmax.value * exp(log(Vertex.drand()) / i);
 226         cur.y = Vertex.drand();
 227         cur.norm = cur.x * cur.x + cur.y * cur.y;
 228         cur.right = right;
 229         curmax.value = cur.X();
 230         Vertex left = Vertex.createPoints(n / 2, curmax, i - 1);
 231         cur.left = left;
 232         return cur;
 233     }
 234
 235     /**
 236     * Builds delaunay triangulation.
 237     **/
 238     Edge buildDelaunayTriangulation(Vertex extra) {
 239         EdgePair retVal = buildDelaunay(extra);
 240         return retVal.getLeft();
 241     }
 242
 243     /**
 244     * Recursive delaunay triangulation procedure
 245     * Contains modifications for axis-switching division.
 246     **/
 247     EdgePair buildDelaunay(Vertex extra) {
 248         EdgePair retval = null;
 249         if (getRight() !is null && getLeft() !is null)  {
 250             // more than three elements; do recursion
 251             Vertex minx = getLow();
 252             Vertex maxx = extra;
 253
 254             EdgePair delright = getRight().buildDelaunay(extra);
 255             EdgePair delleft = getLeft().buildDelaunay(this);
 256
 257             retval = Edge.doMerge(delleft.getLeft(), delleft.getRight(),
 258             delright.getLeft(), delright.getRight());
 259
 260             Edge ldo = retval.getLeft();
 261             while (ldo.orig() != minx) {
 262                 ldo = ldo.rPrev();
 263             }
 264             Edge rdo = retval.getRight();
 265             while (rdo.orig() != maxx) {
 266                 rdo = rdo.lPrev();
 267             }
 268
 269             retval = new EdgePair(ldo, rdo);
 270
 271         } else if (getLeft() is null) {
 272             // two points
 273             Edge a = Edge.makeEdge(this, extra);
 274             retval = new EdgePair(a, a.symmetric());
 275         } else {
 276             // left, !right  three points
 277             // 3 cases: triangles of 2 orientations, and 3 points on a line. */
 278             Vertex s1 = getLeft();
 279             Vertex s2 = this;
 280             Vertex s3 = extra;
 281             Edge a = Edge.makeEdge(s1, s2);
 282             Edge b = Edge.makeEdge(s2, s3);
 283             a.symmetric().splice(b);
 284             Edge c = b.connectLeft(a);
 285             if (s1.ccw(s3, s2)) {
 286                 retval = new EdgePair(c.symmetric(), c);
 287             } else {
 288                 retval = new EdgePair(a, b.symmetric());
 289                 if (s1.ccw(s2, s3))
 290                     c.deleteEdge();
 291             }
 292         }
 293
 294         return retval;
 295     }
 296
 297     /**
 298     * Print the tree
 299     **/
 300     void print() {
 301         Vertex tleft, tright;
 302
 303         printf("X=%f  Y=%f\n", X(), Y());
 304         if (left is null)
 305             printf("NULL\n");
 306         else
 307             left.print();
 308         if (right is null)
 309             printf("NULL\n");
 310         else
 311             right.print();
 312     }
 313
 314     /**
 315     * Traverse down the left child to the end
 316     **/
 317     Vertex getLow() {
 318         Vertex temp;
 319         Vertex tree = this;
 320
 321         while ((temp=tree.getLeft()) !is null)
 322             tree = temp;
 323         return tree;
 324     }
 325
 326     /****************************************************************/
 327     /*    Geometric primitives
 328     ****************************************************************/
 329     bool incircle(Vertex b, Vertex c, Vertex d) {
 330         // incircle, as in the Guibas-Stolfi paper
 331         double adx, ady, bdx, bdy, cdx, cdy, dx, dy, anorm, bnorm, cnorm, dnorm;
 332         double dret;
 333         Vertex loc_a,loc_b,loc_c,loc_d;
 334
 335         int donedx,donedy,donednorm;
 336         loc_d = d;
 337         dx = loc_d.X(); dy = loc_d.Y(); dnorm = loc_d.Norm();
 338         loc_a = this;
 339         adx = loc_a.X() - dx; ady = loc_a.Y() - dy; anorm = loc_a.Norm();
 340         loc_b = b;
 341         bdx = loc_b.X() - dx; bdy = loc_b.Y() - dy; bnorm = loc_b.Norm();
 342         loc_c = c;
 343         cdx = loc_c.X() - dx; cdy = loc_c.Y() - dy; cnorm = loc_c.Norm();
 344         dret =  (anorm - dnorm) * (bdx * cdy - bdy * cdx);
 345         dret += (bnorm - dnorm) * (cdx * ady - cdy * adx);
 346         dret += (cnorm - dnorm) * (adx * bdy - ady * bdx);
 347         return (0.0 < dret) ? true : false;
 348     }
 349
 350     bool ccw(Vertex b, Vertex c) {
 351         // TRUE iff this, B, C form a counterclockwise oriented triangle
 352         double dret;
 353         double xa,ya, xb, yb, xc, yc;
 354         Vertex loc_a, loc_b, loc_c;
 355         int donexa, doneya, donexb, doneyb, donexc, doneyc;
 356
 357         loc_a = this;
 358         xa = loc_a.X();
 359         ya = loc_a.Y();
 360         loc_b = b;
 361         xb = loc_b.X();
 362         yb = loc_b.Y();
 363         loc_c = c;
 364         xc = loc_c.X();
 365         yc = loc_c.Y();
 366         dret = (xa-xc) * (yb-yc) - (xb-xc) * (ya-yc);
 367         return (dret  > 0.0)? true : false;
 368     }
 369
 370     /**
 371     * A routine used by the random number generator
 372     **/
 373     static int mult(int p,int q) {
 374         int p1, p0, q1, q0;
 375         int CONST_m1 = 10000;
 376
 377         p1 = p / CONST_m1; p0 = p % CONST_m1;
 378         q1 = q / CONST_m1; q0 = q % CONST_m1;
 379         return ((p0 * q1 + p1 * q0) % CONST_m1) * CONST_m1 + p0 * q0;
 380     }
 381
 382     /**
 383     * Generate the nth random number
 384     **/
 385     static int skiprand(int seed, int n) {
 386         for (; n != 0; n--)
 387             seed = random(seed);
 388         return seed;
 389     }
 390
 391     static int random(int seed) {
 392         int CONST_b = 31415821;
 393
 394         seed = mult(seed, CONST_b) + 1;
 395         return seed;
 396     }
 397
 398     static double drand() {
 399         double retval = (cast(double)(Vertex.seed = Vertex.random(Vertex.seed))) /
 400                         cast(double)2147483647;
 401         return retval;
 402     }
 403 } // Vertex ends
 404
 405
 406 /**
 407 * A class that represents the quad edge data structure which implements
 408 * the edge algebra as described in the algorithm.
 409 * <p>
 410 * Each edge contains 4 parts, or other edges.  Any edge in the group may
 411 * be accessed using a series of rotate and flip operations.  The 0th
 412 * edge is the canonical representative of the group.
 413 * <p>
 414 * The quad edge class does not contain separate information for vertice
 415 * or faces; a vertex is implicitly defined as a ring of edges (the ring
 416 * is created using the next field).
 417 **/
 418 final class Edge {
 419     /**
 420     * Group of edges that describe the quad edge
 421     **/
 422     Edge quadList[];
 423
 424     /**
 425     * The position of this edge within the quad list
 426     **/
 427     int     listPos;
 428
 429     /**
 430     * The vertex that this edge represents
 431     **/
 432     Vertex  vertex;
 433
 434     /**
 435     * Contains a reference to a connected quad edge
 436     **/
 437     Edge    next;
 438
 439     /**
 440     * Create a new edge which.
 441     **/
 442     this(Vertex v, Edge[] ql, int pos) {
 443         vertex = v;
 444         quadList = ql;
 445         listPos = pos;
 446     }
 447
 448     /**
 449     * Create a new edge which.
 450     **/
 451     this(Edge[] ql, int pos) {
 452         this(null, ql, pos);
 453     }
 454
 455     /**
 456     * Create a string representation of the edge
 457     **/
 458     public char* toCString() {
 459         if (vertex !is null)
 460             return vertex.toCString();
 461         else
 462             return "None";
 463     }
 464
 465     public static Edge makeEdge(Vertex o, Vertex d) {
 466         Edge ql[] = new Edge[4];
 467         ql[0] = new Edge(ql, 0);
 468         ql[1] = new Edge(ql, 1);
 469         ql[2] = new Edge(ql, 2);
 470         ql[3] = new Edge(ql, 3);
 471
 472         ql[0].next = ql[0];
 473         ql[1].next = ql[3];
 474         ql[2].next = ql[2];
 475         ql[3].next = ql[1];
 476
 477         Edge base = ql[0];
 478         base.setOrig(o);
 479         base.setDest(d);
 480         return base;
 481     }
 482
 483     public void setNext(Edge n) {
 484         next = n;
 485     }
 486
 487     /**
 488     * Initialize the data (vertex) for the edge's origin
 489     **/
 490     public void setOrig(Vertex o) {
 491         vertex = o;
 492     }
 493
 494     /**
 495     * Initialize the data (vertex) for the edge's destination
 496     **/
 497     public void setDest(Vertex d) {
 498         symmetric().setOrig(d);
 499     }
 500
 501     Edge oNext() {
 502         return next;
 503     }
 504
 505     Edge oPrev() {
 506         return this.rotate().oNext().rotate();
 507     }
 508
 509     Edge lNext() {
 510         return this.rotateInv().oNext().rotate();
 511     }
 512
 513     Edge lPrev() {
 514         return this.oNext().symmetric();
 515     }
 516
 517     Edge rNext() {
 518         return this.rotate().oNext().rotateInv();
 519     }
 520
 521     Edge rPrev() {
 522         return this.symmetric().oNext();
 523     }
 524
 525     Edge dNext() {
 526         return this.symmetric().oNext().symmetric();
 527     }
 528
 529     Edge dPrev() {
 530         return this.rotateInv().oNext().rotateInv();
 531     }
 532
 533     Vertex orig() {
 534         return vertex;
 535     }
 536
 537     Vertex dest() {
 538         return symmetric().orig();
 539     }
 540
 541     /**
 542     * Return the symmetric of the edge.  The symmetric is the same edge
 543     * with the opposite direction.
 544     * @return the symmetric of the edge
 545     **/
 546     Edge symmetric() {
 547         return quadList[(listPos + 2) % 4];
 548     }
 549
 550     /**
 551     * Return the rotated version of the edge.  The rotated version is a
 552     * 90 degree counterclockwise turn.
 553     * @return the rotated version of the edge
 554     **/
 555     Edge rotate() {
 556         return quadList[(listPos + 1) % 4];
 557     }
 558
 559     /**
 560     * Return the inverse rotated version of the edge.  The inverse
 561     * is a 90 degree clockwise turn.
 562     * @return the inverse rotated edge.
 563     **/
 564     Edge rotateInv() {
 565         return quadList[(listPos + 3) % 4];
 566     }
 567
 568     Edge nextQuadEdge() {
 569         return quadList[(listPos + 1) % 4];
 570     }
 571
 572     Edge connectLeft(Edge b) {
 573         Vertex t1,t2;
 574         Edge ans, lnexta;
 575
 576         t1 = dest();
 577         lnexta = lNext();
 578         t2 = b.orig();
 579         ans = Edge.makeEdge(t1, t2);
 580         ans.splice(lnexta);
 581         ans.symmetric().splice(b);
 582         return ans;
 583     }
 584
 585     Edge connectRight(Edge b) {
 586         Vertex t1, t2;
 587         Edge ans, oprevb, q1;
 588
 589         t1 = dest();
 590         t2 = b.orig();
 591         oprevb = b.oPrev();
 592
 593         ans = Edge.makeEdge(t1, t2);
 594         ans.splice(symmetric());
 595         ans.symmetric().splice(oprevb);
 596         return ans;
 597     }
 598
 599     /****************************************************************/
 600     /*    Quad-edge manipulation primitives
 601     ****************************************************************/
 602     void swapedge() {
 603         Edge a = oPrev();
 604         Edge syme = symmetric();
 605         Edge b = syme.oPrev();
 606         splice(a);
 607         syme.splice(b);
 608         Edge lnexttmp = a.lNext();
 609         splice(lnexttmp);
 610         lnexttmp = b.lNext();
 611         syme.splice(lnexttmp);
 612         Vertex a1 = a.dest();
 613         Vertex b1 = b.dest();
 614         setOrig(a1);
 615         setDest(b1);
 616     }
 617
 618     void splice(Edge b) {
 619         Edge alpha = oNext().rotate();
 620         Edge beta = b.oNext().rotate();
 621         Edge t1 = beta.oNext();
 622         Edge temp = alpha.oNext();
 623         alpha.setNext(t1);
 624         beta.setNext(temp);
 625         temp = oNext();
 626         t1 = b.oNext();
 627         b.setNext(temp);
 628         setNext(t1);
 629     }
 630
 631     bool valid(Edge basel) {
 632         Vertex t1 = basel.orig();
 633         Vertex t3 = basel.dest();
 634         Vertex t2 = dest();
 635         return t1.ccw(t2, t3);
 636     }
 637
 638     void deleteEdge() {
 639         Edge f = oPrev();
 640         splice(f);
 641         f = symmetric().oPrev();
 642         symmetric().splice(f);
 643     }
 644
 645     static EdgePair doMerge(Edge ldo, Edge ldi, Edge rdi, Edge rdo) {
 646         while (true) {
 647             Vertex t3 = rdi.orig();
 648             Vertex t1 = ldi.orig();
 649             Vertex t2 = ldi.dest();
 650
 651             while (t1.ccw(t2, t3)) {
 652                 ldi = ldi.lNext();
 653
 654                 t1=ldi.orig();
 655                 t2=ldi.dest();
 656             }
 657
 658             t2 = rdi.dest();
 659
 660             if (t2.ccw(t3, t1))
 661                 rdi = rdi.rPrev();
 662             else
 663                 break;
 664         }
 665
 666         Edge basel = rdi.symmetric().connectLeft(ldi);
 667
 668         Edge lcand = basel.rPrev();
 669         Edge rcand = basel.oPrev();
 670         Vertex t1 = basel.orig();
 671         Vertex t2 = basel.dest();
 672
 673         if (t1 == rdo.orig())
 674             rdo = basel;
 675         if (t2 == ldo.orig())
 676             ldo = basel.symmetric();
 677
 678         while (true) {
 679             Edge t = lcand.oNext();
 680             if (t.valid(basel)) {
 681                 Vertex v4 = basel.orig();
 682
 683                 Vertex v1 = lcand.dest();
 684                 Vertex v3 = lcand.orig();
 685                 Vertex v2 = t.dest();
 686                 while (v1.incircle(v2,v3,v4)){
 687                     lcand.deleteEdge();
 688                     lcand = t;
 689
 690                     t = lcand.oNext();
 691                     v1 = lcand.dest();
 692                     v3 = lcand.orig();
 693                     v2 = t.dest();
 694                 }
 695             }
 696
 697             t = rcand.oPrev();
 698             if (t.valid(basel)) {
 699                 Vertex v4 = basel.dest();
 700                 Vertex v1 = t.dest();
 701                 Vertex v2 = rcand.dest();
 702                 Vertex v3 = rcand.orig();
 703                 while (v1.incircle(v2, v3, v4)) {
 704                     rcand.deleteEdge();
 705                     rcand = t;
 706                     t = rcand.oPrev();
 707                     v2 = rcand.dest();
 708                     v3 = rcand.orig();
 709                     v1 = t.dest();
 710                 }
 711             }
 712
 713             bool lvalid = lcand.valid(basel);
 714
 715             bool rvalid = rcand.valid(basel);
 716             if ((!lvalid) && (!rvalid))
 717                 return new EdgePair(ldo, rdo);
 718
 719             Vertex v1 = lcand.dest();
 720             Vertex v2 = lcand.orig();
 721             Vertex v3 = rcand.orig();
 722             Vertex v4 = rcand.dest();
 723             if (!lvalid || (rvalid && v1.incircle(v2, v3, v4))) {
 724                 basel = rcand.connectLeft(basel.symmetric());
 725                 rcand = basel.symmetric().lNext();
 726             } else {
 727                 basel = lcand.connectRight(basel).symmetric();
 728                 lcand = basel.rPrev();
 729             }
 730         }
 731
 732         assert(0);
 733     }
 734
 735
 736     /**
 737     * Print the voronoi diagram and its dual, the delaunay triangle for the
 738     * diagram.
 739     **/
 740     void outputVoronoiDiagram() {
 741         Edge nex = this;
 742         //  Plot voronoi diagram edges with one endpoint at infinity.
 743         do {
 744             Vec2 v21 = nex.dest();
 745             Vec2 v22 = nex.orig();
 746             Edge tmp = nex.oNext();
 747             Vec2 v23 = tmp.dest();
 748             Vec2 cvxvec = v21.sub(v22);
 749             Vec2 center = v22.circle_center(v21, v23);
 750
 751
 752             Vec2 vv1 = v22.sum(v22);
 753             Vec2 vv2 = vv1.times(0.5);
 754             Vec2 vv3 = center.sub(vv2);
 755             double ln = 1.0 + vv3.magn();
 756             double d1 = ln / cvxvec.magn();
 757             vv1 = cvxvec.cross();
 758             vv2 = vv1.times(d1) ;
 759             vv3 = center.sum(vv2);
 760             printf("Vedge %s %s\n", center.toCString(), vv3.toCString());
 761             nex = nex.rNext();
 762         } while (nex != this);
 763
 764         // plot delaunay triangle edges and finite VD edges.
 765         Stack!(Edge) edges = new Stack!(Edge);
 766         bool[Edge] seen;
 767         pushRing(edges, seen);
 768         printf("no. of edges = %d\n", edges.length);
 769         while (!edges.empty()) {
 770             Edge edge = edges.pop();
 771             if (edge in seen && seen[edge]) {
 772                 Edge prev = edge;
 773                 nex = edge.oNext();
 774                 do {
 775                     Vertex v1 = prev.orig();
 776                     double d1 = v1.X();
 777                     Vertex v2 = prev.dest();
 778                     double d2 = v2.X();
 779                     if (d1 >= d2) {
 780                         printf("Dedge %s %s\n", v1.toCString(), v2.toCString());
 781                         Edge sprev = prev.symmetric();
 782                         Edge snex = sprev.oNext();
 783                         v1 = prev.orig();
 784                         v2 = prev.dest();
 785                         Vertex v3 = nex.dest();
 786                         Vertex v4 = snex.dest();
 787                         if (v1.ccw(v2, v3) != v1.ccw(v2, v4)) {
 788                             Vec2 v21 = prev.orig();
 789                             Vec2 v22 = prev.dest();
 790                             Vec2 v23 = nex.dest();
 791                             Vec2 vv1 = v21.circle_center(v22, v23);
 792                             v21 = sprev.orig();
 793                             v22 = sprev.dest();
 794                             v23 = snex.dest();
 795                             Vec2 vv2 = v21.circle_center(v22, v23);
 796                             printf("Vedge %s %s\n", vv1.toCString(), vv2.toCString());
 797                         }
 798                     }
 799                     seen[prev] = false;
 800                     prev = nex;
 801                     nex = nex.oNext();
 802                 } while (prev != edge);
 803             }
 804             edge.symmetric().pushRing(edges, seen);
 805         }
 806     }
 807
 808     void pushRing(Stack!(Edge) stack, ref bool[Edge] seen) {
 809         Edge nex = oNext();
 810         while (nex != this) {
 811             if (!(nex in seen)) {
 812                 seen[nex] = true;
 813                 stack.push(nex);
 814             }
 815             nex = nex.oNext();
 816         }
 817     }
 818
 819     void pushNonezeroRing(Stack!(Edge) stack, ref bool[Edge] seen) {
 820         Edge nex = oNext();
 821         while (nex != this) {
 822             if (nex in seen) {
 823                 seen.remove(nex);
 824                 stack.push(nex);
 825             }
 826             nex = nex.oNext();
 827         }
 828     }
 829 }  // Edge ends
 830
 831
 832 /**
 833 * A class that represents an edge pair
 834 **/
 835 final class EdgePair {
 836     Edge left;
 837     Edge right;
 838
 839     this(Edge l, Edge r) {
 840         left = l;
 841         right = r;
 842     }
 843
 844     public Edge getLeft() {
 845         return left;
 846     }
 847
 848     public Edge getRight() {
 849         return right;
 850     }
 851 } // EdgePair ends
 852
 853
 854 struct Voronoi {
 855     /**
 856     * The number of points in the diagram
 857     **/
 858     private static int points = 0;
 859
 860     /**
 861     * Set to true to print informative messages
 862     **/
 863     private static bool printMsgs = false;
 864
 865     /**
 866     * Set to true to print the voronoi diagram and its dual,
 867     * the delaunay diagram
 868     **/
 869     private static bool printResults = false;
 870
 871     /**
 872     * The main routine which creates the points and then performs
 873     * the delaunay triagulation.
 874     * @param args the command line parameters
 875     **/
 876     public static void main(char[][] args) {
 877         parseCmdLine(args);
 878
 879         if (printMsgs)
 880             printf("Getting %d points\n", points);
 881
 882         auto start0 = myclock();
 883         Vertex.seed = 1023;
 884         Vertex extra = Vertex.createPoints(1, new MyDouble(1.0), points);
 885         Vertex point = Vertex.createPoints(points-1, new MyDouble(extra.X()),
 886         points-1);
 887         auto end0 = myclock();
 888
 889         if (printMsgs)
 890             printf("Doing voronoi on %d nodes\n", points);
 891
 892         auto start1 = myclock();
 893         Edge edge = point.buildDelaunayTriangulation(extra);
 894         auto end1 = myclock();
 895
 896         if (printResults)
 897             edge.outputVoronoiDiagram();
 898
 899         if (printMsgs) {
 900             printf("Build time %f\n", end0 - start0);
 901             printf("Compute  time %f\n", end1 - start1);
 902             printf("Total time %f\n", end1 - start0);
 903             printf("Done!\n");
 904         }
 905     }
 906
 907     /**
 908     * Parse the command line options.
 909     * @param args the command line options.
 910     **/
 911     private static final void parseCmdLine(char[][] args) {
 912         int i = 1;
 913
 914         while (i < args.length && args[i][0] == '-') {
 915             char[] arg = args[i++];
 916
 917             if (arg == "-n") {
 918                 if (i < args.length)
 919                     points = toInt(args[i++]);
 920                 else
 921                     throw new Exception("-n requires the number of points");
 922             } else if (arg == "-p") {
 923                 printResults = true;
 924             } else if (arg == "-v") {
 925                 printMsgs = true;
 926             } else if (arg == "-h") {
 927                 usage();
 928             }
 929         }
 930
 931         if (points == 0)
 932             usage();
 933     }
 934
 935     /**
 936     * The usage routine which describes the program options.
 937     **/
 938     private static final void usage() {
 939         fprintf(stderr, "usage: voronoi_d -n <points> [-p] [-m] [-h]\n");
 940         fprintf(stderr, "    -n the number of points in the diagram\n");
 941         fprintf(stderr, "    -p (print detailed results/messages - the voronoi diagram>)\n");
 942         fprintf(stderr, "    -v (print informative message)\n");
 943         fprintf(stderr, "    -h (this message)\n");
 944         exit(0);
 945     }
 946 } // Voronoi ends
 947
 948
 949 void main(char[][] args) {
 950     Voronoi.main(args);
 951 }