10 #if !defined(GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP) 11 #define GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP 1 25 #if !defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION) 26 #define GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 0 28 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 29 #include <boost/serialization/nvp.hpp> 30 #include <boost/serialization/split_member.hpp> 31 #include <boost/serialization/array.hpp> 32 #include <boost/serialization/vector.hpp> 37 # pragma warning (push) 38 # pragma warning (disable: 4127) 105 template <
typename dist_t,
typename pos_t,
class distfun_t>
108 static const int version = 1;
111 static const int maxbucket = (2 + ((4 *
sizeof(dist_t)) /
sizeof(int) >= 2 ?
112 (4 *
sizeof(dist_t)) /
sizeof(
int) : 2));
172 void Initialize(
const std::vector<pos_t>& pts,
const distfun_t& dist,
174 GEOGRAPHICLIB_STATIC_ASSERT(std::numeric_limits<dist_t>::is_signed,
175 "dist_t must be a signed type");
176 if (!( 0 <= bucket && bucket <= maxbucket ))
178 (
"bucket must lie in [0, 2 + 4*sizeof(dist_t)/sizeof(int)]");
179 if (pts.size() > size_t(std::numeric_limits<int>::max()))
182 std::vector<item> ids(pts.size());
183 for (
int k =
int(ids.size()); k--;)
184 ids[k] = std::make_pair(dist_t(0), k);
186 std::vector<Node> tree;
187 init(pts, dist, bucket, tree, ids, cost,
188 0,
int(ids.size()),
int(ids.size()/2));
190 _numpoints = int(pts.size());
193 _cost = cost; _c1 = _k = _cmax = 0;
194 _cmin = std::numeric_limits<int>::max();
255 dist_t
Search(
const std::vector<pos_t>& pts,
const distfun_t& dist,
257 std::vector<int>& ind,
259 dist_t maxdist = std::numeric_limits<dist_t>::max(),
261 bool exhaustive =
true,
262 dist_t tol = 0)
const {
263 std::priority_queue<item> results;
264 if (_numpoints > 0 && _numpoints ==
int(pts.size()) &&
265 k > 0 && maxdist > mindist) {
267 dist_t tau = maxdist;
271 std::priority_queue<item> todo;
272 todo.push(std::make_pair(dist_t(1),
int(_tree.size()) - 1));
274 while (!todo.empty()) {
275 int n = todo.top().second;
276 dist_t d = -todo.top().first;
278 dist_t tau1 = tau - tol;
280 if (!( n >= 0 && tau1 >= d ))
continue;
281 const Node& current = _tree[n];
283 bool exitflag =
false, leaf = current.index < 0;
284 for (
int i = 0; i < (leaf ? _bucket : 1); ++i) {
285 int index = leaf ? current.leaves[i] : current.index;
286 if (index < 0)
break;
287 dst = dist(pts[index], query);
290 if (dst > mindist && dst <= tau) {
291 if (
int(results.size()) == k) results.pop();
292 results.push(std::make_pair(dst, index));
293 if (
int(results.size()) == k) {
295 tau = results.top().first;
309 if (current.index < 0)
continue;
311 for (
int l = 0; l < 2; ++l) {
312 if (current.data.child[l] >= 0 &&
313 dst + current.data.upper[l] >= mindist) {
314 if (dst < current.data.lower[l]) {
315 d = current.data.lower[l] - dst;
317 todo.push(std::make_pair(-d, current.data.child[l]));
318 }
else if (dst > current.data.upper[l]) {
319 d = dst - current.data.upper[l];
321 todo.push(std::make_pair(-d, current.data.child[l]));
323 todo.push(std::make_pair(dist_t(1), current.data.child[l]));
330 _mc += (c - omc) / _k;
331 _sc += (c - omc) * (c - _mc);
332 if (c > _cmax) _cmax = c;
333 if (c < _cmin) _cmin = c;
337 ind.resize(results.size());
339 for (
int i =
int(ind.size()); i--;) {
340 ind[i] = int(results.top().second);
341 if (i == 0) d = results.top().first;
370 void Save(std::ostream& os,
bool bin =
true)
const {
371 int realspec = std::numeric_limits<dist_t>::digits *
372 (std::numeric_limits<dist_t>::is_integer ? -1 : 1);
374 char id[] =
"NearestNeighbor_";
381 buf[4] = int(_tree.size());
383 os.write(reinterpret_cast<const char *>(buf), 6 *
sizeof(
int));
384 for (
int i = 0; i < int(_tree.size()); ++i) {
385 const Node& node = _tree[i];
386 os.write(reinterpret_cast<const char *>(&node.index),
sizeof(
int));
387 if (node.index >= 0) {
388 os.write(reinterpret_cast<const char *>(node.data.lower),
390 os.write(reinterpret_cast<const char *>(node.data.upper),
392 os.write(reinterpret_cast<const char *>(node.data.child),
395 os.write(reinterpret_cast<const char *>(node.leaves),
396 _bucket *
sizeof(
int));
400 std::stringstream ostring;
403 if (!std::numeric_limits<dist_t>::is_integer)
404 ostring.precision(std::numeric_limits<dist_t>::digits10 + 2);
405 ostring << version <<
" " << realspec <<
" " << _bucket <<
" " 406 << _numpoints <<
" " << _tree.size() <<
" " << _cost;
407 for (
int i = 0; i < int(_tree.size()); ++i) {
408 const Node& node = _tree[i];
409 ostring <<
"\n" << node.index;
410 if (node.index >= 0) {
411 for (
int l = 0; l < 2; ++l)
412 ostring <<
" " << node.data.lower[l] <<
" " << node.data.upper[l]
413 <<
" " << node.data.child[l];
415 for (
int l = 0; l < _bucket; ++l)
416 ostring <<
" " << node.leaves[l];
444 void Load(std::istream& is,
bool bin =
true) {
445 int version1, realspec, bucket, numpoints, treesize, cost;
450 if (!(std::strcmp(
id,
"NearestNeighbor_") == 0))
452 is.read(reinterpret_cast<char *>(&version1),
sizeof(
int));
453 is.read(reinterpret_cast<char *>(&realspec),
sizeof(
int));
454 is.read(reinterpret_cast<char *>(&bucket),
sizeof(
int));
455 is.read(reinterpret_cast<char *>(&numpoints),
sizeof(
int));
456 is.read(reinterpret_cast<char *>(&treesize),
sizeof(
int));
457 is.read(reinterpret_cast<char *>(&cost),
sizeof(
int));
459 if (!( is >> version1 >> realspec >> bucket >> numpoints >> treesize
463 if (!( version1 == version ))
465 if (!( realspec == std::numeric_limits<dist_t>::digits *
466 (std::numeric_limits<dist_t>::is_integer ? -1 : 1) ))
468 if (!( 0 <= bucket && bucket <= maxbucket ))
470 if (!( 0 <= treesize && treesize <= numpoints ))
472 std::vector<Node> tree;
473 tree.reserve(treesize);
474 for (
int i = 0; i < treesize; ++i) {
477 is.read(reinterpret_cast<char *>(&node.index),
sizeof(
int));
478 if (node.index >= 0) {
479 is.read(reinterpret_cast<char *>(node.data.lower),
481 is.read(reinterpret_cast<char *>(node.data.upper),
483 is.read(reinterpret_cast<char *>(node.data.child),
486 is.read(reinterpret_cast<char *>(node.leaves),
487 bucket *
sizeof(
int));
488 for (
int l = bucket; l < maxbucket; ++l)
492 if (!( is >> node.index ))
494 if (node.index >= 0) {
495 for (
int l = 0; l < 2; ++l) {
496 if (!( is >> node.data.lower[l] >> node.data.upper[l]
497 >> node.data.child[l] ))
503 for (
int l = 0; l < bucket; ++l) {
504 if (!( is >> node.leaves[l] ))
507 for (
int l = bucket; l < maxbucket; ++l)
511 node.Check(numpoints, treesize, bucket);
512 tree.push_back(node);
515 _numpoints = numpoints;
518 _cost = cost; _c1 = _k = _cmax = 0;
519 _cmin = std::numeric_limits<int>::max();
531 { t.
Save(os,
false);
return os; }
542 { t.
Load(is,
false);
return is; }
576 void Statistics(
int& setupcost,
int& numsearches,
int& searchcost,
577 int& mincost,
int& maxcost,
578 double& mean,
double& sd)
const {
579 setupcost = _cost; numsearches = _k; searchcost = _c1;
580 mincost = _cmin; maxcost = _cmax;
581 mean = _mc; sd = std::sqrt(_sc / (_k - 1));
590 _c1 = _k = _cmax = 0;
591 _cmin = std::numeric_limits<int>::max();
597 typedef std::pair<dist_t, int> item;
602 dist_t lower[2], upper[2];
607 int leaves[maxbucket];
614 for (
int i = 0; i < 2; ++i) {
615 data.lower[i] = data.upper[i] = 0;
621 void Check(
int numpoints,
int treesize,
int bucket)
const {
622 if (!( -1 <= index && index < numpoints ))
625 if (!( -1 <= data.child[0] && data.child[0] < treesize &&
626 -1 <= data.child[1] && data.child[1] < treesize ))
628 if (!( 0 <= data.lower[0] && data.lower[0] <= data.upper[0] &&
629 data.upper[0] <= data.lower[1] &&
630 data.lower[1] <= data.upper[1] ))
636 for (
int l = 0; l < bucket; ++l) {
638 ((l == 0 ? 0 : -1) <= leaves[l] && leaves[l] < numpoints) :
641 start = leaves[l] >= 0;
643 for (
int l = bucket; l < maxbucket; ++l) {
650 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 651 friend class boost::serialization::access;
652 template<
class Archive>
void save(Archive& ar,
const unsigned int)
const {
653 ar & boost::serialization::make_nvp(
"index", index);
655 ar & boost::serialization::make_nvp(
"leaves", leaves);
657 ar & boost::serialization::make_nvp(
"lower", data.lower)
658 & boost::serialization::make_nvp(
"upper", data.upper)
659 & boost::serialization::make_nvp(
"child", data.child);
661 template<
class Archive>
void load(Archive& ar,
const unsigned int) {
662 ar & boost::serialization::make_nvp(
"index", index);
664 ar & boost::serialization::make_nvp(
"leaves", leaves);
666 ar & boost::serialization::make_nvp(
"lower", data.lower)
667 & boost::serialization::make_nvp(
"upper", data.upper)
668 & boost::serialization::make_nvp(
"child", data.child);
670 template<
class Archive>
671 void serialize(Archive& ar,
const unsigned int file_version)
672 { boost::serialization::split_member(ar, *
this, file_version); }
676 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 677 friend class boost::serialization::access;
678 template<
class Archive>
void save(Archive& ar,
const unsigned)
const {
679 int realspec = std::numeric_limits<dist_t>::digits *
680 (std::numeric_limits<dist_t>::is_integer ? -1 : 1);
683 int version1 = version;
684 ar & boost::serialization::make_nvp(
"version", version1)
685 & boost::serialization::make_nvp(
"realspec", realspec)
686 & boost::serialization::make_nvp(
"bucket", _bucket)
687 & boost::serialization::make_nvp(
"numpoints", _numpoints)
688 & boost::serialization::make_nvp(
"cost", _cost)
689 & boost::serialization::make_nvp(
"tree", _tree);
691 template<
class Archive>
void load(Archive& ar,
const unsigned) {
692 int version1, realspec, bucket, numpoints, cost;
693 ar & boost::serialization::make_nvp(
"version", version1);
694 if (version1 != version)
696 std::vector<Node> tree;
697 ar & boost::serialization::make_nvp(
"realspec", realspec);
698 if (!( realspec == std::numeric_limits<dist_t>::digits *
699 (std::numeric_limits<dist_t>::is_integer ? -1 : 1) ))
701 ar & boost::serialization::make_nvp(
"bucket", bucket);
702 if (!( 0 <= bucket && bucket <= maxbucket ))
704 ar & boost::serialization::make_nvp(
"numpoints", numpoints)
705 & boost::serialization::make_nvp(
"cost", cost)
706 & boost::serialization::make_nvp(
"tree", tree);
707 if (!( 0 <=
int(tree.size()) &&
int(tree.size()) <= numpoints ))
709 for (
int i = 0; i < int(tree.size()); ++i)
710 tree[i].Check(numpoints,
int(tree.size()), bucket);
712 _numpoints = numpoints;
715 _cost = cost; _c1 = _k = _cmax = 0;
716 _cmin = std::numeric_limits<int>::max();
718 template<
class Archive>
719 void serialize(Archive& ar,
const unsigned int file_version)
720 { boost::serialization::split_member(ar, *
this, file_version); }
723 int _numpoints, _bucket, _cost;
724 std::vector<Node> _tree;
726 mutable double _mc, _sc;
727 mutable int _c1, _k, _cmin, _cmax;
729 int init(
const std::vector<pos_t>& pts,
const distfun_t& dist,
int bucket,
730 std::vector<Node>& tree, std::vector<item>& ids,
int& cost,
731 int l,
int u,
int vp) {
737 if (u - l > (bucket == 0 ? 1 : bucket)) {
743 int m = (u + l + 1) / 2;
745 for (
int k = l + 1; k < u; ++k) {
746 ids[k].first = dist(pts[ids[l].second], pts[ids[k].second]);
750 std::nth_element(ids.begin() + l + 1, ids.begin() + m, ids.begin() + u);
751 node.index = ids[l].second;
753 typename std::vector<item>::iterator
754 t = std::min_element(ids.begin() + l + 1, ids.begin() + m);
755 node.data.lower[0] = t->first;
756 t = std::max_element(ids.begin() + l + 1, ids.begin() + m);
757 node.data.upper[0] = t->first;
760 node.data.child[0] = init(pts, dist, bucket, tree, ids, cost,
761 l + 1, m,
int(t - ids.begin()));
763 typename std::vector<item>::iterator
764 t = std::max_element(ids.begin() + m, ids.begin() + u);
765 node.data.lower[1] = ids[m].first;
766 node.data.upper[1] = t->first;
768 node.data.child[1] = init(pts, dist, bucket, tree, ids, cost,
769 m, u,
int(t - ids.begin()));
772 node.index = ids[l].second;
777 std::sort(ids.begin() + l, ids.begin() + u);
778 for (
int i = l; i < u; ++i)
779 node.leaves[i-l] = ids[i].second;
780 for (
int i = u - l; i < bucket; ++i)
782 for (
int i = bucket; i < maxbucket; ++i)
787 tree.push_back(node);
788 return int(tree.size()) - 1;
807 template <
typename dist_t,
typename pos_t,
class distfun_t>
815 #if defined(_MSC_VER) 816 # pragma warning (pop) 819 #endif // GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP
void Statistics(int &setupcost, int &numsearches, int &searchcost, int &mincost, int &maxcost, double &mean, double &sd) const
NearestNeighbor(const std::vector< pos_t > &pts, const distfun_t &dist, int bucket=4)
void Load(std::istream &is, bool bin=true)
dist_t Search(const std::vector< pos_t > &pts, const distfun_t &dist, const pos_t &query, std::vector< int > &ind, int k=1, dist_t maxdist=std::numeric_limits< dist_t >::max(), dist_t mindist=-1, bool exhaustive=true, dist_t tol=0) const
Namespace for GeographicLib.
friend std::ostream & operator<<(std::ostream &os, const NearestNeighbor &t)
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
friend std::istream & operator>>(std::istream &is, NearestNeighbor &t)
Exception handling for GeographicLib.
Header for GeographicLib::Constants class.
void Initialize(const std::vector< pos_t > &pts, const distfun_t &dist, int bucket=4)
void swap(NearestNeighbor &t)
Nearest-neighbor calculations.
void ResetStatistics() const
void Save(std::ostream &os, bool bin=true) const