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