1 #ifndef __FASTJET__VORONOI_H__
2 #define __FASTJET__VORONOI_H__
3
4 //FJSTARTHEADER
5 // $Id: Voronoi.hh 4442 2020-05-05 07:50:11Z soyez $
6 //
7 // Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
8 //
9 //----------------------------------------------------------------------
10 // This file is part of FastJet.
11 //
12 // FastJet is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU General Public License as published by
14 // the Free Software Foundation; either version 2 of the License, or
15 // (at your option) any later version.
16 //
17 // The algorithms that underlie FastJet have required considerable
18 // development. They are described in the original FastJet paper,
19 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
20 // FastJet as part of work towards a scientific publication, please
21 // quote the version you use and include a citation to the manual and
22 // optionally also to hep-ph/0512210.
23 //
24 // FastJet is distributed in the hope that it will be useful,
25 // but WITHOUT ANY WARRANTY; without even the implied warranty of
26 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 // GNU General Public License for more details.
28 //
29 // You should have received a copy of the GNU General Public License
30 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
31 //----------------------------------------------------------------------
32 //FJENDHEADER
33
34
35 /*
36 * The author of this software is Steven Fortune.
37 * Copyright (c) 1994 by AT&T Bell Laboratories.
38 * Permission to use, copy, modify, and distribute this software for any
39 * purpose without fee is hereby granted, provided that this entire notice
40 * is included in all copies of any software which is or includes a copy
41 * or modification of this software and in all copies of the supporting
42 * documentation for such software.
43 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
44 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
45 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
46 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
47 */
48
49 /*
50 * This code was originally written by Stephan Fortune in C code. I,
51 * Shane O'Sullivan, have since modified it, encapsulating it in a C++
52 * class and, fixing memory leaks and adding accessors to the Voronoi
53 * Edges. Permission to use, copy, modify, and distribute this
54 * software for any purpose without fee is hereby granted, provided
55 * that this entire notice is included in all copies of any software
56 * which is or includes a copy or modification of this software and in
57 * all copies of the supporting documentation for such software. THIS
58 * SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
59 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
60 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE
61 * MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR
62 * PURPOSE.
63 */
64
65 /*
66 * This code, included in the FastJet distribution, was originally
67 * written by Stephan Fortune in C and adapted to C++ by Shane
68 * O'Sullivan under the terms repported above.
69 *
70 * Below are the list of changes implemented by the FastJet authors:
71 *
72 * 2011-11-14 Gregory Soyez <soyez@fastjet.fr>
73 *
74 * * removed 'plot' and 'triangulate' (were always 0)
75 * * removed unused plot functions (openpl, circle, range,
76 * out_bisector, out_ep, out_vertex, out_site, out_triple)
77 * * removed unused 'VPoint p' in 'intersect'
78 *
79 *
80 * 2011-07-22 Gregory Soyez <soyez@fastjet.fr>
81 *
82 * * replaced Point by VPoint (to avoid any potential conflict
83 * with an already existing class Point in FastJet
84 *
85 *
86 * 2008-04-01 Gregory Soyez <soyez@fastjet.fr>
87 *
88 * * declared ystar volatile in HalfEdge (apparently fixes a bug
89 * related to VD computations with points on a grid)
90 *
91 *
92 * 2007-05-07 Gregory Soyez <soyez@fastjet.fr>
93 *
94 * * put the code in the fastjet namespace
95 *
96 * * replaced float by double
97 *
98 * * generateVoronoi() takes a vector of Point instead of 2
99 * pointers
100 *
101 * * added info about the parent sites to GraphEdge
102 *
103 * * removed condition on minimal distance between sites
104 *
105 */
106
107 #include "fastjet/LimitedWarning.hh"
108 #include <vector>
109 #include <math.h>
110 #include <stdlib.h>
111 #include <string.h>
112
113 #define DELETED -2
114 #define le 0
115 #define re 1
116
117 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
118
119 /**
120 * \if internal_doc
121 * @ingroup internal
122 * \class VPoint
123 * class to handle a 2d point
124 * \endif
125 */
126 class VPoint{
127 public:
128 /// defailt ctor
VPoint()129 VPoint() : x(0.0), y(0.0) {}
130
131 /// ctor with initialisation
VPoint(double _x,double _y)132 VPoint(double _x, double _y) : x(_x), y(_y) {}
133
134 /// addition
operator +(const VPoint & p) const135 inline VPoint operator + (const VPoint &p) const{
136 return VPoint(x+p.x, y+p.y);
137 }
138
139 /// subtraction
operator -(const VPoint & p) const140 inline VPoint operator - (const VPoint &p) const{
141 return VPoint(x-p.x, y-p.y);
142 }
143
144 /// scalar multiplication
operator *(const double t) const145 inline VPoint operator * (const double t) const{
146 return VPoint(x*t, y*t);
147 }
148
149 /// vector coordinates
150 double x,y;
151 };
152
153
154 /// norm of a vector
norm(const VPoint p)155 inline double norm(const VPoint p){
156 return p.x*p.x+p.y*p.y;
157 }
158
159
160 /// 2D vector product
vector_product(const VPoint & p1,const VPoint & p2)161 inline double vector_product(const VPoint &p1, const VPoint &p2){
162 return p1.x*p2.y-p1.y*p2.x;
163 }
164
165
166 /// scalar product
scalar_product(const VPoint & p1,const VPoint & p2)167 inline double scalar_product(const VPoint &p1, const VPoint &p2){
168 return p1.x*p2.x+p1.y*p2.y;
169 }
170
171
172 /**
173 * \if internal_doc
174 * @ingroup internal
175 * \class GraphEdge
176 * handle an edge of the Voronoi Diagram.
177 * \endif
178 */
179 class GraphEdge{
180 public:
181 /// coordinates of the extreme points
182 double x1,y1,x2,y2;
183
184 /// indices of the parent sites that define the edge
185 int point1, point2;
186
187 /// pointer to the next edge
188 GraphEdge* next;
189 };
190
191
192 /**
193 * \if internal_doc
194 * @ingroup internal
195 * \class Site
196 * structure used both for particle sites and for vertices.
197 * \endif
198 */
199 class Site{
200 public:
201 VPoint coord;
202 int sitenbr;
203 int refcnt;
204 };
205
206
207
208 class Freenode{
209 public:
210 Freenode *nextfree;
211 };
212
213
214 class FreeNodeArrayList{
215 public:
216 Freenode* memory;
217 FreeNodeArrayList* next;
218 };
219
220
221 class Freelist{
222 public:
223 Freenode *head;
224 int nodesize;
225 };
226
227 class Edge{
228 public:
229 double a,b,c;
230 Site *ep[2];
231 Site *reg[2];
232 int edgenbr;
233 };
234
235
236 class Halfedge{
237 public:
238 Halfedge *ELleft, *ELright;
239 Edge *ELedge;
240 int ELrefcnt;
241 char ELpm;
242 Site *vertex;
243 volatile double ystar;
244 Halfedge *PQnext;
245 };
246
247 /**
248 * \if internal_doc
249 * @ingroup internal
250 * \class VoronoiDiagramGenerator
251 * Shane O'Sullivan C++ version of Stephan Fortune Voronoi diagram
252 * generator
253 * \endif
254 */
255 class VoronoiDiagramGenerator{
256 public:
257 VoronoiDiagramGenerator();
258 ~VoronoiDiagramGenerator();
259
260 bool generateVoronoi(std::vector<VPoint> *_parent_sites,
261 double minX, double maxX, double minY, double maxY,
262 double minDist=0);
263
resetIterator()264 inline void resetIterator(){
265 iteratorEdges = allEdges;
266 }
267
getNext(GraphEdge ** e)268 bool getNext(GraphEdge **e){
269 if(iteratorEdges == 0)
270 return false;
271
272 *e = iteratorEdges;
273 iteratorEdges = iteratorEdges->next;
274 return true;
275 }
276
277 std::vector<VPoint> *parent_sites;
278 int n_parent_sites;
279
280 private:
281 void cleanup();
282 void cleanupEdges();
283 char *getfree(Freelist *fl);
284 Halfedge *PQfind();
285 int PQempty();
286
287 Halfedge **ELhash;
288 Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
289 Halfedge *HEcreate(Edge *e,int pm);
290
291 VPoint PQ_min();
292 Halfedge *PQextractmin();
293 void freeinit(Freelist *fl,int size);
294 void makefree(Freenode *curr,Freelist *fl);
295 void geominit();
296 void plotinit();
297
298 // GS: removed the unused (always ==0) argument
299 bool voronoi(/*int triangulate*/);
300 void ref(Site *v);
301 void deref(Site *v);
302 void endpoint(Edge *e,int lr,Site * s);
303
304 void ELdelete(Halfedge *he);
305 Halfedge *ELleftbnd(VPoint *p);
306 Halfedge *ELright(Halfedge *he);
307 void makevertex(Site *v);
308
309 void PQinsert(Halfedge *he,Site * v, double offset);
310 void PQdelete(Halfedge *he);
311 bool ELinitialize();
312 void ELinsert(Halfedge *lb, Halfedge *newHe);
313 Halfedge * ELgethash(int b);
314 Halfedge *ELleft(Halfedge *he);
315 Site *leftreg(Halfedge *he);
316 bool PQinitialize();
317 int PQbucket(Halfedge *he);
318 void clip_line(Edge *e);
319 char *myalloc(unsigned n);
320 int right_of(Halfedge *el,VPoint *p);
321
322 Site *rightreg(Halfedge *he);
323 Edge *bisect(Site *s1, Site *s2);
324 double dist(Site *s,Site *t);
325
326 // GS: 'p' is unused and always ==0 (see also comment by
327 // S. O'Sullivan in the source file), so we remove it
328 Site *intersect(Halfedge *el1, Halfedge *el2 /*, VPoint *p=0*/);
329
330 Site *nextone();
331
332 void pushGraphEdge(double x1, double y1, double x2, double y2,
333 Site *s1, Site *s2);
334
335 // Gregory Soyez: unused plotting methods
336 // void openpl();
337 // void circle(double x, double y, double radius);
338 // void range(double minX, double minY, double maxX, double maxY);
339 //
340 // void out_bisector(Edge *e);
341 // void out_ep(Edge *e);
342 // void out_vertex(Site *v);
343 // void out_site(Site *s);
344 //
345 // void out_triple(Site *s1, Site *s2,Site * s3);
346
347 Freelist hfl;
348 Halfedge *ELleftend, *ELrightend;
349 int ELhashsize;
350
351 int sorted, debug;
352 double xmin, xmax, ymin, ymax, deltax, deltay;
353
354 Site *sites;
355 int nsites;
356 int siteidx;
357 int sqrt_nsites;
358 int nvertices;
359 Freelist sfl;
360 Site *bottomsite;
361
362 int nedges;
363 Freelist efl;
364 int PQhashsize;
365 Halfedge *PQhash;
366 int PQcount;
367 int PQmin;
368
369 int ntry, totalsearch;
370 double pxmin, pxmax, pymin, pymax, cradius;
371 int total_alloc;
372
373 double borderMinX, borderMaxX, borderMinY, borderMaxY;
374
375 FreeNodeArrayList* allMemoryList;
376 FreeNodeArrayList* currentMemoryBlock;
377
378 GraphEdge* allEdges;
379 GraphEdge* iteratorEdges;
380
381 double minDistanceBetweenSites;
382
383 static LimitedWarning _warning_degeneracy;
384 };
385
386 int scomp(const void *p1,const void *p2);
387
388
389 FASTJET_END_NAMESPACE
390
391 #endif
392