1 | //FJSTARTHEADER
|
---|
2 | // $Id: DnnPlane.hh 3442 2014-07-24 07:20:49Z salam $
|
---|
3 | //
|
---|
4 | // Copyright (c) 2005-2014, 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
|
---|
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
|
---|
115 | 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)
|
---|
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'
|
---|
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 | (near != NULL) &&
|
---|
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 |
|
---|
227 |
|
---|
228 | // here follow some inline implementations of the simpler of the
|
---|
229 | // functions defined above
|
---|
230 |
|
---|
231 | inline int DnnPlane::NearestNeighbourIndex(const int ii) const {
|
---|
232 | return _supervertex[ii].NNindex;}
|
---|
233 |
|
---|
234 | inline double DnnPlane::NearestNeighbourDistance(const int ii) const {
|
---|
235 | return _supervertex[ii].NNdistance;}
|
---|
236 |
|
---|
237 | inline bool DnnPlane::Valid(const int index) const {
|
---|
238 | if (index >= 0 && index < static_cast<int>(_supervertex.size())) {
|
---|
239 | return (_supervertex[index].vertex != NULL);} else {return false;} }
|
---|
240 |
|
---|
241 | inline EtaPhi DnnPlane::etaphi(const int i) const {
|
---|
242 | Point * p = & (_supervertex[i].vertex->point());
|
---|
243 | return EtaPhi(p->x(),p->y()); }
|
---|
244 |
|
---|
245 | inline double DnnPlane::eta(const int i) const {
|
---|
246 | return _supervertex[i].vertex->point().x(); }
|
---|
247 |
|
---|
248 | inline double DnnPlane::phi(const int i) const {
|
---|
249 | return _supervertex[i].vertex->point().y(); }
|
---|
250 |
|
---|
251 |
|
---|
252 | FASTJET_END_NAMESPACE
|
---|
253 |
|
---|
254 | #endif // __FASTJET_DNNPLANE_HH__
|
---|
255 |
|
---|
256 | #endif // DROP_CGAL
|
---|