33 #include "../nabo/index_heap.h" 44 template<
typename T,
typename CloudType>
49 for (
int i = 0; i < v.size(); ++i)
60 template<
typename T,
typename CloudType>
68 if (elCount & (1 << i))
71 for (
int j = 0; j <= i; ++j)
77 template<
typename T,
typename CloudType>
81 for (
int i = 0; i < indexes.size(); ++i)
82 cloudIndexes.coeffRef(i) = nodes[indexes[i]].index;
86 template<
typename T,
typename CloudType>
89 const size_t count(last - first);
93 nodes[pos] =
Node(first->pos, -1, first->index);
99 Vector mean(Vector::Zero(this->dim));
102 mean /= last - first;
104 Vector var(Vector::Zero(this->dim));
106 var += (it->pos - mean).cwise() * (it->pos - mean);
108 const size_t cutDim = argMax<T>(var);
114 const size_t recurseCount(count-1);
115 const size_t rightCount(recurseCount/2);
116 const size_t leftCount(recurseCount-rightCount);
117 assert(last - rightCount == first + leftCount + 1);
119 nodes[pos] =
Node((first+leftCount)->pos, cutDim, (first+leftCount)->index);
126 buildNodes(first, first + leftCount, childLeft(pos));
127 buildNodes(first + leftCount + 1, last, childRight(pos));
131 nodes[childLeft(pos)] =
Node(first->pos, -1, first->index);
132 nodes[childRight(pos)] =
Node(
Vector(), -2, 0);
136 template<
typename T,
typename CloudType>
139 const Node& node(nodes[pos]);
144 cout <<
"<circle cx=\"" << 100*node.
pos(0) <<
"\" cy=\"" << 100*node.
pos(1) <<
"\" r=\"1\" stroke=\"black\" stroke-width=\"0.2\" fill=\"red\"/>" << endl;
146 cout <<
"pt at\n" << node.
pos << endl;
153 Vector leftMaxValues(maxValues);
154 leftMaxValues[node.
dim] = node.
pos[node.
dim];
156 Vector rightMinValues(minValues);
157 rightMinValues[node.
dim] = node.
pos[node.
dim];
161 cout <<
"<line x1=\"" << 100*rightMinValues(0) <<
"\" y1=\"" << 100*rightMinValues(1) <<
"\" x2=\"" << 100*leftMaxValues(0) <<
"\" y2=\"" << 100*leftMaxValues(1) <<
"\" style=\"stroke:rgb(0,0,0);stroke-width:0.2\"/>" << endl;
163 cout <<
"cut from\n" << rightMinValues <<
"\nto\n" << leftMaxValues << endl;
165 dump(minValues, leftMaxValues, childLeft(pos));
166 dump(rightMinValues, maxValues, childRight(pos));
170 template<
typename T,
typename CloudType>
176 buildPoints.reserve(cloud.cols());
177 for (
int i = 0; i < cloud.cols(); ++i)
179 const Vector& v(cloud.col(i));
187 buildNodes(buildPoints.begin(), buildPoints.end(), 0);
195 template<
typename T,
typename CloudType>
201 template<
typename T,
typename CloudType>
204 typedef priority_queue<SearchElement> Queue;
206 const T maxError(1 + epsilon);
212 statistics.lastQueryVisitCount = 0;
214 while (!queue.empty())
220 if (el.minDist * maxError > heap.
headValue())
227 assert (node.
dim != -2);
230 const T dist(dist2<T>(node.
pos, query));
232 (allowSelfMatch || (dist > numeric_limits<T>::epsilon())))
239 const T offset(query.coeff(node.
dim) - node.
pos.coeff(node.
dim));
240 const T offset2(offset * offset);
264 ++statistics.lastQueryVisitCount;
267 statistics.totalVisitCount += statistics.lastQueryVisitCount;
284 template<
typename T,
typename CloudType>
290 template<
typename T,
typename CloudType>
295 assert(
nodes.size() > 0);
296 assert(
nodes[0].pos.size() == query.size());
300 statistics.lastQueryVisitCount = 0;
302 recurseKnn(query, 0, 0, heap, off, 1 + epsilon, allowSelfMatch);
307 statistics.totalVisitCount += statistics.lastQueryVisitCount;
312 template<
typename T,
typename CloudType>
316 const int cd(node.
dim);
318 ++statistics.lastQueryVisitCount;
323 const T dist(dist2<T, CloudType>(node.
pos, query));
325 (allowSelfMatch || (dist > numeric_limits<T>::epsilon()))
331 const T old_off(off.coeff(cd));
332 const T new_off(query.coeff(cd) - node.
pos.coeff(cd));
336 rd += - old_off*old_off + new_off*new_off;
339 off.coeffRef(cd) = new_off;
341 off.coeffRef(cd) = old_off;
347 rd += - old_off*old_off + new_off*new_off;
350 off.coeffRef(cd) = new_off;
352 off.coeffRef(cd) = old_off;
367 template<
typename T,
typename CloudType>
377 if (elCount & (1 << i))
380 for (
int j = 0; j <= i; ++j)
387 template<
typename T,
typename CloudType>
390 const size_t count(last - first);
394 const int dim = -2-(first->index);
395 assert(pos <
nodes.size());
408 mean /= last - first;
412 var += (it->pos - mean).cwise() * (it->pos - mean);
414 cutDim = argMax<T>(var);
419 cutDim = argMax<T>(maxValues - minValues);
423 const size_t rightCount(count/2);
424 const size_t leftCount(count - rightCount);
425 assert(last - rightCount == first + leftCount);
429 nth_element(first, first + leftCount, last,
CompareDim(cutDim));
432 const T cutVal((first+leftCount)->pos.coeff(cutDim));
438 Vector leftMaxValues(maxValues);
439 leftMaxValues[cutDim] = cutVal;
441 Vector rightMinValues(minValues);
442 rightMinValues[cutDim] = cutVal;
445 buildNodes(first, first + leftCount,
childLeft(pos), minValues, leftMaxValues, balanceVariance);
446 buildNodes(first + leftCount, last,
childRight(pos), rightMinValues, maxValues, balanceVariance);
449 template<
typename T,
typename CloudType>
455 buildPoints.reserve(cloud.cols());
456 for (
int i = 0; i < cloud.cols(); ++i)
458 const Vector& v(cloud.col(i));
471 template<
typename T,
typename CloudType>
476 assert(
nodes.size() > 0);
478 Vector off(Vector::Zero(query.size()));
480 statistics.lastQueryVisitCount = 0;
482 recurseKnn(query, 0, 0, heap, off, 1 + epsilon, allowSelfMatch);
487 statistics.totalVisitCount += statistics.lastQueryVisitCount;
489 return heap.getIndexes();
492 template<
typename T,
typename CloudType>
496 const int cd(node.
dim);
498 ++statistics.lastQueryVisitCount;
504 const int index(-(cd + 2));
505 const T dist(dist2<T>(query,
cloud.col(index)));
507 (allowSelfMatch || (dist > numeric_limits<T>::epsilon()))
513 const T old_off(off.coeff(cd));
514 const T new_off(query.coeff(cd) - node.
cutVal);
518 rd += - old_off*old_off + new_off*new_off;
521 off.coeffRef(cd) = new_off;
523 off.coeffRef(cd) = old_off;
529 rd += - old_off*old_off + new_off*new_off;
532 off.coeffRef(cd) = new_off;
534 off.coeffRef(cd) = old_off;
549 template<
typename T,
typename Heap,
typename CloudType>
552 const size_t count(last - first);
553 const unsigned pos(
nodes.size());
563 const unsigned cutDim = argMax<T>(maxValues - minValues);
564 T cutVal((maxValues(cutDim) + minValues(cutDim))/2);
571 size_t rightStart(0);
572 while (rightStart < count && (first+rightStart)->pos.coeff(cutDim) < cutVal)
578 cutVal = first->pos.coeff(cutDim);
581 else if (rightStart == count)
583 rightStart = count - 1;
584 cutVal = (first + rightStart)->pos.coeff(cutDim);
588 Vector leftMaxValues(maxValues);
589 leftMaxValues[cutDim] = cutVal;
591 Vector rightMinValues(minValues);
592 rightMinValues[cutDim] = cutVal;
595 const size_t rightCount(count - rightStart);
596 const size_t leftCount(count - rightCount);
599 nodes.push_back(
Node(cutDim, cutVal, 0));
602 const unsigned __attribute__ ((unused)) leftChild =
buildNodes(first, first + leftCount, minValues, leftMaxValues);
603 assert(leftChild == pos + 1);
604 const unsigned rightChild =
buildNodes(first + leftCount, last, rightMinValues, maxValues);
607 nodes[pos].rightChild = rightChild;
611 template<
typename T,
typename Heap,
typename CloudType>
617 buildPoints.reserve(cloud.cols());
618 for (
int i = 0; i < cloud.cols(); ++i)
620 const Vector& v(cloud.col(i));
633 template<
typename T,
typename Heap,
typename CloudType>
638 assert(
nodes.size() > 0);
640 Vector off(Vector::Zero(query.size()));
642 statistics.lastQueryVisitCount = 0;
644 recurseKnn(query, 0, 0, heap, off, 1+epsilon, allowSelfMatch);
649 statistics.totalVisitCount += statistics.lastQueryVisitCount;
651 return heap.getIndexes();
654 template<
typename T,
typename Heap,
typename CloudType>
658 assert(
nodes.size() > 0);
660 assert(
nodes.size() > 0);
665 const int colCount(query.cols());
667 for (
int i = 0; i < colCount; ++i)
674 statistics.lastQueryVisitCount = 0;
676 recurseKnn(
q, 0, 0, heap, off, 1+epsilon, allowSelfMatch);
681 result.col(i) = heap.getIndexes();
683 statistics.totalVisitCount += statistics.lastQueryVisitCount;
689 template<
typename T,
typename Heap,
typename CloudType>
697 const unsigned index(node.
ptIndex);
701 const T* qPtr(&query.coeff(0));
702 const T* dPtr(&
cloud.coeff(0, index));
703 const int dim(query.size());
704 for (
int i = 0; i <
dim; ++i)
706 const T diff(*qPtr - *dPtr);
710 if ((dist < heap.headValue()) &&
711 (allowSelfMatch || (dist > numeric_limits<T>::epsilon()))
713 heap.replaceHead(index, dist);
717 const unsigned cd(node.
dim);
718 const T old_off(off.coeff(cd));
719 const T new_off(query.coeff(cd) - node.
cutVal);
723 rd += - old_off*old_off + new_off*new_off;
724 if (rd * maxError < heap.headValue())
726 off.coeffRef(cd) = new_off;
727 recurseKnn(query, n + 1, rd, heap, off, maxError, allowSelfMatch);
728 off.coeffRef(cd) = old_off;
733 recurseKnn(query, n+1, rd, heap, off, maxError, allowSelfMatch);
734 rd += - old_off*old_off + new_off*new_off;
735 if (rd * maxError < heap.headValue())
737 off.coeffRef(cd) = new_off;
739 off.coeffRef(cd) = old_off;
760 template<
typename T,
typename CloudType>
763 const size_t count(last - first);
764 const unsigned pos(
nodes.size());
769 const int dim = -1-(first->index);
775 const int cutDim = argMax<T>(maxValues - minValues);
776 T cutVal((maxValues(cutDim) + minValues(cutDim))/2);
783 size_t rightStart(0);
784 while (rightStart < count && (first+rightStart)->pos.coeff(cutDim) < cutVal)
790 cutVal = first->pos.coeff(cutDim);
793 else if (rightStart == count)
795 rightStart = count - 1;
796 cutVal = (first + rightStart)->pos.coeff(cutDim);
800 Vector leftMaxValues(maxValues);
801 leftMaxValues[cutDim] = cutVal;
803 Vector rightMinValues(minValues);
804 rightMinValues[cutDim] = cutVal;
807 const size_t rightCount(count - rightStart);
808 const size_t leftCount(count - rightCount);
811 nodes.push_back(
Node(cutDim, cutVal, minValues.coeff(cutDim), maxValues.coeff(cutDim)));
814 const unsigned __attribute__ ((unused)) leftChild =
buildNodes(first, first + leftCount, minValues, leftMaxValues);
815 assert(leftChild == pos + 1);
816 const unsigned rightChild =
buildNodes(first + leftCount, last, rightMinValues, maxValues);
819 nodes[pos].rightChild = rightChild;
823 template<
typename T,
typename CloudType>
829 buildPoints.reserve(cloud.cols());
830 for (
int i = 0; i < cloud.cols(); ++i)
832 const Vector& v(cloud.col(i));
845 template<
typename T,
typename CloudType>
850 assert(
nodes.size() > 0);
853 statistics.lastQueryVisitCount = 0;
855 recurseKnn(query, 0, 0, heap, 1+epsilon, allowSelfMatch);
860 statistics.totalVisitCount += statistics.lastQueryVisitCount;
862 return heap.getIndexes();
865 template<
typename T,
typename CloudType>
869 const int cd(node.
dim);
871 ++statistics.lastQueryVisitCount;
875 const int index(-(cd + 1));
876 const T dist(dist2<T>(query,
cloud.col(index)));
878 (allowSelfMatch || (dist > numeric_limits<T>::epsilon()))
884 const T q_val(query.coeff(cd));
885 const T cut_diff(q_val - node.
cutVal);
888 recurseKnn(query, n+1, rd, heap, maxError, allowSelfMatch);
894 rd += cut_diff*cut_diff - box_diff*box_diff;
907 rd += cut_diff*cut_diff - box_diff*box_diff;
910 recurseKnn(query, n + 1, rd, heap, maxError, allowSelfMatch);
NearestNeighbourSearch< T, CloudType >::Vector Vector
size_t childRight(size_t pos) const
KDTreeBalancedPtInNodesPQ(const CloudType &cloud)
KDTreeUnbalancedPtInLeavesExplicitBoundsStack(const CloudType &cloud)
void sort()
sort the entries, from the smallest to the largest
NearestNeighbourSearch< T, CloudType >::IndexVector IndexVector
std::vector< BuildPoint > BuildPoints
NearestNeighbourSearch< T, CloudType >::Index Index
BuildPoints::iterator BuildPointsIt
KDTreeBalancedPtInLeavesStack(const CloudType &cloud, const bool balanceVariance)
Nearest neighbour search interface, templatized on scalar type.
std::vector< BuildPoint > BuildPoints
KDTreeBalancedPtInNodes< T, CloudType >::Node Node
const VT & headValue() const
get the largest value of the heap
unsigned buildNodes(const BuildPointsIt first, const BuildPointsIt last, const Vector minValues, const Vector maxValues)
NearestNeighbourSearch< T, CloudType >::IndexVector IndexVector
NearestNeighbourSearch< T, CloudType >::IndexMatrix IndexMatrix
virtual IndexVector knn(const Vector &query, const Index k, const T epsilon, const unsigned optionFlags)
size_t getTreeSize(size_t size) const
void recurseKnn(const Vector &query, const size_t n, T rd, Heap &heap, Vector &off, const T maxError, const bool allowSelfMatch)
const CloudType & cloud
the reference to the data-point cloud, which must remain valid during the lifetime of the NearestNeig...
virtual IndexVector knn(const Vector &query, const Index k, const T epsilon, const unsigned optionFlags)
unsigned buildNodes(const BuildPointsIt first, const BuildPointsIt last, const Vector minValues, const Vector maxValues)
size_t childRight(size_t pos) const
NearestNeighbourSearch< T, CloudType >::IndexVector IndexVector
size_t childLeft(size_t pos) const
size_t childLeft(size_t pos) const
void recurseKnn(const Vector &query, const size_t n, T rd, Heap &heap, const T maxError, const bool allowSelfMatch)
NearestNeighbourSearch< T, CloudType >::Matrix Matrix
NearestNeighbourSearch< T, CloudType >::Vector Vector
virtual IndexVector knn(const Vector &query, const Index k, const T epsilon, const unsigned optionFlags)
void buildNodes(const BuildPointsIt first, const BuildPointsIt last, const size_t pos)
size_t argMax(const typename NearestNeighbourSearch< T, CloudType >::Vector &v)
Return the index of the maximum value of a vector of positive values.
void replaceHead(const Index index, const Value value)
put value into heap, replace the largest value if full
const Vector maxBound
the high bound of the search space (axis-aligned bounding box)
NearestNeighbourSearch< T, CloudType >::IndexVector IndexVector
balanced-tree implementation of heap
NearestNeighbourSearch< T, CloudType >::Vector Vector
virtual IndexVector knn(const Vector &query, const Index k, const T epsilon, const unsigned optionFlags)
KDTreeBalancedPtInNodes(const CloudType &cloud)
const Vector minBound
the low bound of the search space (axis-aligned bounding box)
NearestNeighbourSearch< T, CloudType >::Index Index
BuildPoints::iterator BuildPointsIt
NearestNeighbourSearch< T, CloudType >::IndexVector IndexVector
NearestNeighbourSearch< T, CloudType >::IndexVector IndexVector
void buildNodes(const BuildPointsIt first, const BuildPointsIt last, const size_t pos, const Vector minValues, const Vector maxValues, const bool balanceVariance)
virtual IndexVector knn(const Vector &query, const Index k, const T epsilon, const unsigned optionFlags)
const Index dim
the dimensionality of the data-point cloud
size_t getTreeSize(size_t size) const
std::vector< BuildPoint > BuildPoints
NearestNeighbourSearch< T, CloudType >::Index Index
void dump(const Vector minValues, const Vector maxValues, const size_t pos) const
BuildPoints::const_iterator BuildPointsCstIt
IndexVector cloudIndexesFromNodesIndexes(const IndexVector &indexes) const
std::vector< BuildPoint > BuildPoints
virtual IndexMatrix knnM(const Matrix &query, const Index k, const T epsilon, const unsigned optionFlags)
CloudType CloudType
a column-major Eigen matrix in which each column is a point; this matrix has dim rows ...
KDTreeBalancedPtInNodesStack(const CloudType &cloud)
NearestNeighbourSearch< T, CloudType >::Index Index
void recurseKnn(const Vector &query, const unsigned n, T rd, Heap &heap, Vector &off, const T maxError, const bool allowSelfMatch)
BuildPoints::const_iterator BuildPointsCstIt
void recurseKnn(const Vector &query, const size_t n, T rd, Heap &heap, Vector &off, const T maxError, const bool allowSelfMatch)
KDTreeUnbalancedPtInLeavesImplicitBoundsStack(const CloudType &cloud)
BuildPoints::iterator BuildPointsIt
NearestNeighbourSearch< T, CloudType >::Vector Vector
BuildPoints::iterator BuildPointsIt