FastJet  3.3.3
Voronoi.hh
1 #ifndef __FASTJET__VORONOI_H__
2 #define __FASTJET__VORONOI_H__
3 
4 //FJSTARTHEADER
5 // $Id: Voronoi.hh 4420 2019-11-29 09:28:20Z soyez $
6 //
7 // Copyright (c) 2005-2019, 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
129  VPoint() : x(0.0), y(0.0) {}
130 
131  /// ctor with initialisation
132  VPoint(double _x, double _y) : x(_x), y(_y) {}
133 
134  /// addition
135  inline VPoint operator + (const VPoint &p) const{
136  return VPoint(x+p.x, y+p.y);
137  }
138 
139  /// subtraction
140  inline VPoint operator - (const VPoint &p) const{
141  return VPoint(x-p.x, y-p.y);
142  }
143 
144  /// scalar multiplication
145  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
155 inline double norm(const VPoint p){
156  return p.x*p.x+p.y*p.y;
157 }
158 
159 
160 /// 2D vector product
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
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 
264  inline void resetIterator(){
265  iteratorEdges = allEdges;
266  }
267 
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
double norm(const VPoint p)
norm of a vector
Definition: Voronoi.hh:155
Selector operator*(const Selector &s1, const Selector &s2)
successive application of 2 selectors
Definition: Selector.cc:559
double scalar_product(const VPoint &p1, const VPoint &p2)
scalar product
Definition: Voronoi.hh:167
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
double vector_product(const VPoint &p1, const VPoint &p2)
2D vector product
Definition: Voronoi.hh:161