FastJet 3.0.2
|
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