[5] | 1 | //STARTHEADER |
---|
| 2 | // $Id: Voronoi.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
| 3 | // |
---|
| 4 | // Copyright (c) 1994 by AT&T Bell Laboratories (see below) |
---|
| 5 | // |
---|
| 6 | // |
---|
| 7 | //---------------------------------------------------------------------- |
---|
| 8 | // This file is included as part of FastJet but was mostly written by |
---|
| 9 | // S. Fortune in C, put into C++ with memory management by S |
---|
| 10 | // O'Sullivan, and with further interface and memory management |
---|
| 11 | // modifications by Gregory Soyez. |
---|
| 12 | // |
---|
| 13 | // Permission to use, copy, modify, and distribute this software for |
---|
| 14 | // any purpose without fee is hereby granted, provided that this |
---|
| 15 | // entire notice is included in all copies of any software which is or |
---|
| 16 | // includes a copy or modification of this software and in all copies |
---|
| 17 | // of the supporting documentation for such software. THIS SOFTWARE IS |
---|
| 18 | // BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED WARRANTY. |
---|
| 19 | // IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY REPRESENTATION |
---|
| 20 | // OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS |
---|
| 21 | // SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. |
---|
| 22 | // |
---|
| 23 | //---------------------------------------------------------------------- |
---|
| 24 | //ENDHEADER |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | /* |
---|
| 28 | * The author of this software is Steven Fortune. |
---|
| 29 | * Copyright (c) 1994 by AT&T Bell Laboratories. |
---|
| 30 | * Permission to use, copy, modify, and distribute this software for any |
---|
| 31 | * purpose without fee is hereby granted, provided that this entire notice |
---|
| 32 | * is included in all copies of any software which is or includes a copy |
---|
| 33 | * or modification of this software and in all copies of the supporting |
---|
| 34 | * documentation for such software. |
---|
| 35 | * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED |
---|
| 36 | * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY |
---|
| 37 | * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY |
---|
| 38 | * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. |
---|
| 39 | */ |
---|
| 40 | |
---|
| 41 | /* |
---|
| 42 | * This code was originally written by Stephan Fortune in C code. I, |
---|
| 43 | * Shane O'Sullivan, have since modified it, encapsulating it in a C++ |
---|
| 44 | * class and, fixing memory leaks and adding accessors to the Voronoi |
---|
| 45 | * Edges. Permission to use, copy, modify, and distribute this |
---|
| 46 | * software for any purpose without fee is hereby granted, provided |
---|
| 47 | * that this entire notice is included in all copies of any software |
---|
| 48 | * which is or includes a copy or modification of this software and in |
---|
| 49 | * all copies of the supporting documentation for such software. THIS |
---|
| 50 | * SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED |
---|
| 51 | * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY |
---|
| 52 | * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE |
---|
| 53 | * MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR |
---|
| 54 | * PURPOSE. |
---|
| 55 | */ |
---|
| 56 | |
---|
| 57 | /* |
---|
| 58 | * This code, included in the FastJet distribution, was originally |
---|
| 59 | * written by Stephan Fortune in C and adapted to C++ by Shane |
---|
| 60 | * O'Sullivan under the terms reported above. |
---|
| 61 | * |
---|
| 62 | * Below are the list of changes implemented by the FastJet authors: |
---|
| 63 | * |
---|
| 64 | * 2011-11-14 Gregory Soyez <soyez@fastjet.fr> |
---|
| 65 | * |
---|
| 66 | * * removed 'plot' and 'triangulate' (were always 0) |
---|
| 67 | * * removed unused plot functions (openpl, circle, range, |
---|
| 68 | * out_bisector, out_ep, out_vertex, out_site, out_triple) |
---|
| 69 | * * removed unused 'VPoint p' in 'intersect' |
---|
| 70 | * |
---|
| 71 | * |
---|
| 72 | * 2011-07-22 Gregory Soyez <soyez@fastjet.fr> |
---|
| 73 | * |
---|
| 74 | * * replaced Point by VPoint (to avoid any potential conflict |
---|
| 75 | * with an already existing class Point in FastJet |
---|
| 76 | * |
---|
| 77 | * |
---|
| 78 | * 2011-06-28 Gregory Soyez <soyez@fastjet.fr> |
---|
| 79 | * |
---|
| 80 | * * added support for situations with degenerate particles (we just |
---|
| 81 | * discard the particles degenerate wiht an already existing |
---|
| 82 | * one which is perfectly sufficient for our needs) |
---|
| 83 | * * in 'VoronoiDiagramGenerator::intersect', improved the numerical |
---|
| 84 | * precision in cases where the 2 parents are nearly degenerate |
---|
| 85 | * |
---|
| 86 | * |
---|
| 87 | * 2011-06-14 Gregory Soyez <soyez@fastjet.fr> |
---|
| 88 | * |
---|
| 89 | * * fixed a potential overflow bug in VoronoiDiagramGenerator::PQbucket |
---|
| 90 | * |
---|
| 91 | * |
---|
| 92 | * 2007-05-07 Gregory Soyez <soyez@fastjet.fr> |
---|
| 93 | * |
---|
| 94 | * * fied a few memory leaks |
---|
| 95 | * |
---|
| 96 | * * put the code in the fastjet namespace |
---|
| 97 | * |
---|
| 98 | * * replaced float by double |
---|
| 99 | * |
---|
| 100 | * * generateVoronoi() takes a vector of Point instead of 2 |
---|
| 101 | * pointers |
---|
| 102 | * |
---|
| 103 | * * added info about the parent sites to GraphEdge (and clip_line) |
---|
| 104 | * |
---|
| 105 | * * removed condition on minimal distance between sites |
---|
| 106 | */ |
---|
| 107 | |
---|
| 108 | #include <stdio.h> |
---|
| 109 | #include "fastjet/internal/Voronoi.hh" |
---|
| 110 | |
---|
| 111 | using namespace std; |
---|
| 112 | |
---|
| 113 | FASTJET_BEGIN_NAMESPACE |
---|
| 114 | |
---|
| 115 | LimitedWarning VoronoiDiagramGenerator::_warning_degeneracy; |
---|
| 116 | |
---|
| 117 | VoronoiDiagramGenerator::VoronoiDiagramGenerator(){ |
---|
| 118 | siteidx = 0; |
---|
| 119 | sites = NULL; |
---|
| 120 | parent_sites = NULL; |
---|
| 121 | |
---|
| 122 | allMemoryList = new FreeNodeArrayList; |
---|
| 123 | allMemoryList->memory = NULL; |
---|
| 124 | allMemoryList->next = NULL; |
---|
| 125 | currentMemoryBlock = allMemoryList; |
---|
| 126 | allEdges = NULL; |
---|
| 127 | iteratorEdges = 0; |
---|
| 128 | minDistanceBetweenSites = 0; |
---|
| 129 | |
---|
| 130 | ELhash = NULL; |
---|
| 131 | PQhash = NULL; |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | VoronoiDiagramGenerator::~VoronoiDiagramGenerator(){ |
---|
| 135 | cleanup(); |
---|
| 136 | cleanupEdges(); |
---|
| 137 | |
---|
| 138 | if (allMemoryList != NULL) |
---|
| 139 | delete allMemoryList; |
---|
| 140 | } |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | |
---|
| 144 | bool VoronoiDiagramGenerator::generateVoronoi(vector<VPoint> *_parent_sites, |
---|
| 145 | double minX, double maxX, |
---|
| 146 | double minY, double maxY, |
---|
| 147 | double minDist){ |
---|
| 148 | cleanup(); |
---|
| 149 | cleanupEdges(); |
---|
| 150 | int i; |
---|
| 151 | double x, y; |
---|
| 152 | |
---|
| 153 | minDistanceBetweenSites = minDist; |
---|
| 154 | |
---|
| 155 | parent_sites = _parent_sites; |
---|
| 156 | |
---|
| 157 | nsites = n_parent_sites = parent_sites->size(); |
---|
| 158 | debug = 1; |
---|
| 159 | sorted = 0; |
---|
| 160 | freeinit(&sfl, sizeof (Site)); |
---|
| 161 | |
---|
| 162 | //sites = (Site *) myalloc(nsites*sizeof( *sites)); |
---|
| 163 | sites = (Site *) myalloc(nsites*sizeof( Site)); |
---|
| 164 | //parent_sites = (Site *) myalloc(nsites*sizeof( Site)); |
---|
| 165 | |
---|
| 166 | if (sites == 0) |
---|
| 167 | return false; |
---|
| 168 | |
---|
| 169 | xmax = xmin = (*parent_sites)[0].x; |
---|
| 170 | ymax = ymin = (*parent_sites)[0].y; |
---|
| 171 | |
---|
| 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; |
---|
| 178 | sites[i].refcnt = 0; |
---|
| 179 | |
---|
| 180 | if (x<xmin) |
---|
| 181 | xmin=x; |
---|
| 182 | else if (x>xmax) |
---|
| 183 | xmax=x; |
---|
| 184 | |
---|
| 185 | if (y<ymin) |
---|
| 186 | ymin=y; |
---|
| 187 | else if (y>ymax) |
---|
| 188 | ymax=y; |
---|
| 189 | } |
---|
| 190 | |
---|
| 191 | qsort(sites, nsites, sizeof (*sites), scomp); |
---|
| 192 | |
---|
| 193 | // Gregory Soyez |
---|
| 194 | // |
---|
| 195 | // check if some of the particles are degenerate to avoid a crash. |
---|
| 196 | // |
---|
| 197 | // At the moment, we work under the assumption that they will be |
---|
| 198 | // "clustered" later on, so we just keep the 1st one and discard the |
---|
| 199 | // others |
---|
| 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){ |
---|
| 203 | offset++; |
---|
| 204 | } else if (offset>0){ |
---|
| 205 | sites[is-offset] = sites[is]; |
---|
| 206 | } |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | if (offset>0){ |
---|
| 210 | nsites-=offset; |
---|
| 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."); |
---|
| 212 | } |
---|
| 213 | |
---|
| 214 | siteidx = 0; |
---|
| 215 | geominit(); |
---|
| 216 | double temp = 0; |
---|
| 217 | if(minX > maxX){ |
---|
| 218 | temp = minX; |
---|
| 219 | minX = maxX; |
---|
| 220 | maxX = temp; |
---|
| 221 | } |
---|
| 222 | if(minY > maxY){ |
---|
| 223 | temp = minY; |
---|
| 224 | minY = maxY; |
---|
| 225 | maxY = temp; |
---|
| 226 | } |
---|
| 227 | borderMinX = minX; |
---|
| 228 | borderMinY = minY; |
---|
| 229 | borderMaxX = maxX; |
---|
| 230 | borderMaxY = maxY; |
---|
| 231 | |
---|
| 232 | siteidx = 0; |
---|
| 233 | voronoi(); |
---|
| 234 | |
---|
| 235 | return true; |
---|
| 236 | } |
---|
| 237 | |
---|
| 238 | bool VoronoiDiagramGenerator::ELinitialize(){ |
---|
| 239 | int i; |
---|
| 240 | freeinit(&hfl, sizeof(Halfedge)); |
---|
| 241 | ELhashsize = 2 * sqrt_nsites; |
---|
| 242 | //ELhash = (Halfedge **) myalloc ( sizeof *ELhash * ELhashsize); |
---|
| 243 | ELhash = (Halfedge **) myalloc ( sizeof(Halfedge*)*ELhashsize); |
---|
| 244 | |
---|
| 245 | if(ELhash == 0) |
---|
| 246 | return false; |
---|
| 247 | |
---|
| 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; |
---|
| 257 | |
---|
| 258 | return true; |
---|
| 259 | } |
---|
| 260 | |
---|
| 261 | |
---|
| 262 | Halfedge* VoronoiDiagramGenerator::HEcreate(Edge *e,int pm){ |
---|
| 263 | Halfedge *answer; |
---|
| 264 | answer = (Halfedge*) getfree(&hfl); |
---|
| 265 | answer->ELedge = e; |
---|
| 266 | answer->ELpm = pm; |
---|
| 267 | answer->PQnext = (Halfedge*) NULL; |
---|
| 268 | answer->vertex = (Site*) NULL; |
---|
| 269 | answer->ELrefcnt = 0; |
---|
| 270 | return answer; |
---|
| 271 | } |
---|
| 272 | |
---|
| 273 | |
---|
| 274 | void VoronoiDiagramGenerator::ELinsert(Halfedge *lb, Halfedge *newHe) |
---|
| 275 | { |
---|
| 276 | newHe->ELleft = lb; |
---|
| 277 | newHe->ELright = lb->ELright; |
---|
| 278 | (lb->ELright)->ELleft = newHe; |
---|
| 279 | lb->ELright = newHe; |
---|
| 280 | } |
---|
| 281 | |
---|
| 282 | |
---|
| 283 | /* Get entry from hash table, pruning any deleted nodes */ |
---|
| 284 | Halfedge* VoronoiDiagramGenerator::ELgethash(int b){ |
---|
| 285 | Halfedge *he; |
---|
| 286 | |
---|
| 287 | if ((b<0) || (b>=ELhashsize)) |
---|
| 288 | return (Halfedge *) NULL; |
---|
| 289 | |
---|
| 290 | he = ELhash[b]; |
---|
| 291 | if ((he==(Halfedge*) NULL) || (he->ELedge!=(Edge*) DELETED)) |
---|
| 292 | return he; |
---|
| 293 | |
---|
| 294 | /* Hash table points to deleted half edge. Patch as necessary. */ |
---|
| 295 | ELhash[b] = (Halfedge*) NULL; |
---|
| 296 | if ((he->ELrefcnt -= 1) == 0) |
---|
| 297 | makefree((Freenode*)he, &hfl); |
---|
| 298 | return (Halfedge*) NULL; |
---|
| 299 | } |
---|
| 300 | |
---|
| 301 | Halfedge * VoronoiDiagramGenerator::ELleftbnd(VPoint *p){ |
---|
| 302 | int i, bucket; |
---|
| 303 | Halfedge *he; |
---|
| 304 | |
---|
| 305 | /* Use hash table to get close to desired halfedge */ |
---|
| 306 | // use the hash function to find the place in the hash map that this |
---|
| 307 | // HalfEdge should be |
---|
| 308 | // Gregory Soyez: the original code was |
---|
| 309 | // |
---|
| 310 | // bucket = (int)((p->x - xmin)/deltax * ELhashsize); |
---|
| 311 | // // make sure that the bucket position in within the range of the |
---|
| 312 | // //hash array |
---|
| 313 | // if (bucket<0) bucket =0; |
---|
| 314 | // if (bucket>=ELhashsize) bucket = ELhashsize - 1; |
---|
| 315 | // |
---|
| 316 | // but this runs the risk of having a overflow which would |
---|
| 317 | // cause bucket to be truncated at 0 instead of ELhashsize - 1 |
---|
| 318 | // (or vice-versa) |
---|
| 319 | // We fix this by performing the test immediately on the double |
---|
| 320 | // We put in a extra bit of margin to be sure the conversion does |
---|
| 321 | // not play dirty tricks on us |
---|
| 322 | |
---|
| 323 | //const double &px = p->x; |
---|
| 324 | if (p->x < xmin){ bucket=0;} |
---|
| 325 | else if (p->x >= xmax){ bucket = ELhashsize - 1;} |
---|
| 326 | else{ |
---|
| 327 | bucket= (int)((p->x - xmin)/deltax * ELhashsize); |
---|
| 328 | if (bucket>=ELhashsize) bucket = ELhashsize - 1; // the lower cut should be robust |
---|
| 329 | } |
---|
| 330 | |
---|
| 331 | he = ELgethash(bucket); |
---|
| 332 | |
---|
| 333 | // if the HE isn't found, search backwards and forwards in the hash |
---|
| 334 | // map for the first non-null entry |
---|
| 335 | if (he == (Halfedge*) NULL){ |
---|
| 336 | for(i=1;1;i++){ |
---|
| 337 | if ((he=ELgethash(bucket-i)) != (Halfedge*) NULL) |
---|
| 338 | break; |
---|
| 339 | if ((he=ELgethash(bucket+i)) != (Halfedge*) NULL) |
---|
| 340 | break; |
---|
| 341 | }; |
---|
| 342 | totalsearch += i; |
---|
| 343 | }; |
---|
| 344 | ntry += 1; |
---|
| 345 | /* Now search linear list of halfedges for the correct one */ |
---|
| 346 | if ((he==ELleftend) || (he != ELrightend && right_of(he,p))){ |
---|
| 347 | do{ |
---|
| 348 | he = he->ELright; |
---|
| 349 | } while (he!=ELrightend && right_of(he,p)); |
---|
| 350 | // keep going right on the list until either the end is reached, |
---|
| 351 | // or you find the 1st edge which the point |
---|
| 352 | he = he->ELleft; //isn't to the right of |
---|
| 353 | } else |
---|
| 354 | // if the point is to the left of the HalfEdge, then search left |
---|
| 355 | // for the HE just to the left of the point |
---|
| 356 | do{ |
---|
| 357 | he = he->ELleft; |
---|
| 358 | } while ((he!=ELleftend )&& (!right_of(he,p))); |
---|
| 359 | |
---|
| 360 | /* Update hash table and reference counts */ |
---|
| 361 | if ((bucket > 0) && (bucket <ELhashsize-1)){ |
---|
| 362 | if(ELhash[bucket] != (Halfedge *) NULL){ |
---|
| 363 | ELhash[bucket]->ELrefcnt -= 1; |
---|
| 364 | } |
---|
| 365 | ELhash[bucket] = he; |
---|
| 366 | ELhash[bucket]->ELrefcnt += 1; |
---|
| 367 | }; |
---|
| 368 | |
---|
| 369 | return he; |
---|
| 370 | } |
---|
| 371 | |
---|
| 372 | |
---|
| 373 | /* This delete routine can't reclaim node, since pointers from hash |
---|
| 374 | table may be present. */ |
---|
| 375 | void VoronoiDiagramGenerator::ELdelete(Halfedge *he){ |
---|
| 376 | (he->ELleft)->ELright = he->ELright; |
---|
| 377 | (he->ELright)->ELleft = he->ELleft; |
---|
| 378 | he->ELedge = (Edge*) DELETED; |
---|
| 379 | } |
---|
| 380 | |
---|
| 381 | |
---|
| 382 | Halfedge* VoronoiDiagramGenerator::ELright(Halfedge *he){ |
---|
| 383 | return he->ELright; |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | |
---|
| 387 | Halfedge* VoronoiDiagramGenerator::ELleft(Halfedge *he){ |
---|
| 388 | return he->ELleft; |
---|
| 389 | } |
---|
| 390 | |
---|
| 391 | |
---|
| 392 | Site * VoronoiDiagramGenerator::leftreg(Halfedge *he){ |
---|
| 393 | if (he->ELedge == (Edge*) NULL) |
---|
| 394 | return(bottomsite); |
---|
| 395 | return (he->ELpm == le) |
---|
| 396 | ? he->ELedge->reg[le] |
---|
| 397 | : he->ELedge->reg[re]; |
---|
| 398 | } |
---|
| 399 | |
---|
| 400 | Site * VoronoiDiagramGenerator::rightreg(Halfedge *he){ |
---|
| 401 | if (he->ELedge == (Edge*) NULL) |
---|
| 402 | // if this halfedge has no edge, return the bottom site (whatever |
---|
| 403 | // that is) |
---|
| 404 | return bottomsite; |
---|
| 405 | |
---|
| 406 | // if the ELpm field is zero, return the site 0 that this edge |
---|
| 407 | // bisects, otherwise return site number 1 |
---|
| 408 | return (he->ELpm == le) |
---|
| 409 | ? he->ELedge->reg[re] |
---|
| 410 | : he->ELedge->reg[le]; |
---|
| 411 | } |
---|
| 412 | |
---|
| 413 | void VoronoiDiagramGenerator::geominit(){ |
---|
| 414 | double sn; |
---|
| 415 | |
---|
| 416 | freeinit(&efl, sizeof(Edge)); |
---|
| 417 | nvertices = 0; |
---|
| 418 | nedges = 0; |
---|
| 419 | sn = (double)nsites+4; |
---|
| 420 | sqrt_nsites = (int)sqrt(sn); |
---|
| 421 | deltay = ymax - ymin; |
---|
| 422 | deltax = xmax - xmin; |
---|
| 423 | } |
---|
| 424 | |
---|
| 425 | |
---|
| 426 | Edge * VoronoiDiagramGenerator::bisect(Site *s1, Site *s2){ |
---|
| 427 | double dx,dy,adx,ady; |
---|
| 428 | Edge *newedge; |
---|
| 429 | |
---|
| 430 | newedge = (Edge*) getfree(&efl); |
---|
| 431 | |
---|
| 432 | newedge->reg[0] = s1; //store the sites that this edge is bisecting |
---|
| 433 | newedge->reg[1] = s2; |
---|
| 434 | ref(s1); |
---|
| 435 | ref(s2); |
---|
| 436 | |
---|
| 437 | // to begin with, there are no endpoints on the bisector - it goes |
---|
| 438 | // to infinity |
---|
| 439 | newedge->ep[0] = (Site*) NULL; |
---|
| 440 | newedge->ep[1] = (Site*) NULL; |
---|
| 441 | |
---|
| 442 | // get the difference in x dist between the sites |
---|
| 443 | dx = s2->coord.x - s1->coord.x; |
---|
| 444 | dy = s2->coord.y - s1->coord.y; |
---|
| 445 | |
---|
| 446 | // make sure that the difference is positive |
---|
| 447 | adx = dx>0 ? dx : -dx; |
---|
| 448 | ady = dy>0 ? dy : -dy; |
---|
| 449 | |
---|
| 450 | // get the slope of the line |
---|
| 451 | newedge->c = (double)(s1->coord.x * dx + s1->coord.y * dy |
---|
| 452 | + (dx*dx + dy*dy)*0.5); |
---|
| 453 | |
---|
| 454 | if (adx>ady){ |
---|
| 455 | //set formula of line, with x fixed to 1 |
---|
| 456 | newedge->a = 1.0; newedge->b = dy/dx; newedge->c /= dx; |
---|
| 457 | } else { |
---|
| 458 | //set formula of line, with y fixed to 1 |
---|
| 459 | newedge->b = 1.0; newedge->a = dx/dy; newedge->c /= dy; |
---|
| 460 | } |
---|
| 461 | |
---|
| 462 | newedge->edgenbr = nedges; |
---|
| 463 | nedges++; |
---|
| 464 | |
---|
| 465 | return newedge; |
---|
| 466 | } |
---|
| 467 | |
---|
| 468 | |
---|
| 469 | // create a new site where the HalfEdges el1 and el2 intersect - note |
---|
| 470 | // that the VPoint in the argument list is not used, don't know why |
---|
| 471 | // it's there |
---|
| 472 | // |
---|
| 473 | // Gregory Soyez: removed the uinused point p |
---|
| 474 | Site* VoronoiDiagramGenerator::intersect(Halfedge *el1, Halfedge *el2 |
---|
| 475 | /*, VPoint *p*/){ |
---|
| 476 | Edge *e1,*e2, *e; |
---|
| 477 | Halfedge *el; |
---|
| 478 | double d, xint, yint; |
---|
| 479 | int right_of_site; |
---|
| 480 | Site *v; |
---|
| 481 | |
---|
| 482 | e1 = el1->ELedge; |
---|
| 483 | e2 = el2->ELedge; |
---|
| 484 | if ((e1==(Edge*)NULL) || (e2==(Edge*)NULL)) |
---|
| 485 | return ((Site*) NULL); |
---|
| 486 | |
---|
| 487 | // if the two edges bisect the same parent, return null |
---|
| 488 | if (e1->reg[1] == e2->reg[1]) |
---|
| 489 | return (Site*) NULL; |
---|
| 490 | |
---|
| 491 | // Gregory Soyez: |
---|
| 492 | // |
---|
| 493 | // if the 2 parents are too close, the intersection is going to be |
---|
| 494 | // computed from the "long edges" of the triangle which could causes |
---|
| 495 | // large rounding errors. In this case, use the bisector of the 2 |
---|
| 496 | // parents to find the interaction point |
---|
| 497 | // |
---|
| 498 | // The following replaces |
---|
| 499 | // d = e1->a * e2->b - e1->b * e2->a; |
---|
| 500 | // if (-1.0e-10<d && d<1.0e-10) |
---|
| 501 | // return (Site*) NULL; |
---|
| 502 | // |
---|
| 503 | // xint = (e1->c*e2->b - e2->c*e1->b)/d; |
---|
| 504 | // yint = (e2->c*e1->a - e1->c*e2->a)/d; |
---|
| 505 | |
---|
| 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; |
---|
| 510 | |
---|
| 511 | if (dx*dx + dy*dy < 1e-14*(dxref*dxref+dyref*dyref)){ |
---|
| 512 | // make sure that the difference is positive |
---|
| 513 | double adx = dx>0 ? dx : -dx; |
---|
| 514 | double ady = dy>0 ? dy : -dy; |
---|
| 515 | |
---|
| 516 | // get the slope of the line |
---|
| 517 | double a,b; |
---|
| 518 | double c = (double)(e1->reg[1]->coord.x * dx + e1->reg[1]->coord.y * dy |
---|
| 519 | + (dx*dx + dy*dy)*0.5); |
---|
| 520 | |
---|
| 521 | if (adx>ady){ |
---|
| 522 | a = 1.0; b = dy/dx; c /= dx; |
---|
| 523 | } else { |
---|
| 524 | b = 1.0; a = dx/dy; c /= dy; |
---|
| 525 | } |
---|
| 526 | |
---|
| 527 | d = e1->a * b - e1->b * a; |
---|
| 528 | if (-1.0e-10<d && d<1.0e-10) { |
---|
| 529 | return (Site*) NULL; |
---|
| 530 | } |
---|
| 531 | |
---|
| 532 | xint = (e1->c*b - c*e1->b)/d; |
---|
| 533 | yint = (c*e1->a - e1->c*a)/d; |
---|
| 534 | |
---|
| 535 | } else { |
---|
| 536 | d = e1->a * e2->b - e1->b * e2->a; |
---|
| 537 | if (-1.0e-10<d && d<1.0e-10) { |
---|
| 538 | return (Site*) NULL; |
---|
| 539 | } |
---|
| 540 | |
---|
| 541 | xint = (e1->c*e2->b - e2->c*e1->b)/d; |
---|
| 542 | yint = (e2->c*e1->a - e1->c*e2->a)/d; |
---|
| 543 | } |
---|
| 544 | // end of Gregory Soyez's modifications |
---|
| 545 | |
---|
| 546 | volatile double local_y1 = e1->reg[1]->coord.y; |
---|
| 547 | volatile double local_y2 = e2->reg[1]->coord.y; |
---|
| 548 | |
---|
| 549 | if( (local_y1 < local_y2) || |
---|
| 550 | ((local_y1 == local_y2) && |
---|
| 551 | (e1->reg[1]->coord.x < e2->reg[1]->coord.x)) ){ |
---|
| 552 | el = el1; |
---|
| 553 | e = e1; |
---|
| 554 | } else { |
---|
| 555 | el = el2; |
---|
| 556 | e = e2; |
---|
| 557 | } |
---|
| 558 | |
---|
| 559 | right_of_site = xint >= e->reg[1]->coord.x; |
---|
| 560 | if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re)) |
---|
| 561 | return (Site*) NULL; |
---|
| 562 | |
---|
| 563 | // create a new site at the point of intersection - this is a new |
---|
| 564 | // vector event waiting to happen |
---|
| 565 | v = (Site*) getfree(&sfl); |
---|
| 566 | v->refcnt = 0; |
---|
| 567 | v->coord.x = xint; |
---|
| 568 | v->coord.y = yint; |
---|
| 569 | return v; |
---|
| 570 | } |
---|
| 571 | |
---|
| 572 | //HERE |
---|
| 573 | |
---|
| 574 | /* returns 1 if p is to right of halfedge e */ |
---|
| 575 | int VoronoiDiagramGenerator::right_of(Halfedge *el,VPoint *p) |
---|
| 576 | { |
---|
| 577 | Edge *e; |
---|
| 578 | Site *topsite; |
---|
| 579 | int right_of_site, above, fast; |
---|
| 580 | double dxp, dyp, dxs, t1, t2, t3, yl; |
---|
| 581 | |
---|
| 582 | e = el->ELedge; |
---|
| 583 | topsite = e->reg[1]; |
---|
| 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); |
---|
| 587 | |
---|
| 588 | if (e->a == 1.0) |
---|
| 589 | { dyp = p->y - topsite->coord.y; |
---|
| 590 | dxp = p->x - topsite->coord.x; |
---|
| 591 | fast = 0; |
---|
| 592 | if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) ) |
---|
| 593 | { above = dyp>= e->b*dxp; |
---|
| 594 | fast = above; |
---|
| 595 | } |
---|
| 596 | else |
---|
| 597 | { above = p->x + p->y*e->b > e-> c; |
---|
| 598 | if(e->b<0.0) above = !above; |
---|
| 599 | if (!above) fast = 1; |
---|
| 600 | }; |
---|
| 601 | if (!fast) |
---|
| 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; |
---|
| 606 | }; |
---|
| 607 | } |
---|
| 608 | else /*e->b==1.0 */ |
---|
| 609 | { yl = e->c - e->a*p->x; |
---|
| 610 | t1 = p->y - yl; |
---|
| 611 | t2 = p->x - topsite->coord.x; |
---|
| 612 | t3 = yl - topsite->coord.y; |
---|
| 613 | above = t1*t1 > t2*t2 + t3*t3; |
---|
| 614 | }; |
---|
| 615 | return (el->ELpm==le ? above : !above); |
---|
| 616 | } |
---|
| 617 | |
---|
| 618 | |
---|
| 619 | void VoronoiDiagramGenerator::endpoint(Edge *e,int lr,Site * s) |
---|
| 620 | { |
---|
| 621 | e->ep[lr] = s; |
---|
| 622 | ref(s); |
---|
| 623 | if(e->ep[re-lr]== (Site *) NULL) |
---|
| 624 | return; |
---|
| 625 | |
---|
| 626 | clip_line(e); |
---|
| 627 | |
---|
| 628 | deref(e->reg[le]); |
---|
| 629 | deref(e->reg[re]); |
---|
| 630 | makefree((Freenode*)e, &efl); |
---|
| 631 | } |
---|
| 632 | |
---|
| 633 | |
---|
| 634 | double VoronoiDiagramGenerator::dist(Site *s,Site *t) |
---|
| 635 | { |
---|
| 636 | double dx,dy; |
---|
| 637 | dx = s->coord.x - t->coord.x; |
---|
| 638 | dy = s->coord.y - t->coord.y; |
---|
| 639 | return (double)(sqrt(dx*dx + dy*dy)); |
---|
| 640 | } |
---|
| 641 | |
---|
| 642 | |
---|
| 643 | void VoronoiDiagramGenerator::makevertex(Site *v) |
---|
| 644 | { |
---|
| 645 | v->sitenbr = nvertices; |
---|
| 646 | nvertices += 1; |
---|
| 647 | //GS unused plot: out_vertex(v); |
---|
| 648 | } |
---|
| 649 | |
---|
| 650 | |
---|
| 651 | void VoronoiDiagramGenerator::deref(Site *v) |
---|
| 652 | { |
---|
| 653 | v->refcnt -= 1; |
---|
| 654 | if (v->refcnt == 0 ) |
---|
| 655 | makefree((Freenode*)v, &sfl); |
---|
| 656 | } |
---|
| 657 | |
---|
| 658 | void VoronoiDiagramGenerator::ref(Site *v) |
---|
| 659 | { |
---|
| 660 | v->refcnt += 1; |
---|
| 661 | } |
---|
| 662 | |
---|
| 663 | //push the HalfEdge into the ordered linked list of vertices |
---|
| 664 | void VoronoiDiagramGenerator::PQinsert(Halfedge *he,Site * v, double offset) |
---|
| 665 | { |
---|
| 666 | Halfedge *last, *next; |
---|
| 667 | |
---|
| 668 | he->vertex = v; |
---|
| 669 | ref(v); |
---|
| 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))) |
---|
| 675 | { |
---|
| 676 | last = next; |
---|
| 677 | }; |
---|
| 678 | he->PQnext = last->PQnext; |
---|
| 679 | last->PQnext = he; |
---|
| 680 | PQcount += 1; |
---|
| 681 | } |
---|
| 682 | |
---|
| 683 | //remove the HalfEdge from the list of vertices |
---|
| 684 | void VoronoiDiagramGenerator::PQdelete(Halfedge *he) |
---|
| 685 | { |
---|
| 686 | Halfedge *last; |
---|
| 687 | |
---|
| 688 | if(he->vertex != (Site *) NULL) |
---|
| 689 | { |
---|
| 690 | last = &PQhash[PQbucket(he)]; |
---|
| 691 | while (last->PQnext != he) |
---|
| 692 | last = last->PQnext; |
---|
| 693 | |
---|
| 694 | last->PQnext = he->PQnext; |
---|
| 695 | PQcount -= 1; |
---|
| 696 | deref(he->vertex); |
---|
| 697 | he->vertex = (Site *) NULL; |
---|
| 698 | }; |
---|
| 699 | } |
---|
| 700 | |
---|
| 701 | int VoronoiDiagramGenerator::PQbucket(Halfedge *he) |
---|
| 702 | { |
---|
| 703 | // Gregory Soyez: the original code was |
---|
| 704 | // |
---|
| 705 | // bucket = (int)((he->ystar - ymin)/deltay * PQhashsize); |
---|
| 706 | // if (bucket<0) bucket = 0; |
---|
| 707 | // if (bucket>=PQhashsize) bucket = PQhashsize-1 ; |
---|
| 708 | // if (bucket < PQmin) PQmin = bucket; |
---|
| 709 | // return(bucket); |
---|
| 710 | // |
---|
| 711 | // but this runs the risk of having a overflow which would |
---|
| 712 | // cause bucket to be truncated at 0 instead of PQhashsize-1 |
---|
| 713 | // (or vice-versa) |
---|
| 714 | // We fix this by performing the test immediately on the double |
---|
| 715 | // We put in a extra bit of margin to be sure the conversion does |
---|
| 716 | // not play dirty tricks on us |
---|
| 717 | |
---|
| 718 | int bucket; |
---|
| 719 | |
---|
| 720 | double hey = he->ystar; |
---|
| 721 | if (hey < ymin){ bucket = 0;} |
---|
| 722 | else if (hey >= ymax){ bucket = PQhashsize-1;} |
---|
| 723 | else { |
---|
| 724 | bucket = (int)((hey - ymin)/deltay * PQhashsize); |
---|
| 725 | if (bucket>=PQhashsize) bucket = PQhashsize-1 ; |
---|
| 726 | } |
---|
| 727 | |
---|
| 728 | if (bucket < PQmin) PQmin = bucket; |
---|
| 729 | return(bucket); |
---|
| 730 | } |
---|
| 731 | |
---|
| 732 | |
---|
| 733 | |
---|
| 734 | int VoronoiDiagramGenerator::PQempty() |
---|
| 735 | { |
---|
| 736 | return(PQcount==0); |
---|
| 737 | } |
---|
| 738 | |
---|
| 739 | |
---|
| 740 | VPoint VoronoiDiagramGenerator::PQ_min() |
---|
| 741 | { |
---|
| 742 | VPoint answer; |
---|
| 743 | |
---|
| 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; |
---|
| 747 | return (answer); |
---|
| 748 | } |
---|
| 749 | |
---|
| 750 | Halfedge * VoronoiDiagramGenerator::PQextractmin() |
---|
| 751 | { |
---|
| 752 | Halfedge *curr; |
---|
| 753 | |
---|
| 754 | curr = PQhash[PQmin].PQnext; |
---|
| 755 | PQhash[PQmin].PQnext = curr->PQnext; |
---|
| 756 | PQcount -= 1; |
---|
| 757 | return(curr); |
---|
| 758 | } |
---|
| 759 | |
---|
| 760 | |
---|
| 761 | bool VoronoiDiagramGenerator::PQinitialize() |
---|
| 762 | { |
---|
| 763 | int i; |
---|
| 764 | |
---|
| 765 | PQcount = 0; |
---|
| 766 | PQmin = 0; |
---|
| 767 | PQhashsize = 4 * sqrt_nsites; |
---|
| 768 | //PQhash = (Halfedge *) myalloc(PQhashsize * sizeof *PQhash); |
---|
| 769 | PQhash = (Halfedge *) myalloc(PQhashsize * sizeof(Halfedge)); |
---|
| 770 | |
---|
| 771 | if(PQhash == 0) |
---|
| 772 | return false; |
---|
| 773 | |
---|
| 774 | for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (Halfedge *)NULL; |
---|
| 775 | |
---|
| 776 | return true; |
---|
| 777 | } |
---|
| 778 | |
---|
| 779 | |
---|
| 780 | void VoronoiDiagramGenerator::freeinit(Freelist *fl,int size) |
---|
| 781 | { |
---|
| 782 | fl->head = (Freenode *) NULL; |
---|
| 783 | fl->nodesize = size; |
---|
| 784 | } |
---|
| 785 | |
---|
| 786 | char * VoronoiDiagramGenerator::getfree(Freelist *fl) |
---|
| 787 | { |
---|
| 788 | int i; |
---|
| 789 | Freenode *t; |
---|
| 790 | |
---|
| 791 | if(fl->head == (Freenode *) NULL) |
---|
| 792 | { |
---|
| 793 | t = (Freenode *) myalloc(sqrt_nsites * fl->nodesize); |
---|
| 794 | |
---|
| 795 | if(t == 0) |
---|
| 796 | return 0; |
---|
| 797 | |
---|
| 798 | currentMemoryBlock->next = new FreeNodeArrayList; |
---|
| 799 | currentMemoryBlock = currentMemoryBlock->next; |
---|
| 800 | currentMemoryBlock->memory = t; |
---|
| 801 | currentMemoryBlock->next = 0; |
---|
| 802 | |
---|
| 803 | for(i=0; i<sqrt_nsites; i+=1) |
---|
| 804 | makefree((Freenode *)((char *)t+i*fl->nodesize), fl); |
---|
| 805 | }; |
---|
| 806 | t = fl->head; |
---|
| 807 | fl->head = (fl->head)->nextfree; |
---|
| 808 | return((char *)t); |
---|
| 809 | } |
---|
| 810 | |
---|
| 811 | |
---|
| 812 | |
---|
| 813 | void VoronoiDiagramGenerator::makefree(Freenode *curr,Freelist *fl) |
---|
| 814 | { |
---|
| 815 | curr->nextfree = fl->head; |
---|
| 816 | fl->head = curr; |
---|
| 817 | } |
---|
| 818 | |
---|
| 819 | void VoronoiDiagramGenerator::cleanup() |
---|
| 820 | { |
---|
| 821 | if(sites != NULL){ |
---|
| 822 | free(sites); |
---|
| 823 | sites = 0; |
---|
| 824 | } |
---|
| 825 | |
---|
| 826 | FreeNodeArrayList* current=NULL, *prev=NULL; |
---|
| 827 | |
---|
| 828 | current = prev = allMemoryList; |
---|
| 829 | |
---|
| 830 | while(current->next!=NULL){ |
---|
| 831 | prev = current; |
---|
| 832 | current = current->next; |
---|
| 833 | free(prev->memory); |
---|
| 834 | delete prev; |
---|
| 835 | prev = NULL; |
---|
| 836 | } |
---|
| 837 | |
---|
| 838 | if(current!=NULL){ |
---|
| 839 | if (current->memory!=NULL){ |
---|
| 840 | free(current->memory); |
---|
| 841 | } |
---|
| 842 | delete current; |
---|
| 843 | } |
---|
| 844 | |
---|
| 845 | allMemoryList = new FreeNodeArrayList; |
---|
| 846 | allMemoryList->next = NULL; |
---|
| 847 | allMemoryList->memory = NULL; |
---|
| 848 | currentMemoryBlock = allMemoryList; |
---|
| 849 | |
---|
| 850 | if (ELhash!=NULL) |
---|
| 851 | free(ELhash); |
---|
| 852 | |
---|
| 853 | if (PQhash!=NULL) |
---|
| 854 | free(PQhash); |
---|
| 855 | } |
---|
| 856 | |
---|
| 857 | void VoronoiDiagramGenerator::cleanupEdges() |
---|
| 858 | { |
---|
| 859 | GraphEdge* geCurrent = NULL, *gePrev = NULL; |
---|
| 860 | geCurrent = gePrev = allEdges; |
---|
| 861 | |
---|
| 862 | while(geCurrent!=NULL){ |
---|
| 863 | gePrev = geCurrent; |
---|
| 864 | geCurrent = geCurrent->next; |
---|
| 865 | delete gePrev; |
---|
| 866 | } |
---|
| 867 | |
---|
| 868 | allEdges = 0; |
---|
| 869 | |
---|
| 870 | } |
---|
| 871 | |
---|
| 872 | void VoronoiDiagramGenerator::pushGraphEdge(double x1, double y1, double x2, double y2, |
---|
| 873 | Site *s1, Site *s2) |
---|
| 874 | { |
---|
| 875 | GraphEdge* newEdge = new GraphEdge; |
---|
| 876 | newEdge->next = allEdges; |
---|
| 877 | allEdges = newEdge; |
---|
| 878 | newEdge->x1 = x1; |
---|
| 879 | newEdge->y1 = y1; |
---|
| 880 | newEdge->x2 = x2; |
---|
| 881 | newEdge->y2 = y2; |
---|
| 882 | |
---|
| 883 | newEdge->point1 = s1->sitenbr; |
---|
| 884 | newEdge->point2 = s2->sitenbr; |
---|
| 885 | } |
---|
| 886 | |
---|
| 887 | |
---|
| 888 | char * VoronoiDiagramGenerator::myalloc(unsigned n) |
---|
| 889 | { |
---|
| 890 | char *t=0; |
---|
| 891 | t=(char*)malloc(n); |
---|
| 892 | total_alloc += n; |
---|
| 893 | return(t); |
---|
| 894 | } |
---|
| 895 | |
---|
| 896 | |
---|
| 897 | // unused plot functions |
---|
| 898 | // |
---|
| 899 | // /* for those who don't have Cherry's plot */ |
---|
| 900 | // /* #include <plot.h> */ |
---|
| 901 | // void VoronoiDiagramGenerator::openpl(){} |
---|
| 902 | // void VoronoiDiagramGenerator::circle(double x, double y, double radius){} |
---|
| 903 | // void VoronoiDiagramGenerator::range(double minX, double minY, double maxX, double maxY){} |
---|
| 904 | // |
---|
| 905 | // |
---|
| 906 | // |
---|
| 907 | // void VoronoiDiagramGenerator::out_bisector(Edge *e) |
---|
| 908 | // { |
---|
| 909 | // |
---|
| 910 | // |
---|
| 911 | // } |
---|
| 912 | // |
---|
| 913 | // |
---|
| 914 | // void VoronoiDiagramGenerator::out_ep(Edge *e) |
---|
| 915 | // { |
---|
| 916 | // |
---|
| 917 | // |
---|
| 918 | // } |
---|
| 919 | // |
---|
| 920 | // void VoronoiDiagramGenerator::out_vertex(Site *v) |
---|
| 921 | // { |
---|
| 922 | // |
---|
| 923 | // } |
---|
| 924 | // |
---|
| 925 | // |
---|
| 926 | // void VoronoiDiagramGenerator::out_site(Site *s) |
---|
| 927 | // { |
---|
| 928 | // // Gregory Soyez: |
---|
| 929 | // // plot was always 0 so the expression below was always false |
---|
| 930 | // // and even if it was not, 'circle' does nothing! |
---|
| 931 | // // |
---|
| 932 | // // if(!triangulate & plot & !debug) |
---|
| 933 | // // circle (s->coord.x, s->coord.y, cradius); |
---|
| 934 | // |
---|
| 935 | // } |
---|
| 936 | // |
---|
| 937 | // |
---|
| 938 | // void VoronoiDiagramGenerator::out_triple(Site *s1, Site *s2,Site * s3) |
---|
| 939 | // { |
---|
| 940 | // |
---|
| 941 | // } |
---|
| 942 | |
---|
| 943 | |
---|
| 944 | |
---|
| 945 | void VoronoiDiagramGenerator::plotinit() |
---|
| 946 | { |
---|
| 947 | double dx,dy,d; |
---|
| 948 | |
---|
| 949 | dy = ymax - ymin; |
---|
| 950 | dx = xmax - xmin; |
---|
| 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); |
---|
| 957 | //GS unused: openpl(); |
---|
| 958 | //GS unused: range(pxmin, pymin, pxmax, pymax); |
---|
| 959 | } |
---|
| 960 | |
---|
| 961 | |
---|
| 962 | void VoronoiDiagramGenerator::clip_line(Edge *e) |
---|
| 963 | { |
---|
| 964 | Site *s1, *s2; |
---|
| 965 | double x1=0,x2=0,y1=0,y2=0; //, temp = 0; |
---|
| 966 | |
---|
| 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; |
---|
| 971 | |
---|
| 972 | //if the distance between the two points this line was created from is less than |
---|
| 973 | //the square root of 2, then ignore it |
---|
| 974 | //TODO improve/remove |
---|
| 975 | //if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites) |
---|
| 976 | // { |
---|
| 977 | // return; |
---|
| 978 | // } |
---|
| 979 | pxmin = borderMinX; |
---|
| 980 | pxmax = borderMaxX; |
---|
| 981 | pymin = borderMinY; |
---|
| 982 | pymax = borderMaxY; |
---|
| 983 | |
---|
| 984 | if(e->a == 1.0 && e ->b >= 0.0) |
---|
| 985 | { |
---|
| 986 | s1 = e->ep[1]; |
---|
| 987 | s2 = e->ep[0]; |
---|
| 988 | } |
---|
| 989 | else |
---|
| 990 | { |
---|
| 991 | s1 = e->ep[0]; |
---|
| 992 | s2 = e->ep[1]; |
---|
| 993 | }; |
---|
| 994 | |
---|
| 995 | if(e->a == 1.0) |
---|
| 996 | { |
---|
| 997 | y1 = pymin; |
---|
| 998 | if (s1!=(Site *)NULL && s1->coord.y > pymin) |
---|
| 999 | { |
---|
| 1000 | y1 = s1->coord.y; |
---|
| 1001 | } |
---|
| 1002 | if(y1>pymax) |
---|
| 1003 | { |
---|
| 1004 | // printf("\nClipped (1) y1 = %f to %f",y1,pymax); |
---|
| 1005 | y1 = pymax; |
---|
| 1006 | //return; |
---|
| 1007 | } |
---|
| 1008 | x1 = e->c - e->b * y1; |
---|
| 1009 | y2 = pymax; |
---|
| 1010 | if (s2!=(Site *)NULL && s2->coord.y < pymax) |
---|
| 1011 | y2 = s2->coord.y; |
---|
| 1012 | |
---|
| 1013 | if(y2<pymin) |
---|
| 1014 | { |
---|
| 1015 | //printf("\nClipped (2) y2 = %f to %f",y2,pymin); |
---|
| 1016 | y2 = pymin; |
---|
| 1017 | //return; |
---|
| 1018 | } |
---|
| 1019 | x2 = (e->c) - (e->b) * y2; |
---|
| 1020 | if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin))) |
---|
| 1021 | { |
---|
| 1022 | //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax); |
---|
| 1023 | return; |
---|
| 1024 | } |
---|
| 1025 | if(x1> pxmax) |
---|
| 1026 | { x1 = pxmax; y1 = (e->c - x1)/e->b;}; |
---|
| 1027 | if(x1<pxmin) |
---|
| 1028 | { x1 = pxmin; y1 = (e->c - x1)/e->b;}; |
---|
| 1029 | if(x2>pxmax) |
---|
| 1030 | { x2 = pxmax; y2 = (e->c - x2)/e->b;}; |
---|
| 1031 | if(x2<pxmin) |
---|
| 1032 | { x2 = pxmin; y2 = (e->c - x2)/e->b;}; |
---|
| 1033 | } |
---|
| 1034 | else |
---|
| 1035 | { |
---|
| 1036 | x1 = pxmin; |
---|
| 1037 | if (s1!=(Site *)NULL && s1->coord.x > pxmin) |
---|
| 1038 | x1 = s1->coord.x; |
---|
| 1039 | if(x1>pxmax) |
---|
| 1040 | { |
---|
| 1041 | //printf("\nClipped (3) x1 = %f to %f",x1,pxmin); |
---|
| 1042 | //return; |
---|
| 1043 | x1 = pxmax; |
---|
| 1044 | } |
---|
| 1045 | y1 = e->c - e->a * x1; |
---|
| 1046 | x2 = pxmax; |
---|
| 1047 | if (s2!=(Site *)NULL && s2->coord.x < pxmax) |
---|
| 1048 | x2 = s2->coord.x; |
---|
| 1049 | if(x2<pxmin) |
---|
| 1050 | { |
---|
| 1051 | //printf("\nClipped (4) x2 = %f to %f",x2,pxmin); |
---|
| 1052 | //return; |
---|
| 1053 | x2 = pxmin; |
---|
| 1054 | } |
---|
| 1055 | y2 = e->c - e->a * x2; |
---|
| 1056 | if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin))) |
---|
| 1057 | { |
---|
| 1058 | //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax); |
---|
| 1059 | return; |
---|
| 1060 | } |
---|
| 1061 | if(y1> pymax) |
---|
| 1062 | { y1 = pymax; x1 = (e->c - y1)/e->a;}; |
---|
| 1063 | if(y1<pymin) |
---|
| 1064 | { y1 = pymin; x1 = (e->c - y1)/e->a;}; |
---|
| 1065 | if(y2>pymax) |
---|
| 1066 | { y2 = pymax; x2 = (e->c - y2)/e->a;}; |
---|
| 1067 | if(y2<pymin) |
---|
| 1068 | { y2 = pymin; x2 = (e->c - y2)/e->a;}; |
---|
| 1069 | }; |
---|
| 1070 | |
---|
| 1071 | //printf("\nPushing line (%f,%f,%f,%f)",x1,y1,x2,y2); |
---|
| 1072 | //fprintf(stdout, "Line with vertices (%f,%f) and (%f,%f)\n", |
---|
| 1073 | // e->reg[0]->coord.x, e->reg[1]->coord.x, e->reg[0]->coord.y, e->reg[1]->coord.y); |
---|
| 1074 | pushGraphEdge(x1,y1,x2,y2,e->reg[0],e->reg[1]); |
---|
| 1075 | } |
---|
| 1076 | |
---|
| 1077 | |
---|
| 1078 | /* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax, |
---|
| 1079 | deltax, deltay (can all be estimates). |
---|
| 1080 | Performance suffers if they are wrong; better to make nsites, |
---|
| 1081 | deltax, and deltay too big than too small. (?) */ |
---|
| 1082 | |
---|
| 1083 | bool VoronoiDiagramGenerator::voronoi() |
---|
| 1084 | { |
---|
| 1085 | Site *newsite, *bot, *top, *temp, *p; |
---|
| 1086 | Site *v; |
---|
| 1087 | VPoint newintstar; |
---|
| 1088 | int pm; |
---|
| 1089 | Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector; |
---|
| 1090 | Edge *e; |
---|
| 1091 | |
---|
| 1092 | PQinitialize(); |
---|
| 1093 | bottomsite = nextone(); |
---|
| 1094 | //GS unused plot: out_site(bottomsite); |
---|
| 1095 | bool retval = ELinitialize(); |
---|
| 1096 | |
---|
| 1097 | if(!retval) |
---|
| 1098 | return false; |
---|
| 1099 | |
---|
| 1100 | newsite = nextone(); |
---|
| 1101 | while(1) |
---|
| 1102 | { |
---|
| 1103 | |
---|
| 1104 | if(!PQempty()) |
---|
| 1105 | newintstar = PQ_min(); |
---|
| 1106 | |
---|
| 1107 | //if the lowest site has a smaller y value than the lowest vector intersection, process the site |
---|
| 1108 | //otherwise process the vector intersection |
---|
| 1109 | if (newsite != (Site *)NULL && (PQempty() || newsite->coord.y < newintstar.y |
---|
| 1110 | || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x))) |
---|
| 1111 | {/* new site is smallest - this is a site event*/ |
---|
| 1112 | //GS unused plot: out_site(newsite); //output the site |
---|
| 1113 | lbnd = ELleftbnd(&(newsite->coord)); //get the first HalfEdge to the LEFT of the new site |
---|
| 1114 | rbnd = ELright(lbnd); //get the first HalfEdge to the RIGHT of the new site |
---|
| 1115 | bot = rightreg(lbnd); //if this halfedge has no edge, , bot = bottom site (whatever that is) |
---|
| 1116 | e = bisect(bot, newsite); //create a new edge that bisects |
---|
| 1117 | bisector = HEcreate(e, le); //create a new HalfEdge, setting its ELpm field to 0 |
---|
| 1118 | ELinsert(lbnd, bisector); //insert this new bisector edge between the left and right vectors in a linked list |
---|
| 1119 | |
---|
| 1120 | 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 |
---|
| 1121 | { |
---|
| 1122 | PQdelete(lbnd); |
---|
| 1123 | PQinsert(lbnd, p, dist(p,newsite)); |
---|
| 1124 | }; |
---|
| 1125 | lbnd = bisector; |
---|
| 1126 | bisector = HEcreate(e, re); //create a new HalfEdge, setting its ELpm field to 1 |
---|
| 1127 | ELinsert(lbnd, bisector); //insert the new HE to the right of the original bisector earlier in the IF stmt |
---|
| 1128 | |
---|
| 1129 | if ((p = intersect(bisector, rbnd)) != (Site *) NULL) //if this new bisector intersects with the |
---|
| 1130 | { |
---|
| 1131 | PQinsert(bisector, p, dist(p,newsite)); //push the HE into the ordered linked list of vertices |
---|
| 1132 | }; |
---|
| 1133 | newsite = nextone(); |
---|
| 1134 | } |
---|
| 1135 | else if (!PQempty()) /* intersection is smallest - this is a vector event */ |
---|
| 1136 | { |
---|
| 1137 | lbnd = PQextractmin(); //pop the HalfEdge with the lowest vector off the ordered list of vectors |
---|
| 1138 | llbnd = ELleft(lbnd); //get the HalfEdge to the left of the above HE |
---|
| 1139 | rbnd = ELright(lbnd); //get the HalfEdge to the right of the above HE |
---|
| 1140 | rrbnd = ELright(rbnd); //get the HalfEdge to the right of the HE to the right of the lowest HE |
---|
| 1141 | bot = leftreg(lbnd); //get the Site to the left of the left HE which it bisects |
---|
| 1142 | top = rightreg(rbnd); //get the Site to the right of the right HE which it bisects |
---|
| 1143 | |
---|
| 1144 | //GS unused plot: out_triple(bot, top, rightreg(lbnd)); //output the triple of sites, stating that a circle goes through them |
---|
| 1145 | |
---|
| 1146 | v = lbnd->vertex; //get the vertex that caused this event |
---|
| 1147 | makevertex(v); //set the vertex number - couldn't do this earlier since we didn't know when it would be processed |
---|
| 1148 | endpoint(lbnd->ELedge,lbnd->ELpm,v); //set the endpoint of the left HalfEdge to be this vector |
---|
| 1149 | endpoint(rbnd->ELedge,rbnd->ELpm,v); //set the endpoint of the right HalfEdge to be this vector |
---|
| 1150 | ELdelete(lbnd); //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map |
---|
| 1151 | PQdelete(rbnd); //remove all vertex events to do with the right HE |
---|
| 1152 | ELdelete(rbnd); //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map |
---|
| 1153 | pm = le; //set the pm variable to zero |
---|
| 1154 | |
---|
| 1155 | if (bot->coord.y > top->coord.y) //if the site to the left of the event is higher than the Site |
---|
| 1156 | { //to the right of it, then swap them and set the 'pm' variable to 1 |
---|
| 1157 | temp = bot; |
---|
| 1158 | bot = top; |
---|
| 1159 | top = temp; |
---|
| 1160 | pm = re; |
---|
| 1161 | } |
---|
| 1162 | e = bisect(bot, top); //create an Edge (or line) that is between the two Sites. This creates |
---|
| 1163 | //the formula of the line, and assigns a line number to it |
---|
| 1164 | bisector = HEcreate(e, pm); //create a HE from the Edge 'e', and make it point to that edge with its ELedge field |
---|
| 1165 | ELinsert(llbnd, bisector); //insert the new bisector to the right of the left HE |
---|
| 1166 | endpoint(e, re-pm, v); //set one endpoint to the new edge to be the vector point 'v'. |
---|
| 1167 | //If the site to the left of this bisector is higher than the right |
---|
| 1168 | //Site, then this endpoint is put in position 0; otherwise in pos 1 |
---|
| 1169 | deref(v); //delete the vector 'v' |
---|
| 1170 | |
---|
| 1171 | //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it |
---|
| 1172 | if((p = intersect(llbnd, bisector)) != (Site *) NULL) |
---|
| 1173 | { |
---|
| 1174 | PQdelete(llbnd); |
---|
| 1175 | PQinsert(llbnd, p, dist(p,bot)); |
---|
| 1176 | }; |
---|
| 1177 | |
---|
| 1178 | //if right HE and the new bisector don't intersect, then reinsert it |
---|
| 1179 | if ((p = intersect(bisector, rrbnd)) != (Site *) NULL) |
---|
| 1180 | { |
---|
| 1181 | PQinsert(bisector, p, dist(p,bot)); |
---|
| 1182 | }; |
---|
| 1183 | } |
---|
| 1184 | else break; |
---|
| 1185 | }; |
---|
| 1186 | |
---|
| 1187 | |
---|
| 1188 | |
---|
| 1189 | |
---|
| 1190 | for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd)) |
---|
| 1191 | { |
---|
| 1192 | e = lbnd->ELedge; |
---|
| 1193 | |
---|
| 1194 | clip_line(e); |
---|
| 1195 | }; |
---|
| 1196 | |
---|
| 1197 | //cleanup(); |
---|
| 1198 | |
---|
| 1199 | return true; |
---|
| 1200 | |
---|
| 1201 | } |
---|
| 1202 | |
---|
| 1203 | |
---|
| 1204 | int scomp(const void *p1,const void *p2) |
---|
| 1205 | { |
---|
| 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); |
---|
| 1211 | return(0); |
---|
| 1212 | } |
---|
| 1213 | |
---|
| 1214 | /* return a single in-storage site */ |
---|
| 1215 | Site * VoronoiDiagramGenerator::nextone() |
---|
| 1216 | { |
---|
| 1217 | Site *s; |
---|
| 1218 | if(siteidx < nsites) |
---|
| 1219 | { |
---|
| 1220 | s = &sites[siteidx]; |
---|
| 1221 | siteidx += 1; |
---|
| 1222 | return(s); |
---|
| 1223 | } |
---|
| 1224 | else |
---|
| 1225 | return( (Site *)NULL); |
---|
| 1226 | } |
---|
| 1227 | |
---|
| 1228 | FASTJET_END_NAMESPACE |
---|