FastJet 3.0.2
Voronoi.cc
00001 //STARTHEADER
00002 // $Id: Voronoi.cc 2767 2011-11-24 21:43:09Z salam $
00003 //
00004 // Copyright (c) 1994 by AT&T Bell Laboratories (see below)
00005 //
00006 //
00007 //----------------------------------------------------------------------
00008 // This file is included as part of FastJet but was mostly written by
00009 // S. Fortune in C, put into C++ with memory management by S
00010 // O'Sullivan, and with further interface and memory management
00011 // modifications by Gregory Soyez.
00012 //
00013 // Permission to use, copy, modify, and distribute this software for
00014 // any purpose without fee is hereby granted, provided that this
00015 // entire notice is included in all copies of any software which is or
00016 // includes a copy or modification of this software and in all copies
00017 // of the supporting documentation for such software. THIS SOFTWARE IS
00018 // BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED WARRANTY.
00019 // IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY REPRESENTATION
00020 // OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS
00021 // SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00022 //
00023 //----------------------------------------------------------------------
00024 //ENDHEADER
00025 
00026 
00027 /*
00028  * The author of this software is Steven Fortune.  
00029  * Copyright (c) 1994 by AT&T Bell Laboratories.
00030  * Permission to use, copy, modify, and distribute this software for any
00031  * purpose without fee is hereby granted, provided that this entire notice
00032  * is included in all copies of any software which is or includes a copy
00033  * or modification of this software and in all copies of the supporting
00034  * documentation for such software.
00035  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00036  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
00037  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
00038  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
00039  */
00040 
00041 /* 
00042  * This code was originally written by Stephan Fortune in C code.  I,
00043  * Shane O'Sullivan, have since modified it, encapsulating it in a C++
00044  * class and, fixing memory leaks and adding accessors to the Voronoi
00045  * Edges.  Permission to use, copy, modify, and distribute this
00046  * software for any purpose without fee is hereby granted, provided
00047  * that this entire notice is included in all copies of any software
00048  * which is or includes a copy or modification of this software and in
00049  * all copies of the supporting documentation for such software.  THIS
00050  * SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
00051  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
00052  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE
00053  * MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR
00054  * PURPOSE.
00055  */
00056 
00057 /*
00058  * This code, included in the FastJet distribution, was originally
00059  * written by Stephan Fortune in C and adapted to C++ by Shane
00060  * O'Sullivan under the terms reported above.
00061  *
00062  * Below are the list of changes implemented by the FastJet authors:
00063  *
00064  * 2011-11-14  Gregory Soyez  <soyez@fastjet.fr>
00065  * 
00066  *      * removed 'plot' and 'triangulate' (were always 0)
00067  *      * removed unused plot functions (openpl, circle, range, 
00068  *        out_bisector, out_ep, out_vertex, out_site, out_triple)
00069  *      * removed unused 'VPoint p' in 'intersect'
00070  *
00071  *
00072  * 2011-07-22  Gregory Soyez  <soyez@fastjet.fr>
00073  * 
00074  *      * replaced Point by VPoint (to avoid any potential conflict
00075  *        with an already existing class Point in FastJet
00076  * 
00077  * 
00078  * 2011-06-28  Gregory Soyez  <soyez@fastjet.fr>
00079  * 
00080  *      * added support for situations with degenerate particles (we just
00081  *        discard the particles degenerate wiht an already existing
00082  *        one which is perfectly sufficient for our needs)
00083  *      * in 'VoronoiDiagramGenerator::intersect', improved the numerical
00084  *        precision in cases where the 2 parents are nearly degenerate
00085  * 
00086  * 
00087  * 2011-06-14  Gregory Soyez  <soyez@fastjet.fr>
00088  * 
00089  *      * fixed a potential overflow bug in VoronoiDiagramGenerator::PQbucket
00090  * 
00091  * 
00092  * 2007-05-07  Gregory Soyez  <soyez@fastjet.fr>
00093  * 
00094  *      * fied a few memory leaks
00095  *
00096  *      * put the code in the fastjet namespace
00097  * 
00098  *      * replaced float by double
00099  * 
00100  *      * generateVoronoi() takes a vector of Point instead of 2
00101  *        pointers
00102  * 
00103  *      * added info about the parent sites to GraphEdge (and clip_line)
00104  * 
00105  *      * removed condition on minimal distance between sites
00106  */
00107 
00108 #include <stdio.h>
00109 #include "fastjet/internal/Voronoi.hh"
00110 
00111 using namespace std;
00112 
00113 FASTJET_BEGIN_NAMESPACE
00114 
00115 LimitedWarning VoronoiDiagramGenerator::_warning_degeneracy;
00116 
00117 VoronoiDiagramGenerator::VoronoiDiagramGenerator(){
00118   siteidx = 0;
00119   sites = NULL;
00120   parent_sites = NULL;
00121 
00122   allMemoryList = new FreeNodeArrayList;
00123   allMemoryList->memory = NULL;
00124   allMemoryList->next = NULL;
00125   currentMemoryBlock = allMemoryList;
00126   allEdges = NULL;
00127   iteratorEdges = 0;
00128   minDistanceBetweenSites = 0;
00129 
00130   ELhash = NULL;
00131   PQhash = NULL;
00132 }
00133 
00134 VoronoiDiagramGenerator::~VoronoiDiagramGenerator(){
00135   cleanup();
00136   cleanupEdges();
00137 
00138   if (allMemoryList != NULL)
00139     delete allMemoryList;
00140 }
00141 
00142 
00143 
00144 bool VoronoiDiagramGenerator::generateVoronoi(vector<VPoint> *_parent_sites,
00145                                               double minX, double maxX, 
00146                                               double minY, double maxY, 
00147                                               double minDist){
00148   cleanup();
00149   cleanupEdges();
00150   int i;
00151   double x, y;
00152 
00153   minDistanceBetweenSites = minDist;
00154 
00155   parent_sites = _parent_sites;
00156 
00157   nsites = n_parent_sites = parent_sites->size();
00158   debug = 1;
00159   sorted = 0; 
00160   freeinit(&sfl, sizeof (Site));
00161                 
00162   //sites = (Site *) myalloc(nsites*sizeof( *sites));
00163   sites = (Site *) myalloc(nsites*sizeof( Site));
00164   //parent_sites = (Site *) myalloc(nsites*sizeof( Site));
00165  
00166   if (sites == 0)
00167     return false;
00168 
00169   xmax = xmin = (*parent_sites)[0].x;
00170   ymax = ymin = (*parent_sites)[0].y;
00171   
00172   for(i=0;i<nsites;i++){
00173     x = (*parent_sites)[i].x;
00174     y = (*parent_sites)[i].y;
00175     sites[i].coord.x = x;
00176     sites[i].coord.y = y;
00177     sites[i].sitenbr = i;
00178     sites[i].refcnt  = 0;
00179     
00180     if (x<xmin)
00181       xmin=x;
00182     else if (x>xmax)
00183       xmax=x;
00184     
00185     if (y<ymin)
00186       ymin=y;
00187     else if (y>ymax)
00188       ymax=y;
00189   }
00190         
00191   qsort(sites, nsites, sizeof (*sites), scomp);
00192 
00193   // Gregory Soyez
00194   // 
00195   // check if some of the particles are degenerate to avoid a crash.
00196   //
00197   // At the moment, we work under the assumption that they will be
00198   // "clustered" later on, so we just keep the 1st one and discard the
00199   // others
00200   unsigned int offset=0;
00201   for (int is=1;is<nsites;is++){
00202     if (sites[is].coord.y==sites[is-1].coord.y && sites[is].coord.x==sites[is-1].coord.x){
00203       offset++;
00204     } else if (offset>0){
00205       sites[is-offset] = sites[is];
00206     }
00207   }
00208 
00209   if (offset>0){
00210     nsites-=offset;
00211     _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.");
00212   }
00213 
00214   siteidx = 0;
00215   geominit();
00216   double temp = 0;
00217   if(minX > maxX){
00218     temp = minX;
00219     minX = maxX;
00220     maxX = temp;
00221   }
00222   if(minY > maxY){
00223     temp = minY;
00224     minY = maxY;
00225     maxY = temp;
00226   }
00227   borderMinX = minX;
00228   borderMinY = minY;
00229   borderMaxX = maxX;
00230   borderMaxY = maxY;
00231         
00232   siteidx = 0;
00233   voronoi();
00234 
00235   return true;
00236 }
00237 
00238 bool VoronoiDiagramGenerator::ELinitialize(){
00239   int i;
00240   freeinit(&hfl, sizeof(Halfedge));
00241   ELhashsize = 2 * sqrt_nsites;
00242   //ELhash = (Halfedge **) myalloc ( sizeof *ELhash * ELhashsize);
00243   ELhash = (Halfedge **) myalloc ( sizeof(Halfedge*)*ELhashsize);
00244   
00245   if(ELhash == 0)
00246     return false;
00247   
00248   for(i=0; i<ELhashsize; i +=1) ELhash[i] = (Halfedge *)NULL;
00249   ELleftend = HEcreate((Edge*) NULL, 0);
00250   ELrightend = HEcreate((Edge*) NULL, 0);
00251   ELleftend->ELleft = (Halfedge*) NULL;
00252   ELleftend->ELright = ELrightend;
00253   ELrightend->ELleft = ELleftend;
00254   ELrightend->ELright = (Halfedge*) NULL;
00255   ELhash[0] = ELleftend;
00256   ELhash[ELhashsize-1] = ELrightend;
00257 
00258   return true;
00259 }
00260 
00261 
00262 Halfedge* VoronoiDiagramGenerator::HEcreate(Edge *e,int pm){
00263   Halfedge *answer;
00264   answer = (Halfedge*) getfree(&hfl);
00265   answer->ELedge = e;
00266   answer->ELpm = pm;
00267   answer->PQnext = (Halfedge*) NULL;
00268   answer->vertex = (Site*) NULL;
00269   answer->ELrefcnt = 0;
00270   return answer;
00271 }
00272 
00273 
00274 void VoronoiDiagramGenerator::ELinsert(Halfedge *lb, Halfedge *newHe)
00275 {
00276   newHe->ELleft = lb;
00277   newHe->ELright = lb->ELright;
00278   (lb->ELright)->ELleft = newHe;
00279   lb->ELright = newHe;
00280 }
00281 
00282 
00283 /* Get entry from hash table, pruning any deleted nodes */
00284 Halfedge* VoronoiDiagramGenerator::ELgethash(int b){
00285   Halfedge *he;
00286         
00287   if ((b<0) || (b>=ELhashsize)) 
00288     return (Halfedge *) NULL;
00289 
00290   he = ELhash[b]; 
00291   if ((he==(Halfedge*) NULL) || (he->ELedge!=(Edge*) DELETED)) 
00292     return he;
00293         
00294   /* Hash table points to deleted half edge.  Patch as necessary. */
00295   ELhash[b] = (Halfedge*) NULL;
00296   if ((he->ELrefcnt -= 1) == 0) 
00297     makefree((Freenode*)he, &hfl);
00298   return (Halfedge*) NULL;
00299 }       
00300 
00301 Halfedge * VoronoiDiagramGenerator::ELleftbnd(VPoint *p){
00302   int i, bucket;
00303   Halfedge *he;
00304         
00305   /* Use hash table to get close to desired halfedge */
00306   // use the hash function to find the place in the hash map that this
00307   // HalfEdge should be
00308   // Gregory Soyez: the original code was 
00309   //
00310   //   bucket = (int)((p->x - xmin)/deltax * ELhashsize);       
00311   //   // make sure that the bucket position in within the range of the
00312   //   //hash array
00313   //   if (bucket<0) bucket =0;
00314   //   if (bucket>=ELhashsize) bucket = ELhashsize - 1;
00315   //
00316   // but this runs the risk of having a overflow which would 
00317   // cause bucket to be truncated at 0 instead of ELhashsize - 1
00318   // (or vice-versa)
00319   // We fix this by performing the test immediately on the double
00320   // We put in a extra bit of margin to be sure the conversion does
00321   // not play dirty tricks on us
00322 
00323   //const double &px = p->x;
00324   if (p->x < xmin){ bucket=0;}
00325   else if (p->x >= xmax){ bucket = ELhashsize - 1;}
00326   else{
00327     bucket= (int)((p->x - xmin)/deltax * ELhashsize);
00328     if (bucket>=ELhashsize) bucket = ELhashsize - 1;  // the lower cut should be robust
00329   }
00330 
00331   he = ELgethash(bucket);
00332 
00333   // if the HE isn't found, search backwards and forwards in the hash
00334   // map for the first non-null entry
00335   if (he == (Halfedge*) NULL){   
00336     for(i=1;1;i++){     
00337       if ((he=ELgethash(bucket-i)) != (Halfedge*) NULL) 
00338         break;
00339       if ((he=ELgethash(bucket+i)) != (Halfedge*) NULL) 
00340         break;
00341     };
00342     totalsearch += i;
00343   };
00344   ntry += 1;
00345   /* Now search linear list of halfedges for the correct one */
00346   if ((he==ELleftend)  || (he != ELrightend && right_of(he,p))){
00347     do{
00348       he = he->ELright;
00349     } while (he!=ELrightend && right_of(he,p));
00350       // keep going right on the list until either the end is reached,
00351       // or you find the 1st edge which the point
00352       he = he->ELleft;  //isn't to the right of
00353   } else 
00354     // if the point is to the left of the HalfEdge, then search left
00355     // for the HE just to the left of the point
00356     do{
00357       he = he->ELleft;
00358     } while ((he!=ELleftend )&& (!right_of(he,p)));
00359                 
00360   /* Update hash table and reference counts */
00361   if ((bucket > 0) && (bucket <ELhashsize-1)){  
00362     if(ELhash[bucket] != (Halfedge *) NULL){
00363       ELhash[bucket]->ELrefcnt -= 1;
00364     }
00365     ELhash[bucket] = he;
00366     ELhash[bucket]->ELrefcnt += 1;
00367   };
00368 
00369   return he;
00370 }
00371 
00372 
00373 /* This delete routine can't reclaim node, since pointers from hash
00374    table may be present.   */
00375 void VoronoiDiagramGenerator::ELdelete(Halfedge *he){
00376   (he->ELleft)->ELright = he->ELright;
00377   (he->ELright)->ELleft = he->ELleft;
00378   he->ELedge = (Edge*) DELETED;
00379 }
00380 
00381 
00382 Halfedge* VoronoiDiagramGenerator::ELright(Halfedge *he){
00383   return he->ELright;
00384 }
00385 
00386 
00387 Halfedge* VoronoiDiagramGenerator::ELleft(Halfedge *he){
00388   return he->ELleft;
00389 }
00390 
00391 
00392 Site * VoronoiDiagramGenerator::leftreg(Halfedge *he){
00393   if (he->ELedge == (Edge*) NULL) 
00394     return(bottomsite);
00395   return (he->ELpm == le) 
00396     ? he->ELedge->reg[le] 
00397     : he->ELedge->reg[re];
00398 }
00399 
00400 Site * VoronoiDiagramGenerator::rightreg(Halfedge *he){
00401   if (he->ELedge == (Edge*) NULL)
00402     // if this halfedge has no edge, return the bottom site (whatever
00403     // that is)
00404     return bottomsite;
00405   
00406   // if the ELpm field is zero, return the site 0 that this edge
00407   // bisects, otherwise return site number 1
00408   return (he->ELpm == le)
00409     ? he->ELedge->reg[re] 
00410     : he->ELedge->reg[le];
00411 }
00412 
00413 void VoronoiDiagramGenerator::geominit(){
00414   double sn;
00415   
00416   freeinit(&efl, sizeof(Edge));
00417   nvertices = 0;
00418   nedges = 0;
00419   sn = (double)nsites+4;
00420   sqrt_nsites = (int)sqrt(sn);
00421   deltay = ymax - ymin;
00422   deltax = xmax - xmin;
00423 }
00424 
00425 
00426 Edge * VoronoiDiagramGenerator::bisect(Site *s1, Site *s2){
00427   double dx,dy,adx,ady;
00428   Edge *newedge;        
00429 
00430   newedge = (Edge*) getfree(&efl);
00431         
00432   newedge->reg[0] = s1; //store the sites that this edge is bisecting
00433   newedge->reg[1] = s2;
00434   ref(s1); 
00435   ref(s2);
00436 
00437   // to begin with, there are no endpoints on the bisector - it goes
00438   // to infinity
00439   newedge->ep[0] = (Site*) NULL;
00440   newedge->ep[1] = (Site*) NULL;
00441   
00442   // get the difference in x dist between the sites
00443   dx = s2->coord.x - s1->coord.x;
00444   dy = s2->coord.y - s1->coord.y;
00445 
00446   // make sure that the difference is positive
00447   adx = dx>0 ? dx : -dx;
00448   ady = dy>0 ? dy : -dy;
00449 
00450   // get the slope of the line
00451   newedge->c = (double)(s1->coord.x * dx + s1->coord.y * dy
00452                         + (dx*dx + dy*dy)*0.5);
00453 
00454   if (adx>ady){ 
00455     //set formula of line, with x fixed to 1
00456     newedge->a = 1.0; newedge->b = dy/dx; newedge->c /= dx;
00457   } else {      
00458     //set formula of line, with y fixed to 1
00459     newedge->b = 1.0; newedge->a = dx/dy; newedge->c /= dy;
00460   }
00461         
00462   newedge->edgenbr = nedges;
00463   nedges++;
00464 
00465   return newedge;
00466 }
00467 
00468 
00469 // create a new site where the HalfEdges el1 and el2 intersect - note
00470 // that the VPoint in the argument list is not used, don't know why
00471 // it's there
00472 //
00473 // Gregory Soyez: removed the uinused point p
00474 Site* VoronoiDiagramGenerator::intersect(Halfedge *el1, Halfedge *el2
00475                                          /*, VPoint *p*/){
00476   Edge *e1,*e2, *e;
00477   Halfedge *el;
00478   double d, xint, yint;
00479   int right_of_site;
00480   Site *v;
00481         
00482   e1 = el1->ELedge;
00483   e2 = el2->ELedge;
00484   if ((e1==(Edge*)NULL) || (e2==(Edge*)NULL)) 
00485     return ((Site*) NULL);
00486 
00487   // if the two edges bisect the same parent, return null
00488   if (e1->reg[1] == e2->reg[1]) 
00489     return (Site*) NULL;
00490 
00491   // Gregory Soyez:
00492   //    
00493   // if the 2 parents are too close, the intersection is going to be
00494   // computed from the "long edges" of the triangle which could causes
00495   // large rounding errors. In this case, use the bisector of the 2
00496   // parents to find the interaction point
00497   // 
00498   // The following replaces 
00499   //   d = e1->a * e2->b - e1->b * e2->a;
00500   //   if (-1.0e-10<d && d<1.0e-10) 
00501   //     return (Site*) NULL;
00502   //    
00503   //   xint = (e1->c*e2->b - e2->c*e1->b)/d;
00504   //   yint = (e2->c*e1->a - e1->c*e2->a)/d;
00505 
00506   double dx = e2->reg[1]->coord.x - e1->reg[1]->coord.x;
00507   double dy = e2->reg[1]->coord.y - e1->reg[1]->coord.y;
00508   double dxref = e1->reg[1]->coord.x - e1->reg[0]->coord.x;
00509   double dyref = e1->reg[1]->coord.y - e1->reg[0]->coord.y;
00510 
00511   if (dx*dx + dy*dy < 1e-14*(dxref*dxref+dyref*dyref)){
00512     // make sure that the difference is positive
00513     double adx = dx>0 ? dx : -dx;
00514     double ady = dy>0 ? dy : -dy;
00515     
00516     // get the slope of the line
00517     double a,b;
00518     double c = (double)(e1->reg[1]->coord.x * dx + e1->reg[1]->coord.y * dy
00519                         + (dx*dx + dy*dy)*0.5);
00520     
00521     if (adx>ady){
00522       a = 1.0; b = dy/dx; c /= dx;
00523     } else {
00524       b = 1.0; a = dx/dy; c /= dy;
00525     }
00526 
00527     d = e1->a * b - e1->b * a;
00528     if (-1.0e-10<d && d<1.0e-10) {
00529       return (Site*) NULL;
00530     }
00531         
00532     xint = (e1->c*b - c*e1->b)/d;
00533     yint = (c*e1->a - e1->c*a)/d;
00534     
00535   } else {      
00536     d = e1->a * e2->b - e1->b * e2->a;
00537     if (-1.0e-10<d && d<1.0e-10) {
00538       return (Site*) NULL;
00539     }
00540         
00541     xint = (e1->c*e2->b - e2->c*e1->b)/d;
00542     yint = (e2->c*e1->a - e1->c*e2->a)/d;
00543   }
00544   // end of Gregory Soyez's modifications
00545 
00546   volatile double local_y1 = e1->reg[1]->coord.y;
00547   volatile double local_y2 = e2->reg[1]->coord.y;
00548         
00549   if( (local_y1 < local_y2) ||
00550       ((local_y1 == local_y2) &&
00551        (e1->reg[1]->coord.x <  e2->reg[1]->coord.x)) ){ 
00552     el = el1; 
00553     e = e1;
00554   } else {      
00555     el = el2; 
00556     e = e2;
00557   }
00558         
00559   right_of_site = xint >= e->reg[1]->coord.x;
00560   if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re)) 
00561     return (Site*) NULL;
00562         
00563   // create a new site at the point of intersection - this is a new
00564   // vector event waiting to happen
00565   v = (Site*) getfree(&sfl);
00566   v->refcnt = 0;
00567   v->coord.x = xint;
00568   v->coord.y = yint;
00569   return v;
00570 }
00571 
00572 //HERE
00573 
00574 /* returns 1 if p is to right of halfedge e */
00575 int VoronoiDiagramGenerator::right_of(Halfedge *el,VPoint *p)
00576 {
00577   Edge *e;
00578   Site *topsite;
00579   int right_of_site, above, fast;
00580   double dxp, dyp, dxs, t1, t2, t3, yl;
00581         
00582   e = el->ELedge;
00583   topsite = e->reg[1];
00584   right_of_site = p->x > topsite->coord.x;
00585   if(right_of_site && el->ELpm == le) return(1);
00586   if(!right_of_site && el->ELpm == re) return (0);
00587         
00588   if (e->a == 1.0)
00589     {   dyp = p->y - topsite->coord.y;
00590     dxp = p->x - topsite->coord.x;
00591     fast = 0;
00592     if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )
00593       { above = dyp>= e->b*dxp; 
00594       fast = above;
00595       }
00596     else 
00597       { above = p->x + p->y*e->b > e-> c;
00598       if(e->b<0.0) above = !above;
00599       if (!above) fast = 1;
00600       };
00601     if (!fast)
00602       { dxs = topsite->coord.x - (e->reg[0])->coord.x;
00603       above = e->b * (dxp*dxp - dyp*dyp) <
00604         dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
00605       if(e->b<0.0) above = !above;
00606       };
00607     }
00608   else  /*e->b==1.0 */
00609     {   yl = e->c - e->a*p->x;
00610     t1 = p->y - yl;
00611     t2 = p->x - topsite->coord.x;
00612     t3 = yl - topsite->coord.y;
00613     above = t1*t1 > t2*t2 + t3*t3;
00614     };
00615   return (el->ELpm==le ? above : !above);
00616 }
00617 
00618 
00619 void VoronoiDiagramGenerator::endpoint(Edge *e,int lr,Site * s)
00620 {
00621   e->ep[lr] = s;
00622   ref(s);
00623   if(e->ep[re-lr]== (Site *) NULL) 
00624     return;
00625 
00626   clip_line(e);
00627 
00628   deref(e->reg[le]);
00629   deref(e->reg[re]);
00630   makefree((Freenode*)e, &efl);
00631 }
00632 
00633 
00634 double VoronoiDiagramGenerator::dist(Site *s,Site *t)
00635 {
00636   double dx,dy;
00637   dx = s->coord.x - t->coord.x;
00638   dy = s->coord.y - t->coord.y;
00639   return (double)(sqrt(dx*dx + dy*dy));
00640 }
00641 
00642 
00643 void VoronoiDiagramGenerator::makevertex(Site *v)
00644 {
00645   v->sitenbr = nvertices;
00646   nvertices += 1;
00647   //GS unused plot: out_vertex(v);
00648 }
00649 
00650 
00651 void VoronoiDiagramGenerator::deref(Site *v)
00652 {
00653   v->refcnt -= 1;
00654   if (v->refcnt == 0 ) 
00655     makefree((Freenode*)v, &sfl);
00656 }
00657 
00658 void VoronoiDiagramGenerator::ref(Site *v)
00659 {
00660   v->refcnt += 1;
00661 }
00662 
00663 //push the HalfEdge into the ordered linked list of vertices
00664 void VoronoiDiagramGenerator::PQinsert(Halfedge *he,Site * v, double offset)
00665 {
00666   Halfedge *last, *next;
00667         
00668   he->vertex = v;
00669   ref(v);
00670   he->ystar = (double)(v->coord.y + offset);
00671   last = &PQhash[PQbucket(he)];
00672   while ((next = last->PQnext) != (Halfedge *) NULL &&
00673          (he->ystar  > next->ystar  ||
00674           (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x)))
00675     {   
00676       last = next;
00677     };
00678   he->PQnext = last->PQnext; 
00679   last->PQnext = he;
00680   PQcount += 1;
00681 }
00682 
00683 //remove the HalfEdge from the list of vertices 
00684 void VoronoiDiagramGenerator::PQdelete(Halfedge *he)
00685 {
00686   Halfedge *last;
00687         
00688   if(he->vertex != (Site *) NULL)
00689     {   
00690       last = &PQhash[PQbucket(he)];
00691       while (last->PQnext != he) 
00692         last = last->PQnext;
00693 
00694       last->PQnext = he->PQnext;
00695       PQcount -= 1;
00696       deref(he->vertex);
00697       he->vertex = (Site *) NULL;
00698     };
00699 }
00700 
00701 int VoronoiDiagramGenerator::PQbucket(Halfedge *he)
00702 {
00703   // Gregory Soyez: the original code was 
00704   //
00705   //   bucket = (int)((he->ystar - ymin)/deltay * PQhashsize);
00706   //   if (bucket<0) bucket = 0;
00707   //   if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
00708   //   if (bucket < PQmin) PQmin = bucket;
00709   //   return(bucket);
00710   //
00711   // but this runs the risk of having a overflow which would 
00712   // cause bucket to be truncated at 0 instead of PQhashsize-1
00713   // (or vice-versa)
00714   // We fix this by performing the test immediately on the double
00715   // We put in a extra bit of margin to be sure the conversion does
00716   // not play dirty tricks on us
00717 
00718   int bucket;
00719         
00720   double hey = he->ystar;
00721   if (hey < ymin){ bucket = 0;}
00722   else if (hey >= ymax){ bucket = PQhashsize-1;}
00723   else {
00724     bucket = (int)((hey - ymin)/deltay * PQhashsize);
00725     if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
00726   }
00727 
00728   if (bucket < PQmin) PQmin = bucket;
00729   return(bucket);
00730 }
00731 
00732 
00733 
00734 int VoronoiDiagramGenerator::PQempty()
00735 {
00736   return(PQcount==0);
00737 }
00738 
00739 
00740 VPoint VoronoiDiagramGenerator::PQ_min()
00741 {
00742   VPoint answer;
00743         
00744   while(PQhash[PQmin].PQnext == (Halfedge *)NULL) {PQmin += 1;};
00745   answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
00746   answer.y = PQhash[PQmin].PQnext->ystar;
00747   return (answer);
00748 }
00749 
00750 Halfedge * VoronoiDiagramGenerator::PQextractmin()
00751 {
00752   Halfedge *curr;
00753         
00754   curr = PQhash[PQmin].PQnext;
00755   PQhash[PQmin].PQnext = curr->PQnext;
00756   PQcount -= 1;
00757   return(curr);
00758 }
00759 
00760 
00761 bool VoronoiDiagramGenerator::PQinitialize()
00762 {
00763   int i; 
00764         
00765   PQcount = 0;
00766   PQmin = 0;
00767   PQhashsize = 4 * sqrt_nsites;
00768   //PQhash = (Halfedge *) myalloc(PQhashsize * sizeof *PQhash);
00769   PQhash = (Halfedge *) myalloc(PQhashsize * sizeof(Halfedge));
00770 
00771   if(PQhash == 0)
00772     return false;
00773 
00774   for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (Halfedge *)NULL;
00775 
00776   return true;
00777 }
00778 
00779 
00780 void VoronoiDiagramGenerator::freeinit(Freelist *fl,int size)
00781 {
00782   fl->head = (Freenode *) NULL;
00783   fl->nodesize = size;
00784 }
00785 
00786 char * VoronoiDiagramGenerator::getfree(Freelist *fl)
00787 {
00788   int i; 
00789   Freenode *t;
00790 
00791   if(fl->head == (Freenode *) NULL)
00792     {   
00793       t =  (Freenode *) myalloc(sqrt_nsites * fl->nodesize);
00794 
00795       if(t == 0)
00796         return 0;
00797                 
00798       currentMemoryBlock->next = new FreeNodeArrayList;
00799       currentMemoryBlock = currentMemoryBlock->next;
00800       currentMemoryBlock->memory = t;
00801       currentMemoryBlock->next = 0;
00802 
00803       for(i=0; i<sqrt_nsites; i+=1)     
00804         makefree((Freenode *)((char *)t+i*fl->nodesize), fl);           
00805     };
00806   t = fl->head;
00807   fl->head = (fl->head)->nextfree;
00808   return((char *)t);
00809 }
00810 
00811 
00812 
00813 void VoronoiDiagramGenerator::makefree(Freenode *curr,Freelist *fl)
00814 {
00815   curr->nextfree = fl->head;
00816   fl->head = curr;
00817 }
00818 
00819 void VoronoiDiagramGenerator::cleanup()
00820 {
00821   if(sites != NULL){
00822     free(sites);
00823     sites = 0;
00824   }
00825 
00826   FreeNodeArrayList* current=NULL, *prev=NULL;
00827 
00828   current = prev = allMemoryList;
00829 
00830   while(current->next!=NULL){
00831     prev = current;
00832     current = current->next;
00833     free(prev->memory);
00834     delete prev;
00835     prev = NULL;
00836   }
00837 
00838   if(current!=NULL){
00839     if (current->memory!=NULL){
00840       free(current->memory);
00841     }
00842     delete current;
00843   }
00844 
00845   allMemoryList = new FreeNodeArrayList;
00846   allMemoryList->next = NULL;
00847   allMemoryList->memory = NULL;
00848   currentMemoryBlock = allMemoryList;
00849 
00850   if (ELhash!=NULL)
00851     free(ELhash);
00852 
00853   if (PQhash!=NULL)
00854     free(PQhash);
00855 }
00856 
00857 void VoronoiDiagramGenerator::cleanupEdges()
00858 {
00859   GraphEdge* geCurrent = NULL, *gePrev = NULL;
00860   geCurrent = gePrev = allEdges;
00861 
00862   while(geCurrent!=NULL){
00863     gePrev = geCurrent;
00864     geCurrent = geCurrent->next;
00865     delete gePrev;
00866   }
00867 
00868   allEdges = 0;
00869 
00870 }
00871 
00872 void VoronoiDiagramGenerator::pushGraphEdge(double x1, double y1, double x2, double y2,
00873                                             Site *s1, Site *s2)
00874 {
00875   GraphEdge* newEdge = new GraphEdge;
00876   newEdge->next = allEdges;
00877   allEdges = newEdge;
00878   newEdge->x1 = x1;
00879   newEdge->y1 = y1;
00880   newEdge->x2 = x2;
00881   newEdge->y2 = y2;
00882 
00883   newEdge->point1 = s1->sitenbr;
00884   newEdge->point2 = s2->sitenbr;
00885 }
00886 
00887 
00888 char * VoronoiDiagramGenerator::myalloc(unsigned n)
00889 {
00890   char *t=0;    
00891   t=(char*)malloc(n);
00892   total_alloc += n;
00893   return(t);
00894 }
00895 
00896 
00897 // unused plot functions
00898 //
00899 // /* for those who don't have Cherry's plot */
00900 // /* #include <plot.h> */
00901 // void VoronoiDiagramGenerator::openpl(){}
00902 // void VoronoiDiagramGenerator::circle(double x, double y, double radius){}
00903 // void VoronoiDiagramGenerator::range(double minX, double minY, double maxX, double maxY){}
00904 // 
00905 // 
00906 // 
00907 // void VoronoiDiagramGenerator::out_bisector(Edge *e)
00908 // {
00909 //      
00910 // 
00911 // }
00912 // 
00913 // 
00914 // void VoronoiDiagramGenerator::out_ep(Edge *e)
00915 // {
00916 //      
00917 //      
00918 // }
00919 // 
00920 // void VoronoiDiagramGenerator::out_vertex(Site *v)
00921 // {
00922 //      
00923 // }
00924 // 
00925 // 
00926 // void VoronoiDiagramGenerator::out_site(Site *s)
00927 // {
00928 //   // Gregory Soyez: 
00929 //   //   plot was always 0 so the expression below was always false
00930 //   //   and even if it was not, 'circle' does nothing!
00931 //   //
00932 //   // if(!triangulate & plot & !debug)
00933 //   //   circle (s->coord.x, s->coord.y, cradius);
00934 //      
00935 // }
00936 // 
00937 // 
00938 // void VoronoiDiagramGenerator::out_triple(Site *s1, Site *s2,Site * s3)
00939 // {
00940 //      
00941 // }
00942 
00943 
00944 
00945 void VoronoiDiagramGenerator::plotinit()
00946 {
00947   double dx,dy,d;
00948         
00949   dy = ymax - ymin;
00950   dx = xmax - xmin;
00951   d = (double)(( dx > dy ? dx : dy) * 1.1);
00952   pxmin = (double)(xmin - (d-dx)/2.0);
00953   pxmax = (double)(xmax + (d-dx)/2.0);
00954   pymin = (double)(ymin - (d-dy)/2.0);
00955   pymax = (double)(ymax + (d-dy)/2.0);
00956   cradius = (double)((pxmax - pxmin)/350.0);
00957   //GS unused: openpl();
00958   //GS unused: range(pxmin, pymin, pxmax, pymax);
00959 }
00960 
00961 
00962 void VoronoiDiagramGenerator::clip_line(Edge *e)
00963 {
00964   Site *s1, *s2;
00965   double x1=0,x2=0,y1=0,y2=0; //, temp = 0;
00966 
00967   x1 = e->reg[0]->coord.x;
00968   x2 = e->reg[1]->coord.x;
00969   y1 = e->reg[0]->coord.y;
00970   y2 = e->reg[1]->coord.y;
00971 
00972   //if the distance between the two points this line was created from is less than 
00973   //the square root of 2, then ignore it
00974   //TODO improve/remove
00975   //if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites)
00976   //  {
00977   //    return;
00978   //  }
00979   pxmin = borderMinX;
00980   pxmax = borderMaxX;
00981   pymin = borderMinY;
00982   pymax = borderMaxY;
00983 
00984   if(e->a == 1.0 && e ->b >= 0.0)
00985     {   
00986       s1 = e->ep[1];
00987       s2 = e->ep[0];
00988     }
00989   else 
00990     {
00991       s1 = e->ep[0];
00992       s2 = e->ep[1];
00993     };
00994         
00995   if(e->a == 1.0)
00996     {
00997       y1 = pymin;
00998       if (s1!=(Site *)NULL && s1->coord.y > pymin)
00999         {
01000           y1 = s1->coord.y;
01001         }
01002       if(y1>pymax) 
01003         {
01004           //    printf("\nClipped (1) y1 = %f to %f",y1,pymax);
01005           y1 = pymax;
01006           //return;
01007         }
01008       x1 = e->c - e->b * y1;
01009       y2 = pymax;
01010       if (s2!=(Site *)NULL && s2->coord.y < pymax) 
01011         y2 = s2->coord.y;
01012 
01013       if(y2<pymin) 
01014         {
01015           //printf("\nClipped (2) y2 = %f to %f",y2,pymin);
01016           y2 = pymin;
01017           //return;
01018         }
01019       x2 = (e->c) - (e->b) * y2;
01020       if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin))) 
01021         {
01022           //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);
01023           return;
01024         }
01025       if(x1> pxmax)
01026         {       x1 = pxmax; y1 = (e->c - x1)/e->b;};
01027       if(x1<pxmin)
01028         {       x1 = pxmin; y1 = (e->c - x1)/e->b;};
01029       if(x2>pxmax)
01030         {       x2 = pxmax; y2 = (e->c - x2)/e->b;};
01031       if(x2<pxmin)
01032         {       x2 = pxmin; y2 = (e->c - x2)/e->b;};
01033     }
01034   else
01035     {
01036       x1 = pxmin;
01037       if (s1!=(Site *)NULL && s1->coord.x > pxmin) 
01038         x1 = s1->coord.x;
01039       if(x1>pxmax) 
01040         {
01041           //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);
01042           //return;
01043           x1 = pxmax;
01044         }
01045       y1 = e->c - e->a * x1;
01046       x2 = pxmax;
01047       if (s2!=(Site *)NULL && s2->coord.x < pxmax) 
01048         x2 = s2->coord.x;
01049       if(x2<pxmin) 
01050         {
01051           //printf("\nClipped (4) x2 = %f to %f",x2,pxmin);
01052           //return;
01053           x2 = pxmin;
01054         }
01055       y2 = e->c - e->a * x2;
01056       if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin))) 
01057         {
01058           //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);
01059           return;
01060         }
01061       if(y1> pymax)
01062         {       y1 = pymax; x1 = (e->c - y1)/e->a;};
01063       if(y1<pymin)
01064         {       y1 = pymin; x1 = (e->c - y1)/e->a;};
01065       if(y2>pymax)
01066         {       y2 = pymax; x2 = (e->c - y2)/e->a;};
01067       if(y2<pymin)
01068         {       y2 = pymin; x2 = (e->c - y2)/e->a;};
01069     };
01070         
01071   //printf("\nPushing line (%f,%f,%f,%f)",x1,y1,x2,y2);
01072   //fprintf(stdout, "Line with vertices (%f,%f) and (%f,%f)\n", 
01073   //    e->reg[0]->coord.x, e->reg[1]->coord.x, e->reg[0]->coord.y, e->reg[1]->coord.y);
01074   pushGraphEdge(x1,y1,x2,y2,e->reg[0],e->reg[1]);
01075 }
01076 
01077 
01078 /* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
01079    deltax, deltay (can all be estimates).
01080    Performance suffers if they are wrong; better to make nsites,
01081    deltax, and deltay too big than too small.  (?) */
01082 
01083 bool VoronoiDiagramGenerator::voronoi()
01084 {
01085   Site *newsite, *bot, *top, *temp, *p;
01086   Site *v;
01087   VPoint newintstar;
01088   int pm;
01089   Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
01090   Edge *e;
01091         
01092   PQinitialize();
01093   bottomsite = nextone();
01094   //GS unused plot: out_site(bottomsite);
01095   bool retval = ELinitialize();
01096 
01097   if(!retval)
01098     return false;
01099         
01100   newsite = nextone();
01101   while(1)
01102     {
01103 
01104       if(!PQempty()) 
01105         newintstar = PQ_min();
01106                 
01107       //if the lowest site has a smaller y value than the lowest vector intersection, process the site
01108       //otherwise process the vector intersection               
01109       if (newsite != (Site *)NULL  && (PQempty() || newsite->coord.y < newintstar.y
01110                                        || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
01111         {/* new site is smallest - this is a site event*/
01112           //GS unused plot: out_site(newsite);                                          //output the site
01113           lbnd = ELleftbnd(&(newsite->coord));                          //get the first HalfEdge to the LEFT of the new site
01114           rbnd = ELright(lbnd);                                         //get the first HalfEdge to the RIGHT of the new site
01115           bot = rightreg(lbnd);                                         //if this halfedge has no edge, , bot = bottom site (whatever that is)
01116           e = bisect(bot, newsite);                                     //create a new edge that bisects 
01117           bisector = HEcreate(e, le);                                   //create a new HalfEdge, setting its ELpm field to 0                    
01118           ELinsert(lbnd, bisector);                                     //insert this new bisector edge between the left and right vectors in a linked list     
01119             
01120           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
01121             {   
01122               PQdelete(lbnd);
01123               PQinsert(lbnd, p, dist(p,newsite));
01124             };
01125           lbnd = bisector;                                              
01126           bisector = HEcreate(e, re);                                   //create a new HalfEdge, setting its ELpm field to 1
01127           ELinsert(lbnd, bisector);                                     //insert the new HE to the right of the original bisector earlier in the IF stmt
01128             
01129           if ((p = intersect(bisector, rbnd)) != (Site *) NULL) //if this new bisector intersects with the
01130             {   
01131               PQinsert(bisector, p, dist(p,newsite));                   //push the HE into the ordered linked list of vertices
01132             };
01133           newsite = nextone();  
01134         }
01135       else if (!PQempty()) /* intersection is smallest - this is a vector event */                      
01136         {       
01137           lbnd = PQextractmin();                                                //pop the HalfEdge with the lowest vector off the ordered list of vectors                               
01138           llbnd = ELleft(lbnd);                                         //get the HalfEdge to the left of the above HE
01139           rbnd = ELright(lbnd);                                         //get the HalfEdge to the right of the above HE
01140           rrbnd = ELright(rbnd);                                                //get the HalfEdge to the right of the HE to the right of the lowest HE 
01141           bot = leftreg(lbnd);                                          //get the Site to the left of the left HE which it bisects
01142           top = rightreg(rbnd);                                         //get the Site to the right of the right HE which it bisects
01143             
01144           //GS unused plot: out_triple(bot, top, rightreg(lbnd));               //output the triple of sites, stating that a circle goes through them
01145             
01146           v = lbnd->vertex;                                             //get the vertex that caused this event
01147           makevertex(v);                                                        //set the vertex number - couldn't do this earlier since we didn't know when it would be processed
01148           endpoint(lbnd->ELedge,lbnd->ELpm,v);  //set the endpoint of the left HalfEdge to be this vector
01149           endpoint(rbnd->ELedge,rbnd->ELpm,v);  //set the endpoint of the right HalfEdge to be this vector
01150           ELdelete(lbnd);                                                       //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map  
01151           PQdelete(rbnd);                                                       //remove all vertex events to do with the  right HE
01152           ELdelete(rbnd);                                                       //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map   
01153           pm = le;                                                              //set the pm variable to zero
01154             
01155           if (bot->coord.y > top->coord.y)              //if the site to the left of the event is higher than the Site
01156             {                                                                           //to the right of it, then swap them and set the 'pm' variable to 1
01157               temp = bot; 
01158               bot = top; 
01159               top = temp; 
01160               pm = re;
01161             }
01162           e = bisect(bot, top);                                 //create an Edge (or line) that is between the two Sites. This creates
01163           //the formula of the line, and assigns a line number to it
01164           bisector = HEcreate(e, pm);                           //create a HE from the Edge 'e', and make it point to that edge with its ELedge field
01165           ELinsert(llbnd, bisector);                            //insert the new bisector to the right of the left HE
01166           endpoint(e, re-pm, v);                                        //set one endpoint to the new edge to be the vector point 'v'.
01167           //If the site to the left of this bisector is higher than the right
01168           //Site, then this endpoint is put in position 0; otherwise in pos 1
01169           deref(v);                                                             //delete the vector 'v'
01170             
01171           //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it 
01172           if((p = intersect(llbnd, bisector)) != (Site *) NULL)
01173             {   
01174               PQdelete(llbnd);
01175               PQinsert(llbnd, p, dist(p,bot));
01176             };
01177             
01178           //if right HE and the new bisector don't intersect, then reinsert it 
01179           if ((p = intersect(bisector, rrbnd)) != (Site *) NULL)
01180             {   
01181               PQinsert(bisector, p, dist(p,bot));
01182             };
01183         }
01184       else break;
01185     };
01186 
01187         
01188 
01189 
01190   for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
01191     {   
01192       e = lbnd->ELedge;
01193 
01194       clip_line(e);
01195     };
01196 
01197   //cleanup();
01198 
01199   return true;
01200         
01201 }
01202 
01203 
01204 int scomp(const void *p1,const void *p2)
01205 {
01206   VPoint *s1 = (VPoint*)p1, *s2=(VPoint*)p2;
01207   if(s1->y < s2->y) return(-1);
01208   if(s1->y > s2->y) return(1);
01209   if(s1->x < s2->x) return(-1);
01210   if(s1->x > s2->x) return(1);
01211   return(0);
01212 }
01213 
01214 /* return a single in-storage site */
01215 Site * VoronoiDiagramGenerator::nextone()
01216 {
01217   Site *s;
01218   if(siteidx < nsites)
01219     {   
01220       s = &sites[siteidx];
01221       siteidx += 1;
01222       return(s);
01223     }
01224   else  
01225     return( (Site *)NULL);
01226 }
01227 
01228 FASTJET_END_NAMESPACE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends