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