diff --git a/Shape_detection/examples/Shape_detection/CMakeLists.txt b/Shape_detection/examples/Shape_detection/CMakeLists.txt index ec2042109ec..89925216167 100644 --- a/Shape_detection/examples/Shape_detection/CMakeLists.txt +++ b/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) diff --git a/Shape_detection/examples/Shape_detection/data/circles.ply b/Shape_detection/examples/Shape_detection/data/circles.ply new file mode 100644 index 00000000000..e8767879066 Binary files /dev/null and b/Shape_detection/examples/Shape_detection/data/circles.ply differ diff --git a/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp b/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp new file mode 100644 index 00000000000..324af538dce --- /dev/null +++ b/Shape_detection/examples/Shape_detection/region_growing_circles_on_point_set_2.cpp @@ -0,0 +1,104 @@ +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +using Kernel = CGAL::Simple_cartesian; +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; +using Point_set_2 = CGAL::Point_set_3; +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 + ; +using Region_type = Shape_detection::Least_squares_circle_fit_region + ; +using Region_growing = CGAL::Shape_detection::Region_growing + ; + +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 + red = points3.add_property_map("red", 0).first, + green = points3.add_property_map("green", 0).first, + blue = points3.add_property_map("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& 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; +} diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h index 1c25058d84c..cdacadf1c5d 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h @@ -20,6 +20,7 @@ #include #include +#include #include #include diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_region.h new file mode 100644 index 00000000000..afd2f2c64c0 --- /dev/null +++ b/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 + +#include + +#include +#include +#include +#include + +#include + +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 +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; + 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( + std::cos( + CGAL::to_double( + (angle_threshold * static_cast(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& 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& 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& 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; + 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 diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h index 94a318e46be..a142d34b789 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_region.h +++ b/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),