From: Luca Heltai Date: Sun, 13 Aug 2017 22:48:56 +0000 (-0600) Subject: Add KDTreeDistance implementation X-Git-Tag: v9.0.0-rc1~1262^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b8dbe68c94025f3b43391017649ab6b6eed82117;p=dealii.git Add KDTreeDistance implementation --- diff --git a/cmake/configure/configure_nanoflann.cmake b/cmake/configure/configure_nanoflann.cmake new file mode 100644 index 0000000000..c6bc02b693 --- /dev/null +++ b/cmake/configure/configure_nanoflann.cmake @@ -0,0 +1,20 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2012 - 2014 by the deal.II authors +## +## This file is part of the deal.II library. +## +## The deal.II library is free software; you can use it, redistribute +## it, and/or modify it under the terms of the GNU Lesser General +## Public License as published by the Free Software Foundation; either +## version 2.1 of the License, or (at your option) any later version. +## The full text of the license can be found in the file LICENSE at +## the top level of the deal.II distribution. +## +## --------------------------------------------------------------------- + +# +# Configuration for the nanoflann library: +# + +CONFIGURE_FEATURE(NANOFLANN) diff --git a/cmake/modules/FindNANOFLANN.cmake b/cmake/modules/FindNANOFLANN.cmake new file mode 100644 index 0000000000..b34d80a11d --- /dev/null +++ b/cmake/modules/FindNANOFLANN.cmake @@ -0,0 +1,37 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2014 - 2015 by the deal.II authors +## +## This file is part of the deal.II library. +## +## The deal.II library is free software; you can use it, redistribute +## it, and/or modify it under the terms of the GNU Lesser General +## Public License as published by the Free Software Foundation; either +## version 2.1 of the License, or (at your option) any later version. +## The full text of the license can be found in the file LICENSE at +## the top level of the deal.II distribution. +## +## --------------------------------------------------------------------- + +# +# Try to find the NANOFLANN library +# +# This module exports +# +# NANOFLANN_INCLUDE_DIRS +# + +SET(NANOFLANN_DIR "" CACHE PATH "An optional hint to a NANOFLANN installation") +SET_IF_EMPTY(NANOFLANN_DIR "$ENV{NANOFLANN_DIR}") + +DEAL_II_FIND_PATH(NANOFLANN_INCLUDE_DIR nanoflann.hpp + HINTS ${NANOFLANN_DIR} + PATH_SUFFIXES include + ) + +DEAL_II_PACKAGE_HANDLE(NANOFLANN + INCLUDE_DIRS REQUIRED NANOFLANN_INCLUDE_DIR + USER_INCLUDE_DIRS REQUIRED NANOFLANN_INCLUDE_DIR + CLEAR + NANOFLANN_INCLUDE_DIR + ) diff --git a/doc/news/changes/major/20170813LucaHeltai b/doc/news/changes/major/20170813LucaHeltai new file mode 100644 index 0000000000..a1ccfab047 --- /dev/null +++ b/doc/news/changes/major/20170813LucaHeltai @@ -0,0 +1,7 @@ +New: A new KDTreeDistance class has been added that interfaces +to the nanoflann library (https://github.com/jlblancoc/nanoflann). +This class can be used to extract nearest neighbour information +on collection of points, query for the closes points to a target +point or all points contained within a given distance. +
+(Luca Heltai, 2017/08/13) diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index 4e84a33888..bbce416294 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -46,6 +46,7 @@ #cmakedefine DEAL_II_WITH_METIS #cmakedefine DEAL_II_WITH_MPI #cmakedefine DEAL_II_WITH_MUPARSER +#cmakedefine DEAL_II_WITH_NANOFLANN #cmakedefine DEAL_II_WITH_NETCDF #cmakedefine DEAL_II_WITH_OPENCASCADE #cmakedefine DEAL_II_WITH_P4EST diff --git a/include/deal.II/numerics/kdtree_distance.h b/include/deal.II/numerics/kdtree_distance.h new file mode 100644 index 0000000000..810018eafa --- /dev/null +++ b/include/deal.II/numerics/kdtree_distance.h @@ -0,0 +1,266 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef _dealii__numerics_kdtree_distance_h +#define _dealii__numerics_kdtree_distance_h + +#include + +# ifdef DEAL_II_WITH_NANOFLANN + +#include + +DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +#include +DEAL_II_ENABLE_EXTRA_DIAGNOSTICS + + +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 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 thefunctionalities + * of the nanoflann library, like searching the n nearest neighbors, or + * searching the points that fall within a raidius of a target point. + */ +template +class KDTreeDistance +{ +public: + /** + * The max leaf parameter is used to decide how many points per leaf + * are used in the kdtree algorithm. + * + * If the points are not passed to this constructor, 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 undefinite 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 + * construction time, or when you call set_points()). If you update + * your points, and do not call again set_points(), your results + * will likely be wrong. + */ + KDTreeDistance(const unsigned int &max_leaf_size=10, + const std::vector > &pts=std::vector >()); + + + /** + * Adaptor class used internally by nanoflann. Class actually stores + * a reference to the vector of points, and generates some helper + * functions for nanoflann. + */ + struct PointCloudAdaptor + { + /** + * A typedef used by nanoflann. + */ + typedef double coord_t; + + + /** + * Reference to the vector of points from which we want to compute + * the distance. + */ + const std::vector > &points; //!< A const ref to the data set origin + + + /** + * The constrcutor needs the data set source. + */ + PointCloudAdaptor(const std::vector > &_points) : points(_points) { } + + + /** + * Return number of points in the data set (required by nanoflann). + */ + inline size_t kdtree_get_point_count() const + { + return points.size(); + } + + + /** + * Return the L2 distance between points + */ + inline coord_t kdtree_distance(const coord_t *p1, const size_t idx_p2,size_t size) const + { + AssertDimension(size, dim); + coord_t res=0.0; + for (size_t d=0; d + bool kdtree_get_bbox(BBOX &) const; + }; + + + /** + * A typedef for the actual KDTree object. + */ + typedef typename nanoflann::KDTreeSingleIndexAdaptor , + PointCloudAdaptor, dim, unsigned int> KDTree; + + + /** + * Store a reference to the passed points. After you called this + * method, you can call the value() method to compute the minimum + * distance between an evaluation point and the collection of points + * you passed to this method, or the get_points_within_ball() and + * the get_closest_points() methods. + * + * Notice that the constructor calls this method internally if you + * 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 + * storing the actual tree internally. If you change your points 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 + */ + void set_points(const std::vector > &pts); + + + /** + * A const accessor to the underlying points. + */ + const Point &operator[](unsigned int i) const; + + + /** + * The size of the vector stored by this class. + */ + unsigned int size() const; + + + /** + * Fill 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. + * + * @param[in] point: the target point + * @param[in] radius: the radius of the ball + * @param[out] mathes: indices and distances of the matching points + * @param[in] sorted: sort the output results in ascending order with respect to distances + * + * @return number of points that are within the ball + */ + unsigned int get_points_within_ball(const Point &target, const double &radius, + std::vector > &matches, + bool sorted=false) const; + + /** + * Fill two vectors with the indices and distances of the closest + * points to the given target point. The vectors are filled with + * indices and distances until there is space in them. You should + * resize them to the number of closest points you wish to get. An + * assertion is thrown if the vectors do not have the same size. + * + * @param[in] target: the target point + * @param[out] indices: indices of the matching points + * @param[out] distances: distances of the matching points + */ + void get_closest_points(const Point &target, + std::vector &indices, + std::vector &distances) const; + +private: + /** + * Max number of points per leaf. + */ + unsigned int max_leaf_size; + + + /** + * A point cloud adaptor, to be filled when set points is called. + */ + std::unique_ptr adaptor; + + + /** + * The actual kdtree. + */ + std::unique_ptr kdtree; +}; + + +//------------ inline functions ------------- + +template +inline +unsigned int KDTreeDistance::size() const +{ + if (adaptor) + return adaptor->points.size(); + else + return 0; +}; + +template +inline const Point & +KDTreeDistance::operator[](unsigned int i) const +{ + AssertIndexRange(i, size()); + return adaptor->points[i]; +} + + +template +template +inline bool +KDTreeDistance::PointCloudAdaptor::kdtree_get_bbox(BBOX &) const +{ + return false; +} + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_WITH_NANO_FLANN +#endif diff --git a/source/numerics/CMakeLists.txt b/source/numerics/CMakeLists.txt index 8de366e960..0d5252104a 100644 --- a/source/numerics/CMakeLists.txt +++ b/source/numerics/CMakeLists.txt @@ -23,6 +23,7 @@ SET(_unity_include_src data_postprocessor.cc dof_output_operator.cc histogram.cc + kdtree_distance.cc matrix_tools_once.cc matrix_tools.cc time_dependent.cc diff --git a/source/numerics/kdtree_distance.cc b/source/numerics/kdtree_distance.cc new file mode 100644 index 0000000000..616a3e847f --- /dev/null +++ b/source/numerics/kdtree_distance.cc @@ -0,0 +1,64 @@ +#include + +#ifdef DEAL_II_WITH_NANOFLANN + +#include + +DEAL_II_NAMESPACE_OPEN + + +template +KDTreeDistance::KDTreeDistance(const unsigned int &max_leaf_size, + const std::vector > &pts) + : max_leaf_size(max_leaf_size) +{ + if (pts.size() > 0) + set_points(pts); +} + + +template +unsigned int KDTreeDistance::get_points_within_ball(const Point ¢er, const double &radius, + std::vector > &matches, + bool sorted) const +{ + Assert(adaptor, ExcNotInitialized()); + Assert(kdtree, ExcInternalError()); + + Assert(radius > 0, + ExcMessage("Radius is expected to be positive.")); + + nanoflann::SearchParams params; + params.sorted = sorted; + return kdtree->radiusSearch(¢er[0], radius, matches, params); +} + +template +void KDTreeDistance::get_closest_points(const Point &target, + std::vector &indices, + std::vector &distances) const +{ + Assert(adaptor, ExcNotInitialized()); + Assert(kdtree, ExcInternalError()); + AssertDimension(indices.size(), distances.size()); + + kdtree->knnSearch(&target[0], indices.size(), &indices[0], &distances[0]); +} + +template +void KDTreeDistance::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->buildIndex(); +} + + +template class KDTreeDistance<1>; +template class KDTreeDistance<2>; +template class KDTreeDistance<3>; + +DEAL_II_NAMESPACE_CLOSE + +#endif