Voronoi.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: Voronoi.cc 621 2007-05-09 10:34:30Z 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 #include <stdio.h>
00058 #include "fastjet/internal/Voronoi.hh"
00059 
00060 FASTJET_BEGIN_NAMESPACE
00061 
00062 VoronoiDiagramGenerator::VoronoiDiagramGenerator(){
00063   siteidx = 0;
00064   sites = NULL;
00065   parent_sites = NULL;
00066 
00067   allMemoryList = new FreeNodeArrayList;
00068   allMemoryList->memory = NULL;
00069   allMemoryList->next = NULL;
00070   currentMemoryBlock = allMemoryList;
00071   allEdges = NULL;
00072   iteratorEdges = 0;
00073   minDistanceBetweenSites = 0;
00074 
00075   ELhash = NULL;
00076   PQhash = NULL;
00077 }
00078 
00079 VoronoiDiagramGenerator::~VoronoiDiagramGenerator(){
00080   cleanup();
00081   cleanupEdges();
00082 
00083   if (allMemoryList != NULL)
00084     delete allMemoryList;
00085 }
00086 
00087 
00088 
00089 bool VoronoiDiagramGenerator::generateVoronoi(vector<Point> *_parent_sites,
00090                                               double minX, double maxX, 
00091                                               double minY, double maxY, 
00092                                               double minDist){
00093   cleanup();
00094   cleanupEdges();
00095   int i;
00096   double x, y;
00097 
00098   minDistanceBetweenSites = minDist;
00099 
00100   parent_sites = _parent_sites;
00101 
00102   nsites = n_parent_sites = parent_sites->size();
00103   plot = 0;
00104   triangulate = 0;      
00105   debug = 1;
00106   sorted = 0; 
00107   freeinit(&sfl, sizeof (Site));
00108                 
00109   //sites = (Site *) myalloc(nsites*sizeof( *sites));
00110   sites = (Site *) myalloc(nsites*sizeof( Site));
00111   //parent_sites = (Site *) myalloc(nsites*sizeof( Site));
00112  
00113   if (sites == 0)
00114     return false;
00115 
00116   xmax = xmin = (*parent_sites)[0].x;
00117   ymax = ymin = (*parent_sites)[0].y;
00118   
00119   for(i=0;i<nsites;i++){
00120     x = (*parent_sites)[i].x;
00121     y = (*parent_sites)[i].y;
00122     sites[i].coord.x = x;
00123     sites[i].coord.y = y;
00124     sites[i].sitenbr = i;
00125     sites[i].refcnt  = 0;
00126     
00127     if (x<xmin)
00128       xmin=x;
00129     else if (x>xmax)
00130       xmax=x;
00131     
00132     if (y<ymin)
00133       ymin=y;
00134     else if (y>ymax)
00135       ymax=y;
00136   }
00137         
00138   qsort(sites, nsites, sizeof (*sites), scomp);
00139         
00140   siteidx = 0;
00141   geominit();
00142   double temp = 0;
00143   if(minX > maxX){
00144     temp = minX;
00145     minX = maxX;
00146     maxX = temp;
00147   }
00148   if(minY > maxY){
00149     temp = minY;
00150     minY = maxY;
00151     maxY = temp;
00152   }
00153   borderMinX = minX;
00154   borderMinY = minY;
00155   borderMaxX = maxX;
00156   borderMaxY = maxY;
00157         
00158   siteidx = 0;
00159   voronoi(triangulate);
00160 
00161   return true;
00162 }
00163 
00164 bool VoronoiDiagramGenerator::ELinitialize(){
00165   int i;
00166   freeinit(&hfl, sizeof(Halfedge));
00167   ELhashsize = 2 * sqrt_nsites;
00168   //ELhash = (Halfedge **) myalloc ( sizeof *ELhash * ELhashsize);
00169   ELhash = (Halfedge **) myalloc ( sizeof(Halfedge*)*ELhashsize);
00170   
00171   if(ELhash == 0)
00172     return false;
00173   
00174   for(i=0; i<ELhashsize; i +=1) ELhash[i] = (Halfedge *)NULL;
00175   ELleftend = HEcreate((Edge*) NULL, 0);
00176   ELrightend = HEcreate((Edge*) NULL, 0);
00177   ELleftend->ELleft = (Halfedge*) NULL;
00178   ELleftend->ELright = ELrightend;
00179   ELrightend->ELleft = ELleftend;
00180   ELrightend->ELright = (Halfedge*) NULL;
00181   ELhash[0] = ELleftend;
00182   ELhash[ELhashsize-1] = ELrightend;
00183 
00184   return true;
00185 }
00186 
00187 
00188 Halfedge* VoronoiDiagramGenerator::HEcreate(Edge *e,int pm){
00189   Halfedge *answer;
00190   answer = (Halfedge*) getfree(&hfl);
00191   answer->ELedge = e;
00192   answer->ELpm = pm;
00193   answer->PQnext = (Halfedge*) NULL;
00194   answer->vertex = (Site*) NULL;
00195   answer->ELrefcnt = 0;
00196   return answer;
00197 }
00198 
00199 
00200 void VoronoiDiagramGenerator::ELinsert(Halfedge *lb, Halfedge *newHe)
00201 {
00202   newHe->ELleft = lb;
00203   newHe->ELright = lb->ELright;
00204   (lb->ELright)->ELleft = newHe;
00205   lb->ELright = newHe;
00206 }
00207 
00208 
00209 /* Get entry from hash table, pruning any deleted nodes */
00210 Halfedge* VoronoiDiagramGenerator::ELgethash(int b){
00211   Halfedge *he;
00212         
00213   if ((b<0) || (b>=ELhashsize)) 
00214     return (Halfedge *) NULL;
00215 
00216   he = ELhash[b]; 
00217   if ((he==(Halfedge*) NULL) || (he->ELedge!=(Edge*) DELETED)) 
00218     return he;
00219         
00220   /* Hash table points to deleted half edge.  Patch as necessary. */
00221   ELhash[b] = (Halfedge*) NULL;
00222   if ((he->ELrefcnt -= 1) == 0) 
00223     makefree((Freenode*)he, &hfl);
00224   return (Halfedge*) NULL;
00225 }       
00226 
00227 Halfedge * VoronoiDiagramGenerator::ELleftbnd(Point *p){
00228   int i, bucket;
00229   Halfedge *he;
00230         
00231   /* Use hash table to get close to desired halfedge */
00232   // use the hash function to find the place in the hash map that this
00233   // HalfEdge should be
00234   bucket = (int)((p->x - xmin)/deltax * ELhashsize);    
00235 
00236   // make sure that the bucket position in within the range of the
00237   //hash array
00238   if (bucket<0) bucket =0;
00239   if (bucket>=ELhashsize) bucket = ELhashsize - 1;
00240 
00241   he = ELgethash(bucket);
00242 
00243   // if the HE isn't found, search backwards and forwards in the hash
00244   // map for the first non-null entry
00245   if (he == (Halfedge*) NULL){   
00246     for(i=1;1;i++){     
00247       if ((he=ELgethash(bucket-i)) != (Halfedge*) NULL) 
00248         break;
00249       if ((he=ELgethash(bucket+i)) != (Halfedge*) NULL) 
00250         break;
00251     };
00252     totalsearch += i;
00253   };
00254   ntry += 1;
00255   /* Now search linear list of halfedges for the correct one */
00256   if ((he==ELleftend)  || (he != ELrightend && right_of(he,p))){
00257     do{
00258       he = he->ELright;
00259     } while (he!=ELrightend && right_of(he,p));
00260       // keep going right on the list until either the end is reached,
00261       // or you find the 1st edge which the point
00262       he = he->ELleft;  //isn't to the right of
00263   } else 
00264     // if the point is to the left of the HalfEdge, then search left
00265     // for the HE just to the left of the point
00266     do{
00267       he = he->ELleft;
00268     } while ((he!=ELleftend )&& (!right_of(he,p)));
00269                 
00270   /* Update hash table and reference counts */
00271   if ((bucket > 0) && (bucket <ELhashsize-1)){  
00272     if(ELhash[bucket] != (Halfedge *) NULL){
00273       ELhash[bucket]->ELrefcnt -= 1;
00274     }
00275     ELhash[bucket] = he;
00276     ELhash[bucket]->ELrefcnt += 1;
00277   };
00278 
00279   return he;
00280 }
00281 
00282 
00283 /* This delete routine can't reclaim node, since pointers from hash
00284    table may be present.   */
00285 void VoronoiDiagramGenerator::ELdelete(Halfedge *he){
00286   (he->ELleft)->ELright = he->ELright;
00287   (he->ELright)->ELleft = he->ELleft;
00288   he->ELedge = (Edge*) DELETED;
00289 }
00290 
00291 
00292 Halfedge* VoronoiDiagramGenerator::ELright(Halfedge *he){
00293   return he->ELright;
00294 }
00295 
00296 
00297 Halfedge* VoronoiDiagramGenerator::ELleft(Halfedge *he){
00298   return he->ELleft;
00299 }
00300 
00301 
00302 Site * VoronoiDiagramGenerator::leftreg(Halfedge *he){
00303   if (he->ELedge == (Edge*) NULL) 
00304     return(bottomsite);
00305   return (he->ELpm == le) 
00306     ? he->ELedge->reg[le] 
00307     : he->ELedge->reg[re];
00308 }
00309 
00310 Site * VoronoiDiagramGenerator::rightreg(Halfedge *he){
00311   if (he->ELedge == (Edge*) NULL)
00312     // if this halfedge has no edge, return the bottom site (whatever
00313     // that is)
00314     return bottomsite;
00315   
00316   // if the ELpm field is zero, return the site 0 that this edge
00317   // bisects, otherwise return site number 1
00318   return (he->ELpm == le)
00319     ? he->ELedge->reg[re] 
00320     : he->ELedge->reg[le];
00321 }
00322 
00323 void VoronoiDiagramGenerator::geominit(){
00324   double sn;
00325   
00326   freeinit(&efl, sizeof(Edge));
00327   nvertices = 0;
00328   nedges = 0;
00329   sn = (double)nsites+4;
00330   sqrt_nsites = (int)sqrt(sn);
00331   deltay = ymax - ymin;
00332   deltax = xmax - xmin;
00333 }
00334 
00335 
00336 Edge * VoronoiDiagramGenerator::bisect(Site *s1, Site *s2){
00337   double dx,dy,adx,ady;
00338   Edge *newedge;        
00339 
00340   newedge = (Edge*) getfree(&efl);
00341         
00342   newedge->reg[0] = s1; //store the sites that this edge is bisecting
00343   newedge->reg[1] = s2;
00344   ref(s1); 
00345   ref(s2);
00346 
00347   // to begin with, there are no endpoints on the bisector - it goes
00348   // to infinity
00349   newedge->ep[0] = (Site*) NULL;
00350   newedge->ep[1] = (Site*) NULL;
00351   
00352   // get the difference in x dist between the sites
00353   dx = s2->coord.x - s1->coord.x;
00354   dy = s2->coord.y - s1->coord.y;
00355 
00356   // make sure that the difference is positive
00357   adx = dx>0 ? dx : -dx;
00358   ady = dy>0 ? dy : -dy;
00359 
00360   // get the slope of the line
00361   newedge->c = (double)(s1->coord.x * dx + s1->coord.y * dy
00362                         + (dx*dx + dy*dy)*0.5);
00363 
00364   if (adx>ady){ 
00365     //set formula of line, with x fixed to 1
00366     newedge->a = 1.0; newedge->b = dy/dx; newedge->c /= dx;
00367   } else {      
00368     //set formula of line, with y fixed to 1
00369     newedge->b = 1.0; newedge->a = dx/dy; newedge->c /= dy;
00370   }
00371         
00372   newedge->edgenbr = nedges;
00373   nedges++;
00374 
00375   return newedge;
00376 }
00377 
00378 
00379 // create a new site where the HalfEdges el1 and el2 intersect - note
00380 // that the Point in the argument list is not used, don't know why
00381 // it's there
00382 Site* VoronoiDiagramGenerator::intersect(Halfedge *el1, Halfedge *el2, Point *p){
00383   Edge *e1,*e2, *e;
00384   Halfedge *el;
00385   double d, xint, yint;
00386   int right_of_site;
00387   Site *v;
00388         
00389   e1 = el1->ELedge;
00390   e2 = el2->ELedge;
00391   if ((e1==(Edge*)NULL) || (e2==(Edge*)NULL)) 
00392     return ((Site*) NULL);
00393 
00394   // if the two edges bisect the same parent, return null
00395   if (e1->reg[1] == e2->reg[1]) 
00396     return (Site*) NULL;
00397         
00398   d = e1->a * e2->b - e1->b * e2->a;
00399   if (-1.0e-10<d && d<1.0e-10) 
00400     return (Site*) NULL;
00401         
00402   xint = (e1->c*e2->b - e2->c*e1->b)/d;
00403   yint = (e2->c*e1->a - e1->c*e2->a)/d;
00404         
00405   if( (e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
00406       ((e1->reg[1]->coord.y == e2->reg[1]->coord.y) &&
00407        (e1->reg[1]->coord.x <  e2->reg[1]->coord.x)) ){ 
00408     el = el1; 
00409     e = e1;
00410   } else {      
00411     el = el2; 
00412     e = e2;
00413   }
00414         
00415   right_of_site = xint >= e->reg[1]->coord.x;
00416   if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re)) 
00417     return (Site*) NULL;
00418         
00419   // create a new site at the point of intersection - this is a new
00420   // vector event waiting to happen
00421   v = (Site*) getfree(&sfl);
00422   v->refcnt = 0;
00423   v->coord.x = xint;
00424   v->coord.y = yint;
00425   return v;
00426 }
00427 
00428 //HERE
00429 
00430 /* returns 1 if p is to right of halfedge e */
00431 int VoronoiDiagramGenerator::right_of(Halfedge *el,Point *p)
00432 {
00433   Edge *e;
00434   Site *topsite;
00435   int right_of_site, above, fast;
00436   double dxp, dyp, dxs, t1, t2, t3, yl;
00437         
00438   e = el->ELedge;
00439   topsite = e->reg[1];
00440   right_of_site = p->x > topsite->coord.x;
00441   if(right_of_site && el->ELpm == le) return(1);
00442   if(!right_of_site && el->ELpm == re) return (0);
00443         
00444   if (e->a == 1.0)
00445     {   dyp = p->y - topsite->coord.y;
00446     dxp = p->x - topsite->coord.x;
00447     fast = 0;
00448     if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )
00449       { above = dyp>= e->b*dxp; 
00450       fast = above;
00451       }
00452     else 
00453       { above = p->x + p->y*e->b > e-> c;
00454       if(e->b<0.0) above = !above;
00455       if (!above) fast = 1;
00456       };
00457     if (!fast)
00458       { dxs = topsite->coord.x - (e->reg[0])->coord.x;
00459       above = e->b * (dxp*dxp - dyp*dyp) <
00460         dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
00461       if(e->b<0.0) above = !above;
00462       };
00463     }
00464   else  /*e->b==1.0 */
00465     {   yl = e->c - e->a*p->x;
00466     t1 = p->y - yl;
00467     t2 = p->x - topsite->coord.x;
00468     t3 = yl - topsite->coord.y;
00469     above = t1*t1 > t2*t2 + t3*t3;
00470     };
00471   return (el->ELpm==le ? above : !above);
00472 }
00473 
00474 
00475 void VoronoiDiagramGenerator::endpoint(Edge *e,int lr,Site * s)
00476 {
00477   e->ep[lr] = s;
00478   ref(s);
00479   if(e->ep[re-lr]== (Site *) NULL) 
00480     return;
00481 
00482   clip_line(e);
00483 
00484   deref(e->reg[le]);
00485   deref(e->reg[re]);
00486   makefree((Freenode*)e, &efl);
00487 }
00488 
00489 
00490 double VoronoiDiagramGenerator::dist(Site *s,Site *t)
00491 {
00492   double dx,dy;
00493   dx = s->coord.x - t->coord.x;
00494   dy = s->coord.y - t->coord.y;
00495   return (double)(sqrt(dx*dx + dy*dy));
00496 }
00497 
00498 
00499 void VoronoiDiagramGenerator::makevertex(Site *v)
00500 {
00501   v->sitenbr = nvertices;
00502   nvertices += 1;
00503   out_vertex(v);
00504 }
00505 
00506 
00507 void VoronoiDiagramGenerator::deref(Site *v)
00508 {
00509   v->refcnt -= 1;
00510   if (v->refcnt == 0 ) 
00511     makefree((Freenode*)v, &sfl);
00512 }
00513 
00514 void VoronoiDiagramGenerator::ref(Site *v)
00515 {
00516   v->refcnt += 1;
00517 }
00518 
00519 //push the HalfEdge into the ordered linked list of vertices
00520 void VoronoiDiagramGenerator::PQinsert(Halfedge *he,Site * v, double offset)
00521 {
00522   Halfedge *last, *next;
00523         
00524   he->vertex = v;
00525   ref(v);
00526   he->ystar = (double)(v->coord.y + offset);
00527   last = &PQhash[PQbucket(he)];
00528   while ((next = last->PQnext) != (Halfedge *) NULL &&
00529          (he->ystar  > next->ystar  ||
00530           (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x)))
00531     {   
00532       last = next;
00533     };
00534   he->PQnext = last->PQnext; 
00535   last->PQnext = he;
00536   PQcount += 1;
00537 }
00538 
00539 //remove the HalfEdge from the list of vertices 
00540 void VoronoiDiagramGenerator::PQdelete(Halfedge *he)
00541 {
00542   Halfedge *last;
00543         
00544   if(he->vertex != (Site *) NULL)
00545     {   
00546       last = &PQhash[PQbucket(he)];
00547       while (last->PQnext != he) 
00548         last = last->PQnext;
00549 
00550       last->PQnext = he->PQnext;
00551       PQcount -= 1;
00552       deref(he->vertex);
00553       he->vertex = (Site *) NULL;
00554     };
00555 }
00556 
00557 int VoronoiDiagramGenerator::PQbucket(Halfedge *he)
00558 {
00559   int bucket;
00560         
00561   bucket = (int)((he->ystar - ymin)/deltay * PQhashsize);
00562   if (bucket<0) bucket = 0;
00563   if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
00564   if (bucket < PQmin) PQmin = bucket;
00565   return(bucket);
00566 }
00567 
00568 
00569 
00570 int VoronoiDiagramGenerator::PQempty()
00571 {
00572   return(PQcount==0);
00573 }
00574 
00575 
00576 Point VoronoiDiagramGenerator::PQ_min()
00577 {
00578   Point answer;
00579         
00580   while(PQhash[PQmin].PQnext == (Halfedge *)NULL) {PQmin += 1;};
00581   answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
00582   answer.y = PQhash[PQmin].PQnext->ystar;
00583   return (answer);
00584 }
00585 
00586 Halfedge * VoronoiDiagramGenerator::PQextractmin()
00587 {
00588   Halfedge *curr;
00589         
00590   curr = PQhash[PQmin].PQnext;
00591   PQhash[PQmin].PQnext = curr->PQnext;
00592   PQcount -= 1;
00593   return(curr);
00594 }
00595 
00596 
00597 bool VoronoiDiagramGenerator::PQinitialize()
00598 {
00599   int i; 
00600         
00601   PQcount = 0;
00602   PQmin = 0;
00603   PQhashsize = 4 * sqrt_nsites;
00604   //PQhash = (Halfedge *) myalloc(PQhashsize * sizeof *PQhash);
00605   PQhash = (Halfedge *) myalloc(PQhashsize * sizeof(Halfedge));
00606 
00607   if(PQhash == 0)
00608     return false;
00609 
00610   for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (Halfedge *)NULL;
00611 
00612   return true;
00613 }
00614 
00615 
00616 void VoronoiDiagramGenerator::freeinit(Freelist *fl,int size)
00617 {
00618   fl->head = (Freenode *) NULL;
00619   fl->nodesize = size;
00620 }
00621 
00622 char * VoronoiDiagramGenerator::getfree(Freelist *fl)
00623 {
00624   int i; 
00625   Freenode *t;
00626 
00627   if(fl->head == (Freenode *) NULL)
00628     {   
00629       t =  (Freenode *) myalloc(sqrt_nsites * fl->nodesize);
00630 
00631       if(t == 0)
00632         return 0;
00633                 
00634       currentMemoryBlock->next = new FreeNodeArrayList;
00635       currentMemoryBlock = currentMemoryBlock->next;
00636       currentMemoryBlock->memory = t;
00637       currentMemoryBlock->next = 0;
00638 
00639       for(i=0; i<sqrt_nsites; i+=1)     
00640         makefree((Freenode *)((char *)t+i*fl->nodesize), fl);           
00641     };
00642   t = fl->head;
00643   fl->head = (fl->head)->nextfree;
00644   return((char *)t);
00645 }
00646 
00647 
00648 
00649 void VoronoiDiagramGenerator::makefree(Freenode *curr,Freelist *fl)
00650 {
00651   curr->nextfree = fl->head;
00652   fl->head = curr;
00653 }
00654 
00655 void VoronoiDiagramGenerator::cleanup()
00656 {
00657   if(sites != NULL){
00658     free(sites);
00659     sites = 0;
00660   }
00661 
00662   FreeNodeArrayList* current=NULL, *prev=NULL;
00663 
00664   current = prev = allMemoryList;
00665 
00666   while(current->next!=NULL){
00667     prev = current;
00668     current = current->next;
00669     free(prev->memory);
00670     delete prev;
00671     prev = NULL;
00672   }
00673 
00674   if(current!=NULL){
00675     if (current->memory!=NULL){
00676       free(current->memory);
00677     }
00678     delete current;
00679   }
00680 
00681   allMemoryList = new FreeNodeArrayList;
00682   allMemoryList->next = NULL;
00683   allMemoryList->memory = NULL;
00684   currentMemoryBlock = allMemoryList;
00685 
00686   if (ELhash!=NULL)
00687     free(ELhash);
00688 
00689   if (PQhash!=NULL)
00690     free(PQhash);
00691 }
00692 
00693 void VoronoiDiagramGenerator::cleanupEdges()
00694 {
00695   GraphEdge* geCurrent = NULL, *gePrev = NULL;
00696   geCurrent = gePrev = allEdges;
00697 
00698   while(geCurrent!=NULL){
00699     gePrev = geCurrent;
00700     geCurrent = geCurrent->next;
00701     delete gePrev;
00702   }
00703 
00704   allEdges = 0;
00705 
00706 }
00707 
00708 void VoronoiDiagramGenerator::pushGraphEdge(double x1, double y1, double x2, double y2,
00709                                             Site *s1, Site *s2)
00710 {
00711   GraphEdge* newEdge = new GraphEdge;
00712   newEdge->next = allEdges;
00713   allEdges = newEdge;
00714   newEdge->x1 = x1;
00715   newEdge->y1 = y1;
00716   newEdge->x2 = x2;
00717   newEdge->y2 = y2;
00718 
00719   newEdge->point1 = s1->sitenbr;
00720   newEdge->point2 = s2->sitenbr;
00721 }
00722 
00723 
00724 char * VoronoiDiagramGenerator::myalloc(unsigned n)
00725 {
00726   char *t=0;    
00727   t=(char*)malloc(n);
00728   total_alloc += n;
00729   return(t);
00730 }
00731 
00732 
00733 /* for those who don't have Cherry's plot */
00734 /* #include <plot.h> */
00735 void VoronoiDiagramGenerator::openpl(){}
00736 void VoronoiDiagramGenerator::circle(double x, double y, double radius){}
00737 void VoronoiDiagramGenerator::range(double minX, double minY, double maxX, double maxY){}
00738 
00739 
00740 
00741 void VoronoiDiagramGenerator::out_bisector(Edge *e)
00742 {
00743         
00744 
00745 }
00746 
00747 
00748 void VoronoiDiagramGenerator::out_ep(Edge *e)
00749 {
00750         
00751         
00752 }
00753 
00754 void VoronoiDiagramGenerator::out_vertex(Site *v)
00755 {
00756         
00757 }
00758 
00759 
00760 void VoronoiDiagramGenerator::out_site(Site *s)
00761 {
00762   if(!triangulate & plot & !debug)
00763     circle (s->coord.x, s->coord.y, cradius);
00764         
00765 }
00766 
00767 
00768 void VoronoiDiagramGenerator::out_triple(Site *s1, Site *s2,Site * s3)
00769 {
00770         
00771 }
00772 
00773 
00774 
00775 void VoronoiDiagramGenerator::plotinit()
00776 {
00777   double dx,dy,d;
00778         
00779   dy = ymax - ymin;
00780   dx = xmax - xmin;
00781   d = (double)(( dx > dy ? dx : dy) * 1.1);
00782   pxmin = (double)(xmin - (d-dx)/2.0);
00783   pxmax = (double)(xmax + (d-dx)/2.0);
00784   pymin = (double)(ymin - (d-dy)/2.0);
00785   pymax = (double)(ymax + (d-dy)/2.0);
00786   cradius = (double)((pxmax - pxmin)/350.0);
00787   openpl();
00788   range(pxmin, pymin, pxmax, pymax);
00789 }
00790 
00791 
00792 void VoronoiDiagramGenerator::clip_line(Edge *e)
00793 {
00794   Site *s1, *s2;
00795   double x1=0,x2=0,y1=0,y2=0; //, temp = 0;
00796 
00797   x1 = e->reg[0]->coord.x;
00798   x2 = e->reg[1]->coord.x;
00799   y1 = e->reg[0]->coord.y;
00800   y2 = e->reg[1]->coord.y;
00801 
00802   //if the distance between the two points this line was created from is less than 
00803   //the square root of 2, then ignore it
00804   //TODO improve/remove
00805   //if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites)
00806   //  {
00807   //    return;
00808   //  }
00809   pxmin = borderMinX;
00810   pxmax = borderMaxX;
00811   pymin = borderMinY;
00812   pymax = borderMaxY;
00813 
00814   if(e->a == 1.0 && e ->b >= 0.0)
00815     {   
00816       s1 = e->ep[1];
00817       s2 = e->ep[0];
00818     }
00819   else 
00820     {
00821       s1 = e->ep[0];
00822       s2 = e->ep[1];
00823     };
00824         
00825   if(e->a == 1.0)
00826     {
00827       y1 = pymin;
00828       if (s1!=(Site *)NULL && s1->coord.y > pymin)
00829         {
00830           y1 = s1->coord.y;
00831         }
00832       if(y1>pymax) 
00833         {
00834           //    printf("\nClipped (1) y1 = %f to %f",y1,pymax);
00835           y1 = pymax;
00836           //return;
00837         }
00838       x1 = e->c - e->b * y1;
00839       y2 = pymax;
00840       if (s2!=(Site *)NULL && s2->coord.y < pymax) 
00841         y2 = s2->coord.y;
00842 
00843       if(y2<pymin) 
00844         {
00845           //printf("\nClipped (2) y2 = %f to %f",y2,pymin);
00846           y2 = pymin;
00847           //return;
00848         }
00849       x2 = (e->c) - (e->b) * y2;
00850       if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin))) 
00851         {
00852           //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);
00853           return;
00854         }
00855       if(x1> pxmax)
00856         {       x1 = pxmax; y1 = (e->c - x1)/e->b;};
00857       if(x1<pxmin)
00858         {       x1 = pxmin; y1 = (e->c - x1)/e->b;};
00859       if(x2>pxmax)
00860         {       x2 = pxmax; y2 = (e->c - x2)/e->b;};
00861       if(x2<pxmin)
00862         {       x2 = pxmin; y2 = (e->c - x2)/e->b;};
00863     }
00864   else
00865     {
00866       x1 = pxmin;
00867       if (s1!=(Site *)NULL && s1->coord.x > pxmin) 
00868         x1 = s1->coord.x;
00869       if(x1>pxmax) 
00870         {
00871           //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);
00872           //return;
00873           x1 = pxmax;
00874         }
00875       y1 = e->c - e->a * x1;
00876       x2 = pxmax;
00877       if (s2!=(Site *)NULL && s2->coord.x < pxmax) 
00878         x2 = s2->coord.x;
00879       if(x2<pxmin) 
00880         {
00881           //printf("\nClipped (4) x2 = %f to %f",x2,pxmin);
00882           //return;
00883           x2 = pxmin;
00884         }
00885       y2 = e->c - e->a * x2;
00886       if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin))) 
00887         {
00888           //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);
00889           return;
00890         }
00891       if(y1> pymax)
00892         {       y1 = pymax; x1 = (e->c - y1)/e->a;};
00893       if(y1<pymin)
00894         {       y1 = pymin; x1 = (e->c - y1)/e->a;};
00895       if(y2>pymax)
00896         {       y2 = pymax; x2 = (e->c - y2)/e->a;};
00897       if(y2<pymin)
00898         {       y2 = pymin; x2 = (e->c - y2)/e->a;};
00899     };
00900         
00901   //printf("\nPushing line (%f,%f,%f,%f)",x1,y1,x2,y2);
00902   //fprintf(stdout, "Line with vertices (%f,%f) and (%f,%f)\n", 
00903   //    e->reg[0]->coord.x, e->reg[1]->coord.x, e->reg[0]->coord.y, e->reg[1]->coord.y);
00904   pushGraphEdge(x1,y1,x2,y2,e->reg[0],e->reg[1]);
00905 }
00906 
00907 
00908 /* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
00909    deltax, deltay (can all be estimates).
00910    Performance suffers if they are wrong; better to make nsites,
00911    deltax, and deltay too big than too small.  (?) */
00912 
00913 bool VoronoiDiagramGenerator::voronoi(int triangulate)
00914 {
00915   Site *newsite, *bot, *top, *temp, *p;
00916   Site *v;
00917   Point newintstar;
00918   int pm;
00919   Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
00920   Edge *e;
00921         
00922   PQinitialize();
00923   bottomsite = nextone();
00924   out_site(bottomsite);
00925   bool retval = ELinitialize();
00926 
00927   if(!retval)
00928     return false;
00929         
00930   newsite = nextone();
00931   while(1)
00932     {
00933 
00934       if(!PQempty()) 
00935         newintstar = PQ_min();
00936                 
00937       //if the lowest site has a smaller y value than the lowest vector intersection, process the site
00938       //otherwise process the vector intersection               
00939 
00940       if (newsite != (Site *)NULL       && (PQempty() || newsite->coord.y < newintstar.y
00941                                                     || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
00942         {/* new site is smallest - this is a site event*/
00943           out_site(newsite);                                            //output the site
00944           lbnd = ELleftbnd(&(newsite->coord));                          //get the first HalfEdge to the LEFT of the new site
00945           rbnd = ELright(lbnd);                                         //get the first HalfEdge to the RIGHT of the new site
00946           bot = rightreg(lbnd);                                         //if this halfedge has no edge, , bot = bottom site (whatever that is)
00947           e = bisect(bot, newsite);                                     //create a new edge that bisects 
00948           bisector = HEcreate(e, le);                                   //create a new HalfEdge, setting its ELpm field to 0                    
00949           ELinsert(lbnd, bisector);                                     //insert this new bisector edge between the left and right vectors in a linked list     
00950 
00951           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
00952             {   
00953               PQdelete(lbnd);
00954               PQinsert(lbnd, p, dist(p,newsite));
00955             };
00956           lbnd = bisector;                                              
00957           bisector = HEcreate(e, re);                                   //create a new HalfEdge, setting its ELpm field to 1
00958           ELinsert(lbnd, bisector);                                     //insert the new HE to the right of the original bisector earlier in the IF stmt
00959 
00960           if ((p = intersect(bisector, rbnd)) != (Site *) NULL) //if this new bisector intersects with the
00961             {   
00962               PQinsert(bisector, p, dist(p,newsite));                   //push the HE into the ordered linked list of vertices
00963             };
00964           newsite = nextone();  
00965         }
00966       else if (!PQempty()) /* intersection is smallest - this is a vector event */                      
00967         {       
00968           lbnd = PQextractmin();                                                //pop the HalfEdge with the lowest vector off the ordered list of vectors                               
00969           llbnd = ELleft(lbnd);                                         //get the HalfEdge to the left of the above HE
00970           rbnd = ELright(lbnd);                                         //get the HalfEdge to the right of the above HE
00971           rrbnd = ELright(rbnd);                                                //get the HalfEdge to the right of the HE to the right of the lowest HE 
00972           bot = leftreg(lbnd);                                          //get the Site to the left of the left HE which it bisects
00973           top = rightreg(rbnd);                                         //get the Site to the right of the right HE which it bisects
00974 
00975           out_triple(bot, top, rightreg(lbnd));         //output the triple of sites, stating that a circle goes through them
00976 
00977           v = lbnd->vertex;                                             //get the vertex that caused this event
00978           makevertex(v);                                                        //set the vertex number - couldn't do this earlier since we didn't know when it would be processed
00979           endpoint(lbnd->ELedge,lbnd->ELpm,v);  //set the endpoint of the left HalfEdge to be this vector
00980           endpoint(rbnd->ELedge,rbnd->ELpm,v);  //set the endpoint of the right HalfEdge to be this vector
00981           ELdelete(lbnd);                                                       //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map  
00982           PQdelete(rbnd);                                                       //remove all vertex events to do with the  right HE
00983           ELdelete(rbnd);                                                       //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map   
00984           pm = le;                                                              //set the pm variable to zero
00985                         
00986           if (bot->coord.y > top->coord.y)              //if the site to the left of the event is higher than the Site
00987             {                                                                           //to the right of it, then swap them and set the 'pm' variable to 1
00988               temp = bot; 
00989               bot = top; 
00990               top = temp; 
00991               pm = re;
00992             }
00993           e = bisect(bot, top);                                 //create an Edge (or line) that is between the two Sites. This creates
00994           //the formula of the line, and assigns a line number to it
00995           bisector = HEcreate(e, pm);                           //create a HE from the Edge 'e', and make it point to that edge with its ELedge field
00996           ELinsert(llbnd, bisector);                            //insert the new bisector to the right of the left HE
00997           endpoint(e, re-pm, v);                                        //set one endpoint to the new edge to be the vector point 'v'.
00998           //If the site to the left of this bisector is higher than the right
00999           //Site, then this endpoint is put in position 0; otherwise in pos 1
01000           deref(v);                                                             //delete the vector 'v'
01001 
01002           //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it 
01003           if((p = intersect(llbnd, bisector)) != (Site *) NULL)
01004             {   
01005               PQdelete(llbnd);
01006               PQinsert(llbnd, p, dist(p,bot));
01007             };
01008 
01009           //if right HE and the new bisector don't intersect, then reinsert it 
01010           if ((p = intersect(bisector, rrbnd)) != (Site *) NULL)
01011             {   
01012               PQinsert(bisector, p, dist(p,bot));
01013             };
01014         }
01015       else break;
01016     };
01017 
01018         
01019 
01020 
01021   for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
01022     {   
01023       e = lbnd->ELedge;
01024 
01025       clip_line(e);
01026     };
01027 
01028   //cleanup();
01029 
01030   return true;
01031         
01032 }
01033 
01034 
01035 int scomp(const void *p1,const void *p2)
01036 {
01037   Point *s1 = (Point*)p1, *s2=(Point*)p2;
01038   if(s1->y < s2->y) return(-1);
01039   if(s1->y > s2->y) return(1);
01040   if(s1->x < s2->x) return(-1);
01041   if(s1->x > s2->x) return(1);
01042   return(0);
01043 }
01044 
01045 /* return a single in-storage site */
01046 Site * VoronoiDiagramGenerator::nextone()
01047 {
01048   Site *s;
01049   if(siteidx < nsites)
01050     {   
01051       s = &sites[siteidx];
01052       siteidx += 1;
01053       return(s);
01054     }
01055   else  
01056     return( (Site *)NULL);
01057 }
01058 
01059 FASTJET_END_NAMESPACE

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2