109 #include "fastjet/internal/Voronoi.hh" 113 FASTJET_BEGIN_NAMESPACE
115 LimitedWarning VoronoiDiagramGenerator::_warning_degeneracy;
117 VoronoiDiagramGenerator::VoronoiDiagramGenerator(){
122 allMemoryList =
new FreeNodeArrayList;
123 allMemoryList->memory = NULL;
124 allMemoryList->next = NULL;
125 currentMemoryBlock = allMemoryList;
128 minDistanceBetweenSites = 0;
134 VoronoiDiagramGenerator::~VoronoiDiagramGenerator(){
138 if (allMemoryList != NULL)
139 delete allMemoryList;
144 bool VoronoiDiagramGenerator::generateVoronoi(vector<VPoint> *_parent_sites,
145 double minX,
double maxX,
146 double minY,
double maxY,
153 minDistanceBetweenSites = minDist;
155 parent_sites = _parent_sites;
157 nsites = n_parent_sites = parent_sites->size();
160 freeinit(&sfl,
sizeof (Site));
163 sites = (Site *) myalloc(nsites*
sizeof( Site));
169 xmax = xmin = (*parent_sites)[0].x;
170 ymax = ymin = (*parent_sites)[0].y;
172 for(i=0;i<nsites;i++){
173 x = (*parent_sites)[i].x;
174 y = (*parent_sites)[i].y;
175 sites[i].coord.x = x;
176 sites[i].coord.y = y;
177 sites[i].sitenbr = i;
191 qsort(sites, nsites,
sizeof (*sites), scomp);
200 unsigned int offset=0;
201 for (
int is=1;is<nsites;is++){
202 if (sites[is].coord.y==sites[is-1].coord.y && sites[is].coord.x==sites[is-1].coord.x){
204 }
else if (offset>0){
205 sites[is-offset] = sites[is];
211 _warning_degeneracy.warn(
"VoronoiDiagramGenerator: two (or more) particles are degenerate in rapidity and azimuth, Voronoi cell assigned to the first of each set of degenerate particles.");
238 bool VoronoiDiagramGenerator::ELinitialize(){
240 freeinit(&hfl,
sizeof(Halfedge));
241 ELhashsize = 2 * sqrt_nsites;
243 ELhash = (Halfedge **) myalloc (
sizeof(Halfedge*)*ELhashsize);
248 for(i=0; i<ELhashsize; i +=1) ELhash[i] = (Halfedge *)NULL;
249 ELleftend = HEcreate((Edge*) NULL, 0);
250 ELrightend = HEcreate((Edge*) NULL, 0);
251 ELleftend->ELleft = (Halfedge*) NULL;
252 ELleftend->ELright = ELrightend;
253 ELrightend->ELleft = ELleftend;
254 ELrightend->ELright = (Halfedge*) NULL;
255 ELhash[0] = ELleftend;
256 ELhash[ELhashsize-1] = ELrightend;
262 Halfedge* VoronoiDiagramGenerator::HEcreate(Edge *e,
int pm){
264 answer = (Halfedge*) getfree(&hfl);
267 answer->PQnext = (Halfedge*) NULL;
268 answer->vertex = (Site*) NULL;
269 answer->ELrefcnt = 0;
274 void VoronoiDiagramGenerator::ELinsert(Halfedge *lb, Halfedge *newHe)
277 newHe->ELright = lb->ELright;
278 (lb->ELright)->ELleft = newHe;
284 Halfedge* VoronoiDiagramGenerator::ELgethash(
int b){
287 if ((b<0) || (b>=ELhashsize))
288 return (Halfedge *) NULL;
291 if ((he==(Halfedge*) NULL) || (he->ELedge!=(Edge*) DELETED))
295 ELhash[b] = (Halfedge*) NULL;
296 if ((he->ELrefcnt -= 1) == 0)
297 makefree((Freenode*)he, &hfl);
298 return (Halfedge*) NULL;
301 Halfedge * VoronoiDiagramGenerator::ELleftbnd(VPoint *p){
324 if (p->x < xmin){ bucket=0;}
325 else if (p->x >= xmax){ bucket = ELhashsize - 1;}
327 bucket= (int)((p->x - xmin)/deltax * ELhashsize);
328 if (bucket>=ELhashsize) bucket = ELhashsize - 1;
331 he = ELgethash(bucket);
335 if (he == (Halfedge*) NULL){
337 if ((he=ELgethash(bucket-i)) != (Halfedge*) NULL)
339 if ((he=ELgethash(bucket+i)) != (Halfedge*) NULL)
346 if ((he==ELleftend) || (he != ELrightend && right_of(he,p))){
349 }
while (he!=ELrightend && right_of(he,p));
358 }
while ((he!=ELleftend )&& (!right_of(he,p)));
361 if ((bucket > 0) && (bucket <ELhashsize-1)){
362 if(ELhash[bucket] != (Halfedge *) NULL){
363 ELhash[bucket]->ELrefcnt -= 1;
366 ELhash[bucket]->ELrefcnt += 1;
375 void VoronoiDiagramGenerator::ELdelete(Halfedge *he){
376 (he->ELleft)->ELright = he->ELright;
377 (he->ELright)->ELleft = he->ELleft;
378 he->ELedge = (Edge*) DELETED;
382 Halfedge* VoronoiDiagramGenerator::ELright(Halfedge *he){
387 Halfedge* VoronoiDiagramGenerator::ELleft(Halfedge *he){
392 Site * VoronoiDiagramGenerator::leftreg(Halfedge *he){
393 if (he->ELedge == (Edge*) NULL)
395 return (he->ELpm == le)
396 ? he->ELedge->reg[le]
397 : he->ELedge->reg[re];
400 Site * VoronoiDiagramGenerator::rightreg(Halfedge *he){
401 if (he->ELedge == (Edge*) NULL)
408 return (he->ELpm == le)
409 ? he->ELedge->reg[re]
410 : he->ELedge->reg[le];
413 void VoronoiDiagramGenerator::geominit(){
416 freeinit(&efl,
sizeof(Edge));
419 sn = (double)nsites+4;
420 sqrt_nsites = (int)sqrt(sn);
421 deltay = ymax - ymin;
422 deltax = xmax - xmin;
426 Edge * VoronoiDiagramGenerator::bisect(Site *s1, Site *s2){
427 double dx,dy,adx,ady;
430 newedge = (Edge*) getfree(&efl);
432 newedge->reg[0] = s1;
433 newedge->reg[1] = s2;
439 newedge->ep[0] = (Site*) NULL;
440 newedge->ep[1] = (Site*) NULL;
443 dx = s2->coord.x - s1->coord.x;
444 dy = s2->coord.y - s1->coord.y;
447 adx = dx>0 ? dx : -dx;
448 ady = dy>0 ? dy : -dy;
451 newedge->c = (double)(s1->coord.x * dx + s1->coord.y * dy
452 + (dx*dx + dy*dy)*0.5);
456 newedge->a = 1.0; newedge->b = dy/dx; newedge->c /= dx;
459 newedge->b = 1.0; newedge->a = dx/dy; newedge->c /= dy;
462 newedge->edgenbr = nedges;
474 Site* VoronoiDiagramGenerator::intersect(Halfedge *el1, Halfedge *el2
478 double d, xint, yint;
484 if ((e1==(Edge*)NULL) || (e2==(Edge*)NULL))
485 return ((Site*) NULL);
488 if (e1->reg[1] == e2->reg[1])
506 double dx = e2->reg[1]->coord.x - e1->reg[1]->coord.x;
507 double dy = e2->reg[1]->coord.y - e1->reg[1]->coord.y;
508 double dxref = e1->reg[1]->coord.x - e1->reg[0]->coord.x;
509 double dyref = e1->reg[1]->coord.y - e1->reg[0]->coord.y;
511 if (dx*dx + dy*dy < 1e-14*(dxref*dxref+dyref*dyref)){
513 double adx = dx>0 ? dx : -dx;
514 double ady = dy>0 ? dy : -dy;
518 double c = (double)(e1->reg[1]->coord.x * dx + e1->reg[1]->coord.y * dy
519 + (dx*dx + dy*dy)*0.5);
522 a = 1.0; b = dy/dx; c /= dx;
524 b = 1.0; a = dx/dy; c /= dy;
527 d = e1->a * b - e1->b * a;
528 if (-1.0e-10<d && d<1.0e-10) {
532 xint = (e1->c*b - c*e1->b)/d;
533 yint = (c*e1->a - e1->c*a)/d;
536 d = e1->a * e2->b - e1->b * e2->a;
537 if (-1.0e-10<d && d<1.0e-10) {
541 xint = (e1->c*e2->b - e2->c*e1->b)/d;
542 yint = (e2->c*e1->a - e1->c*e2->a)/d;
546 volatile double local_y1 = e1->reg[1]->coord.y;
547 volatile double local_y2 = e2->reg[1]->coord.y;
549 if( (local_y1 < local_y2) ||
550 ((local_y1 == local_y2) &&
551 (e1->reg[1]->coord.x < e2->reg[1]->coord.x)) ){
559 right_of_site = xint >= e->reg[1]->coord.x;
560 if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re))
565 v = (Site*) getfree(&sfl);
575 int VoronoiDiagramGenerator::right_of(Halfedge *el,VPoint *p)
579 int right_of_site, above, fast;
580 double dxp, dyp, dxs, t1, t2, t3, yl;
584 right_of_site = p->x > topsite->coord.x;
585 if(right_of_site && el->ELpm == le)
return(1);
586 if(!right_of_site && el->ELpm == re)
return (0);
589 { dyp = p->y - topsite->coord.y;
590 dxp = p->x - topsite->coord.x;
592 if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )
593 { above = dyp>= e->b*dxp;
597 { above = p->x + p->y*e->b > e-> c;
598 if(e->b<0.0) above = !above;
599 if (!above) fast = 1;
602 { dxs = topsite->coord.x - (e->reg[0])->coord.x;
603 above = e->b * (dxp*dxp - dyp*dyp) <
604 dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
605 if(e->b<0.0) above = !above;
609 { yl = e->c - e->a*p->x;
611 t2 = p->x - topsite->coord.x;
612 t3 = yl - topsite->coord.y;
613 above = t1*t1 > t2*t2 + t3*t3;
615 return (el->ELpm==le ? above : !above);
619 void VoronoiDiagramGenerator::endpoint(Edge *e,
int lr,Site * s)
623 if(e->ep[re-lr]== (Site *) NULL)
630 makefree((Freenode*)e, &efl);
634 double VoronoiDiagramGenerator::dist(Site *s,Site *t)
637 dx = s->coord.x - t->coord.x;
638 dy = s->coord.y - t->coord.y;
639 return (
double)(sqrt(dx*dx + dy*dy));
643 void VoronoiDiagramGenerator::makevertex(Site *v)
645 v->sitenbr = nvertices;
651 void VoronoiDiagramGenerator::deref(Site *v)
655 makefree((Freenode*)v, &sfl);
658 void VoronoiDiagramGenerator::ref(Site *v)
664 void VoronoiDiagramGenerator::PQinsert(Halfedge *he,Site * v,
double offset)
666 Halfedge *last, *next;
670 he->ystar = (double)(v->coord.y + offset);
671 last = &PQhash[PQbucket(he)];
672 while ((next = last->PQnext) != (Halfedge *) NULL &&
673 (he->ystar > next->ystar ||
674 (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x)))
678 he->PQnext = last->PQnext;
684 void VoronoiDiagramGenerator::PQdelete(Halfedge *he)
688 if(he->vertex != (Site *) NULL)
690 last = &PQhash[PQbucket(he)];
691 while (last->PQnext != he)
694 last->PQnext = he->PQnext;
697 he->vertex = (Site *) NULL;
701 int VoronoiDiagramGenerator::PQbucket(Halfedge *he)
720 double hey = he->ystar;
721 if (hey < ymin){ bucket = 0;}
722 else if (hey >= ymax){ bucket = PQhashsize-1;}
724 bucket = (int)((hey - ymin)/deltay * PQhashsize);
725 if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
728 if (bucket < PQmin) PQmin = bucket;
734 int VoronoiDiagramGenerator::PQempty()
740 VPoint VoronoiDiagramGenerator::PQ_min()
744 while(PQhash[PQmin].PQnext == (Halfedge *)NULL) {PQmin += 1;};
745 answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
746 answer.y = PQhash[PQmin].PQnext->ystar;
750 Halfedge * VoronoiDiagramGenerator::PQextractmin()
754 curr = PQhash[PQmin].PQnext;
755 PQhash[PQmin].PQnext = curr->PQnext;
761 bool VoronoiDiagramGenerator::PQinitialize()
767 PQhashsize = 4 * sqrt_nsites;
769 PQhash = (Halfedge *) myalloc(PQhashsize *
sizeof(Halfedge));
774 for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (Halfedge *)NULL;
780 void VoronoiDiagramGenerator::freeinit(Freelist *fl,
int size)
782 fl->head = (Freenode *) NULL;
786 char * VoronoiDiagramGenerator::getfree(Freelist *fl)
791 if(fl->head == (Freenode *) NULL)
793 t = (Freenode *) myalloc(sqrt_nsites * fl->nodesize);
798 currentMemoryBlock->next =
new FreeNodeArrayList;
799 currentMemoryBlock = currentMemoryBlock->next;
800 currentMemoryBlock->memory = t;
801 currentMemoryBlock->next = 0;
803 for(i=0; i<sqrt_nsites; i+=1)
804 makefree((Freenode *)((
char *)t+i*fl->nodesize), fl);
807 fl->head = (fl->head)->nextfree;
813 void VoronoiDiagramGenerator::makefree(Freenode *curr,Freelist *fl)
815 curr->nextfree = fl->head;
819 void VoronoiDiagramGenerator::cleanup()
826 FreeNodeArrayList* current=NULL, *prev=NULL;
828 current = prev = allMemoryList;
830 while(current->next!=NULL){
832 current = current->next;
839 if (current->memory!=NULL){
840 free(current->memory);
845 allMemoryList =
new FreeNodeArrayList;
846 allMemoryList->next = NULL;
847 allMemoryList->memory = NULL;
848 currentMemoryBlock = allMemoryList;
857 void VoronoiDiagramGenerator::cleanupEdges()
859 GraphEdge* geCurrent = NULL, *gePrev = NULL;
860 geCurrent = gePrev = allEdges;
862 while(geCurrent!=NULL){
864 geCurrent = geCurrent->next;
872 void VoronoiDiagramGenerator::pushGraphEdge(
double x1,
double y1,
double x2,
double y2,
875 GraphEdge* newEdge =
new GraphEdge;
876 newEdge->next = allEdges;
883 newEdge->point1 = s1->sitenbr;
884 newEdge->point2 = s2->sitenbr;
888 char * VoronoiDiagramGenerator::myalloc(
unsigned n)
945 void VoronoiDiagramGenerator::plotinit()
951 d = (double)(( dx > dy ? dx : dy) * 1.1);
952 pxmin = (double)(xmin - (d-dx)/2.0);
953 pxmax = (double)(xmax + (d-dx)/2.0);
954 pymin = (double)(ymin - (d-dy)/2.0);
955 pymax = (double)(ymax + (d-dy)/2.0);
956 cradius = (double)((pxmax - pxmin)/350.0);
962 void VoronoiDiagramGenerator::clip_line(Edge *e)
965 double x1=0,x2=0,y1=0,y2=0;
967 x1 = e->reg[0]->coord.x;
968 x2 = e->reg[1]->coord.x;
969 y1 = e->reg[0]->coord.y;
970 y2 = e->reg[1]->coord.y;
984 if(e->a == 1.0 && e ->b >= 0.0)
998 if (s1!=(Site *)NULL && s1->coord.y > pymin)
1008 x1 = e->c - e->b * y1;
1010 if (s2!=(Site *)NULL && s2->coord.y < pymax)
1019 x2 = (e->c) - (e->b) * y2;
1020 if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin)))
1026 { x1 = pxmax; y1 = (e->c - x1)/e->b;};
1028 { x1 = pxmin; y1 = (e->c - x1)/e->b;};
1030 { x2 = pxmax; y2 = (e->c - x2)/e->b;};
1032 { x2 = pxmin; y2 = (e->c - x2)/e->b;};
1037 if (s1!=(Site *)NULL && s1->coord.x > pxmin)
1045 y1 = e->c - e->a * x1;
1047 if (s2!=(Site *)NULL && s2->coord.x < pxmax)
1055 y2 = e->c - e->a * x2;
1056 if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin)))
1062 { y1 = pymax; x1 = (e->c - y1)/e->a;};
1064 { y1 = pymin; x1 = (e->c - y1)/e->a;};
1066 { y2 = pymax; x2 = (e->c - y2)/e->a;};
1068 { y2 = pymin; x2 = (e->c - y2)/e->a;};
1074 pushGraphEdge(x1,y1,x2,y2,e->reg[0],e->reg[1]);
1083 bool VoronoiDiagramGenerator::voronoi()
1085 Site *newsite, *bot, *top, *temp, *p;
1089 Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
1093 bottomsite = nextone();
1095 bool retval = ELinitialize();
1100 newsite = nextone();
1105 newintstar = PQ_min();
1109 if (newsite != (Site *)NULL && (PQempty() || newsite->coord.y < newintstar.y
1110 || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
1113 lbnd = ELleftbnd(&(newsite->coord));
1114 rbnd = ELright(lbnd);
1115 bot = rightreg(lbnd);
1116 e = bisect(bot, newsite);
1117 bisector = HEcreate(e, le);
1118 ELinsert(lbnd, bisector);
1120 if ((p = intersect(lbnd, bisector)) != (Site *) NULL)
1123 PQinsert(lbnd, p, dist(p,newsite));
1126 bisector = HEcreate(e, re);
1127 ELinsert(lbnd, bisector);
1129 if ((p = intersect(bisector, rbnd)) != (Site *) NULL)
1131 PQinsert(bisector, p, dist(p,newsite));
1133 newsite = nextone();
1135 else if (!PQempty())
1137 lbnd = PQextractmin();
1138 llbnd = ELleft(lbnd);
1139 rbnd = ELright(lbnd);
1140 rrbnd = ELright(rbnd);
1141 bot = leftreg(lbnd);
1142 top = rightreg(rbnd);
1148 endpoint(lbnd->ELedge,lbnd->ELpm,v);
1149 endpoint(rbnd->ELedge,rbnd->ELpm,v);
1155 if (bot->coord.y > top->coord.y)
1162 e = bisect(bot, top);
1164 bisector = HEcreate(e, pm);
1165 ELinsert(llbnd, bisector);
1166 endpoint(e, re-pm, v);
1172 if((p = intersect(llbnd, bisector)) != (Site *) NULL)
1175 PQinsert(llbnd, p, dist(p,bot));
1179 if ((p = intersect(bisector, rrbnd)) != (Site *) NULL)
1181 PQinsert(bisector, p, dist(p,bot));
1190 for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
1204 int scomp(
const void *p1,
const void *p2)
1206 VPoint *s1 = (VPoint*)p1, *s2=(VPoint*)p2;
1207 if(s1->y < s2->y)
return(-1);
1208 if(s1->y > s2->y)
return(1);
1209 if(s1->x < s2->x)
return(-1);
1210 if(s1->x > s2->x)
return(1);
1215 Site * VoronoiDiagramGenerator::nextone()
1218 if(siteidx < nsites)
1220 s = &sites[siteidx];
1225 return( (Site *)NULL);
1228 FASTJET_END_NAMESPACE