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 |
---|