[1] | 1 | //STARTHEADER |
---|
| 2 | // $Id: DnnPlane.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
| 3 | // |
---|
| 4 | // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use |
---|
| 16 | // FastJet as part of work towards a scientific publication, please |
---|
| 17 | // include a citation to the FastJet paper. |
---|
| 18 | // |
---|
| 19 | // FastJet is distributed in the hope that it will be useful, |
---|
| 20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 22 | // GNU General Public License for more details. |
---|
| 23 | // |
---|
| 24 | // You should have received a copy of the GNU General Public License |
---|
| 25 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>. |
---|
| 26 | //---------------------------------------------------------------------- |
---|
| 27 | //ENDHEADER |
---|
| 28 | |
---|
| 29 | |
---|
| 30 | #ifndef DROP_CGAL // in case we do not have the code for CGAL |
---|
| 31 | |
---|
| 32 | #include<set> |
---|
| 33 | #include<list> |
---|
| 34 | #include "fastjet/internal/DnnPlane.hh" |
---|
| 35 | using namespace std; |
---|
| 36 | |
---|
| 37 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | /// Initialiser from a set of points on an Eta-Phi plane, where both |
---|
| 41 | /// eta and phi can have arbitrary ranges |
---|
| 42 | DnnPlane::DnnPlane(const vector<EtaPhi> & input_points, |
---|
| 43 | const bool & verbose ) { |
---|
| 44 | |
---|
| 45 | _verbose = verbose; |
---|
| 46 | int n = input_points.size(); |
---|
| 47 | |
---|
| 48 | // construct Voronoi diagram in such a way as to get the vertex handles |
---|
| 49 | // and remember to set CGAL info with the index of the vertex |
---|
| 50 | SuperVertex sv; |
---|
| 51 | for (int i = 0; i < n; i++) { |
---|
| 52 | sv.vertex = |
---|
| 53 | _TR.insert(Point(input_points[i].first, input_points[i].second)); |
---|
| 54 | |
---|
| 55 | // we are not up to dealing with coincident vertices, so make |
---|
| 56 | // sure the user knows! |
---|
| 57 | _CrashIfVertexPresent(sv.vertex, i); |
---|
| 58 | |
---|
| 59 | // we need to assicate an index to each vertex -- thus when we get |
---|
| 60 | // a vertex (e.g. as a nearest neighbour) from CGAL, we will be |
---|
| 61 | // able to figure out which particle it corresponded to. |
---|
| 62 | sv.vertex->info() = i; |
---|
| 63 | _supervertex.push_back(sv); |
---|
| 64 | } |
---|
| 65 | |
---|
| 66 | // label infinite vertex info with negative index |
---|
| 67 | _TR.infinite_vertex()->info() = INFINITE_VERTEX; |
---|
| 68 | |
---|
| 69 | // set up the structure that holds nearest distances and neighbours |
---|
| 70 | for (int j = 0; j < n; j++) {_SetNearest(j);} |
---|
| 71 | |
---|
| 72 | } |
---|
| 73 | |
---|
| 74 | |
---|
| 75 | //---------------------------------------------------------------------- |
---|
| 76 | /// Crashes if the given vertex handle already exists. Otherwise |
---|
| 77 | /// it does the bookkeeping for future such tests |
---|
| 78 | void DnnPlane::_CrashIfVertexPresent( |
---|
| 79 | const Vertex_handle & vertex, const int & its_index) { |
---|
| 80 | if (!_crash_on_coincidence) return; |
---|
| 81 | |
---|
| 82 | // vertices that do not have the same geometric position as any |
---|
| 83 | // other vertex so far added have info().val() == NEW_VERTEX -- this |
---|
| 84 | // is ensured by the InitialisedInt class, which forms the "info" |
---|
| 85 | // part of our |
---|
| 86 | // CGAL::Triangulation_vertex_base_with_info_2<InitialisedInt,K> |
---|
| 87 | // |
---|
| 88 | // If the vertex coincides with one that already exists, then |
---|
| 89 | // info().val() it's info().val() will have been updated (in |
---|
| 90 | // DNN:DNN) to be equal to a vertex "index". |
---|
| 91 | if (vertex->info().val() != NEW_VERTEX) { |
---|
| 92 | ostringstream err; |
---|
| 93 | err << "ERROR in DnnPlane::_CrashIfVertexPresent" |
---|
| 94 | <<endl << "Point "<<its_index<<" coincides with point " |
---|
| 95 | <<vertex->info().val() << endl; |
---|
| 96 | throw DnnError(err.str()); |
---|
| 97 | } |
---|
| 98 | } |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | //---------------------------------------------------------------------- |
---|
| 102 | /// remove the points labelled by the vector indices_to_remove, and |
---|
| 103 | /// add the points specified by the vector points_to_add |
---|
| 104 | /// (corresponding indices will be calculated automatically); the |
---|
| 105 | /// idea behind this routine is that the points to be added will |
---|
| 106 | /// somehow be close to the one or other of the points being removed |
---|
| 107 | /// and this can be used by the implementation to provide hints for |
---|
| 108 | /// inserting the new points in whatever structure it is using. In a |
---|
| 109 | /// kt-algorithm the points being added will be a result of a |
---|
| 110 | /// combination of the points to be removed -- hence the proximity |
---|
| 111 | /// is (more or less) guaranteed. |
---|
| 112 | void DnnPlane::RemoveAndAddPoints( |
---|
| 113 | const vector<int> & indices_to_remove, |
---|
| 114 | const vector<EtaPhi> & points_to_add, |
---|
| 115 | vector<int> & indices_added, |
---|
| 116 | vector<int> & indices_of_updated_neighbours) { |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | // build set of UNION of Voronoi neighbours of a pair of nearest |
---|
| 120 | // neighbours |
---|
| 121 | set<int> NeighbourUnion; |
---|
| 122 | // later on it will be convenient to have access to a set (rather |
---|
| 123 | // than vector) of indices being removed |
---|
| 124 | set<int> indices_removed; |
---|
| 125 | |
---|
| 126 | // for each of the indices to be removed add the voronoi neighbourhood to |
---|
| 127 | // the NeighbourUnion set. |
---|
| 128 | for (size_t ir = 0; ir < indices_to_remove.size(); ir++) { |
---|
| 129 | int index = indices_to_remove[ir]; |
---|
| 130 | indices_removed.insert(index); |
---|
| 131 | if (_verbose) cout << " Starting RemoveAndAddPoints" << endl; |
---|
| 132 | if (_verbose) cout << " point " << index << endl; |
---|
| 133 | // have a circulators that will go round the Voronoi neighbours of |
---|
| 134 | // _supervertex[index1].vertex |
---|
| 135 | Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex); |
---|
| 136 | Vertex_circulator done = vc; |
---|
| 137 | do { |
---|
| 138 | // if a neighbouring vertex not the infinite vertex, then add it |
---|
| 139 | // to our union of neighbouring vertices. |
---|
| 140 | if (_verbose) cout << "examining " << vc->info().val() << endl; |
---|
| 141 | if (vc->info().val() != INFINITE_VERTEX) { |
---|
| 142 | // NB: from it=1 onwards occasionally it might already have |
---|
| 143 | // been inserted -- but double insertion still leaves only one |
---|
| 144 | // copy in the set, so there's no problem |
---|
| 145 | NeighbourUnion.insert(vc->info().val()); |
---|
| 146 | if (_verbose) cout << "inserted " << vc->info().val() << endl; |
---|
| 147 | } |
---|
| 148 | } while (++vc != done); |
---|
| 149 | } |
---|
| 150 | |
---|
| 151 | if (_verbose) { |
---|
| 152 | set<int>::iterator it = NeighbourUnion.begin(); |
---|
| 153 | cout << "Union of neighbours of combined points" << endl; |
---|
| 154 | for ( ; it != NeighbourUnion.end(); ++it ) { |
---|
| 155 | cout << *it << endl; |
---|
| 156 | } |
---|
| 157 | } |
---|
| 158 | |
---|
| 159 | // update set, triangulation and supervertex info |
---|
| 160 | for (size_t ir = 0; ir < indices_to_remove.size(); ir++) { |
---|
| 161 | int index = indices_to_remove[ir]; |
---|
| 162 | |
---|
| 163 | // NeighbourUnion should not contain the points to be removed |
---|
| 164 | // (because later we will assume they still exist). |
---|
| 165 | NeighbourUnion.erase(indices_to_remove[ir]); |
---|
| 166 | |
---|
| 167 | // points to be removed should also be eliminated from the |
---|
| 168 | // triangulation and the supervertex structure should be updated |
---|
| 169 | // to reflect the fact that the points are no longer valid. |
---|
| 170 | _TR.remove(_supervertex[index].vertex); |
---|
| 171 | _supervertex[index].vertex = NULL; |
---|
| 172 | } |
---|
| 173 | |
---|
| 174 | // add new point: give a "hint" to the inserter that |
---|
| 175 | // the new point should be added close to old points -- the easiest way |
---|
| 176 | // of getting this is to take a point from the NeighbourUnion AFTER we have |
---|
| 177 | // removed point1, point2, and to get one of its incident faces. |
---|
| 178 | // |
---|
| 179 | // This hinting improves speed by c. 25% for 10^4 points because it |
---|
| 180 | // avoids the need for a costly (sqrt{N}) location step (at least |
---|
| 181 | // with a non-hierarchical triangulation -- with a hierarchical one, |
---|
| 182 | // this step could be done away with, though there will still be a |
---|
| 183 | // cost of O(ln N) to pay. |
---|
| 184 | // |
---|
| 185 | // For some reason inserting the point before the two removals |
---|
| 186 | // slows things down by c. 25%. This importance of the order |
---|
| 187 | // is not understood. |
---|
| 188 | // |
---|
| 189 | // At some point it might be worth trying to select the "nearest" |
---|
| 190 | // of the various points in the neighbour union to avoid large |
---|
| 191 | // steps in cases where we have 0..2pi periodicity and the first member |
---|
| 192 | // of the neighbour union happens to be on the wrong side. |
---|
| 193 | Face_handle face; |
---|
| 194 | if (indices_to_remove.size() > 0) { |
---|
| 195 | // face can only be found if there were points to remove in first place |
---|
| 196 | face = _TR.incident_faces( |
---|
| 197 | _supervertex[*NeighbourUnion.begin()].vertex);} |
---|
| 198 | // make sure the output arrays are empty |
---|
| 199 | indices_added.clear(); |
---|
| 200 | indices_of_updated_neighbours.clear(); |
---|
| 201 | for (size_t ia = 0; ia < points_to_add.size(); ia++) { |
---|
| 202 | SuperVertex sv; |
---|
| 203 | _supervertex.push_back(sv); |
---|
| 204 | int index = _supervertex.size()-1; |
---|
| 205 | indices_added.push_back(index); |
---|
| 206 | |
---|
| 207 | if (indices_to_remove.size() > 0) { |
---|
| 208 | // be careful of using face (for location hinting) only when it exists |
---|
| 209 | _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, |
---|
| 210 | points_to_add[ia].second),face);} |
---|
| 211 | else { |
---|
| 212 | _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first, |
---|
| 213 | points_to_add[ia].second)); |
---|
| 214 | } |
---|
| 215 | // we are not up to dealing with coincident vertices, so make |
---|
| 216 | // sure the user knows! |
---|
| 217 | _CrashIfVertexPresent(_supervertex[index].vertex, index); |
---|
| 218 | _supervertex[index].vertex->info() = index; |
---|
| 219 | |
---|
| 220 | // first find nearest neighbour of "newpoint" (shorthand for |
---|
| 221 | // _supervertex[index].vertex); while we're at it, for each of the |
---|
| 222 | // voronoi neighbours, "D", of newpoint, examine whether newpoint is |
---|
| 223 | // closer to "D" than D's current nearest neighbour -- when this |
---|
| 224 | // occurs, put D into indices_of_updated_neighbours. |
---|
| 225 | // |
---|
| 226 | // manually put newpoint on indices_of_updated_neighbours |
---|
| 227 | indices_of_updated_neighbours.push_back(index); |
---|
| 228 | _SetAndUpdateNearest(index, indices_of_updated_neighbours); |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | // for Voronoi neighbours j of any of the removed points for which |
---|
| 232 | // one of those removed points was the nearest neighbour, |
---|
| 233 | // redetermine the nearest neighbour of j and add j onto the vector |
---|
| 234 | // of indices_of_updated_neighbours. |
---|
| 235 | set<int>::iterator it2 = NeighbourUnion.begin(); |
---|
| 236 | for ( ; it2 != NeighbourUnion.end(); ++it2 ) { |
---|
| 237 | int j = *it2; |
---|
| 238 | // the if avoids the vertex at infinity, which gets a negative index |
---|
| 239 | if( j != INFINITE_VERTEX ) { |
---|
| 240 | // this is where we check if the nearest neighbour of j was one |
---|
| 241 | // of the removed points |
---|
| 242 | if (indices_removed.count(_supervertex[j].NNindex)) { |
---|
| 243 | if (_verbose) cout << "j " << j << endl; |
---|
| 244 | _SetNearest(j); |
---|
| 245 | indices_of_updated_neighbours.push_back(j); |
---|
| 246 | if (_verbose) cout << "NN of " << j << " : " |
---|
| 247 | << _supervertex[j].NNindex |
---|
| 248 | << ", dist = " << _supervertex[j].NNdistance <<endl; |
---|
| 249 | } |
---|
| 250 | } |
---|
| 251 | } |
---|
| 252 | |
---|
| 253 | } |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | //---------------------------------------------------------------------- |
---|
| 257 | /// Determines the index and distance of the nearest neighbour to |
---|
| 258 | /// point j and puts the information into the _supervertex entry for j. |
---|
| 259 | void DnnPlane::_SetNearest (const int & j) { |
---|
| 260 | Vertex_handle current = _supervertex[j].vertex; |
---|
| 261 | Vertex_circulator vc = _TR.incident_vertices(current); |
---|
| 262 | Vertex_circulator done = vc; |
---|
| 263 | double dist; |
---|
| 264 | double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double? |
---|
| 265 | Vertex_handle nearest = _TR.infinite_vertex(); |
---|
| 266 | |
---|
| 267 | // when there is only one finite point left in the triangulation, |
---|
| 268 | // there are no triangles. Presumably this is why voronoi returns |
---|
| 269 | // NULL for the incident vertex circulator. Check if this is |
---|
| 270 | // happening before circulating over it... (Otherwise it crashes |
---|
| 271 | // when looking for neighbours of last point) |
---|
| 272 | if (vc != NULL) do { |
---|
| 273 | if ( vc->info().val() != INFINITE_VERTEX) { |
---|
| 274 | // find distance between j and its Voronoi neighbour (vc) |
---|
| 275 | if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl; |
---|
| 276 | dist = _euclid_distance(current->point(), vc->point()); |
---|
| 277 | // check if j is closer to vc than vc's currently registered |
---|
| 278 | // nearest neighbour (and update things if it is) |
---|
| 279 | if (dist < mindist) { |
---|
| 280 | mindist = dist; nearest = vc; |
---|
| 281 | if (_verbose) cout << "nearer "; |
---|
| 282 | } |
---|
| 283 | if (_verbose) cout << vc->point() << "; "<< dist << endl; |
---|
| 284 | } |
---|
| 285 | } while (++vc != done); // move on to next Voronoi neighbour |
---|
| 286 | // set j's supervertex info about nearest neighbour |
---|
| 287 | _supervertex[j].NNindex = nearest->info().val(); |
---|
| 288 | _supervertex[j].NNdistance = mindist; |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | //---------------------------------------------------------------------- |
---|
| 292 | /// Determines and stores the nearest neighbour of j, and where |
---|
| 293 | /// necessary update the nearest-neighbour info of Voronoi neighbours |
---|
| 294 | /// of j; |
---|
| 295 | /// |
---|
| 296 | /// For each voronoi neighbour D of j if the distance between j and D |
---|
| 297 | /// is less than D's own nearest neighbour, then update the |
---|
| 298 | /// nearest-neighbour info in D; push D's index onto |
---|
| 299 | /// indices_of_updated_neighbours |
---|
| 300 | /// |
---|
| 301 | /// Note that j is NOT pushed onto indices_of_updated_neighbours -- |
---|
| 302 | /// if you want it there, put it there yourself. |
---|
| 303 | /// |
---|
| 304 | /// NB: note that we have _SetAndUpdateNearest as a completely |
---|
| 305 | /// separate routine from _SetNearest because we want to |
---|
| 306 | /// use one single ciruclation over voronoi neighbours to find the |
---|
| 307 | /// nearest neighbour and to update the voronoi neighbours if need |
---|
| 308 | /// be. |
---|
| 309 | void DnnPlane::_SetAndUpdateNearest( |
---|
| 310 | const int & j, |
---|
| 311 | vector<int> & indices_of_updated_neighbours) { |
---|
| 312 | |
---|
| 313 | Vertex_handle current = _supervertex[j].vertex; |
---|
| 314 | Vertex_circulator vc = _TR.incident_vertices(current); |
---|
| 315 | Vertex_circulator done = vc; |
---|
| 316 | double dist; |
---|
| 317 | double mindist = HUGE_DOUBLE; // change this to "HUGE" or max_double? |
---|
| 318 | Vertex_handle nearest = _TR.infinite_vertex(); |
---|
| 319 | |
---|
| 320 | // when there is only one finite point left in the triangulation, |
---|
| 321 | // there are no triangles. Presumably this is why voronoi returns |
---|
| 322 | // NULL for the incident vertex circulator. Check if this is |
---|
| 323 | // happening before circulating over it... (Otherwise it crashes |
---|
| 324 | // when looking for neighbours of last point) |
---|
| 325 | if (vc != NULL) do { |
---|
| 326 | if (vc->info().val() != INFINITE_VERTEX) { |
---|
| 327 | if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl; |
---|
| 328 | // find distance between j and its Voronoi neighbour (vc) |
---|
| 329 | dist = _euclid_distance(current->point(), vc->point()); |
---|
| 330 | // update the mindist if we are closer than anything found so far |
---|
| 331 | if (dist < mindist) { |
---|
| 332 | mindist = dist; nearest = vc; |
---|
| 333 | if (_verbose) cout << "nearer "; |
---|
| 334 | } |
---|
| 335 | // find index corresponding to vc for easy manipulation |
---|
| 336 | int vcindx = vc->info().val(); |
---|
| 337 | if (_verbose) cout << vc->point() << "; "<< dist << endl; |
---|
| 338 | // check if j is closer to vc than vc's currently registered |
---|
| 339 | // nearest neighbour (and update things if it is) |
---|
| 340 | if (dist < _supervertex[vcindx].NNdistance) { |
---|
| 341 | if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl; |
---|
| 342 | _supervertex[vcindx].NNdistance = dist; |
---|
| 343 | _supervertex[vcindx].NNindex = j; |
---|
| 344 | indices_of_updated_neighbours.push_back(vcindx); |
---|
| 345 | } |
---|
| 346 | } |
---|
| 347 | } while (++vc != done); // move on to next Voronoi neighbour |
---|
| 348 | // set j's supervertex info about nearest neighbour |
---|
| 349 | _supervertex[j].NNindex = nearest->info().val(); |
---|
| 350 | _supervertex[j].NNdistance = mindist; |
---|
| 351 | } |
---|
| 352 | |
---|
| 353 | FASTJET_END_NAMESPACE |
---|
| 354 | |
---|
| 355 | #endif // DROP_CGAL |
---|