30 #ifndef DROP_CGAL // in case we do not have the code for CGAL 34 #include "fastjet/internal/DnnPlane.hh" 37 FASTJET_BEGIN_NAMESPACE
42 DnnPlane::DnnPlane(
const vector<EtaPhi> & input_points,
43 const bool & verbose ) {
46 int n = input_points.size();
51 for (
int i = 0; i < n; i++) {
53 _TR.insert(Point(input_points[i].first, input_points[i].second));
57 _CrashIfVertexPresent(sv.vertex, i);
62 sv.vertex->info() = i;
63 _supervertex.push_back(sv);
67 _TR.infinite_vertex()->info() = INFINITE_VERTEX;
70 for (
int j = 0; j < n; j++) {_SetNearest(j);}
78 void DnnPlane::_CrashIfVertexPresent(
79 const Vertex_handle & vertex,
const int & its_index) {
80 if (!_crash_on_coincidence)
return;
91 if (vertex->info().val() != NEW_VERTEX) {
93 err <<
"ERROR in DnnPlane::_CrashIfVertexPresent" 94 <<endl <<
"Point "<<its_index<<
" coincides with point " 95 <<vertex->info().val() << endl;
96 throw DnnError(err.str());
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) {
121 set<int> NeighbourUnion;
124 set<int> indices_removed;
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;
140 if (_verbose) cout <<
"examining " << vc->info().val() << endl;
141 if (vc->info().val() != INFINITE_VERTEX) {
145 NeighbourUnion.insert(vc->info().val());
146 if (_verbose) cout <<
"inserted " << vc->info().val() << endl;
148 }
while (++vc != done);
152 set<int>::iterator it = NeighbourUnion.begin();
153 cout <<
"Union of neighbours of combined points" << endl;
154 for ( ; it != NeighbourUnion.end(); ++it ) {
160 for (
size_t ir = 0; ir < indices_to_remove.size(); ir++) {
161 int index = indices_to_remove[ir];
165 NeighbourUnion.erase(indices_to_remove[ir]);
170 _TR.remove(_supervertex[index].vertex);
171 _supervertex[index].vertex = NULL;
194 if (indices_to_remove.size() > 0) {
196 face = _TR.incident_faces(
197 _supervertex[*NeighbourUnion.begin()].vertex);}
199 indices_added.clear();
200 indices_of_updated_neighbours.clear();
201 for (
size_t ia = 0; ia < points_to_add.size(); ia++) {
203 _supervertex.push_back(sv);
204 int index = _supervertex.size()-1;
205 indices_added.push_back(index);
207 if (indices_to_remove.size() > 0) {
209 _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
210 points_to_add[ia].second),face);}
212 _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
213 points_to_add[ia].second));
217 _CrashIfVertexPresent(_supervertex[index].vertex, index);
218 _supervertex[index].vertex->info() = index;
227 indices_of_updated_neighbours.push_back(index);
228 _SetAndUpdateNearest(index, indices_of_updated_neighbours);
235 set<int>::iterator it2 = NeighbourUnion.begin();
236 for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
239 if( j != INFINITE_VERTEX ) {
242 if (indices_removed.count(_supervertex[j].NNindex)) {
243 if (_verbose) cout <<
"j " << j << endl;
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;
259 void DnnPlane::_SetNearest (
const int & j) {
260 Vertex_handle current = _supervertex[j].vertex;
264 double mindist = HUGE_DOUBLE;
265 Vertex_handle nearest = _TR.infinite_vertex();
273 if ( vc->info().val() != INFINITE_VERTEX) {
275 if (_verbose) cout << current->info().val() <<
" " << vc->info().val() << endl;
276 dist = _euclid_distance(current->point(), vc->point());
279 if (dist < mindist) {
280 mindist = dist; nearest = vc;
281 if (_verbose) cout <<
"nearer ";
283 if (_verbose) cout << vc->point() <<
"; "<< dist << endl;
285 }
while (++vc != done);
287 _supervertex[j].NNindex = nearest->info().val();
288 _supervertex[j].NNdistance = mindist;
309 void DnnPlane::_SetAndUpdateNearest(
311 vector<int> & indices_of_updated_neighbours) {
313 Vertex_handle current = _supervertex[j].vertex;
317 double mindist = HUGE_DOUBLE;
318 Vertex_handle nearest = _TR.infinite_vertex();
326 if (vc->info().val() != INFINITE_VERTEX) {
327 if (_verbose) cout << current->info().val() <<
" " << vc->info().val() << endl;
329 dist = _euclid_distance(current->point(), vc->point());
331 if (dist < mindist) {
332 mindist = dist; nearest = vc;
333 if (_verbose) cout <<
"nearer ";
336 int vcindx = vc->info().val();
337 if (_verbose) cout << vc->point() <<
"; "<< dist << endl;
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);
347 }
while (++vc != done);
349 _supervertex[j].NNindex = nearest->info().val();
350 _supervertex[j].NNdistance = mindist;
353 FASTJET_END_NAMESPACE
Triangulation::Vertex_circulator Vertex_circulator
CGAL Point structure.