98 static const int version = 1;
101 static const int maxbucket =
102 (2 + ((4 *
sizeof(dist_t)) /
sizeof(int) >= 2 ?
103 (4 *
sizeof(dist_t)) /
sizeof(
int) : 2));
163 void Initialize(
const std::vector<pos_t>& pts,
const distfun_t& dist,
165 static_assert(std::numeric_limits<dist_t>::is_signed,
166 "dist_t must be a signed type");
167 if (!( 0 <= bucket && bucket <= maxbucket ))
169 (
"bucket must lie in [0, 2 + 4*sizeof(dist_t)/sizeof(int)]");
170 if (pts.size() > size_t(std::numeric_limits<int>::max()))
173 std::vector<item> ids(pts.size());
174 for (
int k =
int(ids.size()); k--;)
175 ids[k] = std::make_pair(dist_t(0), k);
177 std::vector<Node> tree;
178 init(pts, dist, bucket, tree, ids, cost,
179 0,
int(ids.size()),
int(ids.size()/2));
181 _numpoints = int(pts.size());
184 _cost = cost; _c1 = _k = _cmax = 0;
185 _cmin = std::numeric_limits<int>::max();
251 dist_t
Search(
const std::vector<pos_t>& pts,
const distfun_t& dist,
253 std::vector<int>& ind,
255 dist_t maxdist = std::numeric_limits<dist_t>::max(),
257 bool exhaustive =
true,
258 dist_t tol = 0)
const {
259 if (_numpoints !=
int(pts.size()))
261 std::priority_queue<item> results;
262 if (_numpoints > 0 && k > 0 && maxdist > mindist) {
264 dist_t tau = maxdist;
268 std::priority_queue<item> todo;
269 todo.push(std::make_pair(dist_t(1),
int(_tree.size()) - 1));
271 while (!todo.empty()) {
272 int n = todo.top().second;
273 dist_t d = -todo.top().first;
275 dist_t tau1 = tau - tol;
277 if (!( n >= 0 && tau1 >= d ))
continue;
278 const Node& current = _tree[n];
280 bool exitflag =
false, leaf = current.index < 0;
281 for (
int i = 0; i < (leaf ? _bucket : 1); ++i) {
282 int index = leaf ? current.leaves[i] : current.index;
283 if (index < 0)
break;
284 dst = dist(pts[index], query);
287 if (dst > mindist && dst <= tau) {
288 if (
int(results.size()) == k) results.pop();
289 results.push(std::make_pair(dst, index));
290 if (
int(results.size()) == k) {
292 tau = results.top().first;
306 if (current.index < 0)
continue;
308 for (
int l = 0; l < 2; ++l) {
309 if (current.data.child[l] >= 0 &&
310 dst + current.data.upper[l] >= mindist) {
311 if (dst < current.data.lower[l]) {
312 d = current.data.lower[l] - dst;
314 todo.push(std::make_pair(-d, current.data.child[l]));
315 }
else if (dst > current.data.upper[l]) {
316 d = dst - current.data.upper[l];
318 todo.push(std::make_pair(-d, current.data.child[l]));
320 todo.push(std::make_pair(dist_t(1), current.data.child[l]));
327 _mc += (c - omc) / _k;
328 _sc += (c - omc) * (c - _mc);
329 if (c > _cmax) _cmax = c;
330 if (c < _cmin) _cmin = c;
334 ind.resize(results.size());
336 for (
int i =
int(ind.size()); i--;) {
337 ind[i] = int(results.top().second);
338 if (i == 0) d = results.top().first;
367 void Save(std::ostream& os,
bool bin =
true)
const {
368 int realspec = std::numeric_limits<dist_t>::digits *
369 (std::numeric_limits<dist_t>::is_integer ? -1 : 1);
371 char id[] =
"NearestNeighbor_";
378 buf[4] = int(_tree.size());
380 os.write(
reinterpret_cast<const char *
>(buf), 6 *
sizeof(
int));
381 for (
int i = 0; i < int(_tree.size()); ++i) {
382 const Node& node = _tree[i];
383 os.write(
reinterpret_cast<const char *
>(&node.index),
sizeof(int));
384 if (node.index >= 0) {
385 os.write(
reinterpret_cast<const char *
>(node.data.lower),
387 os.write(
reinterpret_cast<const char *
>(node.data.upper),
389 os.write(
reinterpret_cast<const char *
>(node.data.child),
392 os.write(
reinterpret_cast<const char *
>(node.leaves),
393 _bucket *
sizeof(int));
397 std::stringstream ostring;
400 if (!std::numeric_limits<dist_t>::is_integer) {
401 static const int prec
402 = int(std::ceil(std::numeric_limits<dist_t>::digits *
403 std::log10(2.0) + 1));
404 ostring.precision(prec);
406 ostring << version <<
" " << realspec <<
" " << _bucket <<
" "
407 << _numpoints <<
" " << _tree.size() <<
" " << _cost;
408 for (
int i = 0; i < int(_tree.size()); ++i) {
409 const Node& node = _tree[i];
410 ostring <<
"\n" << node.index;
411 if (node.index >= 0) {
412 for (
int l = 0; l < 2; ++l)
413 ostring <<
" " << node.data.lower[l] <<
" " << node.data.upper[l]
414 <<
" " << node.data.child[l];
416 for (
int l = 0; l < _bucket; ++l)
417 ostring <<
" " << node.leaves[l];
445 void Load(std::istream& is,
bool bin =
true) {
446 int version1, realspec, bucket, numpoints, treesize, cost;
451 if (!(std::strcmp(
id,
"NearestNeighbor_") == 0))
453 is.read(
reinterpret_cast<char *
>(&version1),
sizeof(
int));
454 is.read(
reinterpret_cast<char *
>(&realspec),
sizeof(
int));
455 is.read(
reinterpret_cast<char *
>(&bucket),
sizeof(
int));
456 is.read(
reinterpret_cast<char *
>(&numpoints),
sizeof(
int));
457 is.read(
reinterpret_cast<char *
>(&treesize),
sizeof(
int));
458 is.read(
reinterpret_cast<char *
>(&cost),
sizeof(
int));
460 if (!( is >> version1 >> realspec >> bucket >> numpoints >> treesize
464 if (!( version1 == version ))
466 if (!( realspec == std::numeric_limits<dist_t>::digits *
467 (std::numeric_limits<dist_t>::is_integer ? -1 : 1) ))
469 if (!( 0 <= bucket && bucket <= maxbucket ))
471 if (!( 0 <= treesize && treesize <= numpoints ))
476 std::vector<Node> tree;
477 tree.reserve(treesize);
478 for (
int i = 0; i < treesize; ++i) {
481 is.read(
reinterpret_cast<char *
>(&node.index),
sizeof(int));
482 if (node.index >= 0) {
483 is.read(
reinterpret_cast<char *
>(node.data.lower),
485 is.read(
reinterpret_cast<char *
>(node.data.upper),
487 is.read(
reinterpret_cast<char *
>(node.data.child),
490 is.read(
reinterpret_cast<char *
>(node.leaves),
491 bucket *
sizeof(int));
492 for (
int l = bucket; l < maxbucket; ++l)
496 if (!( is >> node.index ))
498 if (node.index >= 0) {
499 for (
int l = 0; l < 2; ++l) {
500 if (!( is >> node.data.lower[l] >> node.data.upper[l]
501 >> node.data.child[l] ))
507 for (
int l = 0; l < bucket; ++l) {
508 if (!( is >> node.leaves[l] ))
511 for (
int l = bucket; l < maxbucket; ++l)
515 node.Check(numpoints, treesize, bucket);
516 tree.push_back(node);
519 _numpoints = numpoints;
522 _cost = cost; _c1 = _k = _cmax = 0;
523 _cmin = std::numeric_limits<int>::max();
535 { t.
Save(os,
false);
return os; }
546 { t.
Load(is,
false);
return is; }
554 std::swap(_numpoints, t._numpoints);
555 std::swap(_bucket, t._bucket);
556 std::swap(_cost, t._cost);
558 std::swap(_mc, t._mc);
559 std::swap(_sc, t._sc);
560 std::swap(_c1, t._c1);
562 std::swap(_cmin, t._cmin);
563 std::swap(_cmax, t._cmax);
580 void Statistics(
int& setupcost,
int& numsearches,
int& searchcost,
581 int& mincost,
int& maxcost,
582 double& mean,
double& sd)
const {
583 setupcost = _cost; numsearches = _k; searchcost = _c1;
584 mincost = _cmin; maxcost = _cmax;
585 mean = _mc; sd = std::sqrt(_sc / (_k - 1));
594 _c1 = _k = _cmax = 0;
595 _cmin = std::numeric_limits<int>::max();
601 typedef std::pair<dist_t, int> item;
606 dist_t lower[2], upper[2];
611 int leaves[maxbucket];
618 for (
int i = 0; i < 2; ++i) {
619 data.lower[i] = data.upper[i] = 0;
625 void Check(
int numpoints,
int treesize,
int bucket)
const {
626 if (!( -1 <= index && index < numpoints ))
629 if (!( -1 <= data.child[0] && data.child[0] < treesize &&
630 -1 <= data.child[1] && data.child[1] < treesize ))
632 if (!( 0 <= data.lower[0] && data.lower[0] <= data.upper[0] &&
633 data.upper[0] <= data.lower[1] &&
634 data.lower[1] <= data.upper[1] ))
640 for (
int l = 0; l < bucket; ++l) {
642 ((l == 0 ? 0 : -1) <= leaves[l] && leaves[l] < numpoints) :
645 start = leaves[l] >= 0;
647 for (
int l = bucket; l < maxbucket; ++l) {
654#if defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION) && \
655 GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
656 friend class boost::serialization::access;
657 template<
class Archive>
658 void save(Archive& ar,
const unsigned int)
const {
659 ar & boost::serialization::make_nvp(
"index", index);
661 ar & boost::serialization::make_nvp(
"leaves", leaves);
663 ar & boost::serialization::make_nvp(
"lower", data.lower)
664 & boost::serialization::make_nvp(
"upper", data.upper)
665 & boost::serialization::make_nvp(
"child", data.child);
667 template<
class Archive>
668 void load(Archive& ar,
const unsigned int) {
669 ar & boost::serialization::make_nvp(
"index", index);
671 ar & boost::serialization::make_nvp(
"leaves", leaves);
673 ar & boost::serialization::make_nvp(
"lower", data.lower)
674 & boost::serialization::make_nvp(
"upper", data.upper)
675 & boost::serialization::make_nvp(
"child", data.child);
677 template<
class Archive>
678 void serialize(Archive& ar,
const unsigned int file_version)
679 { boost::serialization::split_member(ar, *
this, file_version); }
683#if defined(GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION) && \
684 GEOGRAPHICLIB_HAVE_BOOST_SERIALIZATION
685 friend class boost::serialization::access;
686 template<
class Archive>
void save(Archive& ar,
const unsigned)
const {
687 int realspec = std::numeric_limits<dist_t>::digits *
688 (std::numeric_limits<dist_t>::is_integer ? -1 : 1);
691 int version1 = version;
692 ar & boost::serialization::make_nvp(
"version", version1)
693 & boost::serialization::make_nvp(
"realspec", realspec)
694 & boost::serialization::make_nvp(
"bucket", _bucket)
695 & boost::serialization::make_nvp(
"numpoints", _numpoints)
696 & boost::serialization::make_nvp(
"cost", _cost)
697 & boost::serialization::make_nvp(
"tree", _tree);
699 template<
class Archive>
void load(Archive& ar,
const unsigned) {
700 int version1, realspec, bucket, numpoints, cost;
701 ar & boost::serialization::make_nvp(
"version", version1);
702 if (version1 != version)
704 std::vector<Node> tree;
705 ar & boost::serialization::make_nvp(
"realspec", realspec);
706 if (!( realspec == std::numeric_limits<dist_t>::digits *
707 (std::numeric_limits<dist_t>::is_integer ? -1 : 1) ))
708 throw
GeographicLib::GeographicErr(
"Different dist_t types");
709 ar & boost::serialization::make_nvp(
"bucket", bucket);
710 if (!( 0 <= bucket && bucket <= maxbucket ))
712 ar & boost::serialization::make_nvp(
"numpoints", numpoints)
713 & boost::serialization::make_nvp(
"cost", cost)
714 & boost::serialization::make_nvp(
"tree", tree);
715 if (!( 0 <=
int(tree.size()) &&
int(tree.size()) <= numpoints ))
718 for (
int i = 0; i < int(tree.size()); ++i)
719 tree[i].Check(numpoints,
int(tree.size()), bucket);
721 _numpoints = numpoints;
724 _cost = cost; _c1 = _k = _cmax = 0;
725 _cmin = std::numeric_limits<int>::max();
727 template<
class Archive>
728 void serialize(Archive& ar,
const unsigned int file_version)
729 { boost::serialization::split_member(ar, *
this, file_version); }
732 int _numpoints, _bucket, _cost;
733 std::vector<Node> _tree;
735 mutable double _mc, _sc;
736 mutable int _c1, _k, _cmin, _cmax;
738 int init(
const std::vector<pos_t>& pts,
const distfun_t& dist,
int bucket,
739 std::vector<Node>& tree, std::vector<item>& ids,
int& cost,
740 int l,
int u,
int vp) {
746 if (u - l > (bucket == 0 ? 1 : bucket)) {
750 std::swap(ids[l], ids[i]);
752 int m = (u + l + 1) / 2;
754 for (
int k = l + 1; k < u; ++k) {
755 ids[k].first = dist(pts[ids[l].second], pts[ids[k].second]);
759 std::nth_element(ids.begin() + l + 1,
762 node.index = ids[l].second;
764 typename std::vector<item>::iterator
765 t = std::min_element(ids.begin() + l + 1, ids.begin() + m);
766 node.data.lower[0] = t->first;
767 t = std::max_element(ids.begin() + l + 1, ids.begin() + m);
768 node.data.upper[0] = t->first;
771 node.data.child[0] = init(pts, dist, bucket, tree, ids, cost,
772 l + 1, m,
int(t - ids.begin()));
774 typename std::vector<item>::iterator
775 t = std::max_element(ids.begin() + m, ids.begin() + u);
776 node.data.lower[1] = ids[m].first;
777 node.data.upper[1] = t->first;
779 node.data.child[1] = init(pts, dist, bucket, tree, ids, cost,
780 m, u,
int(t - ids.begin()));
783 node.index = ids[l].second;
788 std::sort(ids.begin() + l, ids.begin() + u);
789 for (
int i = l; i < u; ++i)
790 node.leaves[i-l] = ids[i].second;
791 for (
int i = u - l; i < bucket; ++i)
793 for (
int i = bucket; i < maxbucket; ++i)
798 tree.push_back(node);
799 return int(tree.size()) - 1;