00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00033
00034 #include<set>
00035 #include<list>
00036 #include "fastjet/internal/DnnPlane.hh"
00037 using namespace std;
00038
00039 FASTJET_BEGIN_NAMESPACE
00040
00041
00044 DnnPlane::DnnPlane(const vector<EtaPhi> & input_points,
00045 const bool & verbose ) {
00046
00047 _verbose = verbose;
00048 int n = input_points.size();
00049
00050
00051
00052 SuperVertex sv;
00053 for (int i = 0; i < n; i++) {
00054 sv.vertex =
00055 _TR.insert(Point(input_points[i].first, input_points[i].second));
00056
00057
00058
00059 _CrashIfVertexPresent(sv.vertex, i);
00060
00061
00062
00063
00064 sv.vertex->info() = i;
00065 _supervertex.push_back(sv);
00066 }
00067
00068
00069 _TR.infinite_vertex()->info() = INFINITE_VERTEX;
00070
00071
00072 for (int j = 0; j < n; j++) {_SetNearest(j);}
00073
00074 }
00075
00076
00077
00080 void DnnPlane::_CrashIfVertexPresent(
00081 const Vertex_handle & vertex, const int & its_index) {
00082 if (!_crash_on_coincidence) return;
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 if (vertex->info().val() != NEW_VERTEX) {
00094 ostringstream err;
00095 err << "ERROR in DnnPlane::_CrashIfVertexPresent"
00096 <<endl << "Point "<<its_index<<" coincides with point "
00097 <<vertex->info().val() << endl;
00098 throw DnnError(err.str());
00099 }
00100 }
00101
00102
00103
00114 void DnnPlane::RemoveAndAddPoints(
00115 const vector<int> & indices_to_remove,
00116 const vector<EtaPhi> & points_to_add,
00117 vector<int> & indices_added,
00118 vector<int> & indices_of_updated_neighbours) {
00119
00120
00121
00122
00123 set<int> NeighbourUnion;
00124
00125
00126 set<int> indices_removed;
00127
00128
00129
00130 for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
00131 int index = indices_to_remove[ir];
00132 indices_removed.insert(index);
00133 if (_verbose) cout << " Starting RemoveAndAddPoints" << endl;
00134 if (_verbose) cout << " point " << index << endl;
00135
00136
00137 Vertex_circulator vc = _TR.incident_vertices(_supervertex[index].vertex);
00138 Vertex_circulator done = vc;
00139 do {
00140
00141
00142 if (_verbose) cout << "examining " << vc->info().val() << endl;
00143 if (vc->info().val() != INFINITE_VERTEX) {
00144
00145
00146
00147 NeighbourUnion.insert(vc->info().val());
00148 if (_verbose) cout << "inserted " << vc->info().val() << endl;
00149 }
00150 } while (++vc != done);
00151 }
00152
00153 if (_verbose) {
00154 set<int>::iterator it = NeighbourUnion.begin();
00155 cout << "Union of neighbours of combined points" << endl;
00156 for ( ; it != NeighbourUnion.end(); ++it ) {
00157 cout << *it << endl;
00158 }
00159 }
00160
00161
00162 for (size_t ir = 0; ir < indices_to_remove.size(); ir++) {
00163 int index = indices_to_remove[ir];
00164
00165
00166
00167 NeighbourUnion.erase(indices_to_remove[ir]);
00168
00169
00170
00171
00172 _TR.remove(_supervertex[index].vertex);
00173 _supervertex[index].vertex = NULL;
00174 }
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 Face_handle face;
00196 if (indices_to_remove.size() > 0) {
00197
00198 face = _TR.incident_faces(
00199 _supervertex[*NeighbourUnion.begin()].vertex);}
00200
00201 indices_added.clear();
00202 indices_of_updated_neighbours.clear();
00203 for (size_t ia = 0; ia < points_to_add.size(); ia++) {
00204 SuperVertex sv;
00205 _supervertex.push_back(sv);
00206 int index = _supervertex.size()-1;
00207 indices_added.push_back(index);
00208
00209 if (indices_to_remove.size() > 0) {
00210
00211 _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
00212 points_to_add[ia].second),face);}
00213 else {
00214 _supervertex[index].vertex = _TR.insert(Point(points_to_add[ia].first,
00215 points_to_add[ia].second));
00216 }
00217
00218
00219 _CrashIfVertexPresent(_supervertex[index].vertex, index);
00220 _supervertex[index].vertex->info() = index;
00221
00222
00223
00224
00225
00226
00227
00228
00229 indices_of_updated_neighbours.push_back(index);
00230 _SetAndUpdateNearest(index, indices_of_updated_neighbours);
00231 }
00232
00233
00234
00235
00236
00237 set<int>::iterator it2 = NeighbourUnion.begin();
00238 for ( ; it2 != NeighbourUnion.end(); ++it2 ) {
00239 int j = *it2;
00240
00241 if( j != INFINITE_VERTEX ) {
00242
00243
00244 if (indices_removed.count(_supervertex[j].NNindex)) {
00245 if (_verbose) cout << "j " << j << endl;
00246 _SetNearest(j);
00247 indices_of_updated_neighbours.push_back(j);
00248 if (_verbose) cout << "NN of " << j << " : "
00249 << _supervertex[j].NNindex
00250 << ", dist = " << _supervertex[j].NNdistance <<endl;
00251 }
00252 }
00253 }
00254
00255 }
00256
00257
00258
00261 void DnnPlane::_SetNearest (const int & j) {
00262 Vertex_handle current = _supervertex[j].vertex;
00263 Vertex_circulator vc = _TR.incident_vertices(current);
00264 Vertex_circulator done = vc;
00265 double dist;
00266 double mindist = HUGE_DOUBLE;
00267 Vertex_handle nearest = _TR.infinite_vertex();
00268
00269
00270
00271
00272
00273
00274 if (vc != NULL) do {
00275 if ( vc->info().val() != INFINITE_VERTEX) {
00276
00277 if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
00278 dist = _euclid_distance(current->point(), vc->point());
00279
00280
00281 if (dist < mindist) {
00282 mindist = dist; nearest = vc;
00283 if (_verbose) cout << "nearer ";
00284 }
00285 if (_verbose) cout << vc->point() << "; "<< dist << endl;
00286 }
00287 } while (++vc != done);
00288
00289 _supervertex[j].NNindex = nearest->info().val();
00290 _supervertex[j].NNdistance = mindist;
00291 }
00292
00293
00311 void DnnPlane::_SetAndUpdateNearest(
00312 const int & j,
00313 vector<int> & indices_of_updated_neighbours) {
00314
00315 Vertex_handle current = _supervertex[j].vertex;
00316 Vertex_circulator vc = _TR.incident_vertices(current);
00317 Vertex_circulator done = vc;
00318 double dist;
00319 double mindist = HUGE_DOUBLE;
00320 Vertex_handle nearest = _TR.infinite_vertex();
00321
00322
00323
00324
00325
00326
00327 if (vc != NULL) do {
00328 if (vc->info().val() != INFINITE_VERTEX) {
00329 if (_verbose) cout << current->info().val() << " " << vc->info().val() << endl;
00330
00331 dist = _euclid_distance(current->point(), vc->point());
00332
00333 if (dist < mindist) {
00334 mindist = dist; nearest = vc;
00335 if (_verbose) cout << "nearer ";
00336 }
00337
00338 int vcindx = vc->info().val();
00339 if (_verbose) cout << vc->point() << "; "<< dist << endl;
00340
00341
00342 if (dist < _supervertex[vcindx].NNdistance) {
00343 if (_verbose) cout << vcindx << "'s NN becomes " << current->info().val() << endl;
00344 _supervertex[vcindx].NNdistance = dist;
00345 _supervertex[vcindx].NNindex = j;
00346 indices_of_updated_neighbours.push_back(vcindx);
00347 }
00348 }
00349 } while (++vc != done);
00350
00351 _supervertex[j].NNindex = nearest->info().val();
00352 _supervertex[j].NNdistance = mindist;
00353 }
00354
00355 FASTJET_END_NAMESPACE
00356
00357 #endif // DROP_CGAL