From 7e2b26cfe5849a29402992a8703440f8ba3f798b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 22 Aug 2017 20:41:38 -0600 Subject: [PATCH] Minor code updates to the kdtree implementation. --- include/deal.II/numerics/kdtree.h | 121 +++++++++++++++++------------- source/numerics/kdtree.cc | 37 ++++++--- 2 files changed, 95 insertions(+), 63 deletions(-) diff --git a/include/deal.II/numerics/kdtree.h b/include/deal.II/numerics/kdtree.h index d87cde2156..ea8a15be66 100644 --- a/include/deal.II/numerics/kdtree.h +++ b/include/deal.II/numerics/kdtree.h @@ -34,30 +34,31 @@ DEAL_II_NAMESPACE_OPEN /** * 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. @@ -67,28 +68,33 @@ class KDTree { 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 > &pts=std::vector >()); + KDTree(const unsigned int &max_leaf_size = 10, + const std::vector > &pts = std::vector >()); /** @@ -115,7 +121,7 @@ public: * The constructor needs the vector of points from which we want to build * the tree. */ - PointCloudAdaptor(const std::vector > &_points); + PointCloudAdaptor (const std::vector > &_points); /** @@ -127,13 +133,16 @@ public: /** * 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; /** @@ -144,15 +153,17 @@ public: * expected dimensionality (e.g. 2 or 3 for point clouds). */ template - bool kdtree_get_bbox(BBOX &) const; + bool kdtree_get_bbox (BBOX &) const; }; /** * A typedef for the actual KDTree object. */ - typedef typename nanoflann::KDTreeSingleIndexAdaptor , - PointCloudAdaptor, dim, unsigned int> NanoFlannKDTree; + typedef + typename nanoflann::KDTreeSingleIndexAdaptor , + PointCloudAdaptor, dim, unsigned int> + NanoFlannKDTree; /** @@ -163,7 +174,7 @@ public: * 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 @@ -171,56 +182,57 @@ public: * 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 > &pts); + void set_points (const std::vector > &pts); /** - * A const accessor to the underlying points. + * A const accessor to the @p i'th one among the underlying points. */ - const Point &operator[](unsigned int i) const; + const Point &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 > get_points_within_ball(const Point &target, - const double &radius, - bool sorted=false) const; + std::vector > + get_points_within_ball (const Point &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 > get_closest_points(const Point &target, - const unsigned int n_points) const; + std::vector > + get_closest_points (const Point &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; /** @@ -253,7 +265,7 @@ unsigned int KDTree::size() const template inline const Point & -KDTree::operator[](unsigned int i) const +KDTree::operator[] (const unsigned int i) const { AssertIndexRange(i, size()); return adaptor->points[i]; @@ -263,7 +275,8 @@ KDTree::operator[](unsigned int i) const template KDTree::PointCloudAdaptor::PointCloudAdaptor(const std::vector > &_points) - : points(_points) + : + points(_points) {} @@ -290,7 +303,7 @@ KDTree::PointCloudAdaptor::kdtree_get_pt(const size_t idx, int d) const template template inline bool -KDTree::PointCloudAdaptor::kdtree_get_bbox(BBOX &) const +KDTree::PointCloudAdaptor::kdtree_get_bbox (BBOX &) const { return false; } @@ -299,12 +312,14 @@ KDTree::PointCloudAdaptor::kdtree_get_bbox(BBOX &) const template inline double -KDTree::PointCloudAdaptor::kdtree_distance(const double *p1, const size_t idx_p2,size_t size) const +KDTree::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 -KDTree::KDTree(const unsigned int &max_leaf_size, +KDTree::KDTree(const unsigned int &max_leaf_size, const std::vector > &pts) - : max_leaf_size(max_leaf_size) + : + max_leaf_size(max_leaf_size) { if (pts.size() > 0) set_points(pts); } + template -std::vector > KDTree::get_points_within_ball(const Point ¢er, - const double &radius, - bool sorted) const +std::vector > +KDTree::get_points_within_ball(const Point ¢er, + const double &radius, + bool sorted) const { - std::vector > matches; Assert(adaptor, ExcNotInitialized()); Assert(kdtree, ExcInternalError()); @@ -31,32 +33,47 @@ std::vector > KDTree::get_points_within_bal nanoflann::SearchParams params; params.sorted = sorted; + + std::vector > matches; kdtree->radiusSearch(¢er[0], radius, matches, params); + return matches; } + + template -std::vector > KDTree::get_closest_points(const Point &target, - const unsigned int n_points) const +std::vector > +KDTree::get_closest_points(const Point &target, + const unsigned int n_points) const { Assert(adaptor, ExcNotInitialized()); Assert(kdtree, ExcInternalError()); + + // get the information out of nanoflann std::vector indices(n_points); std::vector distances(n_points); - std::vector > matches(n_points); kdtree->knnSearch(&target[0], n_points, &indices[0], &distances[0]); + + // convert it to the format we want to return + std::vector > matches(n_points); for (unsigned int i=0; i void KDTree::set_points(const std::vector > &pts) { Assert(pts.size() > 0, ExcMessage("Expecting a non zero set of points.")); adaptor = std_cxx14::make_unique(pts); - kdtree = std_cxx14::make_unique(dim, *adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size)); + kdtree = std_cxx14::make_unique(dim, + *adaptor, + nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size)); kdtree->buildIndex(); } -- 2.39.5