/**
* A wrapper for the nanoflann library, used to compute the distance from a
* collection of points, and to efficiently return nearest neighbors to a
- * target point. This function uses nanoflann to efficiently partition the
- * space in a k-dimensional tree. The cost of each query is then roughly of
- * order log(n), where n is the number of points stored in this class.
+ * target point. This class uses nanoflann to efficiently partition the
+ * space in a $k$-dimensional tree. The cost of each query is then roughly of
+ * order $\log(n)$, where $n$ is the number of points stored in this class.
*
* The wrapper provides methods that give access to some of the functionalities
- * of the nanoflann library, like searching the n nearest neighbors, or
+ * of the nanoflann library, like searching the $p$ nearest neighbors of
+ * a given point, or
* searching the points that fall within a radius of a target point.
*
* @quotation
* From wikipedia (https://en.wikipedia.org/wiki/K-d_tree):
*
- * A k-d tree is a binary tree in which every node is a k-dimensional point.
+ * A k-d tree is a binary tree in which every node is a $k$-dimensional point.
* Every non-leaf node can be thought of as implicitly generating a splitting
* hyperplane that divides the space into two parts, known as half-spaces.
* Points to the left of this hyperplane are represented by the left subtree of
* that node and points right of the hyperplane are represented by the right
* subtree. The hyperplane direction is chosen in the following way: every node
- * in the tree is associated with one of the k-dimensions, with the hyperplane
+ * in the tree is associated with one of the $k$-dimensions, with the hyperplane
* perpendicular to that dimension's axis. So, for example, if for a particular
* split the "x" axis is chosen, all points in the subtree with a smaller "x"
* value than the node will appear in the left subtree and all points with
* larger "x" value will be in the right subtree. In such a case, the
- * hyperplane would be set by the x-value of the point, and its normal would be
- * the unit x-axis.
+ * hyperplane would be set by the $x$-value of the point, and its normal would be
+ * the unit $x$-axis.
* @endquotation
*
* @author Luca Heltai, 2017.
{
public:
/**
- * The max leaf parameter is used to decide how many points per leaf
+ * Constructor.
+ *
+ * @param[in] max_leaf_size A number denoting how many points per leaf
* are used in the kdtree algorithm.
*
- * If the points are not passed to this constructor, then you have
+ * @param[in] pts A vector of points that are to be represented by
+ * the current object. If no points are passed to this constructor
+ * (or if the default value of the argument is used), then you have
* to pass them later to this object by calling the set_points()
* method.
*
* Access to any of the methods without first passing a reference to
* a vector of points will result in an exception. Only a reference
* to the points is stored, so you should make sure that the life of
- * the the vector you pass is longer than the life of this class, or
- * you'll get undefined behaviour.
+ * the vector you pass is longer than the life of this class, or
+ * you will get undefined behaviour.
*
- * If you update the vector of points in someway, remember to call
- * again the set_points() method. The tree and the index are
- * constructed only once, when you pass the points (either at
+ * @warning If you change the contents of the vector of points that you
+ * passed either to the constructor or to set_points(), remember to call
+ * the set_points() method again. The tree and the index are
+ * constructed only once when you pass the points (either at
* construction time, or when you call set_points()). If you update
- * your points, and do not call again set_points(), your results
+ * your points, and do not call set_points() again, then all following results
* will likely be wrong.
*/
- KDTree(const unsigned int &max_leaf_size=10,
- const std::vector<Point<dim> > &pts=std::vector<Point<dim> >());
+ KDTree(const unsigned int &max_leaf_size = 10,
+ const std::vector<Point<dim> > &pts = std::vector<Point<dim> >());
/**
* The constructor needs the vector of points from which we want to build
* the tree.
*/
- PointCloudAdaptor(const std::vector<Point<dim> > &_points);
+ PointCloudAdaptor (const std::vector<Point<dim> > &_points);
/**
/**
* Return the L2 distance between points
*/
- coord_t kdtree_distance(const coord_t *p1, const size_t idx_p2,size_t size) const;
+ coord_t kdtree_distance (const coord_t *p1,
+ const size_t idx_p2,
+ const size_t size) const;
/**
* Return the d-th component of the idx-th point in the class.
*/
- coord_t kdtree_get_pt(const size_t idx, int d) const;
+ coord_t kdtree_get_pt (const size_t idx,
+ const int d) const;
/**
* expected dimensionality (e.g. 2 or 3 for point clouds).
*/
template <class BBOX>
- bool kdtree_get_bbox(BBOX &) const;
+ bool kdtree_get_bbox (BBOX &) const;
};
/**
* A typedef for the actual KDTree object.
*/
- typedef typename nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloudAdaptor> ,
- PointCloudAdaptor, dim, unsigned int> NanoFlannKDTree;
+ typedef
+ typename nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloudAdaptor> ,
+ PointCloudAdaptor, dim, unsigned int>
+ NanoFlannKDTree;
/**
* the get_closest_points() methods.
*
* Notice that the constructor calls this method internally if you
- * pass it a non empty vector of points.
+ * pass it a non-empty vector of points.
*
* Whenever your points change, you should call this method again,
* since this is the method responsible for building the index and
* don't call again this method, any function you call later will
* happily return wrong values without you noticing.
*
- * @param pts: a collection of points
+ * @param[in] pts A collection of points
*/
- void set_points(const std::vector<Point<dim> > &pts);
+ void set_points (const std::vector<Point<dim> > &pts);
/**
- * A const accessor to the underlying points.
+ * A const accessor to the @p i'th one among the underlying points.
*/
- const Point<dim> &operator[](unsigned int i) const;
+ const Point<dim> &operator[] (const unsigned int i) const;
/**
- * The size of the vector stored by this class.
+ * The number of points currently stored by this class.
*/
unsigned int size() const;
/**
- * Fill a vector with the indices and the distance of the points
+ * Fill and return a vector with the indices and the distance of the points
* that are at distance less than or equal to the given radius from
- * the target point. Consider preallocating the size of the return
- * vector if you have a wild guess of how many should be there.
+ * the target point.
*
- * @param[in] point: the target point
- * @param[in] radius: the radius of the ball
- * @param[in] sorted: sort the output results in ascending order with respect to distances
+ * @param[in] point The target point
+ * @param[in] radius The radius of the ball
+ * @param[in] sorted Sort the output results in ascending order with respect to distances
*
- * @return vector of indices and distances of the matching points
+ * @return A vector of indices and distances of the matching points
*/
- std::vector<std::pair<unsigned int, double> > get_points_within_ball(const Point<dim> &target,
- const double &radius,
- bool sorted=false) const;
+ std::vector<std::pair<unsigned int, double> >
+ get_points_within_ball (const Point<dim> &target,
+ const double &radius,
+ const bool sorted=false) const;
/**
- * Fill a vectors with the indices and distances of the closest @p n_points
+ * Fill and return a vector with the indices and distances of the closest @p n_points
* points to the given target point.
*
- * @param[in] target: the target point
- * @param[in] n_points: the number of requested points
+ * @param[in] target The target point
+ * @param[in] n_points The number of requested points
*
- * @return a vector of pairs of indices and distances of the matching points
+ * @return A vector of pairs of indices and distances of the matching points
*/
- std::vector<std::pair<unsigned int, double> > get_closest_points(const Point<dim> &target,
- const unsigned int n_points) const;
+ std::vector<std::pair<unsigned int, double> >
+ get_closest_points (const Point<dim> &target,
+ const unsigned int n_points) const;
private:
/**
- * Max number of points per leaf.
+ * Max number of points per leaf as set in the constructor.
*/
- unsigned int max_leaf_size;
+ const unsigned int max_leaf_size;
/**
template <int dim>
inline const Point<dim> &
-KDTree<dim>::operator[](unsigned int i) const
+KDTree<dim>::operator[] (const unsigned int i) const
{
AssertIndexRange(i, size());
return adaptor->points[i];
template <int dim>
KDTree<dim>::PointCloudAdaptor::PointCloudAdaptor(const std::vector<Point<dim> > &_points)
- : points(_points)
+ :
+ points(_points)
{}
template <int dim>
template <class BBOX>
inline bool
-KDTree<dim>::PointCloudAdaptor::kdtree_get_bbox(BBOX &) const
+KDTree<dim>::PointCloudAdaptor::kdtree_get_bbox (BBOX &) const
{
return false;
}
template <int dim>
inline double
-KDTree<dim>::PointCloudAdaptor::kdtree_distance(const double *p1, const size_t idx_p2,size_t size) const
+KDTree<dim>::PointCloudAdaptor::kdtree_distance (const double *p1,
+ const size_t idx_p2,
+ const size_t size) const
{
AssertDimension(size, dim);
double res=0.0;
for (size_t d=0; d<size; ++d)
- res += (p1[d]-points[idx_p2][d])*(p1[d]-points[idx_p2][d]);
+ res += (p1[d]-points[idx_p2][d]) * (p1[d]-points[idx_p2][d]);
return std::sqrt(res);
}
#endif
template <int dim>
-KDTree<dim>::KDTree(const unsigned int &max_leaf_size,
+KDTree<dim>::KDTree(const unsigned int &max_leaf_size,
const std::vector<Point<dim> > &pts)
- : max_leaf_size(max_leaf_size)
+ :
+ max_leaf_size(max_leaf_size)
{
if (pts.size() > 0)
set_points(pts);
}
+
template <int dim>
-std::vector<std::pair<unsigned int, double> > KDTree<dim>::get_points_within_ball(const Point<dim> ¢er,
- const double &radius,
- bool sorted) const
+std::vector<std::pair<unsigned int, double> >
+KDTree<dim>::get_points_within_ball(const Point<dim> ¢er,
+ const double &radius,
+ bool sorted) const
{
- std::vector<std::pair<unsigned int, double> > matches;
Assert(adaptor, ExcNotInitialized());
Assert(kdtree, ExcInternalError());
nanoflann::SearchParams params;
params.sorted = sorted;
+
+ std::vector<std::pair<unsigned int, double> > matches;
kdtree->radiusSearch(¢er[0], radius, matches, params);
+
return matches;
}
+
+
template <int dim>
-std::vector<std::pair<unsigned int, double> > KDTree<dim>::get_closest_points(const Point<dim> &target,
- const unsigned int n_points) const
+std::vector<std::pair<unsigned int, double> >
+KDTree<dim>::get_closest_points(const Point<dim> &target,
+ const unsigned int n_points) const
{
Assert(adaptor, ExcNotInitialized());
Assert(kdtree, ExcInternalError());
+
+ // get the information out of nanoflann
std::vector<unsigned int> indices(n_points);
std::vector<double> distances(n_points);
- std::vector<std::pair<unsigned int, double> > matches(n_points);
kdtree->knnSearch(&target[0], n_points, &indices[0], &distances[0]);
+
+ // convert it to the format we want to return
+ std::vector<std::pair<unsigned int, double> > matches(n_points);
for (unsigned int i=0; i<n_points; ++i)
matches[i] = std::make_pair(indices[i], distances[i]);
+
return matches;
}
+
+
template <int dim>
void KDTree<dim>::set_points(const std::vector<Point<dim> > &pts)
{
Assert(pts.size() > 0, ExcMessage("Expecting a non zero set of points."));
adaptor = std_cxx14::make_unique<PointCloudAdaptor>(pts);
- kdtree = std_cxx14::make_unique<NanoFlannKDTree>(dim, *adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
+ kdtree = std_cxx14::make_unique<NanoFlannKDTree>(dim,
+ *adaptor,
+ nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
kdtree->buildIndex();
}