nominatim/voronoi/VoronoiDiagramGenerator.cpp

   1 /*\r
   2 * The author of this software is Steven Fortune.  Copyright (c) 1994 by AT&T\r
   3 * Bell Laboratories.\r
   4 * Permission to use, copy, modify, and distribute this software for any\r
   5 * purpose without fee is hereby granted, provided that this entire notice\r
   6 * is included in all copies of any software which is or includes a copy\r
   7 * or modification of this software and in all copies of the supporting\r
   8 * documentation for such software.\r
   9 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED\r
  10 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY\r
  11 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY\r
  12 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.\r
  13 */\r
  14 \r
  15 /* \r
  16 * This code was originally written by Stephan Fortune in C code.  I, Shane O'Sullivan, \r
  17 * have since modified it, encapsulating it in a C++ class and, fixing memory leaks and \r
  18 * adding accessors to the Voronoi Edges.\r
  19 * Permission to use, copy, modify, and distribute this software for any\r
  20 * purpose without fee is hereby granted, provided that this entire notice\r
  21 * is included in all copies of any software which is or includes a copy\r
  22 * or modification of this software and in all copies of the supporting\r
  23 * documentation for such software.\r
  24 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED\r
  25 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY\r
  26 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY\r
  27 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.\r
  28 */\r
  29 \r
  30 #include "VoronoiDiagramGenerator.h"\r
  31 #include <stdio.h>\r
  32 #include <sys/mman.h>\r
  33 \r
  34 VoronoiDiagramGenerator::VoronoiDiagramGenerator()\r
  35 {\r
  36         siteidx = 0;\r
  37         sites = 0;\r
  38 \r
  39         allMemoryList = new FreeNodeArrayList;\r
  40         allMemoryList->memory = 0;\r
  41         allMemoryList->next = 0;\r
  42         currentMemoryBlock = allMemoryList;\r
  43         allEdges = 0;\r
  44         iteratorEdges = 0;\r
  45         minDistanceBetweenSites = 0;\r
  46 }\r
  47 \r
  48 VoronoiDiagramGenerator::~VoronoiDiagramGenerator()\r
  49 {\r
  50         cleanup();\r
  51         cleanupEdges();\r
  52 \r
  53         if(allMemoryList != 0)\r
  54                 delete allMemoryList;\r
  55 }\r
  56 \r
  57 bool VoronoiDiagramGenerator::generateVoronoi(struct SourcePoint* srcPoints, int numPoints, float minX, float maxX, float minY, float maxY, float minDist)\r
  58 {\r
  59         cleanup();\r
  60         cleanupEdges();\r
  61         int i;\r
  62 \r
  63         minDistanceBetweenSites = minDist;\r
  64 \r
  65         nsites = numPoints;\r
  66         plot = 0;\r
  67         triangulate = 0;\r
  68         debug = 1;\r
  69         sorted = 0;\r
  70         freeinit(&sfl, sizeof (Site));\r
  71 \r
  72         sites = (struct Site *) myalloc(nsites * sizeof(*sites));\r
  73         polygons = (struct Polygon *) myalloc(nsites * sizeof(*polygons));\r
  74 \r
  75         if(sites == 0) return false;\r
  76 \r
  77         xmin = srcPoints[0].x;\r
  78         ymin = srcPoints[0].y;\r
  79         xmax = srcPoints[0].x;\r
  80         ymax = srcPoints[0].y;\r
  81 \r
  82         for(i = 0; i < nsites; i++)\r
  83         {\r
  84                 sites[i].coord.x = srcPoints[i].x;\r
  85                 sites[i].coord.y = srcPoints[i].y;\r
  86                 sites[i].weight = srcPoints[i].weight;\r
  87                 sites[i].sitenbr = i;\r
  88                 sites[i].refcnt = 0; // prevent reuse?\r
  89 \r
  90                 if(sites[i].coord.x < xmin)\r
  91                         xmin = sites[i].coord.x;\r
  92                 else if(sites[i].coord.x > xmax)\r
  93                         xmax = sites[i].coord.x;\r
  94 \r
  95                 if(sites[i].coord.y < ymin)\r
  96                         ymin = sites[i].coord.y;\r
  97                 else if(sites[i].coord.y > ymax)\r
  98                         ymax = sites[i].coord.y;\r
  99 \r
 100                 polygons[i].coord.x = sites[i].coord.x;\r
 101                 polygons[i].coord.y = sites[i].coord.y;\r
 102                 polygons[i].numpoints = 0;\r
 103                 polygons[i].pointlist = NULL;\r
 104                 polygons[i].boundary = 0;\r
 105 \r
 106                 //printf("\n%lf %lf\n", sites[i].coord.x, sites[i].coord.y);\r
 107         }\r
 108 \r
 109         qsort(sites, nsites, sizeof (*sites), scomp);\r
 110 \r
 111         siteidx = 0;\r
 112         geominit();\r
 113         float temp = 0;\r
 114         if(minX > maxX)\r
 115         {\r
 116                 temp = minX;\r
 117                 minX = maxX;\r
 118                 maxX = temp;\r
 119         }\r
 120         if(minY > maxY)\r
 121         {\r
 122                 temp = minY;\r
 123                 minY = maxY;\r
 124                 maxY = temp;\r
 125         }\r
 126         borderMinX = minX;\r
 127         borderMinY = minY;\r
 128         borderMaxX = maxX;\r
 129         borderMaxY = maxY;\r
 130 \r
 131         corners[0].x = borderMinX;\r
 132         corners[0].y = borderMinY;\r
 133         corners[1].x = borderMinX;\r
 134         corners[1].y = borderMaxY;\r
 135         corners[2].x = borderMaxX;\r
 136         corners[2].y = borderMaxY;\r
 137         corners[3].x = borderMaxX;\r
 138         corners[3].y = borderMinY;\r
 139 \r
 140         siteidx = 0;\r
 141         voronoi(triangulate);\r
 142 \r
 143         return true;\r
 144 }\r
 145 \r
 146 bool VoronoiDiagramGenerator::ELinitialize()\r
 147 {\r
 148         int i;\r
 149         freeinit(&hfl, sizeof **ELhash);\r
 150         ELhashsize = 2 * sqrt_nsites;\r
 151         ELhash = (struct Halfedge **) myalloc ( sizeof *ELhash * ELhashsize);\r
 152 \r
 153         if(ELhash == 0)\r
 154                 return false;\r
 155 \r
 156         for(i=0; i<ELhashsize; i +=1) ELhash[i] = (struct Halfedge *)NULL;\r
 157         ELleftend = HEcreate( (struct Edge *)NULL, 0);\r
 158         ELrightend = HEcreate( (struct Edge *)NULL, 0);\r
 159         ELleftend -> ELleft = (struct Halfedge *)NULL;\r
 160         ELleftend -> ELright = ELrightend;\r
 161         ELrightend -> ELleft = ELleftend;\r
 162         ELrightend -> ELright = (struct Halfedge *)NULL;\r
 163         ELhash[0] = ELleftend;\r
 164         ELhash[ELhashsize-1] = ELrightend;\r
 165 \r
 166         return true;\r
 167 }\r
 168 \r
 169 \r
 170 struct Halfedge* VoronoiDiagramGenerator::HEcreate(struct Edge *e,int pm)\r
 171 {\r
 172         struct Halfedge *answer;\r
 173         answer = (struct Halfedge *) getfree(&hfl);\r
 174         answer -> ELedge = e;\r
 175         answer -> ELpm = pm;\r
 176         answer -> PQnext = (struct Halfedge *) NULL;\r
 177         answer -> vertex = (struct Site *) NULL;\r
 178         answer -> ELrefcnt = 0;\r
 179         return(answer);\r
 180 }\r
 181 \r
 182 \r
 183 void VoronoiDiagramGenerator::ELinsert(struct   Halfedge *lb, struct Halfedge *newHe)\r
 184 {\r
 185         newHe -> ELleft = lb;\r
 186         newHe -> ELright = lb -> ELright;\r
 187         (lb -> ELright) -> ELleft = newHe;\r
 188         lb -> ELright = newHe;\r
 189 }\r
 190 \r
 191 /* Get entry from hash table, pruning any deleted nodes */\r
 192 struct Halfedge * VoronoiDiagramGenerator::ELgethash(int b)\r
 193 {\r
 194         struct Halfedge *he;\r
 195 \r
 196         if(b<0 || b>=ELhashsize)\r
 197                 return((struct Halfedge *) NULL);\r
 198         he = ELhash[b];\r
 199         if (he == (struct Halfedge *) NULL || he->ELedge != (struct Edge *) DELETED )\r
 200                 return (he);\r
 201 \r
 202         /* Hash table points to deleted half edge.  Patch as necessary. */\r
 203         ELhash[b] = (struct Halfedge *) NULL;\r
 204         if ((he -> ELrefcnt -= 1) == 0) \r
 205                 makefree((Freenode*)he, &hfl);\r
 206         return ((struct Halfedge *) NULL);\r
 207 }\r
 208 \r
 209 struct Halfedge * VoronoiDiagramGenerator::ELleftbnd(struct Point *p)\r
 210 {\r
 211         int i, bucket;\r
 212         struct Halfedge *he;\r
 213 \r
 214         /* Use hash table to get close to desired halfedge */\r
 215         bucket = (int)((p->x - xmin)/deltax * ELhashsize);      //use the hash function to find the place in the hash map that this HalfEdge should be\r
 216 \r
 217         if(bucket<0) bucket =0;                                 //make sure that the bucket position in within the range of the hash array\r
 218         if(bucket>=ELhashsize) bucket = ELhashsize - 1;\r
 219 \r
 220         he = ELgethash(bucket);\r
 221         if(he == (struct Halfedge *) NULL)                      //if the HE isn't found, search backwards and forwards in the hash map for the first non-null entry\r
 222         {\r
 223                 for(i=1; 1 ; i += 1)\r
 224                 {\r
 225                         if ((he=ELgethash(bucket-i)) != (struct Halfedge *) NULL)\r
 226                                 break;\r
 227                         if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL)\r
 228                                 break;\r
 229                 };\r
 230                 totalsearch += i;\r
 231         };\r
 232         ntry += 1;\r
 233         /* Now search linear list of halfedges for the correct one */\r
 234         if (he==ELleftend  || (he != ELrightend && right_of(he,p)))\r
 235         {\r
 236                 do \r
 237                 {\r
 238                         he = he -> ELright;\r
 239                 } while (he!=ELrightend && right_of(he,p));     //keep going right on the list until either the end is reached, or you find the 1st edge which the point\r
 240                 he = he -> ELleft;                              //isn't to the right of\r
 241         }\r
 242         else                                                    //if the point is to the left of the HalfEdge, then search left for the HE just to the left of the point\r
 243                 do\r
 244                 {\r
 245                         he = he -> ELleft;\r
 246                 } while (he!=ELleftend && !right_of(he,p));\r
 247 \r
 248         /* Update hash table and reference counts */\r
 249         if(bucket > 0 && bucket <ELhashsize-1)\r
 250         {\r
 251                 if(ELhash[bucket] != (struct Halfedge *) NULL) \r
 252                 {\r
 253                         ELhash[bucket] -> ELrefcnt -= 1;\r
 254                 }\r
 255                 ELhash[bucket] = he;\r
 256                 ELhash[bucket] -> ELrefcnt += 1;\r
 257         };\r
 258         return (he);\r
 259 }\r
 260 \r
 261 \r
 262 /* This delete routine can't reclaim node, since pointers from hash\r
 263 table may be present.   */\r
 264 void VoronoiDiagramGenerator::ELdelete(struct Halfedge *he)\r
 265 {\r
 266         (he -> ELleft) -> ELright = he -> ELright;\r
 267         (he -> ELright) -> ELleft = he -> ELleft;\r
 268         he -> ELedge = (struct Edge *)DELETED;\r
 269 }\r
 270 \r
 271 \r
 272 struct Halfedge * VoronoiDiagramGenerator::ELright(struct Halfedge *he)\r
 273 {\r
 274         return (he -> ELright);\r
 275 }\r
 276 \r
 277 struct Halfedge * VoronoiDiagramGenerator::ELleft(struct Halfedge *he)\r
 278 {\r
 279         return (he -> ELleft);\r
 280 }\r
 281 \r
 282 \r
 283 struct Site * VoronoiDiagramGenerator::leftreg(struct Halfedge *he)\r
 284 {\r
 285         if(he -> ELedge == (struct Edge *)NULL)\r
 286                 return(bottomsite);\r
 287         return( he -> ELpm == le ?\r
 288                 he -> ELedge -> reg[le] : he -> ELedge -> reg[re]);\r
 289 }\r
 290 \r
 291 struct Site * VoronoiDiagramGenerator::rightreg(struct Halfedge *he)\r
 292 {\r
 293         if(he -> ELedge == (struct Edge *)NULL) //if this halfedge has no edge, return the bottom site (whatever that is)\r
 294                 return(bottomsite);\r
 295 \r
 296         //if the ELpm field is zero, return the site 0 that this edge bisects, otherwise return site number 1\r
 297         return( he -> ELpm == le ? he -> ELedge -> reg[re] : he -> ELedge -> reg[le]);\r
 298 }\r
 299 \r
 300 void VoronoiDiagramGenerator::geominit()\r
 301 {\r
 302         float sn;\r
 303 \r
 304         freeinit(&efl, sizeof(Edge));\r
 305         nvertices = 0;\r
 306         nedges = 0;\r
 307         sn = (float)nsites+4;\r
 308         sqrt_nsites = (int)sqrt(sn);\r
 309         deltay = ymax - ymin;\r
 310         deltax = xmax - xmin;\r
 311 }\r
 312 \r
 313 \r
 314 struct Edge * VoronoiDiagramGenerator::bisect(struct Site *s1,struct    Site *s2)\r
 315 {\r
 316         float dx,dy,adx,ady;\r
 317         struct Edge *newedge;\r
 318 \r
 319         newedge = (struct Edge *) getfree(&efl);\r
 320 \r
 321         newedge -> reg[0] = s1; //store the sites that this edge is bisecting\r
 322         newedge -> reg[1] = s2;\r
 323         ref(s1);\r
 324         ref(s2);\r
 325         newedge -> ep[0] = (struct Site *) NULL; //to begin with, there are no endpoints on the bisector - it goes to infinity\r
 326         newedge -> ep[1] = (struct Site *) NULL;\r
 327 \r
 328         dx = s2->coord.x - s1->coord.x;                 //get the difference in x dist between the sites\r
 329         dy = s2->coord.y - s1->coord.y;\r
 330         adx = dx>0 ? dx : -dx;                                  //make sure that the difference in positive\r
 331         ady = dy>0 ? dy : -dy;\r
 332         newedge -> c = (float)(s1->coord.x * dx + s1->coord.y * dy + (dx*dx + dy*dy)*0.5);//get the slope of the line\r
 333 \r
 334         if (adx>ady)\r
 335         {\r
 336                 newedge -> a = 1.0; newedge -> b = dy/dx; newedge -> c /= dx;//set formula of line, with x fixed to 1\r
 337         }\r
 338         else\r
 339         {\r
 340                 newedge -> b = 1.0; newedge -> a = dx/dy; newedge -> c /= dy;//set formula of line, with y fixed to 1\r
 341         };\r
 342 \r
 343         newedge -> edgenbr = nedges;\r
 344 \r
 345         //printf("\nbisect(%d) ((%f,%f) and (%f,%f)",nedges,s1->coord.x,s1->coord.y,s2->coord.x,s2->coord.y);\r
 346 \r
 347         nedges += 1;\r
 348         return(newedge);\r
 349 }\r
 350 \r
 351 //create a new site where the HalfEdges el1 and el2 intersect - note that the Point in the argument list is not used, don't know why it's there\r
 352 struct Site * VoronoiDiagramGenerator::intersect(struct Halfedge *el1, struct Halfedge *el2, struct Point *p)\r
 353 {\r
 354         struct  Edge *e1,*e2, *e;\r
 355         struct  Halfedge *el;\r
 356         float d, xint, yint;\r
 357         int right_of_site;\r
 358         struct Site *v;\r
 359 \r
 360         e1 = el1 -> ELedge;\r
 361         e2 = el2 -> ELedge;\r
 362         if(e1 == (struct Edge*)NULL || e2 == (struct Edge*)NULL)\r
 363                 return ((struct Site *) NULL);\r
 364 \r
 365         //if the two edges bisect the same parent, return null\r
 366         if (e1->reg[1] == e2->reg[1])\r
 367                 return ((struct Site *) NULL);\r
 368 \r
 369         d = e1->a * e2->b - e1->b * e2->a;\r
 370         if (-1.0e-10<d && d<1.0e-10)\r
 371                 return ((struct Site *) NULL);\r
 372 \r
 373         xint = (e1->c*e2->b - e2->c*e1->b)/d;\r
 374         yint = (e2->c*e1->a - e1->c*e2->a)/d;\r
 375 \r
 376         if( (e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||\r
 377                 (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&\r
 378                 e1->reg[1]->coord.x < e2->reg[1]->coord.x) )\r
 379         {\r
 380                 el = el1;\r
 381                 e = e1;\r
 382         }\r
 383         else\r
 384         {\r
 385                 el = el2;\r
 386                 e = e2;\r
 387         };\r
 388 \r
 389         right_of_site = xint >= e -> reg[1] -> coord.x;\r
 390         if ((right_of_site && el -> ELpm == le) || (!right_of_site && el -> ELpm == re))\r
 391                 return ((struct Site *) NULL);\r
 392 \r
 393         //create a new site at the point of intersection - this is a new vector event waiting to happen\r
 394         v = (struct Site *) getfree(&sfl);\r
 395         v -> refcnt = 0;\r
 396         v -> coord.x = xint;\r
 397         v -> coord.y = yint;\r
 398         return(v);\r
 399 }\r
 400 \r
 401 /* returns 1 if p is to right of halfedge e */\r
 402 int VoronoiDiagramGenerator::right_of(struct Halfedge *el,struct Point *p)\r
 403 {\r
 404         struct Edge *e;\r
 405         struct Site *topsite;\r
 406         int right_of_site, above, fast;\r
 407         float dxp, dyp, dxs, t1, t2, t3, yl;\r
 408 \r
 409         e = el -> ELedge;\r
 410         topsite = e -> reg[1];\r
 411         right_of_site = p -> x > topsite -> coord.x;\r
 412         if(right_of_site && el -> ELpm == le) return(1);\r
 413         if(!right_of_site && el -> ELpm == re) return (0);\r
 414         if (e->a == 1.0)\r
 415         {\r
 416                 dyp = p->y - topsite->coord.y;\r
 417                 dxp = p->x - topsite->coord.x;\r
 418                 fast = 0;\r
 419                 if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )\r
 420                 {\r
 421                         above = dyp>= e->b*dxp;\r
 422                         fast = above;\r
 423                 }\r
 424                 else\r
 425                 {\r
 426                         above = p->x + p->y*e->b > e-> c;\r
 427                         if(e->b<0.0) above = !above;\r
 428                         if (!above) fast = 1;\r
 429                 }\r
 430 \r
 431                 if (!fast)\r
 432                 {\r
 433                         dxs = topsite->coord.x - (e->reg[0])->coord.x;\r
 434                         above = e->b * (dxp*dxp - dyp*dyp) <\r
 435                         dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);\r
 436                         if(e->b<0.0) above = !above;\r
 437                 }\r
 438         }\r
 439         else  /*e->b==1.0 */\r
 440         {\r
 441                 yl = e->c - e->a*p->x;\r
 442                 t1 = p->y - yl;\r
 443                 t2 = p->x - topsite->coord.x;\r
 444                 t3 = yl - topsite->coord.y;\r
 445                 above = t1*t1 > t2*t2 + t3*t3;\r
 446         }\r
 447         return (el->ELpm==le ? above : !above);\r
 448 }\r
 449 \r
 450 \r
 451 void VoronoiDiagramGenerator::endpoint(struct Edge *e,int lr,struct Site * s)\r
 452 {\r
 453         e -> ep[lr] = s;\r
 454         ref(s);\r
 455         return;\r
 456 \r
 457         if(e -> ep[re-lr]== (struct Site *) NULL)\r
 458                 return;\r
 459 \r
 460         clip_line(e);\r
 461 \r
 462         deref(e->reg[le]);\r
 463         deref(e->reg[re]);\r
 464         makefree((Freenode*)e, &efl);\r
 465 }\r
 466 \r
 467 void VoronoiDiagramGenerator::endpoint(struct Edge *e1,int lr,struct Site * s, struct Edge *e2, struct Edge *e3)\r
 468 {\r
 469         e1 -> ep[lr] = s;\r
 470         ref(s);\r
 471 \r
 472         s->coordout.x = s->coord.x;\r
 473         s->coordout.y = s->coord.y;\r
 474 \r
 475         if(e1 -> ep[le] != (struct Site *) NULL && e1 -> ep[re] != (struct Site *) NULL)\r
 476         {\r
 477                 clip_line(e1);\r
 478                 deref(e1->reg[le]);\r
 479                 deref(e1->reg[re]);\r
 480                 makefree((Freenode*)e1, &efl);\r
 481         }\r
 482 \r
 483         if(e2 -> ep[le] != (struct Site *) NULL && e2 -> ep[re] != (struct Site *) NULL)\r
 484         {\r
 485                 clip_line(e2);\r
 486                 deref(e2->reg[le]);\r
 487                 deref(e2->reg[re]);\r
 488                 makefree((Freenode*)e2, &efl);\r
 489         }\r
 490 \r
 491         if(e3 -> ep[le] != (struct Site *) NULL && e3 -> ep[re] != (struct Site *) NULL)\r
 492         {\r
 493                 clip_line(e3);\r
 494                 deref(e3->reg[le]);\r
 495                 deref(e3->reg[re]);\r
 496                 makefree((Freenode*)e3, &efl);\r
 497         }\r
 498 \r
 499         return; \r
 500 }\r
 501 \r
 502 \r
 503 float VoronoiDiagramGenerator::dist(struct Site *s,struct Site *t)\r
 504 {\r
 505         float dx,dy;\r
 506         dx = s->coord.x - t->coord.x;\r
 507         dy = s->coord.y - t->coord.y;\r
 508         return (float)(sqrt(dx*dx + dy*dy));\r
 509 }\r
 510 \r
 511 \r
 512 void VoronoiDiagramGenerator::makevertex(struct Site *v)\r
 513 {\r
 514         v -> sitenbr = nvertices;\r
 515         nvertices += 1;\r
 516         out_vertex(v);\r
 517 }\r
 518 \r
 519 \r
 520 void VoronoiDiagramGenerator::deref(struct Site *v)\r
 521 {\r
 522         v -> refcnt -= 1;\r
 523         if (v -> refcnt == 0 ) \r
 524                 makefree((Freenode*)v, &sfl);\r
 525 }\r
 526 \r
 527 void VoronoiDiagramGenerator::ref(struct Site *v)\r
 528 {\r
 529         v -> refcnt += 1;\r
 530 }\r
 531 \r
 532 //push the HalfEdge into the ordered linked list of vertices\r
 533 void VoronoiDiagramGenerator::PQinsert(struct Halfedge *he,struct Site * v, float offset)\r
 534 {\r
 535         struct Halfedge *last, *next;\r
 536 \r
 537         he -> vertex = v;\r
 538         ref(v);\r
 539         he -> ystar = (float)(v -> coord.y + offset);\r
 540         last = &PQhash[PQbucket(he)];\r
 541         while ((next = last -> PQnext) != (struct Halfedge *) NULL &&\r
 542                 (he -> ystar  > next -> ystar  ||\r
 543                 (he -> ystar == next -> ystar && v -> coord.x > next->vertex->coord.x)))\r
 544         {\r
 545                 last = next;\r
 546         }\r
 547         he -> PQnext = last -> PQnext;\r
 548         last -> PQnext = he;\r
 549         PQcount += 1;\r
 550 }\r
 551 \r
 552 //remove the HalfEdge from the list of vertices\r
 553 void VoronoiDiagramGenerator::PQdelete(struct Halfedge *he)\r
 554 {\r
 555         struct Halfedge *last;\r
 556 \r
 557         if(he -> vertex != (struct Site *) NULL)\r
 558         {\r
 559                 last = &PQhash[PQbucket(he)];\r
 560                 while (last -> PQnext != he)\r
 561                         last = last -> PQnext;\r
 562 \r
 563                 last -> PQnext = he -> PQnext;\r
 564                 PQcount -= 1;\r
 565                 deref(he -> vertex);\r
 566                 he -> vertex = (struct Site *) NULL;\r
 567         };\r
 568 }\r
 569 \r
 570 int VoronoiDiagramGenerator::PQbucket(struct Halfedge *he)\r
 571 {\r
 572         int bucket;\r
 573 \r
 574         bucket = (int)((he->ystar - ymin)/deltay * PQhashsize);\r
 575         if (bucket<0) bucket = 0;\r
 576         if (bucket>=PQhashsize) bucket = PQhashsize-1 ;\r
 577         if (bucket < PQmin) PQmin = bucket;\r
 578         return(bucket);\r
 579 }\r
 580 \r
 581 int VoronoiDiagramGenerator::PQempty()\r
 582 {\r
 583         return(PQcount==0);\r
 584 }\r
 585 \r
 586 \r
 587 struct Point VoronoiDiagramGenerator::PQ_min()\r
 588 {\r
 589         struct Point answer;\r
 590 \r
 591         while(PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {PQmin += 1;};\r
 592         answer.x = PQhash[PQmin].PQnext -> vertex -> coord.x;\r
 593         answer.y = PQhash[PQmin].PQnext -> ystar;\r
 594         return (answer);\r
 595 }\r
 596 \r
 597 struct Halfedge * VoronoiDiagramGenerator::PQextractmin()\r
 598 {\r
 599         struct Halfedge *curr;\r
 600 \r
 601         curr = PQhash[PQmin].PQnext;\r
 602         PQhash[PQmin].PQnext = curr -> PQnext;\r
 603         PQcount -= 1;\r
 604         return(curr);\r
 605 }\r
 606 \r
 607 \r
 608 bool VoronoiDiagramGenerator::PQinitialize()\r
 609 {\r
 610         int i;\r
 611 \r
 612         PQcount = 0;\r
 613         PQmin = 0;\r
 614         PQhashsize = 4 * sqrt_nsites;\r
 615         PQhash = (struct Halfedge *) myalloc(PQhashsize * sizeof *PQhash);\r
 616 \r
 617         if(PQhash == 0)\r
 618                 return false;\r
 619 \r
 620         for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (struct Halfedge *)NULL;\r
 621 \r
 622         return true;\r
 623 }\r
 624 \r
 625 \r
 626 void VoronoiDiagramGenerator::freeinit(struct Freelist *fl,int size)\r
 627 {\r
 628         fl -> head = (struct Freenode *) NULL;\r
 629         fl -> nodesize = size;\r
 630 }\r
 631 \r
 632 char * VoronoiDiagramGenerator::getfree(struct Freelist *fl)\r
 633 {\r
 634         int i;\r
 635         struct Freenode *t;\r
 636 \r
 637         if(fl->head == (struct Freenode *) NULL)\r
 638         {\r
 639                 t =  (struct Freenode *) myalloc(sqrt_nsites * fl->nodesize);\r
 640 \r
 641                 if(t == 0)\r
 642                         return 0;\r
 643 \r
 644                 currentMemoryBlock->next = new FreeNodeArrayList;\r
 645                 currentMemoryBlock = currentMemoryBlock->next;\r
 646                 currentMemoryBlock->memory = t;\r
 647                 currentMemoryBlock->next = 0;\r
 648 \r
 649                 for(i=0; i<sqrt_nsites; i+=1)\r
 650                         makefree((struct Freenode *)((char *)t+i*fl->nodesize), fl);\r
 651         };\r
 652         t = fl -> head;\r
 653         fl -> head = (fl -> head) -> nextfree;\r
 654         return((char *)t);\r
 655 }\r
 656 \r
 657 \r
 658 \r
 659 void VoronoiDiagramGenerator::makefree(struct Freenode *curr,struct Freelist *fl)\r
 660 {\r
 661         curr -> nextfree = fl -> head;\r
 662         fl -> head = curr;\r
 663 }\r
 664 \r
 665 void VoronoiDiagramGenerator::cleanup()\r
 666 {\r
 667         if(sites != 0)\r
 668         {\r
 669                 free(sites);\r
 670                 sites = 0;\r
 671         }\r
 672 \r
 673         FreeNodeArrayList* current=0, *prev = 0;\r
 674 \r
 675         current = prev = allMemoryList;\r
 676 \r
 677         while(current->next != 0)\r
 678         {\r
 679                 prev = current;\r
 680                 current = current->next;\r
 681                 free(prev->memory);\r
 682                 delete prev;\r
 683                 prev = 0;\r
 684         }\r
 685 \r
 686         if(current != 0 && current->memory != 0)\r
 687         {\r
 688                 free(current->memory);\r
 689                 delete current;\r
 690         }\r
 691 \r
 692         allMemoryList = new FreeNodeArrayList;\r
 693         allMemoryList->next = 0;\r
 694         allMemoryList->memory = 0;\r
 695         currentMemoryBlock = allMemoryList;\r
 696 }\r
 697 \r
 698 void VoronoiDiagramGenerator::cleanupEdges()\r
 699 {\r
 700         GraphEdge* geCurrent = 0, *gePrev = 0;\r
 701         geCurrent = gePrev = allEdges;\r
 702 \r
 703         while(geCurrent != 0 && geCurrent->next != 0)\r
 704         {\r
 705                 gePrev = geCurrent;\r
 706                 geCurrent = geCurrent->next;\r
 707                 delete gePrev;\r
 708         }\r
 709 \r
 710         allEdges = 0;\r
 711 \r
 712 }\r
 713 \r
 714 void VoronoiDiagramGenerator::pushGraphEdge(float x1, float y1, float x2, float y2)\r
 715 {\r
 716         GraphEdge* newEdge = new GraphEdge;\r
 717         newEdge->next = allEdges;\r
 718         allEdges = newEdge;\r
 719         newEdge->x1 = x1;\r
 720         newEdge->y1 = y1;\r
 721         newEdge->x2 = x2;\r
 722         newEdge->y2 = y2;\r
 723 }\r
 724 \r
 725 \r
 726 char * VoronoiDiagramGenerator::myalloc(unsigned n)\r
 727 {\r
 728         char *t=0;\r
 729         t=(char*)malloc(n);\r
 730         total_alloc += n;\r
 731         return(t);\r
 732 }\r
 733 \r
 734 \r
 735 /* for those who don't have Cherry's plot */\r
 736 /* #include <plot.h> */\r
 737 void VoronoiDiagramGenerator::openpl(){}\r
 738 void VoronoiDiagramGenerator::line(float x1, float y1, float x2, float y2)\r
 739 {\r
 740         pushGraphEdge(x1,y1,x2,y2);\r
 741 \r
 742 }\r
 743 void VoronoiDiagramGenerator::circle(float x, float y, float radius){}\r
 744 void VoronoiDiagramGenerator::range(float minX, float minY, float maxX, float maxY){}\r
 745 \r
 746 \r
 747 \r
 748 void VoronoiDiagramGenerator::out_bisector(struct Edge *e)\r
 749 {\r
 750 \r
 751 }\r
 752 \r
 753 \r
 754 void VoronoiDiagramGenerator::out_ep(struct Edge *e)\r
 755 {\r
 756 \r
 757 }\r
 758 \r
 759 void VoronoiDiagramGenerator::out_vertex(struct Site *v)\r
 760 {\r
 761 \r
 762 }\r
 763 \r
 764 \r
 765 void VoronoiDiagramGenerator::out_site(struct Site *s)\r
 766 {\r
 767         if(!triangulate & plot & !debug)\r
 768                 circle (s->coord.x, s->coord.y, cradius);\r
 769 }\r
 770 \r
 771 \r
 772 void VoronoiDiagramGenerator::out_triple(struct Site *s1, struct Site *s2,struct Site * s3)\r
 773 {\r
 774 \r
 775 }\r
 776 \r
 777 \r
 778 \r
 779 void VoronoiDiagramGenerator::plotinit()\r
 780 {\r
 781         float dx,dy,d;\r
 782 \r
 783         dy = ymax - ymin;\r
 784         dx = xmax - xmin;\r
 785         d = (float)(( dx > dy ? dx : dy) * 1.1);\r
 786         pxmin = (float)(xmin - (d-dx)/2.0);\r
 787         pxmax = (float)(xmax + (d-dx)/2.0);\r
 788         pymin = (float)(ymin - (d-dy)/2.0);\r
 789         pymax = (float)(ymax + (d-dy)/2.0);\r
 790         cradius = (float)((pxmax - pxmin)/350.0);\r
 791         openpl();\r
 792         range(pxmin, pymin, pxmax, pymax);\r
 793 }\r
 794 \r
 795 void VoronoiDiagramGenerator::pushpoint(int sitenbr, double x, double y, int boundary)\r
 796 {\r
 797         Polygon *s;\r
 798 \r
 799         s = &polygons[sitenbr];\r
 800 \r
 801         if (s->numpoints == 0)\r
 802         {\r
 803                 s->pointlist = (PolygonPoint *)malloc(sizeof(struct PolygonPoint)*(s->numpoints+10));\r
 804                 if (!s->pointlist)\r
 805                 {\r
 806                         printf("Out of mem\n");\r
 807                 }\r
 808         }\r
 809         else if (s->numpoints % 10 == 0)\r
 810         {\r
 811                 s->pointlist = (PolygonPoint *)realloc(s->pointlist, sizeof(struct PolygonPoint)*(s->numpoints+10));\r
 812                 if (!s->pointlist)\r
 813                 {\r
 814                         printf("Out of remem\n");\r
 815                 }\r
 816         }\r
 817         s->pointlist[s->numpoints].coord.x = x;\r
 818         s->pointlist[s->numpoints].coord.y = y;\r
 819         s->pointlist[s->numpoints].angle = atan2f(x-s->coord.x, y-s->coord.y);\r
 820         s->pointlist[s->numpoints].boundary = boundary;\r
 821 \r
 822         if (boundary) s->boundary = 1;\r
 823 \r
 824         //printf("point #%d in %d (%lf, %lf) [%d] (%lf, %lf): %lf\n", s->numpoints, sitenbr, s->coord.x, s->coord.y, boundary, x, y, (s->pointlist[s->numpoints].angle/M_PI)*180);\r
 825 \r
 826         s->numpoints++;\r
 827 }\r
 828 \r
 829 int VoronoiDiagramGenerator::ccw( Point p0, Point p1, Point p2 )\r
 830 {\r
 831         double dx1, dx2, dy1, dy2;\r
 832 \r
 833         dx1 = p1.x - p0.x; dy1 = p1.y - p0.y;\r
 834         dx2 = p2.x - p0.x; dy2 = p2.y - p0.y;\r
 835 \r
 836         if (dx1*dy2 > dy1*dx2)\r
 837                 return +1;\r
 838         if (dx1*dy2 < dy1*dx2)\r
 839                 return -1;\r
 840         if ((dx1*dx2 < 0) || (dy1*dy2 < 0))\r
 841                 return -1;\r
 842         if ((dx1*dx1 + dy1*dy1) < (dx2*dx2 + dy2*dy2))\r
 843                 return +1;\r
 844         return 0;\r
 845 }\r
 846 \r
 847 void VoronoiDiagramGenerator::getSitePoints(int sitenbr, int* numpoints, PolygonPoint** pS)\r
 848 {\r
 849         int i, j, c, any, centrevalue, cornerinpolygon[4];\r
 850 \r
 851         if (polygons[sitenbr].numpoints == 0)\r
 852         {\r
 853                 for(c = 0; c < 4; c++)\r
 854                 {\r
 855                         pushpoint(sitenbr, corners[c].x, corners[c].y, 0);\r
 856                 }\r
 857         }\r
 858 \r
 859         qsort(polygons[sitenbr].pointlist, polygons[sitenbr].numpoints, sizeof(PolygonPoint), anglecomp);\r
 860 \r
 861         if (polygons[sitenbr].boundary)\r
 862         {\r
 863 //              printf("\nsite %d is boundary intersection\n", sitenbr);\r
 864 \r
 865                 for(c = 0; c < 4; c++) cornerinpolygon[c] = 1;\r
 866 \r
 867                 for(i = 0; i < polygons[sitenbr].numpoints; i++)\r
 868                 {\r
 869 //                      printf("Point (%lf,%lf) %d\n", polygons[sitenbr].pointlist[i].coord.x,polygons[sitenbr].pointlist[i].coord.y,polygons[sitenbr].pointlist[i].boundary);\r
 870                         j = i > 0?i-1:polygons[sitenbr].numpoints-1;\r
 871                         if (    (!polygons[sitenbr].pointlist[i].boundary || !polygons[sitenbr].pointlist[j].boundary) &&\r
 872                                 (polygons[sitenbr].pointlist[i].coord.x != polygons[sitenbr].pointlist[j].coord.x ||\r
 873                                 polygons[sitenbr].pointlist[i].coord.y != polygons[sitenbr].pointlist[j].coord.y))\r
 874                         {\r
 875 //                              printf("line side test (%lf,%lf) => (%lf,%lf)\n",polygons[sitenbr].pointlist[i].coord.x,polygons[sitenbr].pointlist[i].coord.y,polygons[sitenbr].pointlist[j].coord.x,polygons[sitenbr].pointlist[j].coord.y);\r
 876                                 any = 0;\r
 877                                 centrevalue = ccw(polygons[sitenbr].pointlist[i].coord, polygons[sitenbr].pointlist[j].coord, polygons[sitenbr].coord);\r
 878 //printf(" test against centre (%lf,%lf) %d\n", polygons[sitenbr].coord.x, polygons[sitenbr].coord.y, centrevalue);\r
 879                                 for(c = 0; c < 4; c++)\r
 880                                 {\r
 881                                         if (cornerinpolygon[c])\r
 882                                         {\r
 883 \r
 884 //printf(" test against corner (%lf,%lf) %d\n", corners[c].x, corners[c].y, ccw(polygons[sitenbr].pointlist[i].coord, polygons[sitenbr].pointlist[j].coord, corners[c]));\r
 885 \r
 886                                                 if (centrevalue == ccw(polygons[sitenbr].pointlist[i].coord, polygons[sitenbr].pointlist[j].coord, corners[c]))\r
 887                                                 {\r
 888                                                         any = 1;\r
 889                                                 }\r
 890                                                 else\r
 891                                                 {\r
 892                                                         cornerinpolygon[c] = 0;\r
 893                                                 }\r
 894                                         }\r
 895                                 }\r
 896                                 if (!any) break;\r
 897                         }\r
 898                 }\r
 899                 if (any)\r
 900                 {\r
 901                         for(c = 0; c < 4; c++)\r
 902                         {\r
 903                                 if (cornerinpolygon[c])\r
 904                                 {\r
 905 //                                      printf("adding corger (%lf,%lf) to %d\n", corners[c].x, corners[c].y, sitenbr);\r
 906                                         pushpoint(sitenbr, corners[c].x, corners[c].y, 0);\r
 907                                 }\r
 908                         }\r
 909                 }\r
 910                 qsort(polygons[sitenbr].pointlist, polygons[sitenbr].numpoints, sizeof(PolygonPoint), anglecomp);\r
 911 \r
 912                 polygons[sitenbr].boundary = 0;\r
 913         }\r
 914 \r
 915         *numpoints = polygons[sitenbr].numpoints;\r
 916         *pS = polygons[sitenbr].pointlist;\r
 917 }\r
 918 \r
 919 \r
 920 void VoronoiDiagramGenerator::clip_line(struct Edge *e)\r
 921 {\r
 922         struct Site *s1, *s2, *ts1, *ts2;\r
 923         float x1=0,x2=0,y1=0,y2=0, temp = 0;\r
 924         int boundary1 = 0, boundary2 = 0;\r
 925 \r
 926 \r
 927         x1 = e->reg[0]->coord.x;\r
 928         x2 = e->reg[1]->coord.x;\r
 929         y1 = e->reg[0]->coord.y;\r
 930         y2 = e->reg[1]->coord.y;\r
 931 \r
 932         if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) == 0)\r
 933         {\r
 934                 return;\r
 935         }\r
 936 \r
 937         pxmin = borderMinX;\r
 938         pxmax = borderMaxX;\r
 939         pymin = borderMinY;\r
 940         pymax = borderMaxY;\r
 941 \r
 942         if(e -> a == 1.0 && e ->b >= 0.0)\r
 943         {\r
 944                 s1 = e -> ep[1];\r
 945                 s2 = e -> ep[0];\r
 946 \r
 947                 ts1 = e -> reg[1];\r
 948                 ts2 = e -> reg[0];\r
 949         }\r
 950         else\r
 951         {\r
 952                 s1 = e -> ep[0];\r
 953                 s2 = e -> ep[1];\r
 954 \r
 955                 ts1 = e -> reg[0];\r
 956                 ts2 = e -> reg[1];\r
 957 \r
 958         };\r
 959 \r
 960         if(e -> a == 1.0)\r
 961         {\r
 962                 if (    s1!=(struct Site *)NULL\r
 963                         && s1->coordout.y > pymin && s1->coordout.y < pymax\r
 964                         && s1->coordout.x > pxmin && s1->coordout.x < pxmax)\r
 965                 {\r
 966                         x1 = s1->coordout.x;\r
 967                         y1 = s1->coordout.y;\r
 968                 }\r
 969                 else\r
 970                 {\r
 971                         boundary1 = 1;\r
 972                         y1 = pymin;\r
 973                         if (s1!=(struct Site *)NULL && s1->coord.y > pymin)\r
 974                         {\r
 975                                 y1 = s1->coord.y;\r
 976                         }\r
 977                         if(y1>pymax)\r
 978                         {\r
 979                                 y1 = pymax;\r
 980                         }\r
 981                         x1 = e -> c - e -> b * y1;\r
 982                 }\r
 983 \r
 984                 if (    s2!=(struct Site *)NULL\r
 985                         && s2->coordout.y > pymin && s2->coordout.y < pymax\r
 986                         && s2->coordout.x > pxmin && s2->coordout.x < pxmax)\r
 987                 {\r
 988                         x2 = s2->coordout.x;\r
 989                         y2 = s2->coordout.y;\r
 990                 }\r
 991                 else\r
 992                 {\r
 993                         boundary2 = 1;\r
 994                         y2 = pymax;\r
 995                         if (s2!=(struct Site *)NULL && s2->coord.y < pymax)\r
 996                                 y2 = s2->coord.y;\r
 997                         if(y2<pymin)\r
 998                         {\r
 999                                 y2 = pymin;\r
1000                         }\r
1001                         x2 = (e->c) - (e->b) * y2;\r
1002                 }\r
1003 \r
1004                 if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin)))\r
1005                 {\r
1006                         // Line completely outside clipbox\r
1007                         //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);\r
1008                         return;\r
1009                 }\r
1010                 if(x1 > pxmax)\r
1011                 {\r
1012                         x1 = pxmax;\r
1013                         y1 = (e -> c - x1)/e -> b;\r
1014                 }\r
1015                 if(x1 < pxmin)\r
1016                 {\r
1017                         x1 = pxmin;\r
1018                         y1 = (e -> c - x1)/e -> b;\r
1019                 }\r
1020                 if(x2 > pxmax)\r
1021                 {\r
1022                         x2 = pxmax;\r
1023                         y2 = (e -> c - x2)/e -> b;\r
1024                 }\r
1025                 if(x2 < pxmin)\r
1026                 {\r
1027                         x2 = pxmin;\r
1028                         y2 = (e -> c - x2)/e -> b;\r
1029                 }\r
1030         }\r
1031         else\r
1032         {\r
1033                 if (    s1!=(struct Site *)NULL\r
1034                         && s1->coordout.y > pymin && s1->coordout.y < pymax\r
1035                         && s1->coordout.x > pxmin && s1->coordout.x < pxmax)\r
1036                 {\r
1037                         x1 = s1->coordout.x;\r
1038                         y1 = s1->coordout.y;\r
1039                 }\r
1040                 else\r
1041                 {\r
1042                         boundary1 = 1;\r
1043                         x1 = pxmin;\r
1044                         if (s1!=(struct Site *)NULL && s1->coord.x > pxmin)\r
1045                                 x1 = s1->coord.x;\r
1046                         if(x1>pxmax) \r
1047                         {\r
1048                                 //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);\r
1049                                 //return;\r
1050                                 x1 = pxmax;\r
1051                         }\r
1052                         y1 = e -> c - e -> a * x1;\r
1053                 }\r
1054 \r
1055                 if (    s2!=(struct Site *)NULL\r
1056                         && s2->coordout.y > pymin && s2->coordout.y < pymax\r
1057                         && s2->coordout.x > pxmin && s2->coordout.x < pxmax)\r
1058                 {\r
1059                         x2 = s2->coordout.x;\r
1060                         y2 = s2->coordout.y;\r
1061                 }\r
1062                 else\r
1063                 {\r
1064                         boundary2 = 1;\r
1065                         x2 = pxmax;\r
1066                         if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)\r
1067                                 x2 = s2->coord.x;\r
1068                         if(x2<pxmin)\r
1069                         {\r
1070                                 //printf("\nClipped (4) x2 = %f to %f",x2,pxmin);\r
1071                                 //return;\r
1072                                 x2 = pxmin;\r
1073                         }\r
1074                         y2 = e -> c - e -> a * x2;\r
1075                 }\r
1076 \r
1077                 if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin)))\r
1078                 {\r
1079                         //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);\r
1080                         return;\r
1081                 }\r
1082                 if(y1 > pymax)\r
1083                 {       y1 = pymax; x1 = (e -> c - y1)/e -> a;};\r
1084                 if(y1 < pymin)\r
1085                 {       y1 = pymin; x1 = (e -> c - y1)/e -> a;};\r
1086                 if(y2 > pymax)\r
1087                 {       y2 = pymax; x2 = (e -> c - y2)/e -> a;};\r
1088                 if(y2 < pymin)\r
1089                 {       y2 = pymin; x2 = (e -> c - y2)/e -> a;};\r
1090         };\r
1091 \r
1092         if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) == 0)\r
1093         {\r
1094                 return;\r
1095         }\r
1096 \r
1097         pushpoint(ts1->sitenbr, x1, y1, boundary1);\r
1098         pushpoint(ts2->sitenbr, x2, y2, boundary2);\r
1099         pushpoint(ts1->sitenbr, x2, y2, boundary2);\r
1100         pushpoint(ts2->sitenbr, x1, y1, boundary1);\r
1101 }\r
1102 \r
1103 \r
1104 /* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,\r
1105 deltax, deltay (can all be estimates).\r
1106 Performance suffers if they are wrong; better to make nsites,\r
1107 deltax, and deltay too big than too small.  (?) */\r
1108 \r
1109 bool VoronoiDiagramGenerator::voronoi(int triangulate)\r
1110 {\r
1111         struct Site *newsite, *bot, *top, *temp, *p;\r
1112         struct Site *v;\r
1113         struct Point newintstar;\r
1114         int pm;\r
1115         struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;\r
1116         struct Edge *e, *e2, *e3;\r
1117 \r
1118         PQinitialize();\r
1119         bottomsite = nextone();\r
1120         out_site(bottomsite);\r
1121         bool retval = ELinitialize();\r
1122 \r
1123         if(!retval)\r
1124                 return false;\r
1125 \r
1126         newsite = nextone();\r
1127         while(1)\r
1128         {\r
1129 \r
1130                 if(!PQempty())\r
1131                         newintstar = PQ_min();\r
1132 \r
1133                 //if the lowest site has a smaller y value than the lowest vector intersection, process the site\r
1134                 //otherwise process the vector intersection             \r
1135 \r
1136                 if (newsite != (struct Site *)NULL      && (PQempty() || newsite -> coord.y < newintstar.y\r
1137                         || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))\r
1138                 {/* new site is smallest - this is a site event*/\r
1139                         out_site(newsite);                                              //output the site\r
1140                         lbnd = ELleftbnd(&(newsite->coord));                            //get the first HalfEdge to the LEFT of the new site\r
1141                         rbnd = ELright(lbnd);                                           //get the first HalfEdge to the RIGHT of the new site\r
1142                         bot = rightreg(lbnd);                                           //if this halfedge has no edge, , bot = bottom site (whatever that is)\r
1143                         e = bisect(bot, newsite);                                       //create a new edge that bisects \r
1144                         bisector = HEcreate(e, le);                                     //create a new HalfEdge, setting its ELpm field to 0                    \r
1145                         ELinsert(lbnd, bisector);                                       //insert this new bisector edge between the left and right vectors in a linked list     \r
1146 \r
1147                         if ((p = intersect(lbnd, bisector)) != (struct Site *) NULL)    //if the new bisector intersects with the left edge, remove the left edge's vertex, and put in the new one\r
1148                         {       \r
1149                                 PQdelete(lbnd);\r
1150                                 PQinsert(lbnd, p, dist(p,newsite));\r
1151                         };\r
1152                         lbnd = bisector;                                                \r
1153                         bisector = HEcreate(e, re);                                     //create a new HalfEdge, setting its ELpm field to 1\r
1154                         ELinsert(lbnd, bisector);                                       //insert the new HE to the right of the original bisector earlier in the IF stmt\r
1155 \r
1156                         if ((p = intersect(bisector, rbnd)) != (struct Site *) NULL)    //if this new bisector intersects with the\r
1157                         {\r
1158                                 PQinsert(bisector, p, dist(p,newsite));                 //push the HE into the ordered linked list of vertices\r
1159                         };\r
1160                         newsite = nextone();\r
1161                 }\r
1162                 else if (!PQempty()) /* intersection is smallest - this is a vector event */\r
1163                 {\r
1164                         lbnd = PQextractmin();                                          //pop the HalfEdge with the lowest vector off the ordered list of vectors                               \r
1165                         llbnd = ELleft(lbnd);                                           //get the HalfEdge to the left of the above HE\r
1166                         rbnd = ELright(lbnd);                                           //get the HalfEdge to the right of the above HE\r
1167                         rrbnd = ELright(rbnd);                                          //get the HalfEdge to the right of the HE to the right of the lowest HE \r
1168                         bot = leftreg(lbnd);                                            //get the Site to the left of the left HE which it bisects\r
1169                         top = rightreg(rbnd);                                           //get the Site to the right of the right HE which it bisects\r
1170 \r
1171                         out_triple(bot, top, rightreg(lbnd));           //output the triple of sites, stating that a circle goes through them\r
1172 \r
1173                         v = lbnd->vertex;                                               //get the vertex that caused this event\r
1174                         makevertex(v);                                                  //set the vertex number - couldn't do this earlier since we didn't know when it would be processed\r
1175                         e2 = lbnd->ELedge;\r
1176                         e3 = rbnd->ELedge;\r
1177                         endpoint(lbnd->ELedge,lbnd->ELpm,v);    //set the endpoint of the left HalfEdge to be this vector\r
1178                         endpoint(rbnd->ELedge,rbnd->ELpm,v);    //set the endpoint of the right HalfEdge to be this vector\r
1179                         ELdelete(lbnd);                                                 //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map  \r
1180                         PQdelete(rbnd);                                                 //remove all vertex events to do with the  right HE\r
1181                         ELdelete(rbnd);                                                 //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map   \r
1182                         pm = le;                                                                //set the pm variable to zero\r
1183 \r
1184                         if (bot->coord.y > top->coord.y)                //if the site to the left of the event is higher than the Site\r
1185                         {                                                                               //to the right of it, then swap them and set the 'pm' variable to 1\r
1186                                 temp = bot;\r
1187                                 bot = top;\r
1188                                 top = temp;\r
1189                                 pm = re;\r
1190                         }\r
1191                         e = bisect(bot, top);                                   //create an Edge (or line) that is between the two Sites. This creates\r
1192                                                                                                         //the formula of the line, and assigns a line number to it\r
1193                         bisector = HEcreate(e, pm);                             //create a HE from the Edge 'e', and make it point to that edge with its ELedge field\r
1194                         ELinsert(llbnd, bisector);                              //insert the new bisector to the right of the left HE\r
1195                         endpoint(e, re-pm, v, e2, e3);                                  //set one endpoint to the new edge to be the vector point 'v'.\r
1196                                                                                                         //If the site to the left of this bisector is higher than the right\r
1197                                                                                                         //Site, then this endpoint is put in position 0; otherwise in pos 1\r
1198                         deref(v);                                                               //delete the vector 'v'\r
1199 \r
1200                         //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it\r
1201                         if((p = intersect(llbnd, bisector)) != (struct Site *) NULL)\r
1202                         {\r
1203                                 PQdelete(llbnd);\r
1204                                 PQinsert(llbnd, p, dist(p,bot));\r
1205                         };\r
1206 \r
1207                         //if right HE and the new bisector don't intersect, then reinsert it\r
1208                         if ((p = intersect(bisector, rrbnd)) != (struct Site *) NULL)\r
1209                         {\r
1210                                 PQinsert(bisector, p, dist(p,bot));\r
1211                         };\r
1212                 }\r
1213                 else break;\r
1214         };\r
1215 \r
1216         for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))\r
1217         {\r
1218                 e = lbnd -> ELedge;\r
1219 \r
1220                 clip_line(e);\r
1221         };\r
1222 \r
1223         cleanup();\r
1224 \r
1225         return true;\r
1226 }\r
1227 \r
1228 \r
1229 int scomp(const void *p1,const void *p2)\r
1230 {\r
1231         struct Point *s1 = (Point*)p1, *s2=(Point*)p2;\r
1232         if(s1 -> y < s2 -> y) return(-1);\r
1233         if(s1 -> y > s2 -> y) return(1);\r
1234         if(s1 -> x < s2 -> x) return(-1);\r
1235         if(s1 -> x > s2 -> x) return(1);\r
1236         return(0);\r
1237 }\r
1238 \r
1239 int spcomp(const void *p1,const void *p2)\r
1240 {\r
1241         struct SourcePoint *s1 = (SourcePoint*)p1, *s2=(SourcePoint*)p2;\r
1242         if(s1 -> y < s2 -> y) return(-1);\r
1243         if(s1 -> y > s2 -> y) return(1);\r
1244         if(s1 -> x < s2 -> x) return(-1);\r
1245         if(s1 -> x > s2 -> x) return(1);\r
1246         return(0);\r
1247 }\r
1248 \r
1249 int anglecomp(const void * p1, const void * p2)\r
1250 {\r
1251         PolygonPoint * s1 = (PolygonPoint *)p1 ;\r
1252         PolygonPoint * s2 = (PolygonPoint *)p2 ;\r
1253 \r
1254         if (s1->angle < s2->angle) {\r
1255                 return (-1) ;\r
1256         }\r
1257         if (s1->angle > s2->angle) {\r
1258                 return (1) ;\r
1259         }\r
1260         return (0) ;\r
1261 }\r
1262 \r
1263 /* return a single in-storage site */\r
1264 struct Site * VoronoiDiagramGenerator::nextone()\r
1265 {\r
1266         struct Site *s;\r
1267         if(siteidx < nsites)\r
1268         {\r
1269                 s = &sites[siteidx];\r
1270                 siteidx += 1;\r
1271                 return(s);\r
1272         }\r
1273         else\r
1274                 return( (struct Site *)NULL);\r
1275 }\r
1276 \r