micro/voronoi.d

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