// Copyright (c) 2017 GeometryFactory Sarl (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). // // $URL: https://github.com/CGAL/cgal/blob/v5.2/Classification/include/CGAL/Classification/Point_set_neighborhood.h $ // $Id: Point_set_neighborhood.h c7cd9cb 2020-09-30T08:36:15+02:00 Simon Giraudot // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // Author(s) : Simon Giraudot #ifndef CGAL_CLASSIFICATION_POINT_SET_NEIGHBORHOOD_H #define CGAL_CLASSIFICATION_POINT_SET_NEIGHBORHOOD_H #include #include #include #include #include #include #include #include #include #include namespace CGAL { namespace Classification { /*! \ingroup PkgClassificationPointSet \brief Class that precomputes spatial searching structures for an input point set and gives access to the local neighborhood of a point as a set of indices. It allows the user to generate models of `NeighborQuery` based on a fixed range neighborhood or on a fixed K number of neighbors. In addition, the spatial searching structures can be computed on a simplified version of the point set to allow for neighbor queries at a higher scale. \tparam GeomTraits is a model of \cgal Kernel. \tparam PointRange model of `ConstRange`. Its iterator type is `RandomAccessIterator` and its value type is the key type of `PointMap`. \tparam PointMap model of `ReadablePropertyMap` whose key type is the value type of the iterator of `PointRange` and value type is `GeomTraits::Point_3`. */ template class Point_set_neighborhood { using FT = typename GeomTraits::FT; using Point = typename GeomTraits::Point_3; class My_point_property_map{ const PointRange* input; PointMap point_map; public: using value_type = typename boost::property_traits::value_type; using reference = typename boost::property_traits::reference; using key_type = std::uint32_t; using category = typename boost::property_traits::category; My_point_property_map () { } My_point_property_map (const PointRange *input, PointMap point_map) : input (input), point_map (point_map) { } friend reference get (const My_point_property_map& ppmap, key_type i) { return get(ppmap.point_map, *(ppmap.input->begin()+std::size_t(i))); } }; using Search_traits_base = Search_traits_3; using Search_traits = Search_traits_adapter; using Splitter = Sliding_midpoint; using Distance = Distance_adapter >; using Tree = Kd_tree; using Sphere = Fuzzy_sphere; using Knn = Orthogonal_k_neighbor_search; std::shared_ptr m_tree; Distance m_distance; public: /*! Functor that computes the neighborhood of an input point with a fixed number of neighbors. \cgalModels CGAL::Classification::NeighborQuery \sa Point_set_neighborhood */ class K_neighbor_query { public: using value_type = typename Point_set_neighborhood::Point; ///< private: const Point_set_neighborhood& neighborhood; unsigned int k; public: /*! \brief constructs a K neighbor query object. \param neighborhood point set neighborhood object. \param k number of neighbors per query. */ K_neighbor_query (const Point_set_neighborhood& neighborhood, unsigned int k) : neighborhood (neighborhood), k(k) { } /// \cond SKIP_IN_MANUAL template OutputIterator operator() (const value_type& query, OutputIterator output) const { neighborhood.k_neighbors (query, k, output); return output; } /// \endcond }; /*! Functor that computes the neighborhood of an input point defined as the points lying in a sphere of fixed radius centered at the input point. \cgalModels CGAL::Classification::NeighborQuery \sa Point_set_neighborhood */ class Sphere_neighbor_query { public: using value_type = typename Point_set_neighborhood::Point; ///< private: const Point_set_neighborhood& neighborhood; float radius; public: /*! \brief constructs a range neighbor query object. \param neighborhood point set neighborhood object. \param radius radius of the neighbor query sphere. */ Sphere_neighbor_query (const Point_set_neighborhood& neighborhood, float radius) : neighborhood (neighborhood), radius(radius) { } /// \cond SKIP_IN_MANUAL template OutputIterator operator() (const value_type& query, OutputIterator output) const { neighborhood.sphere_neighbors (query, radius, output); return output; } /// \endcond }; /// \cond SKIP_IN_MANUAL friend class K_neighbor_query; friend class Sphere_neighbor_query; Point_set_neighborhood () { } /// \endcond /// \name Constructors /// @{ /*! \brief constructs a neighborhood object based on the input range. \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`. If no tag is provided, `Parallel_if_available_tag` is used. \param input point range. \param point_map property map to access the input points. */ template Point_set_neighborhood (const PointRange& input, PointMap point_map, const ConcurrencyTag&) : m_tree (nullptr) { init (input, point_map); } /// \cond SKIP_IN_MANUAL Point_set_neighborhood (const PointRange& input, PointMap point_map) : m_tree (nullptr) { init (input, point_map); } template void init (const PointRange& input, PointMap point_map) { My_point_property_map pmap (&input, point_map); m_tree = std::make_shared (boost::counting_iterator (0), boost::counting_iterator (std::uint32_t(input.size())), Splitter(), Search_traits (pmap)); m_distance = Distance (pmap); m_tree->template build(); } /// \endcond /*! \brief constructs a simplified neighborhood object based on the input range. This method first computes a simplified version of the input point set by voxelization: a 3D grid is defined and for each subset present in one cell, only the point closest to the centroid of this subset is used. \tparam ConcurrencyTag enables sequential versus parallel algorithm. Possible values are `Sequential_tag`, `Parallel_tag`, and `Parallel_if_available_tag`. If no tag is provided, `Parallel_if_available_tag` is used. \param input input range. \param point_map property map to access the input points. \param voxel_size size of the cells of the 3D grid used for simplification. */ template Point_set_neighborhood (const PointRange& input, PointMap point_map, float voxel_size, const ConcurrencyTag&) : m_tree (nullptr) { init (input, point_map, voxel_size); } /// \cond SKIP_IN_MANUAL Point_set_neighborhood (const PointRange& input, PointMap point_map, float voxel_size) { init (input, point_map, voxel_size); } template void init (const PointRange& input, PointMap point_map, float voxel_size) { // First, simplify std::vector indices; My_point_property_map pmap (&input, point_map); voxelize_point_set(input.size(), indices, pmap, voxel_size); m_tree = std::make_shared (indices.begin(), indices.end(), Splitter(), Search_traits (pmap)); m_distance = Distance (pmap); m_tree->template build(); } /// \endcond /// @} /// \name Queries /// @{ /*! \brief returns a neighbor query object with fixed number of neighbors `k`. */ K_neighbor_query k_neighbor_query (const unsigned int k) const { return K_neighbor_query (*this, k); } /*! \brief returns a neighbor query object with fixed radius `radius`. */ Sphere_neighbor_query sphere_neighbor_query (const float radius) const { return Sphere_neighbor_query (*this, radius); } /// @} private: template void sphere_neighbors (const Point& query, const FT radius_neighbors, OutputIterator output) const { CGAL_assertion (m_tree != nullptr); Sphere fs (query, radius_neighbors, 0, m_tree->traits()); m_tree->search (output, fs); } template void k_neighbors (const Point& query, const unsigned int k, OutputIterator output) const { CGAL_assertion (m_tree != nullptr); Knn search (*m_tree, query, k, 0, true, m_distance); for (typename Knn::iterator it = search.begin(); it != search.end(); ++ it) *(output ++) = it->first; } template void voxelize_point_set (std::size_t nb_pts, std::vector& indices, Map point_map, float voxel_size) { std::map > grid; for (std::uint32_t i = 0; i < nb_pts; ++ i) { const Point& p = get(point_map, i); Point ref (std::floor(p.x() / voxel_size), std::floor(p.y() / voxel_size), std::floor(p.z() / voxel_size)); typename std::map >::iterator it; boost::tie (it, boost::tuples::ignore) = grid.insert (std::make_pair (ref, std::vector())); it->second.push_back (i); } for (typename std::map >::iterator it = grid.begin(); it != grid.end(); ++ it) { const std::vector& pts = it->second; Point centroid = CGAL::centroid (CGAL::make_transform_iterator_from_property_map (pts.begin(), point_map), CGAL::make_transform_iterator_from_property_map (pts.end(), point_map)); std::uint32_t chosen = 0; float min_dist = (std::numeric_limits::max)(); for (std::size_t i = 0; i < pts.size(); ++ i) { float dist = float(CGAL::squared_distance(get(point_map, pts[i]), centroid)); if (dist < min_dist) { min_dist = dist; chosen = pts[i]; } } indices.push_back (chosen); } } }; } } #endif // CGAL_CLASSIFICATION_POINT_SET_POINT_SET_NEIGHBORHOOD_H