Add region growing on circles for point_set_2

Shape_detection/examples/Shape_detection/CMakeLists.txt

@@ -23,6 +23,7 @@ if(TARGET CGAL::Eigen_support)
   create_single_source_cgal_program("region_growing_on_point_set_2.cpp")
   create_single_source_cgal_program("region_growing_on_point_set_3.cpp")
   create_single_source_cgal_program("region_growing_spheres_on_point_set_3.cpp")
+  create_single_source_cgal_program("region_growing_circles_on_point_set_2.cpp")
   create_single_source_cgal_program("region_growing_on_polygon_mesh.cpp")
   create_single_source_cgal_program("region_growing_with_custom_classes.cpp")
 
@@ -36,6 +37,7 @@ if(TARGET CGAL::Eigen_support)
     region_growing_on_point_set_2
     region_growing_on_point_set_3
     region_growing_spheres_on_point_set_3
+    region_growing_circles_on_point_set_2
     region_growing_on_polygon_mesh
     region_growing_with_custom_classes
     shape_detection_basic_deprecated)

Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp

@@ -0,0 +1,104 @@
+#include <CGAL/Real_timer.h>
+#include <CGAL/Random.h>
+#include <CGAL/Simple_cartesian.h>
+
+#include <CGAL/Point_set_3.h>
+#include <CGAL/Point_set_3/IO.h>
+
+#include <CGAL/Shape_detection/Region_growing/Region_growing.h>
+#include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
+
+#include <boost/function_output_iterator.hpp>
+
+#include <fstream>
+
+using Kernel = CGAL::Simple_cartesian<double>;
+using Point_3 = Kernel::Point_3;
+using Vector_3 = Kernel::Vector_3;
+using Point_2 = Kernel::Point_2;
+using Vector_2 = Kernel::Vector_2;
+
+using Point_set_3 = CGAL::Point_set_3<Point_3, Vector_3>;
+using Point_set_2 = CGAL::Point_set_3<Point_2, Vector_2>;
+using Point_map = Point_set_2::Point_map;
+using Normal_map = Point_set_2::Vector_map;
+
+namespace Shape_detection = CGAL::Shape_detection::Point_set;
+
+using Neighbor_query = Shape_detection::K_neighbor_query
+  <Kernel, Point_set_2, Point_map>;
+using Region_type = Shape_detection::Least_squares_circle_fit_region
+  <Kernel, Point_set_2, Point_map, Normal_map>;
+using Region_growing = CGAL::Shape_detection::Region_growing
+  <Point_set_2, Neighbor_query, Region_type>;
+
+int main (int argc, char** argv)
+{
+  std::ifstream ifile (argc > 1 ? argv[1] : "data/circle.ply");
+  Point_set_3 points3;
+  ifile >> points3;
+
+  std::cerr << points3.size() << " points read" << std::endl;
+
+  // Input should have normals
+  assert (points.has_normal_map());
+
+  Point_set_2 points;
+  points.add_normal_map();
+  for (Point_set_3::Index idx : points3)
+  {
+    const Point_3& p = points3.point(idx);
+    const Vector_3& n = points3.normal(idx);
+    points.insert (Point_2 (p.x(), p.y()), Vector_2 (n.x(), n.y()));
+  }
+
+  // Default parameters for data/circles.ply
+  const std::size_t k = 12;
+  const double tolerance = 0.01;
+  const double max_angle = 10.;
+  const std::size_t min_region_size = 20;
+
+  Neighbor_query neighbor_query(points, k, points.point_map());
+  Region_type region_type(points, tolerance, max_angle, min_region_size,
+                          points.point_map(), points.normal_map());
+  Region_growing region_growing(points, neighbor_query, region_type);
+
+  // Add maps to get colored output
+  Point_set_3::Property_map<unsigned char>
+    red = points3.add_property_map<unsigned char>("red", 0).first,
+    green = points3.add_property_map<unsigned char>("green", 0).first,
+    blue = points3.add_property_map<unsigned char>("blue", 0).first;
+
+  CGAL::Random random;
+
+  std::size_t nb_circles = 0;
+  CGAL::Real_timer timer;
+  timer.start();
+  region_growing.detect
+    (boost::make_function_output_iterator
+     ([&](const std::vector<std::size_t>& region)
+      {
+        // Assign a random color to each region
+        unsigned char r = (unsigned char)(random.get_int(64, 192));
+        unsigned char g = (unsigned char)(random.get_int(64, 192));
+        unsigned char b = (unsigned char)(random.get_int(64, 192));
+        for (const std::size_t& idx : region)
+        {
+          red[idx] = r;
+          green[idx] = g;
+          blue[idx] = b;
+        }
+        ++ nb_circles;
+      }));
+  timer.stop();
+
+  std::cerr << nb_circles << " circles detected in "
+            << timer.time() << " seconds" << std::endl;
+
+  // Save in colored_spheres.ply
+  std::ofstream out ("colored_circles.ply");
+  CGAL::set_binary_mode (out);
+  out << points3;
+
+  return EXIT_SUCCESS;
+}

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h

@@ -20,6 +20,7 @@
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Sphere_neighbor_query.h>
 
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_line_fit_region.h>
+#include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_plane_fit_region.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h>
 

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h

@@ -0,0 +1,340 @@
+// Copyright (c) 2020 GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Simon Giraudot
+//
+
+#ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_REGION_H
+#define CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_REGION_H
+
+#include <CGAL/license/Shape_detection.h>
+
+#include <vector>
+
+#include <CGAL/assertions.h>
+#include <CGAL/number_utils.h>
+#include <CGAL/Cartesian_converter.h>
+#include <CGAL/Eigen_diagonalize_traits.h>
+
+#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
+
+namespace CGAL {
+namespace Shape_detection {
+namespace Point_set {
+
+/*!
+  \ingroup PkgShapeDetectionRGOnPoints
+
+  \brief Region type based on the quality of the least squares circle
+  fit applied to 2D points.
+
+  This class fits a circle to chunks of points in a 2D point set and
+  controls the quality of this fit.  If all quality conditions are
+  satisfied, the chunk is accepted as a valid region, otherwise
+  rejected.
+
+  \tparam GeomTraits
+  must be a model of `Kernel`.
+
+  \tparam InputRange
+  must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`.
+
+  \tparam PointMap
+  must be an `LvaluePropertyMap` whose key type is the value type of the input
+  range and value type is `Kernel::Point_3`.
+
+  \tparam NormalMap
+  must be an `LvaluePropertyMap` whose key type is the value type of the input
+  range and value type is `Kernel::Vector_3`.
+
+  \cgalModels `RegionType`
+*/
+template<typename GeomTraits,
+         typename InputRange,
+         typename PointMap,
+         typename NormalMap>
+class Least_squares_circle_fit_region
+{
+
+public:
+
+  /// \name Types
+  /// @{
+
+  /// \cond SKIP_IN_MANUAL
+  using Traits = GeomTraits;
+  using Input_range = InputRange;
+  using Point_map = PointMap;
+  using Normal_map = NormalMap;
+  /// \endcond
+
+  /// Number type.
+  typedef typename GeomTraits::FT FT;
+
+  /// \cond SKIP_IN_MANUAL
+  using Point_2 = typename Traits::Point_2;
+  using Vector_2 = typename Traits::Vector_2;
+
+  using Squared_distance_2 = typename Traits::Compute_squared_distance_2;
+
+  using Get_sqrt = internal::Get_sqrt<Traits>;
+  using Sqrt = typename Get_sqrt::Sqrt;
+
+private:
+
+  const Input_range& m_input_range;
+
+  const FT m_distance_threshold;
+  const FT m_normal_threshold;
+  const std::size_t m_min_region_size;
+
+  const Point_map m_point_map;
+  const Normal_map m_normal_map;
+
+  const Squared_distance_2 m_squared_distance_2;
+  const Sqrt m_sqrt;
+
+  Point_2 m_center;
+  FT m_radius;
+
+public:
+
+  /// \endcond
+
+  /// @}
+
+  /// \name Initialization
+  /// @{
+
+  /*!
+    \brief initializes all internal data structures.
+
+    \param input_range an instance of `InputRange` with 2D points and
+    corresponding 2D normal vectors
+
+    \param distance_threshold the maximum distance from a point to a
+    circle. %Default is 1.
+
+    \param angle_threshold the maximum accepted angle in degrees
+    between the normal of a point and the radius of a circle. %Default
+    is 25 degrees.
+
+    \param min_region_size the minimum number of 2D points a region
+    must have. %Default is 3.
+
+    \param point_map an instance of `PointMap` that maps an item from
+    `input_range` to `Kernel::Point_3`
+
+    \param normal_map an instance of `NormalMap` that maps an item
+    from `input_range` to `Kernel::Vector_3`
+
+    \param traits an instance of `GeomTraits`.
+
+    \pre `input_range.size() > 0`
+    \pre `distance_threshold >= 0`
+    \pre `angle_threshold >= 0 && angle_threshold <= 90`
+    \pre `min_region_size > 0`
+  */
+  Least_squares_circle_fit_region  (const InputRange& input_range,
+                                    const FT distance_threshold = FT(1),
+                                    const FT angle_threshold = FT(25),
+                                    const std::size_t min_region_size = 3,
+                                    const PointMap point_map = PointMap(),
+                                    const NormalMap normal_map = NormalMap(),
+                                    const GeomTraits traits = GeomTraits())
+    : m_input_range(input_range)
+    , m_distance_threshold(distance_threshold)
+    , m_normal_threshold(static_cast<FT>(
+                           std::cos(
+                             CGAL::to_double(
+                               (angle_threshold * static_cast<FT>(CGAL_PI)) / FT(180))))),
+    m_min_region_size(min_region_size),
+    m_point_map(point_map),
+    m_normal_map(normal_map),
+    m_squared_distance_2(traits.compute_squared_distance_2_object()),
+    m_sqrt(Get_sqrt::sqrt_object(traits))
+  {
+    CGAL_precondition(input_range.size() > 0);
+    CGAL_precondition(distance_threshold >= FT(0));
+    CGAL_precondition(angle_threshold >= FT(0) && angle_threshold <= FT(90));
+    CGAL_precondition(min_region_size > 0);
+  }
+
+  /// @}
+
+  /// \name Access
+  /// @{
+
+  /*!
+    \brief implements `RegionType::is_part_of_region()`.
+
+    This function controls if a point with the index `query_index` is
+    within the `distance_threshold` from the corresponding circle and
+    if the angle between its normal and the circle radius is within
+    the `angle_threshold`.  If both conditions are satisfied, it
+    returns `true`, otherwise `false`.
+
+    \param query_index
+    index of the query point
+
+    The first and third parameters are not used in this implementation.
+
+    \return Boolean `true` or `false`
+
+    \pre `query_index >= 0 && query_index < input_range.size()`
+  */
+  bool is_part_of_region (const std::size_t,
+                          const std::size_t query_index,
+                          const std::vector<std::size_t>& indices) const
+  {
+    CGAL_precondition(query_index < m_input_range.size());
+
+    // First, we need to integrate at least 6 points so that the
+    // computed circle means something
+    if (indices.size() < 6)
+      return true;
+
+    const auto& key = *(m_input_range.begin() + query_index);
+    const Point_2& query_point = get(m_point_map, key);
+    Vector_2 normal = get(m_normal_map, key);
+
+    FT distance_to_center = m_sqrt(m_squared_distance_2 (query_point, m_center));
+    FT distance_to_circle = CGAL::abs (distance_to_center - m_radius);
+
+    if (distance_to_circle > m_distance_threshold)
+      return false;
+
+    normal = normal / m_sqrt (normal * normal);
+
+    Vector_2 ray (m_center, query_point);
+    ray = ray / m_sqrt (ray * ray);
+
+    if (CGAL::abs (normal * ray) < m_normal_threshold)
+      return false;
+
+    return true;
+  }
+
+  /*!
+    \brief implements `RegionType::is_valid_region()`.
+
+    This function controls if the `region` contains at least
+    `min_region_size` points.
+
+    \param region
+    indices of points included in the region
+
+    \return Boolean `true` or `false`
+  */
+  inline bool is_valid_region(const std::vector<std::size_t>& region) const
+  {
+    return (region.size() >= m_min_region_size);
+  }
+
+  /*!
+    \brief implements `RegionType::update()`.
+
+    This function fits the least squares circle to all points from the `region`.
+
+    \param region
+    indices of points included in the region
+
+    \pre `region.size() > 0`
+  */
+  void update(const std::vector<std::size_t>& region)
+  {
+    CGAL_precondition(region.size() > 0);
+
+    auto unary_function
+      = [&](const std::size_t& idx) -> const Point_2&
+        {
+          return get (m_point_map, *(m_input_range.begin() + idx));
+        };
+
+    // Use bbox to compute diagonalization with smaller coordinates to
+    // avoid loss of precision when inverting large coordinates
+    CGAL::Bbox_2 bbox = CGAL::bbox_2
+      (boost::make_transform_iterator
+       (region.begin(), unary_function),
+       boost::make_transform_iterator
+       (region.end(), unary_function));
+
+    using Diagonalize_traits = Eigen_diagonalize_traits<FT, 4>;
+    using Covariance_matrix = typename Diagonalize_traits::Covariance_matrix;
+    using Vector = typename Diagonalize_traits::Vector;
+    using Matrix = typename Diagonalize_traits::Matrix;
+
+    // Circle least squares fitting,
+    // Circle of center (a,b) and radius R
+    // Ri = sqrt((xi - a)^2 + (yi - b)^2)
+    // Minimize Sum(Ri^2 - R^2)^2
+    // -> Minimize Sum(xi^2 + yi^2 − 2 a*xi − 2 b*yi + a^2 + b^2 − R^2)^2
+    // let B=-2a ; C=-2b; D= a^2 + b^2 - R^2
+    // let ri = x^2 + y^2
+    // -> Minimize Sum(D + B*xi + C*yi + ri)^2
+    // -> Minimize Sum(1 + B/D*xi + C/D*yi + ri/D)^2
+    // -> system of linear equations
+    // -> diagonalize matrix
+    // -> center coordinates = -0.5 * eigenvector(1) / eigenvector(3) ; -0.5 * eigenvector(2) / eigenvector(3)
+    Covariance_matrix A
+      = { FT(0), FT(0), FT(0), FT(0), FT(0),
+          FT(0), FT(0), FT(0), FT(0), FT(0) };
+
+    A[0] = region.size();
+    for (const std::size_t& idx : region)
+    {
+      const auto& key = *(m_input_range.begin() + idx);
+      const Point_2& p = get(m_point_map, key);
+      FT x = p.x() - bbox.xmin();
+      FT y = p.y() - bbox.ymin();
+      FT r = x*x + y*y;
+      A[1] += x;
+      A[2] += y;
+      A[3] += r;
+      A[4] += x * x;
+      A[5] += x * y;
+      A[6] += x * r;
+      A[7] += y * y;
+      A[8] += y * r;
+      A[9] += r * r;
+    }
+
+    Vector eigenvalues = { FT(0), FT(0), FT(0), FT(0) };
+    Matrix eigenvectors = { FT(0), FT(0), FT(0), FT(0),
+                            FT(0), FT(0), FT(0), FT(0),
+                            FT(0), FT(0), FT(0), FT(0),
+                            FT(0), FT(0), FT(0), FT(0) };
+
+    Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix
+      (A, eigenvalues, eigenvectors);
+
+    m_center = Point_2 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[3]),
+                        bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[3]));
+
+    m_radius = FT(0);
+    for (const std::size_t& idx : region)
+    {
+      const auto& key = *(m_input_range.begin() + idx);
+      const Point_2& p = get(m_point_map, key);
+      m_radius += m_sqrt (m_squared_distance_2 (p, m_center));
+    }
+
+    m_radius /= region.size();
+  }
+
+  /// @}
+
+};
+
+} // namespace Point_set
+} // namespace Shape_detection
+} // namespace CGAL
+
+#endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_REGION_H

Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h

@@ -264,6 +264,9 @@ public:
         {
           return get (m_point_map, *(m_input_range.begin() + idx));
         };
+
+    // Use bbox to compute diagonalization with smaller coordinates to
+    // avoid loss of precision when inverting large coordinates
     CGAL::Bbox_3 bbox = CGAL::bbox_3
       (boost::make_transform_iterator
        (region.begin(), unary_function),
@@ -276,6 +279,7 @@ public:
     using Matrix = typename Diagonalize_traits::Matrix;
 
     // Sphere least squares fitting
+    // (see Least_square_circle_fit_region for details about computation)
     Covariance_matrix A
       = { FT(0), FT(0), FT(0), FT(0), FT(0),
           FT(0), FT(0), FT(0), FT(0), FT(0),