1 //FJSTARTHEADER
2 // $Id: DnnPlane.hh 4442 2020-05-05 07:50:11Z soyez $
3 //
4 // Copyright (c) 2005-2020, 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 #ifndef DROP_CGAL // in case we do not have the code for CGAL
33
34 #ifndef __FASTJET_DNNPLANE_HH__
35 #define __FASTJET_DNNPLANE_HH__
36
37 #include "fastjet/internal/Triangulation.hh"
38 #include "fastjet/internal/DynamicNearestNeighbours.hh"
39
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41
42
43 /// \if internal_doc
44 /// @ingroup internal
45 /// \class DnnPlane
46 /// class derived from DynamicNearestNeighbours that provides an
47 /// implementation for the Euclidean plane
48 ///
49 /// This class that uses CGAL Delaunay triangulation for most of the
50 /// work (it allows for easy and efficient removal and addition of
51 /// points and circulation over a point's neighbours). The treatment
52 /// of coincident points is not supported by CGAL and is implemented
53 /// according to the method specified in
54 /// issue-tracker/2012-02-CGAL-coincident/METHOD
55 /// \endif
56 class DnnPlane : public DynamicNearestNeighbours {
57 public:
58 /// empty initaliser
DnnPlane()59 DnnPlane() {}
60
61 /// Initialiser from a set of points on an Eta-Phi plane, where both
62 /// eta and phi can have arbitrary ranges
63 DnnPlane(const std::vector<EtaPhi> &, const bool & verbose = false );
64
65
66 /// Returns the index of the nearest neighbour of point labelled
67 /// by ii (assumes ii is valid)
68 int NearestNeighbourIndex(const int ii) const ;
69
70 /// Returns the distance to the nearest neighbour of point labelled
71 /// by index ii (assumes ii is valid)
72 double NearestNeighbourDistance(const int ii) const ;
73
74 /// Returns true iff the given index corresponds to a point that
75 /// exists in the DNN structure (meaning that it has been added, and
76 /// not removed in the meantime)
77 bool Valid(const int index) const;
78
79 void RemoveAndAddPoints(const std::vector<int> & indices_to_remove,
80 const std::vector<EtaPhi> & points_to_add,
81 std::vector<int> & indices_added,
82 std::vector<int> & indices_of_updated_neighbours);
83
84 /// returns the EtaPhi of point with index i.
85 EtaPhi etaphi(const int i) const;
86 /// returns the eta point with index i.
87 double eta(const int i) const;
88 /// returns the phi point with index i.
89 double phi(const int i) const;
90
91 private:
92
93 /// Structure containing a vertex_handle and cached information on
94 /// the nearest neighbour.
95 struct SuperVertex {
96 Vertex_handle vertex; // NULL indicates inexistence...
97 double NNdistance;
98 int NNindex;
99 int coincidence; // ==vertex->info.val() if no coincidence
100 // points to the coinciding SV in case of coincidence
101 // later on for cylinder put a second vertex?
102 };
103
104 std::vector<SuperVertex> _supervertex;
105 //set<Vertex_handle> _vertex_set;
106 bool _verbose;
107
108 //static const bool _crash_on_coincidence = true;
109 static const bool _crash_on_coincidence = false;
110
111 Triangulation _TR; /// CGAL object for dealing with triangulations
112
113 /// calculates and returns the euclidean distance between points p1
114 /// and p2
_euclid_distance(const Point & p1,const Point & p2) const115 inline double _euclid_distance(const Point& p1, const Point& p2) const {
116 double distx= p1.x()-p2.x();
117 double disty= p1.y()-p2.y();
118 return distx*distx+disty*disty;
119 }
120
121 //----------------------------------------------------------------------
122 /// Determines the index and distance of the nearest neighbour to
123 /// point j and puts the information into the _supervertex entry for j
124 void _SetNearest(const int j);
125
126 //----------------------------------------------------------------------
127 /// Determines and stores the nearest neighbour of j.
128 ///
129 /// For each voronoi neighbour D of j if the distance between j and D
130 /// is less than D's own nearest neighbour, then update the
131 /// nearest-neighbour info in D; push D's index onto
132 /// indices_of_updated_neighbours
133 ///
134 /// Note that j is NOT pushed onto indices_of_updated_neighbours --
135 /// if you want it there, put it there yourself.
136 void _SetAndUpdateNearest(const int j,
137 std::vector<int> & indices_of_updated_neighbours);
138
139 /// given a vertex_handle returned by CGAL on insertion of a new
140 /// points, returns the coinciding vertex's value if it turns out
141 /// that it corresponds to a vertex that we already knew about
142 /// (usually because two points coincide)
143 int _CheckIfVertexPresent(const Vertex_handle & vertex,
144 const int its_index);
145
146 //----------------------------------------------------------------------
147 /// if the distance between 'pref' and 'candidate' is smaller (or
148 /// equal) than the one between 'pref' and 'near', return true and
149 /// set 'mindist' to that distance. Note that it is assumed that
150 /// 'mindist' is the euclidian distance between 'pref' and 'near'
151 ///
152 /// Note that the 'near' point is passed through its vertex rather
153 /// than as a point. This allows us to handle cases where we have no min
154 /// yet (near is the infinite vertex)
_is_closer_to(const Point & pref,const Point & candidate,const Vertex_handle & near,double & dist,double & mindist)155 inline bool _is_closer_to(const Point &pref,
156 const Point &candidate,
157 const Vertex_handle &near,
158 double & dist,
159 double & mindist){
160 dist = _euclid_distance(pref, candidate);
161 return _is_closer_to_with_hint(pref, candidate, near, dist, mindist);
162 }
163
164 /// same as '_is_closer_to' except that 'dist' already contains the
165 /// distance between 'pref' and 'candidate'
_is_closer_to_with_hint(const Point & pref,const Point & candidate,const Vertex_handle & near,const double & dist,double & mindist)166 inline bool _is_closer_to_with_hint(const Point &pref,
167 const Point &candidate,
168 const Vertex_handle &near,
169 const double & dist,
170 double & mindist){
171
172 // check if 'dist', the pre-computed distance between 'candidate'
173 // and 'pref' is smaller than the distance between 'pref' and its
174 // currently registered nearest neighbour 'near' (and update
175 // things if it is)
176 //
177 // Interestingly enough, it has to be pointed out that the use of
178 // 'abs' instead of 'std::abs' returns wrong results (apparently
179 // ints without any compiler warning)
180 //
181 // The (near != NULL) test is there for one single reason: when
182 // checking that a newly inserted point is not closer than a
183 // previous NN, if that distance comparison involves a "nearly
184 // degenerate" distance we need to access near->point. But
185 // sometimes, in the course of RemoveAndAddPoints, its previous NN
186 // has been deleted and its vertex (corresponding to 'near') set
187 // to NULL. This is not a problem as all points having a deleted
188 // point as NN will have their NN explicitly recomputed at the end
189 // of RemoveAndAddPoints so here we should just make sure there is
190 // no crash... that's done by checking (near != NULL)
191 if ((std::abs(dist-mindist)<DISTANCE_FOR_CGAL_CHECKS) &&
192 _is_not_null(near) && // GPS cf. note below about != NULL with clang/CGAL
193 (_euclid_distance(candidate, near->point())<DISTANCE_FOR_CGAL_CHECKS)){
194 // we're in a situation where there might be a rounding issue,
195 // use CGAL's distance computation to get it right
196 //
197 // Note that in the test right above,
198 // (abs(dist-mindist)<1e-12) guarantees that the current
199 // nearest point is not the infinite vertex and thus
200 // nearest->point() is not ill-defined
201 if (_verbose) std::cout << "using CGAL's distance ordering" << std::endl;
202 if (CGAL::compare_distance_to_point(pref, candidate, near->point())!=CGAL::LARGER){
203 mindist = dist;
204 return true;
205 }
206 } else if (dist <= mindist) {
207 // Note that the use of a <= in the above expression (instead of
208 // a strict ordering <) is important in one case: when checking
209 // if a new point is the new NN of one of the points in its
210 // neighbourhood, in case of distances being ==, we are sure
211 // that 'candidate' is in a cell adjacent to 'pref' while it may
212 // no longer be the case for 'near'
213 mindist = dist;
214 return true;
215 }
216
217 return false;
218 }
219
220 /// if a distance between a point and 2 others is smaller than this
221 /// and the distance between the two points is also smaller than this
222 /// then use CGAL to compare the distances.
223 static const double DISTANCE_FOR_CGAL_CHECKS;
224
225
226 /// small routine to check if if a vertex handle is not null.
227 ///
228 /// This is centralised here because of the need for a not-very
229 /// readable workaround for certain CGAL/clang++ compiler
230 /// combinations.
_is_not_null(const Vertex_handle & vh)231 inline static bool _is_not_null(const Vertex_handle & vh) {
232 // as found in issue 2015-07-clang61+CGAL, the form
233 // vh != NULL
234 // appears to be broken with CGAL-3.6 and clang-3.6.0/.1
235 //
236 // The not very "readable"
237 // vh.operator->()
238 // does the job instead
239 return vh.operator->();
240 }
241
242 };
243
244
245 // here follow some inline implementations of the simpler of the
246 // functions defined above
247
NearestNeighbourIndex(const int ii) const248 inline int DnnPlane::NearestNeighbourIndex(const int ii) const {
249 return _supervertex[ii].NNindex;}
250
NearestNeighbourDistance(const int ii) const251 inline double DnnPlane::NearestNeighbourDistance(const int ii) const {
252 return _supervertex[ii].NNdistance;}
253
Valid(const int index) const254 inline bool DnnPlane::Valid(const int index) const {
255 if (index >= 0 && index < static_cast<int>(_supervertex.size())) {
256 return _is_not_null(_supervertex[index].vertex);
257 }
258 else {
259 return false;
260 }
261 }
262
263
etaphi(const int i) const264 inline EtaPhi DnnPlane::etaphi(const int i) const {
265 Point * p = & (_supervertex[i].vertex->point());
266 return EtaPhi(p->x(),p->y()); }
267
eta(const int i) const268 inline double DnnPlane::eta(const int i) const {
269 return _supervertex[i].vertex->point().x(); }
270
phi(const int i) const271 inline double DnnPlane::phi(const int i) const {
272 return _supervertex[i].vertex->point().y(); }
273
274
275 FASTJET_END_NAMESPACE
276
277 #endif // __FASTJET_DNNPLANE_HH__
278
279 #endif // DROP_CGAL
280