]> git.openstreetmap.org Git - nominatim.git/blob - nominatim/voronoi/VoronoiDiagramGenerator.cpp
merge in the postcode and is_in terms into the address vector for search
[nominatim.git] / 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