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) 108 template <
typename real,
typename position,
class distance>
111 static const int version = 1;
114 static const int maxbucket = (2 + ((4 *
sizeof(
real)) /
sizeof(int) >= 2 ?
115 (4 *
sizeof(
real)) /
sizeof(
int) : 2));
175 void Initialize(
const std::vector<position>& pts,
const distance& dist,
177 GEOGRAPHICLIB_STATIC_ASSERT(std::numeric_limits<real>::is_signed,
178 "real must be a signed type");
179 if (!( 0 <= bucket && bucket <= maxbucket ))
181 (
"bucket must lie in [0, 2 + 4*sizeof(real)/sizeof(int)]");
182 if (pts.size() > size_t(std::numeric_limits<int>::max()))
185 std::vector<item> ids(pts.size());
186 for (
int k =
int(ids.size()); k--;)
187 ids[k] = std::make_pair(
real(0), k);
189 std::vector<Node> tree;
190 init(pts, dist, bucket, tree, ids, cost,
191 0,
int(ids.size()),
int(ids.size()/2));
193 _numpoints = int(pts.size());
196 _cost = cost; _c1 = _k = _cmax = 0;
197 _cmin = std::numeric_limits<int>::max();
247 real Search(
const std::vector<position>& pts,
const distance& dist,
248 const position& query,
249 std::vector<int>& ind,
251 real maxdist = std::numeric_limits<real>::max(),
253 bool exhaustive =
true,
254 real tol = 0)
const {
255 std::priority_queue<item> results;
256 if (_numpoints > 0 && _numpoints ==
int(pts.size()) &&
257 k > 0 && maxdist > mindist) {
263 std::priority_queue<item> todo;
264 todo.push(std::make_pair(
real(1),
int(_tree.size()) - 1));
266 while (!todo.empty()) {
267 int n = todo.top().second;
268 real d = -todo.top().first;
270 real tau1 = tau - tol;
272 if (!( n >= 0 && tau1 >= d ))
continue;
273 const Node& current = _tree[n];
275 bool exitflag =
false, leaf = current.index < 0;
276 for (
int i = 0; i < (leaf ? _bucket : 1); ++i) {
277 int index = leaf ? current.leaves[i] : current.index;
278 if (index < 0)
break;
279 dst = dist(pts[index], query);
282 if (dst > mindist && dst <= tau) {
283 if (
int(results.size()) == k) results.pop();
284 results.push(std::make_pair(dst, index));
285 if (
int(results.size()) == k) {
287 tau = results.top().first;
301 if (current.index < 0)
continue;
303 for (
int l = 0; l < 2; ++l) {
304 if (current.data.child[l] >= 0 &&
305 dst + current.data.upper[l] >= mindist) {
306 if (dst < current.data.lower[l]) {
307 d = current.data.lower[l] - dst;
309 todo.push(std::make_pair(-d, current.data.child[l]));
310 }
else if (dst > current.data.upper[l]) {
311 d = dst - current.data.upper[l];
313 todo.push(std::make_pair(-d, current.data.child[l]));
315 todo.push(std::make_pair(
real(1), current.data.child[l]));
322 _mc += (c - omc) / _k;
323 _sc += (c - omc) * (c - _mc);
324 if (c > _cmax) _cmax = c;
325 if (c < _cmin) _cmin = c;
329 ind.resize(results.size());
331 for (
int i =
int(ind.size()); i--;) {
332 ind[i] = int(results.top().second);
333 if (i == 0) d = results.top().first;
362 void Save(std::ostream& os,
bool bin =
true)
const {
363 int realspec = std::numeric_limits<real>::digits *
364 (std::numeric_limits<real>::is_integer ? -1 : 1);
366 char id[] =
"NearestNeighbor_";
373 buf[4] = int(_tree.size());
375 os.write(reinterpret_cast<const char *>(buf), 6 *
sizeof(
int));
376 for (
int i = 0; i < int(_tree.size()); ++i) {
377 const Node& node = _tree[i];
378 os.write(reinterpret_cast<const char *>(&node.index),
sizeof(
int));
379 if (node.index >= 0) {
380 os.write(reinterpret_cast<const char *>(node.data.lower),
382 os.write(reinterpret_cast<const char *>(node.data.upper),
384 os.write(reinterpret_cast<const char *>(node.data.child),
387 os.write(reinterpret_cast<const char *>(node.leaves),
388 _bucket *
sizeof(
int));
392 std::stringstream ostring;
396 ostring.precision(std::numeric_limits<real>::digits10 + 2);
397 ostring << version <<
" " << realspec <<
" " << _bucket <<
" " 398 << _numpoints <<
" " << _tree.size() <<
" " << _cost;
399 for (
int i = 0; i < int(_tree.size()); ++i) {
400 const Node& node = _tree[i];
401 ostring <<
"\n" << node.index;
402 if (node.index >= 0) {
403 for (
int l = 0; l < 2; ++l)
404 ostring <<
" " << node.data.lower[l] <<
" " << node.data.upper[l]
405 <<
" " << node.data.child[l];
407 for (
int l = 0; l < _bucket; ++l)
408 ostring <<
" " << node.leaves[l];
435 void Load(std::istream& is,
bool bin =
true) {
436 int version1, realspec, bucket, numpoints, treesize, cost;
441 if (!(std::strcmp(
id,
"NearestNeighbor_") == 0))
443 is.read(reinterpret_cast<char *>(&version1),
sizeof(
int));
444 is.read(reinterpret_cast<char *>(&realspec),
sizeof(
int));
445 is.read(reinterpret_cast<char *>(&bucket),
sizeof(
int));
446 is.read(reinterpret_cast<char *>(&numpoints),
sizeof(
int));
447 is.read(reinterpret_cast<char *>(&treesize),
sizeof(
int));
448 is.read(reinterpret_cast<char *>(&cost),
sizeof(
int));
450 if (!( is >> version1 >> realspec >> bucket >> numpoints >> treesize
454 if (!( version1 == version ))
456 if (!( realspec == std::numeric_limits<real>::digits *
457 (std::numeric_limits<real>::is_integer ? -1 : 1) ))
459 if (!( 0 <= bucket && bucket <= maxbucket ))
461 if (!( 0 <= treesize && treesize <= numpoints ))
463 std::vector<Node> tree;
464 tree.reserve(treesize);
465 for (
int i = 0; i < treesize; ++i) {
468 is.read(reinterpret_cast<char *>(&node.index),
sizeof(
int));
469 if (node.index >= 0) {
470 is.read(reinterpret_cast<char *>(node.data.lower),
472 is.read(reinterpret_cast<char *>(node.data.upper),
474 is.read(reinterpret_cast<char *>(node.data.child),
477 is.read(reinterpret_cast<char *>(node.leaves),
478 bucket *
sizeof(
int));
479 for (
int l = bucket; l < maxbucket; ++l)
483 if (!( is >> node.index ))
485 if (node.index >= 0) {
486 for (
int l = 0; l < 2; ++l) {
487 if (!( is >> node.data.lower[l] >> node.data.upper[l]
488 >> node.data.child[l] ))
494 for (
int l = 0; l < bucket; ++l) {
495 if (!( is >> node.leaves[l] ))
498 for (
int l = bucket; l < maxbucket; ++l)
502 node.Check(numpoints, treesize, bucket);
503 tree.push_back(node);
506 _numpoints = numpoints;
509 _cost = cost; _c1 = _k = _cmax = 0;
510 _cmin = std::numeric_limits<int>::max();
527 void Statistics(
int& setupcost,
int& numsearches,
int& searchcost,
528 int& mincost,
int& maxcost,
529 double& mean,
double& sd)
const {
530 setupcost = _cost; numsearches = _k; searchcost = _c1;
531 mincost = _cmin; maxcost = _cmax;
532 mean = _mc; sd = std::sqrt(_sc / (_k - 1));
541 _c1 = _k = _cmax = 0;
542 _cmin = std::numeric_limits<int>::max();
548 typedef std::pair<real, int> item;
557 int leaves[maxbucket];
564 for (
int i = 0; i < 2; ++i) {
565 data.lower[i] = data.upper[i] = 0;
571 void Check(
int numpoints,
int treesize,
int bucket)
const {
572 if (!( -1 <= index && index << numpoints ))
575 if (!( -1 <= data.child[0] && data.child[0] < treesize &&
576 -1 <= data.child[1] && data.child[1] < treesize ))
578 if (!( 0 <= data.lower[0] && data.lower[0] <= data.upper[0] &&
579 data.upper[0] <= data.lower[1] &&
580 data.lower[1] <= data.upper[1] ))
586 for (
int l = 0; l < bucket; ++l) {
588 ((l == 0 ? 0 : -1) <= leaves[l] && leaves[l] < numpoints) :
591 start = leaves[l] >= 0;
593 for (
int l = bucket; l < maxbucket; ++l) {
600 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 601 friend class boost::serialization::access;
602 template<
class Archive>
void save(Archive& ar,
const unsigned int)
const {
603 ar & boost::serialization::make_nvp(
"index", index);
605 ar & boost::serialization::make_nvp(
"leaves", leaves);
607 ar & boost::serialization::make_nvp(
"lower", data.lower)
608 & boost::serialization::make_nvp(
"upper", data.upper)
609 & boost::serialization::make_nvp(
"child", data.child);
611 template<
class Archive>
void load(Archive& ar,
const unsigned int) {
612 ar & boost::serialization::make_nvp(
"index", index);
614 ar & boost::serialization::make_nvp(
"leaves", leaves);
616 ar & boost::serialization::make_nvp(
"lower", data.lower)
617 & boost::serialization::make_nvp(
"upper", data.upper)
618 & boost::serialization::make_nvp(
"child", data.child);
620 template<
class Archive>
621 void serialize(Archive& ar,
const unsigned int file_version)
622 { boost::serialization::split_member(ar, *
this, file_version); }
625 #if GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION 626 friend class boost::serialization::access;
627 template<
class Archive>
void save(Archive& ar,
const unsigned)
const {
628 int realspec = std::numeric_limits<real>::digits *
629 (std::numeric_limits<real>::is_integer ? -1 : 1);
632 int version1 = version;
633 ar & boost::serialization::make_nvp(
"version", version1)
634 & boost::serialization::make_nvp(
"realspec", realspec)
635 & boost::serialization::make_nvp(
"bucket", _bucket)
636 & boost::serialization::make_nvp(
"numpoints", _numpoints)
637 & boost::serialization::make_nvp(
"cost", _cost)
638 & boost::serialization::make_nvp(
"tree", _tree);
640 template<
class Archive>
void load(Archive& ar,
const unsigned) {
641 int version1, realspec, bucket, numpoints, cost;
642 ar & boost::serialization::make_nvp(
"version", version1);
643 if (version1 != version)
645 std::vector<Node> tree;
646 ar & boost::serialization::make_nvp(
"realspec", realspec);
647 if (!( realspec == std::numeric_limits<real>::digits *
648 (std::numeric_limits<real>::is_integer ? -1 : 1) ))
650 ar & boost::serialization::make_nvp(
"bucket", bucket);
651 if (!( 0 <= bucket && bucket <= maxbucket ))
653 ar & boost::serialization::make_nvp(
"numpoints", numpoints)
654 & boost::serialization::make_nvp(
"cost", cost)
655 & boost::serialization::make_nvp(
"tree", tree);
656 if (!( 0 <=
int(tree.size()) &&
int(tree.size()) <= numpoints ))
658 for (
int i = 0; i < int(tree.size()); ++i)
659 tree[i].Check(numpoints,
int(tree.size()), bucket);
661 _numpoints = numpoints;
664 _cost = cost; _c1 = _k = _cmax = 0;
665 _cmin = std::numeric_limits<int>::max();
667 template<
class Archive>
668 void serialize(Archive& ar,
const unsigned int file_version)
669 { boost::serialization::split_member(ar, *
this, file_version); }
672 int _numpoints, _bucket, _cost;
673 std::vector<Node> _tree;
675 mutable double _mc, _sc;
676 mutable int _c1, _k, _cmin, _cmax;
678 int init(
const std::vector<position>& pts,
const distance& dist,
int bucket,
679 std::vector<Node>& tree, std::vector<item>& ids,
int& cost,
680 int l,
int u,
int vp) {
686 if (u - l > (bucket == 0 ? 1 : bucket)) {
690 std::swap(ids[l], ids[i]);
692 int m = (u + l + 1) / 2;
694 for (
int k = l + 1; k < u; ++k) {
695 ids[k].first = dist(pts[ids[l].second], pts[ids[k].second]);
699 std::nth_element(ids.begin() + l + 1, ids.begin() + m, ids.begin() + u);
700 node.index = ids[l].second;
702 typename std::vector<item>::iterator
703 t = std::min_element(ids.begin() + l + 1, ids.begin() + m);
704 node.data.lower[0] = t->first;
705 t = std::max_element(ids.begin() + l + 1, ids.begin() + m);
706 node.data.upper[0] = t->first;
709 node.data.child[0] = init(pts, dist, bucket, tree, ids, cost,
710 l + 1, m,
int(t - ids.begin()));
712 typename std::vector<item>::iterator
713 t = std::max_element(ids.begin() + m, ids.begin() + u);
714 node.data.lower[1] = ids[m].first;
715 node.data.upper[1] = t->first;
717 node.data.child[1] = init(pts, dist, bucket, tree, ids, cost,
718 m, u,
int(t - ids.begin()));
721 node.index = ids[l].second;
726 std::sort(ids.begin() + l, ids.begin() + u);
727 for (
int i = l; i < u; ++i)
728 node.leaves[i-l] = ids[i].second;
729 for (
int i = u - l; i < bucket; ++i)
731 for (
int i = bucket; i < maxbucket; ++i)
736 tree.push_back(node);
737 return int(tree.size()) - 1;
744 #if defined(_MSC_VER) 745 # pragma warning (pop) 748 #endif // GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP void ResetStatistics() const
GeographicLib::Math::real real
void Load(std::istream &is, bool bin=true)
void Statistics(int &setupcost, int &numsearches, int &searchcost, int &mincost, int &maxcost, double &mean, double &sd) const
NearestNeighbor(const std::vector< position > &pts, const distance &dist, int bucket=4)
void Initialize(const std::vector< position > &pts, const distance &dist, int bucket=4)
void Save(std::ostream &os, bool bin=true) const
Namespace for GeographicLib.
Exception handling for GeographicLib.
Header for GeographicLib::Constants class.
real Search(const std::vector< position > &pts, const distance &dist, const position &query, std::vector< int > &ind, int k=1, real maxdist=std::numeric_limits< real >::max(), real mindist=-1, bool exhaustive=true, real tol=0) const
Nearest-neighbor calculations.