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