parsed_kdtree_distance.h

Go to the documentation of this file.

//-----------------------------------------------------------
//
//    Copyright (C) 2016 by the deal2lkit authors
//
//    This file is part of the deal2lkit library.
//
//    The deal2lkit 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 deal2lkit distribution.
//
//-----------------------------------------------------------

#ifndef _d2k_parsed_kdtree_distance_h
#define _d2k_parsed_kdtree_distance_h

#include <deal2lkit/config.h>
#include <deal2lkit/utilities.h>
#include <deal.II/base/function.h>
#include <deal2lkit/parameter_acceptor.h>

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#include <external/nanoflann.h>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

using namespace dealii;

D2K_NAMESPACE_OPEN

template<int dim>
class ParsedKDTreeDistance : public ParameterAcceptor, public Function<dim>
{
public:
  ParsedKDTreeDistance(const std::string &name="",
                       const unsigned int &max_leaf_size=10,
                       const std::vector<Point<dim> > &pts=std::vector<Point<dim> >());


  struct PointCloudAdaptor
  {
    typedef double coord_t;


    const std::vector<Point<dim> > &points; 


    PointCloudAdaptor(const std::vector<Point<dim> > &_points) : points(_points) { }


    inline size_t kdtree_get_point_count() const
    {
      return points.size();
    }


    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);
    }


    inline coord_t kdtree_get_pt(const size_t idx, int d) const
    {
      AssertIndexRange(d,dim);
      return points[idx][d];
    }


    template <class BBOX>
    bool kdtree_get_bbox(BBOX &) const
    {
      return false;
    }
  };


  typedef typename nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloudAdaptor> ,
          PointCloudAdaptor, dim, unsigned int>  KDTree;

  virtual void declare_parameters(ParameterHandler &prm);


  virtual double value(const Point<dim> &p, const unsigned int component=0) const;


  void set_points(const std::vector<Point<dim> > &pts);


  inline const Point<dim> &operator[](unsigned int i) const
  {
    AssertIndexRange(i, size());
    return adaptor->points[i];
  }


  inline unsigned int size() const
  {
    if (adaptor)
      return adaptor->points.size();
    else
      return 0;
  };


  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;

  void get_closest_points(const Point<dim> &target,
                          std::vector<unsigned int> &indices,
                          std::vector<double> &distances) const;

private:
  unsigned int max_leaf_size;


  shared_ptr<PointCloudAdaptor> adaptor;


  shared_ptr<KDTree> kdtree;
};

D2K_NAMESPACE_CLOSE

#endif