fastjet/internal/Voronoi.hh

0001 #ifndef __FASTJET__VORONOI_H__
0002 #define __FASTJET__VORONOI_H__
0003
0004 //FJSTARTHEADER
0005 // $Id$
0006 //
0007 // Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
0008 //
0009 //----------------------------------------------------------------------
0010 // This file is part of FastJet.
0011 //
0012 //  FastJet is free software; you can redistribute it and/or modify
0013 //  it under the terms of the GNU General Public License as published by
0014 //  the Free Software Foundation; either version 2 of the License, or
0015 //  (at your option) any later version.
0016 //
0017 //  The algorithms that underlie FastJet have required considerable
0018 //  development. They are described in the original FastJet paper,
0019 //  hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
0020 //  FastJet as part of work towards a scientific publication, please
0021 //  quote the version you use and include a citation to the manual and
0022 //  optionally also to hep-ph/0512210.
0023 //
0024 //  FastJet is distributed in the hope that it will be useful,
0025 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
0026 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0027 //  GNU General Public License for more details.
0028 //
0029 //  You should have received a copy of the GNU General Public License
0030 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
0031 //----------------------------------------------------------------------
0032 //FJENDHEADER
0033
0034
0035 /*
0036 * The author of this software is Steven Fortune.
0037 * Copyright (c) 1994 by AT&T Bell Laboratories.
0038 * Permission to use, copy, modify, and distribute this software for any
0039 * purpose without fee is hereby granted, provided that this entire notice
0040 * is included in all copies of any software which is or includes a copy
0041 * or modification of this software and in all copies of the supporting
0042 * documentation for such software.
0043 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
0044 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
0045 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
0046 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
0047 */
0048
0049 /*
0050 * This code was originally written by Stephan Fortune in C code.  I,
0051 * Shane O'Sullivan, have since modified it, encapsulating it in a C++
0052 * class and, fixing memory leaks and adding accessors to the Voronoi
0053 * Edges.  Permission to use, copy, modify, and distribute this
0054 * software for any purpose without fee is hereby granted, provided
0055 * that this entire notice is included in all copies of any software
0056 * which is or includes a copy or modification of this software and in
0057 * all copies of the supporting documentation for such software.  THIS
0058 * SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
0059 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
0060 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE
0061 * MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR
0062 * PURPOSE.
0063 */
0064
0065 /*
0066  * This code, included in the FastJet distribution, was originally
0067  * written by Stephan Fortune in C and adapted to C++ by Shane
0068  * O'Sullivan under the terms repported above.
0069  *
0070  * Below are the list of changes implemented by the FastJet authors:
0071  *
0072  * 2011-11-14  Gregory Soyez  <soyez@fastjet.fr>
0073  *
0074  *      * removed 'plot' and 'triangulate' (were always 0)
0075  *      * removed unused plot functions (openpl, circle, range,
0076  *        out_bisector, out_ep, out_vertex, out_site, out_triple)
0077  *      * removed unused 'VPoint p' in 'intersect'
0078  *
0079  *
0080  * 2011-07-22  Gregory Soyez  <soyez@fastjet.fr>
0081  *
0082  *      * replaced Point by VPoint (to avoid any potential conflict
0083  *        with an already existing class Point in FastJet
0084  *
0085  *
0086  * 2008-04-01  Gregory Soyez  <soyez@fastjet.fr>
0087  *
0088  *      * declared ystar volatile in HalfEdge (apparently fixes a bug
0089  *        related to VD computations with points on a grid)
0090  *
0091  *
0092  * 2007-05-07  Gregory Soyez  <soyez@fastjet.fr>
0093  *
0094  *      * put the code in the fastjet namespace
0095  *
0096  *      * replaced float by double
0097  *
0098  *      * generateVoronoi() takes a vector of Point instead of 2
0099  *        pointers
0100  *
0101  *      * added info about the parent sites to GraphEdge
0102  *
0103  *      * removed condition on minimal distance between sites
0104  *
0105  */
0106
0107 #include "fastjet/LimitedWarning.hh"
0108 #include <vector>
0109 #include <math.h>
0110 #include <stdlib.h>
0111 #include <string.h>
0112
0113 #define DELETED -2
0114 #define le 0
0115 #define re 1
0116
0117 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
0118
0119 /**
0120  * \if internal_doc
0121  * @ingroup internal
0122  * \class VPoint
0123  * class to handle a 2d point
0124  * \endif
0125  */
0126 class VPoint{
0127 public:
0128   /// defailt ctor
0129   VPoint() : x(0.0), y(0.0) {}
0130
0131   /// ctor with initialisation
0132   VPoint(double _x, double _y) : x(_x), y(_y) {}
0133
0134   /// addition
0135   inline VPoint operator + (const VPoint &p) const{
0136     return VPoint(x+p.x, y+p.y);
0137   }
0138
0139   /// subtraction
0140   inline VPoint operator - (const VPoint &p) const{
0141     return VPoint(x-p.x, y-p.y);
0142   }
0143
0144   /// scalar multiplication
0145   inline VPoint operator * (const double t) const{
0146     return VPoint(x*t, y*t);
0147   }
0148
0149   /// vector coordinates
0150   double x,y;
0151 };
0152
0153
0154 /// norm of a vector
0155 inline double norm(const VPoint p){
0156   return p.x*p.x+p.y*p.y;
0157 }
0158
0159
0160 /// 2D vector product
0161 inline double vector_product(const VPoint &p1, const VPoint &p2){
0162   return p1.x*p2.y-p1.y*p2.x;
0163 }
0164
0165
0166 /// scalar product
0167 inline double scalar_product(const VPoint &p1, const VPoint &p2){
0168   return p1.x*p2.x+p1.y*p2.y;
0169 }
0170
0171
0172 /**
0173  * \if internal_doc
0174  * @ingroup internal
0175  * \class GraphEdge
0176  * handle an edge of the Voronoi Diagram.
0177  * \endif
0178  */
0179 class GraphEdge{
0180 public:
0181   /// coordinates of the extreme points
0182   double x1,y1,x2,y2;
0183
0184   /// indices of the parent sites that define the edge
0185   int point1, point2;
0186
0187   /// pointer to the next edge
0188   GraphEdge* next;
0189 };
0190
0191
0192 /**
0193  * \if internal_doc
0194  * @ingroup internal
0195  * \class Site
0196  * structure used both for particle sites and for vertices.
0197  * \endif
0198  */
0199 class Site{
0200  public:
0201   VPoint    coord;
0202   int sitenbr;
0203   int refcnt;
0204 };
0205
0206
0207
0208 class Freenode{
0209 public:
0210   Freenode *nextfree;
0211 };
0212
0213
0214 class FreeNodeArrayList{
0215 public:
0216   Freenode* memory;
0217   FreeNodeArrayList* next;
0218 };
0219
0220
0221 class Freelist{
0222 public:
0223   Freenode *head;
0224   int nodesize;
0225 };
0226
0227 class Edge{
0228 public:
0229   double a,b,c;
0230   Site *ep[2];
0231   Site *reg[2];
0232   int edgenbr;
0233 };
0234
0235
0236 class Halfedge{
0237 public:
0238   Halfedge *ELleft, *ELright;
0239   Edge *ELedge;
0240   int ELrefcnt;
0241   char ELpm;
0242   Site *vertex;
0243   volatile double ystar;
0244   Halfedge *PQnext;
0245 };
0246
0247 /**
0248  * \if internal_doc
0249  * @ingroup internal
0250  * \class VoronoiDiagramGenerator
0251  * Shane O'Sullivan C++ version of Stephan Fortune Voronoi diagram
0252  * generator
0253  * \endif
0254  */
0255 class VoronoiDiagramGenerator{
0256 public:
0257   VoronoiDiagramGenerator();
0258   ~VoronoiDiagramGenerator();
0259
0260   bool generateVoronoi(std::vector<VPoint> *_parent_sites,
0261                double minX, double maxX, double minY, double maxY,
0262                double minDist=0);
0263
0264   inline void resetIterator(){
0265     iteratorEdges = allEdges;
0266   }
0267
0268   bool getNext(GraphEdge **e){
0269     if(iteratorEdges == 0)
0270       return false;
0271
0272     *e = iteratorEdges;
0273     iteratorEdges = iteratorEdges->next;
0274     return true;
0275   }
0276
0277   std::vector<VPoint> *parent_sites;
0278   int n_parent_sites;
0279
0280 private:
0281   void cleanup();
0282   void cleanupEdges();
0283   char *getfree(Freelist *fl);
0284   Halfedge *PQfind();
0285   int PQempty();
0286
0287   Halfedge **ELhash;
0288   Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
0289   Halfedge *HEcreate(Edge *e,int pm);
0290
0291   VPoint PQ_min();
0292   Halfedge *PQextractmin();
0293   void freeinit(Freelist *fl,int size);
0294   void makefree(Freenode *curr,Freelist *fl);
0295   void geominit();
0296   void plotinit();
0297
0298   // GS: removed the unused (always ==0) argument
0299   bool voronoi(/*int triangulate*/);
0300   void ref(Site *v);
0301   void deref(Site *v);
0302   void endpoint(Edge *e,int lr,Site * s);
0303
0304   void ELdelete(Halfedge *he);
0305   Halfedge *ELleftbnd(VPoint *p);
0306   Halfedge *ELright(Halfedge *he);
0307   void makevertex(Site *v);
0308
0309   void PQinsert(Halfedge *he,Site * v, double offset);
0310   void PQdelete(Halfedge *he);
0311   bool ELinitialize();
0312   void ELinsert(Halfedge *lb, Halfedge *newHe);
0313   Halfedge * ELgethash(int b);
0314   Halfedge *ELleft(Halfedge *he);
0315   Site *leftreg(Halfedge *he);
0316   bool PQinitialize();
0317   int PQbucket(Halfedge *he);
0318   void clip_line(Edge *e);
0319   char *myalloc(unsigned n);
0320   int right_of(Halfedge *el,VPoint *p);
0321
0322   Site *rightreg(Halfedge *he);
0323   Edge *bisect(Site *s1, Site *s2);
0324   double dist(Site *s,Site *t);
0325
0326   // GS: 'p' is unused and always ==0 (see also comment by
0327   //     S. O'Sullivan in the source file), so we remove it
0328   Site *intersect(Halfedge *el1, Halfedge *el2 /*, VPoint *p=0*/);
0329
0330   Site *nextone();
0331
0332   void pushGraphEdge(double x1, double y1, double x2, double y2,
0333              Site *s1, Site *s2);
0334
0335   // Gregory Soyez: unused plotting methods
0336   // void openpl();
0337   // void circle(double x, double y, double radius);
0338   // void range(double minX, double minY, double maxX, double maxY);
0339   //
0340   // void out_bisector(Edge *e);
0341   // void out_ep(Edge *e);
0342   // void out_vertex(Site *v);
0343   // void out_site(Site *s);
0344   //
0345   // void out_triple(Site *s1, Site *s2,Site * s3);
0346
0347   Freelist hfl;
0348   Halfedge *ELleftend, *ELrightend;
0349   int ELhashsize;
0350
0351   int sorted, debug;
0352   double xmin, xmax, ymin, ymax, deltax, deltay;
0353
0354   Site *sites;
0355   int nsites;
0356   int siteidx;
0357   int sqrt_nsites;
0358   int nvertices;
0359   Freelist sfl;
0360   Site *bottomsite;
0361
0362   int nedges;
0363   Freelist efl;
0364   int PQhashsize;
0365   Halfedge *PQhash;
0366   int PQcount;
0367   int PQmin;
0368
0369   int ntry, totalsearch;
0370   double pxmin, pxmax, pymin, pymax, cradius;
0371   int total_alloc;
0372
0373   double borderMinX, borderMaxX, borderMinY, borderMaxY;
0374
0375   FreeNodeArrayList* allMemoryList;
0376   FreeNodeArrayList* currentMemoryBlock;
0377
0378   GraphEdge* allEdges;
0379   GraphEdge* iteratorEdges;
0380
0381   double minDistanceBetweenSites;
0382
0383   static LimitedWarning _warning_degeneracy;
0384 };
0385
0386 int scomp(const void *p1,const void *p2);
0387
0388
0389 FASTJET_END_NAMESPACE
0390
0391 #endif