diff --git a/BGL/include/CGAL/boost/graph/parameters_interface.h b/BGL/include/CGAL/boost/graph/parameters_interface.h index 06298b806f3..6ab2cb936a3 100644 --- a/BGL/include/CGAL/boost/graph/parameters_interface.h +++ b/BGL/include/CGAL/boost/graph/parameters_interface.h @@ -225,6 +225,8 @@ CGAL_add_named_parameter(sphere_radius_t, sphere_radius, sphere_radius) CGAL_add_named_parameter(k_neighbors_t, k_neighbors, k_neighbors) CGAL_add_named_parameter(segment_t, segment_map, segment_map) CGAL_add_named_parameter(cosine_value_t, cosine_value, cosine_value) +CGAL_add_named_parameter(minimum_radius_t, minimum_radius, minimum_radius) +CGAL_add_named_parameter(maximum_radius_t, maximum_radius, maximum_radius) // List of named parameters used in Shape_regularization package CGAL_add_named_parameter(minimum_length_t, minimum_length, minimum_length) 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 index 731ab9890d6..2ca5c92f180 100644 --- 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 @@ -56,27 +56,30 @@ int main (int argc, char** argv) // Default parameters for data/circles.ply const std::size_t k = 12; - const double tolerance = 0.01; + const double max_distance = 0.01; const double max_angle = 10.; const std::size_t min_region_size = 20; - // No constraint on radius - const double min_radius = 0.; - const double max_radius = std::numeric_limits::infinity(); - Neighbor_query neighbor_query( points, CGAL::parameters:: k_neighbors(k). point_map(points.point_map())); - Region_type region_type(points, tolerance, max_angle, min_region_size, - min_radius, max_radius, - points.point_map(), points.normal_map()); + Region_type region_type( + points, CGAL::parameters:: + maximum_distance(max_distance). + maximum_angle(max_angle). + minimum_region_size(min_region_size). + point_map(points.point_map()). + normal_map(points.normal_map())); // Sort indices - Sorting sorting(points, neighbor_query, points.point_map()); + Sorting sorting( + points, neighbor_query, CGAL::parameters:: + point_map(points.point_map())); sorting.sort(); - Region_growing region_growing(points, neighbor_query, region_type, sorting.seed_map()); + Region_growing region_growing( + points, neighbor_query, region_type, sorting.seed_map()); // Add maps to get colored output Point_set_3::Property_map 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 index d6df53529ab..7612ababa26 100644 --- 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 @@ -16,294 +16,401 @@ #include -#include -#include - -#include -#include -#include - +// Internal includes. #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 FT m_min_radius; - const FT m_max_radius; - - 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. + \ingroup PkgShapeDetectionRGOnPoints - \param input_range an instance of `InputRange` with 2D points and - corresponding 2D normal vectors + \brief Region type based on the quality of the least squares circle + fit applied to 2D points. - \param distance_threshold the maximum distance from a point to a - circle. %Default is 1. + 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. - \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. + \tparam GeomTraits + a model of `Kernel` - \param min_region_size the minimum number of 2D points a region - must have. %Default is 3. + \tparam InputRange + a model of `ConstRange` whose iterator type is `RandomAccessIterator` - \param minimum_radius the radius below which an estimated circle - is considered as invalid and discarded. %Default is 0 (no limit). + \tparam PointMap + a model of `ReadablePropertyMap` whose key type is the value type of the input + range and value type is `Kernel::Point_2` - \param maximum_radius the radius above which an estimated circle - is considered as invalid and discarded. %Default is infinity (no - limit). + \tparam NormalMap + a model of `ReadablePropertyMap` whose key type is the value type of the input + range and value type is `Kernel::Vector_2` - \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` + \cgalModels `RegionType` */ - 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 FT minimum_radius = FT(0), - const FT maximum_radius = std::numeric_limits::infinity(), - 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_min_radius (minimum_radius) - , m_max_radius (maximum_radius) - , 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); - CGAL_precondition(minimum_radius >= 0); - CGAL_precondition(maximum_radius > minimum_radius); - } + template< + typename GeomTraits, + typename InputRange, + typename PointMap, + typename NormalMap> + class Least_squares_circle_fit_region { - /// @} + public: - /// \name Access - /// @{ + /// \name Types + /// @{ - /*! - \brief implements `RegionType::is_part_of_region()`. + /// \cond SKIP_IN_MANUAL + using Traits = GeomTraits; + using Input_range = InputRange; + using Point_map = PointMap; + using Normal_map = NormalMap; + /// \endcond - 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`. + /// Number type. + typedef typename GeomTraits::FT FT; - \param query_index - index of the query point + /// @} - \param indices indices of the inliers of the region + private: + using Point_2 = typename Traits::Point_2; + using Vector_2 = typename Traits::Vector_2; - The first parameter is not used in this implementation. + using Squared_distance_2 = typename Traits::Compute_squared_distance_2; + using Get_sqrt = internal::Get_sqrt; + using Sqrt = typename Get_sqrt::Sqrt; - \return Boolean `true` or `false` + public: - \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()); + /// @} + + /// \name Initialization + /// @{ + + /*! + \brief initializes all internal data structures. + + \tparam NamedParameters + a sequence of \ref bgl_namedparameters "Named Parameters" + + \param input_range + an instance of `InputRange` with 2D points and + corresponding 2D normal vectors + + \param np + a sequence of \ref bgl_namedparameters "Named Parameters" + among the ones listed below + + \cgalNamedParamsBegin + \cgalParamNBegin{maximum_distance} + \cgalParamDescription{the maximum distance from a point to a circle} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{1} + \cgalParamNEnd + \cgalParamNBegin{maximum_angle} + \cgalParamDescription{the maximum angle in degrees between + the normal of a point and the radius of a circle} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{25 degrees} + \cgalParamNEnd + \cgalParamNBegin{cosine_value} + \cgalParamDescription{the cos value computed as `cos(maximum_angle * PI / 180)`, + this parameter can be used instead of the `maximum_angle`} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{`cos(25 * PI / 180)`} + \cgalParamNEnd + \cgalParamNBegin{minimum_region_size} + \cgalParamDescription{the minimum number of 2D points a region must have} + \cgalParamType{`std::size_t`} + \cgalParamDefault{3} + \cgalParamNEnd + \cgalParamNBegin{minimum_radius} + \cgalParamDescription{the radius below which an estimated circle + is considered as invalid and discarded} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{0, no limit} + \cgalParamNEnd + \cgalParamNBegin{maximum_radius} + \cgalParamDescription{the radius above which an estimated circle + is considered as invalid and discarded.} + \cgalParamType{`GeomTraits::FT`} + \cgalParamDefault{+infinity, no limit} + \cgalParamNEnd + \cgalParamNBegin{point_map} + \cgalParamDescription{an instance of `PointMap` that maps an item from `input_range` + to `Kernel::Point_2`} + \cgalParamDefault{`PointMap()`} + \cgalParamNEnd + \cgalParamNBegin{normal_map} + \cgalParamDescription{an instance of `NormalMap` that maps an item from `input_range` + to `Kernel::Vector_2`} + \cgalParamDefault{`NormalMap()`} + \cgalParamNEnd + \cgalParamNBegin{geom_traits} + \cgalParamDescription{an instance of `GeomTraits`} + \cgalParamDefault{`GeomTraits()`} + \cgalParamNEnd + \cgalNamedParamsEnd + + \pre `input_range.size() > 0` + \pre `maximum_distance >= 0` + \pre `maximum_angle >= 0 && maximum_angle <= 90` + \pre `cosine_value >= 0 && cosine_value <= 1` + \pre `minimum_region_size > 0` + \pre `minimum_radius >= 0` + \pre `maximum_radius >= minimum_radius` + */ + template + Least_squares_circle_fit_region( + const InputRange& input_range, + const NamedParameters& np) : + m_input_range(input_range), + m_point_map(parameters::choose_parameter(parameters::get_parameter( + np, internal_np::point_map), PointMap())), + m_normal_map(parameters::choose_parameter(parameters::get_parameter( + np, internal_np::normal_map), NormalMap())), + m_traits(parameters::choose_parameter(parameters::get_parameter( + np, internal_np::geom_traits), GeomTraits())), + m_squared_distance_2(m_traits.compute_squared_distance_2_object()), + m_sqrt(Get_sqrt::sqrt_object(m_traits)) { + + CGAL_precondition(input_range.size() > 0); + const FT max_distance = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::maximum_distance), FT(1)); + CGAL_precondition(max_distance >= FT(0)); + m_distance_threshold = max_distance; + + const FT max_angle = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::maximum_angle), FT(25)); + CGAL_precondition(max_angle >= FT(0) && max_angle <= FT(90)); + + m_min_region_size = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::minimum_region_size), 3); + CGAL_precondition(m_min_region_size > 0); + + const FT default_cos_value = static_cast(std::cos(CGAL::to_double( + (max_angle * static_cast(CGAL_PI)) / FT(180)))); + const FT cos_value = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::cosine_value), default_cos_value); + CGAL_precondition(cos_value >= FT(0) && cos_value <= FT(1)); + m_cos_value_threshold = cos_value; + + const FT m_min_radius = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::minimum_radius), FT(0)); + CGAL_precondition(m_min_radius >= FT(0)); + + const FT m_max_radius = parameters::choose_parameter( + parameters::get_parameter(np, internal_np::maximum_radius), FT(std::numeric_limits::max())); + CGAL_precondition(m_max_radius >= m_min_radius); + } + + /*! + \brief initializes all internal data structures. + + \deprecated This constructor is deprecated since the version 5.5 of \cgal. + + \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 minimum_radius + the radius below which an estimated circle is considered as invalid + and discarded. %Default is 0 (no limit). + + \param maximum_radius + the radius above which an estimated circle is considered as invalid + and discarded. %Default is infinity (no limit). + + \param point_map + an instance of `PointMap` that maps an item from `input_range` to `Kernel::Point_2` + + \param normal_map + an instance of `NormalMap` that maps an item from `input_range` to `Kernel::Vector_2` + + \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` + \pre `minimum_radius >= 0` + \pre `maximum_radius >= minimum_radius` + */ + CGAL_DEPRECATED_MSG("This constructor is deprecated since the version 5.5 of CGAL!") + 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 FT minimum_radius = FT(0), + const FT maximum_radius = FT(std::numeric_limits::max()), + const PointMap point_map = PointMap(), + const NormalMap normal_map = NormalMap(), + const GeomTraits traits = GeomTraits()) : + Least_squares_circle_fit_region( + input_range, CGAL::parameters:: + maximum_distance(distance_threshold). + maximum_angle(angle_threshold). + minimum_region_size(min_region_size). + minimum_radius(minimum_radius). + maximum_radius(maximum_radius). + point_map(point_map). + normal_map(normal_map). + geom_traits(traits)) + { } + + /// @} + + /// \name Access + /// @{ + + /*! + \brief implements `RegionType::is_part_of_region()`. + + This function controls if a point with the index `query_index` is within + the `maximum_distance` from the corresponding circle and if the angle between + its normal and the circle radius is within the `maximum_angle`. If both conditions + are satisfied, it returns `true`, otherwise `false`. + + \param query_index + index of the query point + + \param indices indices of the inliers of the region + + The first parameter is not used in this implementation. + + \return Boolean `true` or `false` + + \pre `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; + + // TODO: Why do we get so many nan in this class? + if (std::isnan(m_radius)) + return false; + + // If radius is out of bound, nothing fits, early ending + if (m_radius < m_min_radius || m_radius > m_max_radius) + return false; + + 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); + + const FT sq_dist = m_squared_distance_2(query_point, m_center); + if (std::isnan(sq_dist)) return false; + FT distance_to_center = m_sqrt (sq_dist); + FT distance_to_circle = CGAL::abs (distance_to_center - m_radius); + + if (distance_to_circle > m_distance_threshold) + return false; + + const FT sq_norm = normal * normal; + if (std::isnan(sq_norm)) return false; + normal = normal / m_sqrt (sq_norm); + + Vector_2 ray (m_center, query_point); + const FT sq_ray = ray * ray; + if (std::isnan(sq_ray)) return false; + ray = ray / m_sqrt (sq_ray); + + if (CGAL::abs (normal * ray) < m_cos_value_threshold) + return false; - // First, we need to integrate at least 6 points so that the - // computed circle means something - if (indices.size() < 6) return true; + } - // TODO: Why do we get so many nan in this class? - if (std::isnan(m_radius)) - return false; + /*! + \brief implements `RegionType::is_valid_region()`. - // If radius is out of bound, nothing fits, early ending - if (m_radius < m_min_radius || m_radius > m_max_radius) - return false; + This function controls if the estimated radius is between `minimum_radius` + and `maximum_radius` and if the `region` contains at least `min_region_size` points. - 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); + \param region + indices of points included in the region - const FT sq_dist = m_squared_distance_2(query_point, m_center); - if (std::isnan(sq_dist)) return false; - FT distance_to_center = m_sqrt (sq_dist); - FT distance_to_circle = CGAL::abs (distance_to_center - m_radius); + \return Boolean `true` or `false` + */ + inline bool is_valid_region(const std::vector& region) const { + return ( + (m_min_radius <= m_radius && m_radius <= m_max_radius) && + (region.size() >= m_min_region_size) + ); + } - if (distance_to_circle > m_distance_threshold) - return false; + /*! + \brief implements `RegionType::update()`. - const FT sq_norm = normal * normal; - if (std::isnan(sq_norm)) return false; - normal = normal / m_sqrt (sq_norm); + This function fits the least squares circle to all points from the `region`. - Vector_2 ray (m_center, query_point); - const FT sq_ray = ray * ray; - if (std::isnan(sq_ray)) return false; - ray = ray / m_sqrt (sq_ray); + \param region + indices of points included in the region - if (CGAL::abs (normal * ray) < m_normal_threshold) - return false; + \return Boolean `true` if the circle fitting succeeded and `false` otherwise - return true; - } + \pre `region.size() > 0` + */ + bool update(const std::vector& region) { - /*! - \brief implements `RegionType::is_valid_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)); + }; - This function controls if the estimated radius is between - `minimum_radius` and `maximum_radius` and if the `region` contains - at least `min_region_size` points. + internal::circle_fit + (make_range(boost::make_transform_iterator + (region.begin(), unary_function), + boost::make_transform_iterator + (region.end(), unary_function)), + m_sqrt, m_squared_distance_2, m_center, m_radius); + return true; + } - \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 ((m_min_radius <= m_radius && m_radius <= m_max_radius) - && (region.size() >= m_min_region_size)); - } + private: + const Input_range& m_input_range; + const Point_map m_point_map; + const Normal_map m_normal_map; + const Traits m_traits; - /*! - \brief implements `RegionType::update()`. + FT m_distance_threshold; + FT m_cos_value_threshold; + std::size_t m_min_region_size; + FT m_min_radius; + FT m_max_radius; - This function fits the least squares circle to all points from the `region`. + const Squared_distance_2 m_squared_distance_2; + const Sqrt m_sqrt; - \param region - indices of points included in the region - - \pre `region.size() > 0` - */ - bool 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)); - }; - - internal::circle_fit - (make_range(boost::make_transform_iterator - (region.begin(), unary_function), - boost::make_transform_iterator - (region.end(), unary_function)), - m_sqrt, m_squared_distance_2, m_center, m_radius); - return true; - } - - /// @} - -}; + Point_2 m_center; + FT m_radius; + }; } // namespace Point_set } // namespace Shape_detection diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h index 4b60f1b2de0..64b7f191727 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_circle_fit_sorting.h @@ -1,4 +1,4 @@ -// Copyright (c) 2018 INRIA Sophia-Antipolis (France). +// Copyright (c) 2020 INRIA Sophia-Antipolis (France). // All rights reserved. // // This file is part of CGAL (www.cgal.org). @@ -8,7 +8,7 @@ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // // -// Author(s) : Florent Lafarge, Simon Giraudot, Thien Hoang, Dmitry Anisimov +// Author(s) : Simon Giraudot // #ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_CIRCLE_FIT_SORTING_H @@ -16,201 +16,233 @@ #include -// STL includes. -#include -#include - -// CGAL includes. -#include -#include -#include - // Internal includes. -#include #include -#include namespace CGAL { namespace Shape_detection { namespace Point_set { -/*! - \ingroup PkgShapeDetectionRGOnPoints + /*! + \ingroup PkgShapeDetectionRGOnPoints - \brief Sorting of 2D points with respect to the local circle fit quality. + \brief Sorting of 2D points with respect to the local circle fit quality. - Indices of 2D input points are sorted with respect to the quality of the - least squares circle fit applied to the neighboring points of each point. + Indices of 2D input points are sorted with respect to the quality of the + least squares circle fit applied to the neighboring points of each point. - \tparam GeomTraits - must be a model of `Kernel`. + \tparam GeomTraits + a model of `Kernel` - \tparam InputRange - must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`. + \tparam InputRange + a model of `ConstRange` whose iterator type is `RandomAccessIterator` - \tparam NeighborQuery - must be a model of `NeighborQuery`. + \tparam NeighborQuery + a model of `NeighborQuery` - \tparam PointMap - must be an `LvaluePropertyMap` whose key type is the value type of the input - range and value type is `Kernel::Point_2`. -*/ -template< + \tparam PointMap + a model of `ReadablePropertyMap` whose key type is the value type of the input + range and value type is `Kernel::Point_2` + */ + template< typename GeomTraits, typename InputRange, typename NeighborQuery, typename PointMap> -class Least_squares_circle_fit_sorting { + class Least_squares_circle_fit_sorting { -public: + public: + /// \name Types + /// @{ - /// \name Types - /// @{ + /// \cond SKIP_IN_MANUAL + using Traits = GeomTraits; + using Input_range = InputRange; + using Neighbor_query = NeighborQuery; + using Point_map = PointMap; + using Seed_map = internal::Seed_property_map; + /// \endcond - /// \cond SKIP_IN_MANUAL - using Traits = GeomTraits; - using Input_range = InputRange; - using Neighbor_query = NeighborQuery; - using Point_map = PointMap; - using Seed_map = internal::Seed_property_map; - /// \endcond + #ifdef DOXYGEN_RUNNING + /*! + a model of `ReadablePropertyMap` whose key and value type is `std::size_t`. + This map provides an access to the ordered indices of input points. + */ + typedef unspecified_type Seed_map; + #endif -#ifdef DOXYGEN_RUNNING - /*! - an `LvaluePropertyMap` whose key and value type is `std::size_t`. - This map provides an access to the ordered indices of input points. - */ - typedef unspecified_type Seed_map; -#endif + /// @} - /// @} + private: + // using FT = typename Traits::FT; + // using Compare_scores = internal::Compare_scores; - /// \name Initialization - /// @{ + using Local_traits = Exact_predicates_inexact_constructions_kernel; + using Local_FT = typename Local_traits::FT; + using Local_point_2 = typename Local_traits::Point_2; + using To_local_converter = Cartesian_converter; + using Compare_scores = internal::Compare_scores; - /*! - \brief initializes all internal data structures. + public: + /// \name Initialization + /// @{ - \param input_range - an instance of `InputRange` with 2D points + /*! + \brief initializes all internal data structures. - \param neighbor_query - an instance of `NeighborQuery` that is used internally to - access point's neighbors + \tparam NamedParameters + a sequence of \ref bgl_namedparameters "Named Parameters" - \param point_map - an instance of `PointMap` that maps an item from `input_range` - to `Kernel::Point_2` + \param input_range + an instance of `InputRange` with 2D points - \pre `input_range.size() > 0` - */ - Least_squares_circle_fit_sorting( - const InputRange& input_range, - NeighborQuery& neighbor_query, - const PointMap point_map = PointMap()) : + \param neighbor_query + an instance of `NeighborQuery` that is used internally to + access point's neighbors + + \param np + a sequence of \ref bgl_namedparameters "Named Parameters" + among the ones listed below + + \cgalNamedParamsBegin + \cgalParamNBegin{point_map} + \cgalParamDescription{an instance of `PointMap` that maps an item from `input_range` + to `Kernel::Point_2`} + \cgalParamDefault{`PointMap()`} + \cgalParamNEnd + \cgalParamNBegin{geom_traits} + \cgalParamDescription{an instance of `GeomTraits`} + \cgalParamDefault{`GeomTraits()`} + \cgalParamNEnd + \cgalNamedParamsEnd + + \pre `input_range.size() > 0` + */ + template + Least_squares_circle_fit_sorting( + const InputRange& input_range, + NeighborQuery& neighbor_query, + const NamedParameters& np) : m_input_range(input_range), m_neighbor_query(neighbor_query), - m_point_map(point_map), + m_point_map(parameters::choose_parameter(parameters::get_parameter( + np, internal_np::point_map), PointMap())), + m_traits(parameters::choose_parameter(parameters::get_parameter( + np, internal_np::geom_traits), GeomTraits())), m_to_local_converter() { - CGAL_precondition(input_range.size() > 0); - - m_order.resize(m_input_range.size()); - for (std::size_t i = 0; i < m_input_range.size(); ++i) - m_order[i] = i; - m_scores.resize(m_input_range.size()); - } - - /// @} - - /// \name Sorting - /// @{ - - /*! - \brief sorts indices of input points. - */ - void sort() { - - compute_scores(); - CGAL_postcondition(m_scores.size() > 0); - - Compare_scores cmp(m_scores); - std::sort(m_order.begin(), m_order.end(), cmp); - } - - /// @} - - /// \name Access - /// @{ - - /*! - \brief returns an instance of `Seed_map` to access the ordered indices - of input points. - */ - Seed_map seed_map() { - return Seed_map(m_order); - } - - /// @} - -private: - - // Types. - using Local_traits = Exact_predicates_inexact_constructions_kernel; - using Local_FT = typename Local_traits::FT; - using Local_point_2 = typename Local_traits::Point_2; - using To_local_converter = Cartesian_converter; - using Compare_scores = internal::Compare_scores; - - // Functions. - void compute_scores() { - - std::vector neighbors; - std::vector points; - - typename internal::Get_sqrt::Sqrt sqrt; - typename Local_traits::Compute_squared_distance_2 squared_distance_2; - - for (std::size_t i = 0; i < m_input_range.size(); ++i) { - - neighbors.clear(); - m_neighbor_query(i, neighbors); - neighbors.push_back(i); - - points.clear(); - for (std::size_t j = 0; j < neighbors.size(); ++j) { - CGAL_precondition(neighbors[j] < m_input_range.size()); - - const auto& key = *(m_input_range.begin() + neighbors[j]); - points.push_back(m_to_local_converter(get(m_point_map, key))); - } - CGAL_postcondition(points.size() == neighbors.size()); - - Local_point_2 fitted_center; - Local_FT fitted_radius; - - if (internal::circle_fit (points, sqrt, squared_distance_2, fitted_center, fitted_radius)) - { - // Score is min squared distance to sphere - m_scores[i] = Local_FT(0); - for (const Local_point_2& p : points) - m_scores[i] += abs (sqrt(squared_distance_2(p, fitted_center)) - fitted_radius); - } - else - m_scores[i] = Local_FT(std::numeric_limits::infinity()); + CGAL_precondition(input_range.size() > 0); + m_order.resize(m_input_range.size()); + std::iota(m_order.begin(), m_order.end(), 0); + m_scores.resize(m_input_range.size()); } - } - // Fields. - const Input_range& m_input_range; - Neighbor_query& m_neighbor_query; - const Point_map m_point_map; + /*! + \brief initializes all internal data structures. - std::vector m_order; - std::vector m_scores; + \deprecated This constructor is deprecated since the version 5.5 of \cgal. - const To_local_converter m_to_local_converter; -}; + \param input_range + an instance of `InputRange` with 2D points + + \param neighbor_query + an instance of `NeighborQuery` that is used internally to + access point's neighbors + + \param point_map + an instance of `PointMap` that maps an item from `input_range` + to `Kernel::Point_2` + + \pre `input_range.size() > 0` + */ + CGAL_DEPRECATED_MSG("This constructor is deprecated since the version 5.5 of CGAL!") + Least_squares_circle_fit_sorting( + const InputRange& input_range, + NeighborQuery& neighbor_query, + const PointMap point_map = PointMap()) : + Least_squares_circle_fit_sorting( + input_range, neighbor_query, CGAL::parameters:: + point_map(point_map)) + { } + + /// @} + + /// \name Sorting + /// @{ + + /*! + \brief sorts indices of input points. + */ + void sort() { + + compute_scores(); + CGAL_precondition(m_scores.size() > 0); + Compare_scores cmp(m_scores); + std::sort(m_order.begin(), m_order.end(), cmp); + } + + /// @} + + /// \name Access + /// @{ + + /*! + \brief returns an instance of `Seed_map` to access the ordered indices + of input points. + */ + Seed_map seed_map() { + return Seed_map(m_order); + } + + /// @} + + private: + const Input_range& m_input_range; + Neighbor_query& m_neighbor_query; + const Point_map m_point_map; + const Traits m_traits; + std::vector m_order; + std::vector m_scores; + const To_local_converter m_to_local_converter; + + void compute_scores() { + + std::vector neighbors; + std::vector points; + + typename internal::Get_sqrt::Sqrt sqrt; + typename Local_traits::Compute_squared_distance_2 squared_distance_2; + + for (std::size_t i = 0; i < m_input_range.size(); ++i) { + + neighbors.clear(); + m_neighbor_query(i, neighbors); + neighbors.push_back(i); + + points.clear(); + for (std::size_t j = 0; j < neighbors.size(); ++j) { + CGAL_precondition(neighbors[j] < m_input_range.size()); + + const auto& key = *(m_input_range.begin() + neighbors[j]); + points.push_back(m_to_local_converter(get(m_point_map, key))); + } + CGAL_postcondition(points.size() == neighbors.size()); + + Local_point_2 fitted_center; + Local_FT fitted_radius; + + if (internal::circle_fit (points, sqrt, squared_distance_2, fitted_center, fitted_radius)) + { + // Score is min squared distance to sphere + m_scores[i] = Local_FT(0); + for (const Local_point_2& p : points) + m_scores[i] += abs (sqrt(squared_distance_2(p, fitted_center)) - fitted_radius); + } + else + m_scores[i] = Local_FT(std::numeric_limits::infinity()); + } + } + }; } // namespace Point_set } // namespace Shape_detection diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h index 406dc93b27d..d69e2a04410 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_region.h @@ -24,7 +24,6 @@ #include #include -#include namespace CGAL { namespace Shape_detection { diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h index 7c01bc7a752..b05a2f0011c 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_cylinder_fit_sorting.h @@ -28,7 +28,6 @@ // Internal includes. #include #include -#include namespace CGAL { namespace Shape_detection { 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 3e430d18cf7..43cfde3c5fb 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 @@ -24,7 +24,6 @@ #include #include -#include namespace CGAL { namespace Shape_detection { diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h index 988a4e07d74..a28f131be38 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Least_squares_sphere_fit_sorting.h @@ -28,7 +28,6 @@ // Internal includes. #include #include -#include namespace CGAL { namespace Shape_detection { diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/fitting.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/fitting.h deleted file mode 100644 index b89ac3a0455..00000000000 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/fitting.h +++ /dev/null @@ -1,265 +0,0 @@ -// Copyright (c) 2021 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_INTERNAL_FITTING_H -#define CGAL_SHAPE_DETECTION_REGION_GROWING_INTERNAL_FITTING_H - -#include - -#include -#include - -namespace CGAL { -namespace Shape_detection { -namespace internal { - -template -bool circle_fit (const PointRange& points, - const Sqrt& sqrt, - const Squared_distance_2& squared_distance_2, - Point_2& center, FT& radius) -{ - 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; - - // Use bbox to compute diagonalization with smaller coordinates to - // avoid loss of precision when inverting large coordinates - Bbox_2 bbox = bbox_2 (points.begin(), points.end()); - - // 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 - // NB x y r - // xx xy xr - // yy yr - // rr - // - // -> 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] = static_cast(points.size()); - for (const Point_2& p : points) - { - 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); - - // Perfect line case, no circle can be fitted - if (eigenvectors[3] == 0) - return false; - - center = Point_2 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[3]), - bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[3])); - - radius = FT(0); - for (const Point_2& p : points) - radius += sqrt (squared_distance_2 (p, center)); - radius /= points.size(); - - return true; -} - -template -bool sphere_fit (const PointRange& points, - const Sqrt& sqrt, - const Squared_distance_3& squared_distance_3, - Point_3& center, FT& radius) -{ - 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; - - // Use bbox to compute diagonalization with smaller coordinates to - // avoid loss of precision when inverting large coordinates - Bbox_3 bbox = bbox_3 (points.begin(), points.end()); - - // 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), - FT(0), FT(0), FT(0), FT(0), FT(0) }; - - A[0] = static_cast(points.size()); - for (const Point_3& p : points) - { - FT x = p.x() - bbox.xmin(); - FT y = p.y() - bbox.ymin(); - FT z = p.z() - bbox.zmin(); - FT r = x*x + y*y + z*z; - A[1] += x; - A[2] += y; - A[3] += z; - A[4] += r; - A[5] += x * x; - A[6] += x * y; - A[7] += x * z; - A[8] += x * r; - A[9] += y * y; - A[10] += y * z; - A[11] += y * r; - A[12] += z * z; - A[13] += z * r; - A[14] += r * r; - } - - Vector eigenvalues = { FT(0), 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), 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); - - // Perfect plane case, no sphere can be fitted - if (eigenvectors[4] == 0) - return false; - - center = Point_3 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[4]), - bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[4]), - bbox.zmin() - FT(0.5) * (eigenvectors[3] / eigenvectors[4])); - - radius = FT(0); - for (const Point_3& p : points) - radius += sqrt (squared_distance_3 (p, center)); - radius /= points.size(); - - return true; -} - -template -bool cylinder_fit (const PointRange& points, - PointMap point_map, NormalMap normal_map, - const Sqrt& sqrt, - const Squared_distance_3& squared_distance_3, - Line_3& axis, FT& radius) -{ - using Point_3 = typename boost::property_traits::value_type; - using Vector_3 = typename boost::property_traits::value_type; - - const Point_3& ref = get(point_map, *(points.begin())); - - Vector_3 mean_axis = CGAL::NULL_VECTOR; - std::size_t nb = 0; - radius = FT(0); - Point_3 point_on_axis = ORIGIN; - for (std::size_t i = 0; i < points.size() - 1; ++ i) - { - Vector_3 v0 = get(normal_map, *std::next(points.begin(), i)); - v0 = v0 / sqrt(v0*v0); - Vector_3 v1 = get(normal_map, *std::next(points.begin(), i + 1)); - v1 = v1 / sqrt(v1*v1); - Vector_3 axis = cross_product (v0, v1); - if (sqrt(axis.squared_length()) < (FT)(0.01)) - continue; - axis = axis / sqrt(axis * axis); - - const Point_3& p0 = get(point_map, *std::next(points.begin(), i)); - const Point_3& p1 = get(point_map, *std::next(points.begin(), i + 1)); - - Vector_3 xdir = v0 - axis * (v0 * axis); - xdir = xdir / sqrt (xdir * xdir); - - Vector_3 ydir = cross_product (axis, xdir); - ydir = ydir / sqrt (ydir * ydir); - - FT v1x = v1 * ydir; - FT v1y = -v1 * xdir; - - Vector_3 d (p0, p1); - - FT ox = xdir * d; - FT oy = ydir * d; - FT ldist = v1x * ox + v1y * oy; - - FT radius = ldist / v1x; - - Point_3 point = p0 + xdir * radius; - Line_3 line (point, axis); - point = line.projection(ref); - - point_on_axis = barycenter (point_on_axis, static_cast(nb), point, FT(1)); - - radius += abs(radius); - - if (nb != 0 && axis * mean_axis < 0) - axis = -axis; - - mean_axis = mean_axis + axis; - ++ nb; - } - - if (nb == 0) - return false; - - mean_axis = mean_axis / sqrt(mean_axis * mean_axis); - axis = Line_3 (point_on_axis, mean_axis); - radius /= nb; - - radius = FT(0); - for (const auto& p : points) - { - const Point_3& p0 = get(point_map, p); - radius += sqrt(squared_distance_3(p0, axis)); - } - radius /= points.size(); - - return true; -} - -} // internal -} // namespace Shape_detection -} // namespace CGAL - - -#endif // CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_INTERNAL_FITTING_H diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/region_growing_traits.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/region_growing_traits.h index 5396d2e797f..7c7e38a1983 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/region_growing_traits.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/region_growing_traits.h @@ -16,15 +16,6 @@ #include -// STL includes. -#include -#include -#include - -// CGAL includes. -#include -#include - // Internal includes. #include diff --git a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h index 30037f04338..ef83fb846ba 100644 --- a/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h +++ b/Shape_detection/include/CGAL/Shape_detection/Region_growing/internal/utils.h @@ -19,6 +19,7 @@ // STL includes. #include #include +#include #include #include #include @@ -345,6 +346,240 @@ namespace internal { return std::make_pair(plane, static_cast(score)); } +template +bool circle_fit (const PointRange& points, + const Sqrt& sqrt, + const Squared_distance_2& squared_distance_2, + Point_2& center, FT& radius) +{ + 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; + + // Use bbox to compute diagonalization with smaller coordinates to + // avoid loss of precision when inverting large coordinates + Bbox_2 bbox = bbox_2 (points.begin(), points.end()); + + // 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 + // NB x y r + // xx xy xr + // yy yr + // rr + // + // -> 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] = static_cast(points.size()); + for (const Point_2& p : points) + { + 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); + + // Perfect line case, no circle can be fitted + if (eigenvectors[3] == 0) + return false; + + center = Point_2 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[3]), + bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[3])); + + radius = FT(0); + for (const Point_2& p : points) + radius += sqrt (squared_distance_2 (p, center)); + radius /= points.size(); + + return true; +} + +template +bool sphere_fit (const PointRange& points, + const Sqrt& sqrt, + const Squared_distance_3& squared_distance_3, + Point_3& center, FT& radius) +{ + 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; + + // Use bbox to compute diagonalization with smaller coordinates to + // avoid loss of precision when inverting large coordinates + Bbox_3 bbox = bbox_3 (points.begin(), points.end()); + + // 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), + FT(0), FT(0), FT(0), FT(0), FT(0) }; + + A[0] = static_cast(points.size()); + for (const Point_3& p : points) + { + FT x = p.x() - bbox.xmin(); + FT y = p.y() - bbox.ymin(); + FT z = p.z() - bbox.zmin(); + FT r = x*x + y*y + z*z; + A[1] += x; + A[2] += y; + A[3] += z; + A[4] += r; + A[5] += x * x; + A[6] += x * y; + A[7] += x * z; + A[8] += x * r; + A[9] += y * y; + A[10] += y * z; + A[11] += y * r; + A[12] += z * z; + A[13] += z * r; + A[14] += r * r; + } + + Vector eigenvalues = { FT(0), 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), 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); + + // Perfect plane case, no sphere can be fitted + if (eigenvectors[4] == 0) + return false; + + center = Point_3 (bbox.xmin() - FT(0.5) * (eigenvectors[1] / eigenvectors[4]), + bbox.ymin() - FT(0.5) * (eigenvectors[2] / eigenvectors[4]), + bbox.zmin() - FT(0.5) * (eigenvectors[3] / eigenvectors[4])); + + radius = FT(0); + for (const Point_3& p : points) + radius += sqrt (squared_distance_3 (p, center)); + radius /= points.size(); + + return true; +} + +template +bool cylinder_fit (const PointRange& points, + PointMap point_map, NormalMap normal_map, + const Sqrt& sqrt, + const Squared_distance_3& squared_distance_3, + Line_3& axis, FT& radius) +{ + using Point_3 = typename boost::property_traits::value_type; + using Vector_3 = typename boost::property_traits::value_type; + + const Point_3& ref = get(point_map, *(points.begin())); + + Vector_3 mean_axis = CGAL::NULL_VECTOR; + std::size_t nb = 0; + radius = FT(0); + Point_3 point_on_axis = ORIGIN; + for (std::size_t i = 0; i < points.size() - 1; ++ i) + { + Vector_3 v0 = get(normal_map, *std::next(points.begin(), i)); + v0 = v0 / sqrt(v0*v0); + Vector_3 v1 = get(normal_map, *std::next(points.begin(), i + 1)); + v1 = v1 / sqrt(v1*v1); + Vector_3 axis = cross_product (v0, v1); + if (sqrt(axis.squared_length()) < (FT)(0.01)) + continue; + axis = axis / sqrt(axis * axis); + + const Point_3& p0 = get(point_map, *std::next(points.begin(), i)); + const Point_3& p1 = get(point_map, *std::next(points.begin(), i + 1)); + + Vector_3 xdir = v0 - axis * (v0 * axis); + xdir = xdir / sqrt (xdir * xdir); + + Vector_3 ydir = cross_product (axis, xdir); + ydir = ydir / sqrt (ydir * ydir); + + FT v1x = v1 * ydir; + FT v1y = -v1 * xdir; + + Vector_3 d (p0, p1); + + FT ox = xdir * d; + FT oy = ydir * d; + FT ldist = v1x * ox + v1y * oy; + + FT radius = ldist / v1x; + + Point_3 point = p0 + xdir * radius; + Line_3 line (point, axis); + point = line.projection(ref); + + point_on_axis = barycenter (point_on_axis, static_cast(nb), point, FT(1)); + + radius += abs(radius); + + if (nb != 0 && axis * mean_axis < 0) + axis = -axis; + + mean_axis = mean_axis + axis; + ++ nb; + } + + if (nb == 0) + return false; + + mean_axis = mean_axis / sqrt(mean_axis * mean_axis); + axis = Line_3 (point_on_axis, mean_axis); + radius /= nb; + + radius = FT(0); + for (const auto& p : points) + { + const Point_3& p0 = get(point_map, p); + radius += sqrt(squared_distance_3(p0, axis)); + } + radius /= points.size(); + + return true; +} + } // namespace internal } // namespace Shape_detection } // namespace CGAL