]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add KDTreeDistance implementation
authorLuca Heltai <luca.heltai@sissa.it>
Sun, 13 Aug 2017 22:48:56 +0000 (16:48 -0600)
committerLuca Heltai <luca.heltai@sissa.it>
Sun, 13 Aug 2017 23:31:45 +0000 (17:31 -0600)
cmake/configure/configure_nanoflann.cmake [new file with mode: 0644]
cmake/modules/FindNANOFLANN.cmake [new file with mode: 0644]
doc/news/changes/major/20170813LucaHeltai [new file with mode: 0644]
include/deal.II/base/config.h.in
include/deal.II/numerics/kdtree_distance.h [new file with mode: 0644]
source/numerics/CMakeLists.txt
source/numerics/kdtree_distance.cc [new file with mode: 0644]

diff --git a/cmake/configure/configure_nanoflann.cmake b/cmake/configure/configure_nanoflann.cmake
new file mode 100644 (file)
index 0000000..c6bc02b
--- /dev/null
@@ -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 (file)
index 0000000..b34d80a
--- /dev/null
@@ -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 (file)
index 0000000..a1ccfab
--- /dev/null
@@ -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.
+<br>
+(Luca Heltai, 2017/08/13)
index 4e84a3388888b1839878f60fefeffb675f15dc1a..bbce416294304731d613249c3da47ffec96fe03a 100644 (file)
@@ -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 (file)
index 0000000..810018e
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#  ifdef DEAL_II_WITH_NANOFLANN
+
+#include <deal.II/base/point.h>
+
+DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+#include <nanoflann.hpp>
+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<int dim>
+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<Point<dim> > &pts=std::vector<Point<dim> >());
+
+
+  /**
+   * 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<Point<dim> > &points; //!< A const ref to the data set origin
+
+
+    /**
+     * The constrcutor needs the data set source.
+     */
+    PointCloudAdaptor(const std::vector<Point<dim> > &_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<size; ++d)
+        res += (p1[d]-points[idx_p2][d])*(p1[d]-points[idx_p2][d]);
+      return std::sqrt(res);
+    }
+
+
+    /**
+     * Return the dim'th component of the idx'th point in the class.
+     */
+    inline coord_t kdtree_get_pt(const size_t idx, int d) const
+    {
+      AssertIndexRange(d,dim);
+      return points[idx][d];
+    }
+
+
+    /**
+     * Optional bounding-box computation: return false to default to a
+     * standard bbox computation loop.  Return true if the BBOX was
+     * already computed by the class and returned in "bb" so it can be
+     * avoided to redo it again.  Look at bb.size() to find out the
+     * expected dimensionality (e.g. 2 or 3 for point clouds).
+     */
+    template <class BBOX>
+    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>  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<Point<dim> > &pts);
+
+
+  /**
+   * A const accessor to the underlying points.
+   */
+  const Point<dim> &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<dim> &target, const double &radius,
+                                      std::vector<std::pair<unsigned int, double> > &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<dim> &target,
+                          std::vector<unsigned int> &indices,
+                          std::vector<double> &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<PointCloudAdaptor> adaptor;
+
+
+  /**
+   * The actual kdtree.
+   */
+  std::unique_ptr<KDTree> kdtree;
+};
+
+
+//------------ inline functions -------------
+
+template<int dim>
+inline
+unsigned int KDTreeDistance<dim>::size() const
+{
+  if (adaptor)
+    return adaptor->points.size();
+  else
+    return 0;
+};
+
+template<int dim>
+inline const Point<dim> &
+KDTreeDistance<dim>::operator[](unsigned int i) const
+{
+  AssertIndexRange(i, size());
+  return adaptor->points[i];
+}
+
+
+template<int dim>
+template <class BBOX>
+inline bool
+KDTreeDistance<dim>::PointCloudAdaptor::kdtree_get_bbox(BBOX &) const
+{
+  return false;
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif // DEAL_II_WITH_NANO_FLANN
+#endif
index 8de366e960883882feda85c9d04ad604c5edcab3..0d5252104a318112bc550acbb89f9cb9ac4125fc 100644 (file)
@@ -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 (file)
index 0000000..616a3e8
--- /dev/null
@@ -0,0 +1,64 @@
+#include <deal.II/numerics/kdtree_distance.h>
+
+#ifdef DEAL_II_WITH_NANOFLANN
+
+#include<deal.II/base/std_cxx14/memory.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template<int dim>
+KDTreeDistance<dim>::KDTreeDistance(const unsigned int &max_leaf_size,
+                                    const std::vector<Point<dim> > &pts)
+  : max_leaf_size(max_leaf_size)
+{
+  if (pts.size() > 0)
+    set_points(pts);
+}
+
+
+template<int dim>
+unsigned int KDTreeDistance<dim>::get_points_within_ball(const Point<dim> &center, const double &radius,
+                                                         std::vector<std::pair<unsigned int, double> > &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(&center[0], radius, matches, params);
+}
+
+template<int dim>
+void KDTreeDistance<dim>::get_closest_points(const Point<dim> &target,
+                                             std::vector<unsigned int> &indices,
+                                             std::vector<double> &distances) const
+{
+  Assert(adaptor, ExcNotInitialized());
+  Assert(kdtree, ExcInternalError());
+  AssertDimension(indices.size(), distances.size());
+
+  kdtree->knnSearch(&target[0], indices.size(), &indices[0], &distances[0]);
+}
+
+template<int dim>
+void KDTreeDistance<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<KDTree>(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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.