00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
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
00110 sites = (Site *) myalloc(nsites*sizeof( Site));
00111
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
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
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
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
00232
00233
00234 bucket = (int)((p->x - xmin)/deltax * ELhashsize);
00235
00236
00237
00238 if (bucket<0) bucket =0;
00239 if (bucket>=ELhashsize) bucket = ELhashsize - 1;
00240
00241 he = ELgethash(bucket);
00242
00243
00244
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
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
00261
00262 he = he->ELleft;
00263 } else
00264
00265
00266 do{
00267 he = he->ELleft;
00268 } while ((he!=ELleftend )&& (!right_of(he,p)));
00269
00270
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
00284
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
00313
00314 return bottomsite;
00315
00316
00317
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;
00343 newedge->reg[1] = s2;
00344 ref(s1);
00345 ref(s2);
00346
00347
00348
00349 newedge->ep[0] = (Site*) NULL;
00350 newedge->ep[1] = (Site*) NULL;
00351
00352
00353 dx = s2->coord.x - s1->coord.x;
00354 dy = s2->coord.y - s1->coord.y;
00355
00356
00357 adx = dx>0 ? dx : -dx;
00358 ady = dy>0 ? dy : -dy;
00359
00360
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
00366 newedge->a = 1.0; newedge->b = dy/dx; newedge->c /= dx;
00367 } else {
00368
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
00380
00381
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
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
00420
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
00429
00430
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
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
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
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
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
00734
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;
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
00803
00804
00805
00806
00807
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
00835 y1 = pymax;
00836
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
00846 y2 = pymin;
00847
00848 }
00849 x2 = (e->c) - (e->b) * y2;
00850 if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin)))
00851 {
00852
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
00872
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
00882
00883 x2 = pxmin;
00884 }
00885 y2 = e->c - e->a * x2;
00886 if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin)))
00887 {
00888
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
00902
00903
00904 pushGraphEdge(x1,y1,x2,y2,e->reg[0],e->reg[1]);
00905 }
00906
00907
00908
00909
00910
00911
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
00938
00939
00940 if (newsite != (Site *)NULL && (PQempty() || newsite->coord.y < newintstar.y
00941 || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
00942 {
00943 out_site(newsite);
00944 lbnd = ELleftbnd(&(newsite->coord));
00945 rbnd = ELright(lbnd);
00946 bot = rightreg(lbnd);
00947 e = bisect(bot, newsite);
00948 bisector = HEcreate(e, le);
00949 ELinsert(lbnd, bisector);
00950
00951 if ((p = intersect(lbnd, bisector)) != (Site *) NULL)
00952 {
00953 PQdelete(lbnd);
00954 PQinsert(lbnd, p, dist(p,newsite));
00955 };
00956 lbnd = bisector;
00957 bisector = HEcreate(e, re);
00958 ELinsert(lbnd, bisector);
00959
00960 if ((p = intersect(bisector, rbnd)) != (Site *) NULL)
00961 {
00962 PQinsert(bisector, p, dist(p,newsite));
00963 };
00964 newsite = nextone();
00965 }
00966 else if (!PQempty())
00967 {
00968 lbnd = PQextractmin();
00969 llbnd = ELleft(lbnd);
00970 rbnd = ELright(lbnd);
00971 rrbnd = ELright(rbnd);
00972 bot = leftreg(lbnd);
00973 top = rightreg(rbnd);
00974
00975 out_triple(bot, top, rightreg(lbnd));
00976
00977 v = lbnd->vertex;
00978 makevertex(v);
00979 endpoint(lbnd->ELedge,lbnd->ELpm,v);
00980 endpoint(rbnd->ELedge,rbnd->ELpm,v);
00981 ELdelete(lbnd);
00982 PQdelete(rbnd);
00983 ELdelete(rbnd);
00984 pm = le;
00985
00986 if (bot->coord.y > top->coord.y)
00987 {
00988 temp = bot;
00989 bot = top;
00990 top = temp;
00991 pm = re;
00992 }
00993 e = bisect(bot, top);
00994
00995 bisector = HEcreate(e, pm);
00996 ELinsert(llbnd, bisector);
00997 endpoint(e, re-pm, v);
00998
00999
01000 deref(v);
01001
01002
01003 if((p = intersect(llbnd, bisector)) != (Site *) NULL)
01004 {
01005 PQdelete(llbnd);
01006 PQinsert(llbnd, p, dist(p,bot));
01007 };
01008
01009
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
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
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