10 #if !defined(GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP)
11 #define GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP 1
25 #if defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION) && \
26 GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
27 #include <boost/serialization/nvp.hpp>
28 #include <boost/serialization/split_member.hpp>
29 #include <boost/serialization/array.hpp>
30 #include <boost/serialization/vector.hpp>
35 # pragma warning (push)
36 # pragma warning (disable: 4127)
103 template <
typename dist_t,
typename pos_t,
class distfun_t>
106 static const int version = 1;
109 static const int maxbucket =
110 (2 + ((4 *
sizeof(dist_t)) /
sizeof(int) >= 2 ?
111 (4 *
sizeof(dist_t)) /
sizeof(
int) : 2));
171 void Initialize(
const std::vector<pos_t>& pts,
const distfun_t& dist,
173 GEOGRAPHICLIB_STATIC_ASSERT(std::numeric_limits<dist_t>::is_signed,
174 "dist_t must be a signed type");
175 if (!( 0 <= bucket && bucket <= maxbucket ))
177 (
"bucket must lie in [0, 2 + 4*sizeof(dist_t)/sizeof(int)]");
178 if (pts.size() > size_t(std::numeric_limits<int>::max()))
181 std::vector<item> ids(pts.size());
182 for (
int k =
int(ids.size()); k--;)
183 ids[k] = std::make_pair(dist_t(0), k);
185 std::vector<Node> tree;
186 init(pts, dist, bucket, tree, ids, cost,
187 0,
int(ids.size()),
int(ids.size()/2));
189 _numpoints = int(pts.size());
192 _cost = cost; _c1 = _k = _cmax = 0;
193 _cmin = std::numeric_limits<int>::max();
259 dist_t
Search(
const std::vector<pos_t>& pts,
const distfun_t& dist,
261 std::vector<int>& ind,
263 dist_t maxdist = std::numeric_limits<dist_t>::max(),
265 bool exhaustive =
true,
266 dist_t tol = 0)
const {
267 if (_numpoints !=
int(pts.size()))
269 std::priority_queue<item> results;
270 if (_numpoints > 0 && k > 0 && maxdist > mindist) {
272 dist_t tau = maxdist;
276 std::priority_queue<item> todo;
277 todo.push(std::make_pair(dist_t(1),
int(_tree.size()) - 1));
279 while (!todo.empty()) {
280 int n = todo.top().second;
281 dist_t d = -todo.top().first;
283 dist_t tau1 = tau - tol;
285 if (!( n >= 0 && tau1 >= d ))
continue;
286 const Node& current = _tree[n];
288 bool exitflag =
false, leaf = current.index < 0;
289 for (
int i = 0; i < (leaf ? _bucket : 1); ++i) {
290 int index = leaf ? current.leaves[i] : current.index;
291 if (index < 0)
break;
292 dst = dist(pts[index], query);
295 if (dst > mindist && dst <= tau) {
296 if (
int(results.size()) == k) results.pop();
297 results.push(std::make_pair(dst, index));
298 if (
int(results.size()) == k) {
300 tau = results.top().first;
314 if (current.index < 0)
continue;
316 for (
int l = 0; l < 2; ++l) {
317 if (current.data.child[l] >= 0 &&
318 dst + current.data.upper[l] >= mindist) {
319 if (dst < current.data.lower[l]) {
320 d = current.data.lower[l] - dst;
322 todo.push(std::make_pair(-d, current.data.child[l]));
323 }
else if (dst > current.data.upper[l]) {
324 d = dst - current.data.upper[l];
326 todo.push(std::make_pair(-d, current.data.child[l]));
328 todo.push(std::make_pair(dist_t(1), current.data.child[l]));
335 _mc += (c - omc) / _k;
336 _sc += (c - omc) * (c - _mc);
337 if (c > _cmax) _cmax = c;
338 if (c < _cmin) _cmin = c;
342 ind.resize(results.size());
344 for (
int i =
int(ind.size()); i--;) {
345 ind[i] = int(results.top().second);
346 if (i == 0) d = results.top().first;
375 void Save(std::ostream& os,
bool bin =
true)
const {
376 int realspec = std::numeric_limits<dist_t>::digits *
377 (std::numeric_limits<dist_t>::is_integer ? -1 : 1);
379 char id[] =
"NearestNeighbor_";
386 buf[4] = int(_tree.size());
388 os.write(reinterpret_cast<const char *>(buf), 6 *
sizeof(
int));
389 for (
int i = 0; i < int(_tree.size()); ++i) {
390 const Node& node = _tree[i];
391 os.write(reinterpret_cast<const char *>(&node.index),
sizeof(
int));
392 if (node.index >= 0) {
393 os.write(reinterpret_cast<const char *>(node.data.lower),
395 os.write(reinterpret_cast<const char *>(node.data.upper),
397 os.write(reinterpret_cast<const char *>(node.data.child),
400 os.write(reinterpret_cast<const char *>(node.leaves),
401 _bucket *
sizeof(
int));
405 std::stringstream ostring;
408 if (!std::numeric_limits<dist_t>::is_integer) {
409 static const int prec
410 = int(std::ceil(std::numeric_limits<dist_t>::digits *
411 std::log10(2.0) + 1));
412 ostring.precision(prec);
414 ostring << version <<
" " << realspec <<
" " << _bucket <<
" "
415 << _numpoints <<
" " << _tree.size() <<
" " << _cost;
416 for (
int i = 0; i < int(_tree.size()); ++i) {
417 const Node& node = _tree[i];
418 ostring <<
"\n" << node.index;
419 if (node.index >= 0) {
420 for (
int l = 0; l < 2; ++l)
421 ostring <<
" " << node.data.lower[l] <<
" " << node.data.upper[l]
422 <<
" " << node.data.child[l];
424 for (
int l = 0; l < _bucket; ++l)
425 ostring <<
" " << node.leaves[l];
453 void Load(std::istream& is,
bool bin =
true) {
454 int version1, realspec, bucket, numpoints, treesize, cost;
459 if (!(std::strcmp(
id,
"NearestNeighbor_") == 0))
461 is.read(reinterpret_cast<char *>(&version1),
sizeof(
int));
462 is.read(reinterpret_cast<char *>(&realspec),
sizeof(
int));
463 is.read(reinterpret_cast<char *>(&bucket),
sizeof(
int));
464 is.read(reinterpret_cast<char *>(&numpoints),
sizeof(
int));
465 is.read(reinterpret_cast<char *>(&treesize),
sizeof(
int));
466 is.read(reinterpret_cast<char *>(&cost),
sizeof(
int));
468 if (!( is >> version1 >> realspec >> bucket >> numpoints >> treesize
472 if (!( version1 == version ))
474 if (!( realspec == std::numeric_limits<dist_t>::digits *
475 (std::numeric_limits<dist_t>::is_integer ? -1 : 1) ))
477 if (!( 0 <= bucket && bucket <= maxbucket ))
479 if (!( 0 <= treesize && treesize <= numpoints ))
484 std::vector<Node> tree;
485 tree.reserve(treesize);
486 for (
int i = 0; i < treesize; ++i) {
489 is.read(reinterpret_cast<char *>(&node.index),
sizeof(
int));
490 if (node.index >= 0) {
491 is.read(reinterpret_cast<char *>(node.data.lower),
493 is.read(reinterpret_cast<char *>(node.data.upper),
495 is.read(reinterpret_cast<char *>(node.data.child),
498 is.read(reinterpret_cast<char *>(node.leaves),
499 bucket *
sizeof(
int));
500 for (
int l = bucket; l < maxbucket; ++l)
504 if (!( is >> node.index ))
506 if (node.index >= 0) {
507 for (
int l = 0; l < 2; ++l) {
508 if (!( is >> node.data.lower[l] >> node.data.upper[l]
509 >> node.data.child[l] ))
515 for (
int l = 0; l < bucket; ++l) {
516 if (!( is >> node.leaves[l] ))
519 for (
int l = bucket; l < maxbucket; ++l)
523 node.Check(numpoints, treesize, bucket);
524 tree.push_back(node);
527 _numpoints = numpoints;
530 _cost = cost; _c1 = _k = _cmax = 0;
531 _cmin = std::numeric_limits<int>::max();
543 { t.
Save(os,
false);
return os; }
554 { t.
Load(is,
false);
return is; }
588 void Statistics(
int& setupcost,
int& numsearches,
int& searchcost,
589 int& mincost,
int& maxcost,
590 double& mean,
double& sd)
const {
591 setupcost = _cost; numsearches = _k; searchcost = _c1;
592 mincost = _cmin; maxcost = _cmax;
593 mean = _mc; sd = std::sqrt(_sc / (_k - 1));
602 _c1 = _k = _cmax = 0;
603 _cmin = std::numeric_limits<int>::max();
609 typedef std::pair<dist_t, int> item;
614 dist_t lower[2], upper[2];
619 int leaves[maxbucket];
626 for (
int i = 0; i < 2; ++i) {
627 data.lower[i] = data.upper[i] = 0;
633 void Check(
int numpoints,
int treesize,
int bucket)
const {
634 if (!( -1 <= index && index < numpoints ))
637 if (!( -1 <= data.child[0] && data.child[0] < treesize &&
638 -1 <= data.child[1] && data.child[1] < treesize ))
640 if (!( 0 <= data.lower[0] && data.lower[0] <= data.upper[0] &&
641 data.upper[0] <= data.lower[1] &&
642 data.lower[1] <= data.upper[1] ))
648 for (
int l = 0; l < bucket; ++l) {
650 ((l == 0 ? 0 : -1) <= leaves[l] && leaves[l] < numpoints) :
653 start = leaves[l] >= 0;
655 for (
int l = bucket; l < maxbucket; ++l) {
662 #if defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION) && \
663 GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
664 friend class boost::serialization::access;
665 template<
class Archive>
666 void save(Archive& ar,
const unsigned int)
const {
667 ar & boost::serialization::make_nvp(
"index", index);
669 ar & boost::serialization::make_nvp(
"leaves", leaves);
671 ar & boost::serialization::make_nvp(
"lower", data.lower)
672 & boost::serialization::make_nvp(
"upper", data.upper)
673 & boost::serialization::make_nvp(
"child", data.child);
675 template<
class Archive>
676 void load(Archive& ar,
const unsigned int) {
677 ar & boost::serialization::make_nvp(
"index", index);
679 ar & boost::serialization::make_nvp(
"leaves", leaves);
681 ar & boost::serialization::make_nvp(
"lower", data.lower)
682 & boost::serialization::make_nvp(
"upper", data.upper)
683 & boost::serialization::make_nvp(
"child", data.child);
685 template<
class Archive>
686 void serialize(Archive& ar,
const unsigned int file_version)
687 { boost::serialization::split_member(ar, *
this, file_version); }
691 #if defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION) && \
692 GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
693 friend class boost::serialization::access;
694 template<
class Archive>
void save(Archive& ar,
const unsigned)
const {
695 int realspec = std::numeric_limits<dist_t>::digits *
696 (std::numeric_limits<dist_t>::is_integer ? -1 : 1);
699 int version1 = version;
700 ar & boost::serialization::make_nvp(
"version", version1)
701 & boost::serialization::make_nvp(
"realspec", realspec)
702 & boost::serialization::make_nvp(
"bucket", _bucket)
703 & boost::serialization::make_nvp(
"numpoints", _numpoints)
704 & boost::serialization::make_nvp(
"cost", _cost)
705 & boost::serialization::make_nvp(
"tree", _tree);
707 template<
class Archive>
void load(Archive& ar,
const unsigned) {
708 int version1, realspec, bucket, numpoints, cost;
709 ar & boost::serialization::make_nvp(
"version", version1);
710 if (version1 != version)
712 std::vector<Node> tree;
713 ar & boost::serialization::make_nvp(
"realspec", realspec);
714 if (!( realspec == std::numeric_limits<dist_t>::digits *
715 (std::numeric_limits<dist_t>::is_integer ? -1 : 1) ))
717 ar & boost::serialization::make_nvp(
"bucket", bucket);
718 if (!( 0 <= bucket && bucket <= maxbucket ))
720 ar & boost::serialization::make_nvp(
"numpoints", numpoints)
721 & boost::serialization::make_nvp(
"cost", cost)
722 & boost::serialization::make_nvp(
"tree", tree);
723 if (!( 0 <=
int(tree.size()) &&
int(tree.size()) <= numpoints ))
726 for (
int i = 0; i < int(tree.size()); ++i)
727 tree[i].Check(numpoints,
int(tree.size()), bucket);
729 _numpoints = numpoints;
732 _cost = cost; _c1 = _k = _cmax = 0;
733 _cmin = std::numeric_limits<int>::max();
735 template<
class Archive>
736 void serialize(Archive& ar,
const unsigned int file_version)
737 { boost::serialization::split_member(ar, *
this, file_version); }
740 int _numpoints, _bucket, _cost;
741 std::vector<Node> _tree;
743 mutable double _mc, _sc;
744 mutable int _c1, _k, _cmin, _cmax;
746 int init(
const std::vector<pos_t>& pts,
const distfun_t& dist,
int bucket,
747 std::vector<Node>& tree, std::vector<item>& ids,
int& cost,
748 int l,
int u,
int vp) {
754 if (u - l > (bucket == 0 ? 1 : bucket)) {
760 int m = (u + l + 1) / 2;
762 for (
int k = l + 1; k < u; ++k) {
763 ids[k].first = dist(pts[ids[l].second], pts[ids[k].second]);
767 std::nth_element(ids.begin() + l + 1,
770 node.index = ids[l].second;
772 typename std::vector<item>::iterator
773 t = std::min_element(ids.begin() + l + 1, ids.begin() + m);
774 node.data.lower[0] = t->first;
775 t = std::max_element(ids.begin() + l + 1, ids.begin() + m);
776 node.data.upper[0] = t->first;
779 node.data.child[0] = init(pts, dist, bucket, tree, ids, cost,
780 l + 1, m,
int(t - ids.begin()));
782 typename std::vector<item>::iterator
783 t = std::max_element(ids.begin() + m, ids.begin() + u);
784 node.data.lower[1] = ids[m].first;
785 node.data.upper[1] = t->first;
787 node.data.child[1] = init(pts, dist, bucket, tree, ids, cost,
788 m, u,
int(t - ids.begin()));
791 node.index = ids[l].second;
796 std::sort(ids.begin() + l, ids.begin() + u);
797 for (
int i = l; i < u; ++i)
798 node.leaves[i-l] = ids[i].second;
799 for (
int i = u - l; i < bucket; ++i)
801 for (
int i = bucket; i < maxbucket; ++i)
806 tree.push_back(node);
807 return int(tree.size()) - 1;
826 template <
typename dist_t,
typename pos_t,
class distfun_t>
834 #if defined(_MSC_VER)
835 # pragma warning (pop)
838 #endif // GEOGRAPHICLIB_NEARESTNEIGHBOR_HPP