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