ClosestPair2D.cc
Go to the documentation of this file.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 #include "fastjet/internal/ClosestPair2D.hh"
00032
00033 #include<limits>
00034 #include<iostream>
00035 #include<iomanip>
00036 #include<algorithm>
00037
00038 FASTJET_BEGIN_NAMESPACE
00039
00040 const unsigned int huge_unsigned = 4294967295U;
00041 const unsigned int twopow31 = 2147483648U;
00042
00043 using namespace std;
00044
00045
00047 void ClosestPair2D::_point2shuffle(Point & point, Shuffle & shuffle,
00048 unsigned int shift) {
00049
00050 Coord2D renorm_point = (point.coord - _left_corner)/_range;
00051
00052
00053 assert(renorm_point.x >=0);
00054 assert(renorm_point.x <=1);
00055 assert(renorm_point.y >=0);
00056 assert(renorm_point.y <=1);
00057
00058 shuffle.x = static_cast<unsigned int>(twopow31 * renorm_point.x) + shift;
00059 shuffle.y = static_cast<unsigned int>(twopow31 * renorm_point.y) + shift;
00060 shuffle.point = &point;
00061 }
00062
00063
00065 bool ClosestPair2D::Shuffle::operator<(const Shuffle & q) const {
00066
00067 if (floor_ln2_less(x ^ q.x, y ^ q.y)) {
00068
00069 return (y < q.y);
00070 } else {
00071
00072 return (x < q.x);
00073 }
00074 }
00075
00076
00077
00078
00079 void ClosestPair2D::_initialize(const std::vector<Coord2D> & positions,
00080 const Coord2D & left_corner,
00081 const Coord2D & right_corner,
00082 unsigned int max_size) {
00083
00084 unsigned int n_positions = positions.size();
00085 assert(max_size >= n_positions);
00086
00087
00088
00089
00090 _points.resize(max_size);
00091
00092
00093 for (unsigned int i = n_positions; i < max_size; i++) {
00094 _available_points.push(&(_points[i]));
00095 }
00096
00097 _left_corner = left_corner;
00098 _range = max((right_corner.x - left_corner.x),
00099 (right_corner.y - left_corner.y));
00100
00101
00102
00103 vector<Shuffle> shuffles(n_positions);
00104 for (unsigned int i = 0; i < n_positions; i++) {
00105
00106 _points[i].coord = positions[i];
00107 _points[i].neighbour_dist2 = numeric_limits<double>::max();
00108 _points[i].review_flag = 0;
00109
00110
00111 _point2shuffle(_points[i], shuffles[i], 0);
00112 }
00113
00114
00115 for (unsigned ishift = 0; ishift < _nshift; ishift++) {
00116
00117
00118 _shifts[ishift] = static_cast<unsigned int>(((twopow31*1.0)*ishift)/_nshift);
00119 if (ishift == 0) {_rel_shifts[ishift] = 0;}
00120 else {_rel_shifts[ishift] = _shifts[ishift] - _shifts[ishift-1];}
00121 }
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 _cp_search_range = 30;
00132 _points_under_review.reserve(_nshift * _cp_search_range);
00133
00134
00135 for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
00136
00137
00138 if (ishift > 0) {
00139 unsigned rel_shift = _rel_shifts[ishift];
00140 for (unsigned int i = 0; i < shuffles.size(); i++) {
00141 shuffles[i] += rel_shift; }
00142 }
00143
00144
00145 sort(shuffles.begin(), shuffles.end());
00146
00147
00148 _trees[ishift] = auto_ptr<Tree>(new Tree(shuffles, max_size));
00149
00150
00151 circulator circ = _trees[ishift]->somewhere(), start=circ;
00152
00153 unsigned int CP_range = min(_cp_search_range, n_positions-1);
00154 do {
00155 Point * this_point = circ->point;
00156
00157 this_point->circ[ishift] = circ;
00158
00159 circulator other = circ;
00160 for (unsigned i=0; i < CP_range; i++) {
00161 ++other;
00162 double dist2 = this_point->distance2(*other->point);
00163 if (dist2 < this_point->neighbour_dist2) {
00164 this_point->neighbour_dist2 = dist2;
00165 this_point->neighbour = other->point;
00166 }
00167 }
00168 } while (++circ != start);
00169
00170 }
00171
00172
00173 vector<double> mindists2(n_positions);
00174 for (unsigned int i = 0; i < n_positions; i++) {
00175 mindists2[i] = _points[i].neighbour_dist2;}
00176
00177 _heap = auto_ptr<MinHeap>(new MinHeap(mindists2, max_size));
00178 }
00179
00180
00181
00182 void ClosestPair2D::closest_pair(unsigned int & ID1, unsigned int & ID2,
00183 double & distance2) const {
00184 ID1 = _heap->minloc();
00185 ID2 = _ID(_points[ID1].neighbour);
00186 distance2 = _points[ID1].neighbour_dist2;
00187 if (ID1 > ID2) swap(ID1,ID2);
00188 }
00189
00190
00191
00192 inline void ClosestPair2D::_add_label(Point * point, unsigned int review_flag) {
00193
00194
00195
00196 if (point->review_flag == 0) _points_under_review.push_back(point);
00197
00198
00199 point->review_flag |= review_flag;
00200 }
00201
00202
00203 inline void ClosestPair2D::_set_label(Point * point, unsigned int review_flag) {
00204
00205
00206
00207 if (point->review_flag == 0) _points_under_review.push_back(point);
00208
00209
00210 point->review_flag = review_flag;
00211 }
00212
00213
00214 void ClosestPair2D::remove(unsigned int ID) {
00215
00216
00217
00218 Point * point_to_remove = & (_points[ID]);
00219
00220
00221 _remove_from_search_tree(point_to_remove);
00222
00223
00224
00225 _deal_with_points_to_review();
00226 }
00227
00228
00229
00230 void ClosestPair2D::_remove_from_search_tree(Point * point_to_remove) {
00231
00232
00233
00234
00235 _available_points.push(point_to_remove);
00236
00237
00238 _set_label(point_to_remove, _remove_heap_entry);
00239
00240
00241
00242
00243
00244 unsigned int CP_range = min(_cp_search_range, size()-1);
00245
00246
00247 for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
00248
00249
00250 circulator removed_circ = point_to_remove->circ[ishift];
00251 circulator right_end = removed_circ.next();
00252
00253 _trees[ishift]->remove(removed_circ);
00254
00255
00256 circulator left_end = right_end, orig_right_end = right_end;
00257 for (unsigned int i = 0; i < CP_range; i++) {left_end--;}
00258
00259 if (size()-1 < _cp_search_range) {
00260
00261
00262
00263
00264 left_end--; right_end--;
00265 }
00266
00267
00268
00269
00270 do {
00271 Point * left_point = left_end->point;
00272
00273
00274
00275 if (left_point->neighbour == point_to_remove) {
00276
00277 _add_label(left_point, _review_neighbour);
00278 } else {
00279
00280 double dist2 = left_point->distance2(*right_end->point);
00281 if (dist2 < left_point->neighbour_dist2) {
00282 left_point->neighbour = right_end->point;
00283 left_point->neighbour_dist2 = dist2;
00284
00285 _add_label(left_point, _review_heap_entry);
00286 }
00287 }
00288 ++right_end;
00289 } while (++left_end != orig_right_end);
00290 }
00291
00292 }
00293
00294
00295
00296 void ClosestPair2D::_deal_with_points_to_review() {
00297
00298
00299
00300 unsigned int CP_range = min(_cp_search_range, size()-1);
00301
00302
00303
00304 while(_points_under_review.size() > 0) {
00305
00306 Point * this_point = _points_under_review.back();
00307
00308 _points_under_review.pop_back();
00309
00310 if (this_point->review_flag & _remove_heap_entry) {
00311
00312 assert(!(this_point->review_flag ^ _remove_heap_entry));
00313 _heap->remove(_ID(this_point));
00314 }
00315
00316 else {
00317 if (this_point->review_flag & _review_neighbour) {
00318 this_point->neighbour_dist2 = numeric_limits<double>::max();
00319
00320 for (unsigned int ishift = 0; ishift < _nshift; ishift++) {
00321 circulator other = this_point->circ[ishift];
00322
00323 for (unsigned i=0; i < CP_range; i++) {
00324 ++other;
00325 double dist2 = this_point->distance2(*other->point);
00326 if (dist2 < this_point->neighbour_dist2) {
00327 this_point->neighbour_dist2 = dist2;
00328 this_point->neighbour = other->point;
00329 }
00330 }
00331 }
00332 }
00333
00334
00335 _heap->update(_ID(this_point), this_point->neighbour_dist2);
00336 }
00337
00338
00339 this_point->review_flag = 0;
00340
00341 }
00342
00343 }
00344
00345
00346 unsigned int ClosestPair2D::insert(const Coord2D & new_coord) {
00347
00348
00349 assert(_available_points.size() > 0);
00350 Point * new_point = _available_points.top();
00351 _available_points.pop();
00352
00353
00354 new_point->coord = new_coord;
00355
00356
00357 _insert_into_search_tree(new_point);
00358
00359
00360
00361 _deal_with_points_to_review();
00362
00363
00364 return _ID(new_point);
00365 }
00366
00367
00368 unsigned int ClosestPair2D::replace(unsigned int ID1, unsigned int ID2,
00369 const Coord2D & position) {
00370
00371
00372 Point * point_to_remove = & (_points[ID1]);
00373 _remove_from_search_tree(point_to_remove);
00374
00375 point_to_remove = & (_points[ID2]);
00376 _remove_from_search_tree(point_to_remove);
00377
00378
00379
00380 Point * new_point = _available_points.top();
00381 _available_points.pop();
00382
00383
00384 new_point->coord = position;
00385
00386
00387 _insert_into_search_tree(new_point);
00388
00389
00390
00391 _deal_with_points_to_review();
00392
00393
00394 return _ID(new_point);
00395
00396 }
00397
00398
00399
00400 void ClosestPair2D::replace_many(
00401 const std::vector<unsigned int> & IDs_to_remove,
00402 const std::vector<Coord2D> & new_positions,
00403 std::vector<unsigned int> & new_IDs) {
00404
00405
00406 for (unsigned int i = 0; i < IDs_to_remove.size(); i++) {
00407 _remove_from_search_tree(& (_points[IDs_to_remove[i]]));
00408 }
00409
00410
00411 new_IDs.resize(0);
00412 for (unsigned int i = 0; i < new_positions.size(); i++) {
00413 Point * new_point = _available_points.top();
00414 _available_points.pop();
00415
00416 new_point->coord = new_positions[i];
00417
00418 _insert_into_search_tree(new_point);
00419
00420 new_IDs.push_back(_ID(new_point));
00421 }
00422
00423
00424
00425 _deal_with_points_to_review();
00426
00427 }
00428
00429
00430
00431 void ClosestPair2D::_insert_into_search_tree(Point * new_point) {
00432
00433
00434 _set_label(new_point, _review_heap_entry);
00435
00436
00437 new_point->neighbour_dist2 = numeric_limits<double>::max();
00438
00439
00440 unsigned int CP_range = min(_cp_search_range, size()-1);
00441
00442 for (unsigned ishift = 0; ishift < _nshift; ishift++) {
00443
00444 Shuffle new_shuffle;
00445 _point2shuffle(*new_point, new_shuffle, _shifts[ishift]);
00446
00447
00448 circulator new_circ = _trees[ishift]->insert(new_shuffle);
00449 new_point->circ[ishift] = new_circ;
00450
00451
00452
00453 circulator right_edge = new_circ; right_edge++;
00454 circulator left_edge = new_circ;
00455 for (unsigned int i = 0; i < CP_range; i++) {left_edge--;}
00456
00457
00458 do {
00459 Point * left_point = left_edge->point;
00460 Point * right_point = right_edge->point;
00461
00462
00463
00464 double new_dist2 = left_point->distance2(*new_point);
00465 if (new_dist2 < left_point->neighbour_dist2) {
00466 left_point->neighbour_dist2 = new_dist2;
00467 left_point->neighbour = new_point;
00468 _add_label(left_point, _review_heap_entry);
00469 }
00470
00471
00472
00473 new_dist2 = new_point->distance2(*right_point);
00474 if (new_dist2 < new_point->neighbour_dist2) {
00475 new_point->neighbour_dist2 = new_dist2;
00476 new_point->neighbour = right_point;
00477 }
00478
00479
00480
00481
00482
00483
00484
00485 if (left_point->neighbour == right_point) {
00486 _add_label(left_point, _review_neighbour);
00487 }
00488
00489
00490 right_edge++;
00491 } while (++left_edge != new_circ);
00492 }
00493 }
00494
00495 FASTJET_END_NAMESPACE
00496