adding np to sphere classes

This commit is contained in:
Dmitry Anisimov 2021-09-22 15:27:01 +02:00
parent 0fbc3ea1bb
commit 2c16e6f9bf
4 changed files with 575 additions and 433 deletions

View File

@ -42,22 +42,24 @@ int main (int argc, char** argv)
// Default parameters for data/spheres.ply // Default parameters for data/spheres.ply
const std::size_t k = 12; const std::size_t k = 12;
const double tolerance = 0.01; const double max_distance = 0.01;
const double max_angle = 10.; const double max_angle = 10.;
const std::size_t min_region_size = 50; const std::size_t min_region_size = 50;
// No constraint on radius
const double min_radius = 0.;
const double max_radius = std::numeric_limits<double>::infinity();
Neighbor_query neighbor_query( Neighbor_query neighbor_query(
points, CGAL::parameters:: points, CGAL::parameters::
k_neighbors(k). k_neighbors(k).
point_map(points.point_map())); point_map(points.point_map()));
Region_type region_type(points, tolerance, max_angle, min_region_size, Region_type region_type(
min_radius, max_radius, points, CGAL::parameters::
points.point_map(), points.normal_map()); maximum_distance(max_distance).
Region_growing region_growing(points, neighbor_query, region_type); maximum_angle(max_angle).
minimum_region_size(min_region_size).
point_map(points.point_map()).
normal_map(points.normal_map()));
Region_growing region_growing(
points, neighbor_query, region_type);
// Add maps to get colored output // Add maps to get colored output
Point_set::Property_map<unsigned char> Point_set::Property_map<unsigned char>

View File

@ -83,9 +83,6 @@ namespace Point_set {
using Sqrt = typename Get_sqrt::Sqrt; using Sqrt = typename Get_sqrt::Sqrt;
public: public:
/// @}
/// \name Initialization /// \name Initialization
/// @{ /// @{
@ -202,7 +199,8 @@ namespace Point_set {
CGAL_precondition(m_min_radius >= FT(0)); CGAL_precondition(m_min_radius >= FT(0));
const FT m_max_radius = parameters::choose_parameter( const FT m_max_radius = parameters::choose_parameter(
parameters::get_parameter(np, internal_np::maximum_radius), FT(std::numeric_limits<double>::max())); parameters::get_parameter(np, internal_np::maximum_radius),
FT(std::numeric_limits<double>::max()));
CGAL_precondition(m_max_radius >= m_min_radius); CGAL_precondition(m_max_radius >= m_min_radius);
} }
@ -288,7 +286,8 @@ namespace Point_set {
\param query_index \param query_index
index of the query point index of the query point
\param indices indices of the inliers of the region \param indices
indices of the inliers of the region
The first parameter is not used in this implementation. The first parameter is not used in this implementation.

View File

@ -16,299 +16,407 @@
#include <CGAL/license/Shape_detection.h> #include <CGAL/license/Shape_detection.h>
#include <cmath> // Internal includes.
#include <vector>
#include <CGAL/assertions.h>
#include <CGAL/number_utils.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/Shape_detection/Region_growing/internal/utils.h> #include <CGAL/Shape_detection/Region_growing/internal/utils.h>
namespace CGAL { namespace CGAL {
namespace Shape_detection { namespace Shape_detection {
namespace Point_set { namespace Point_set {
/*!
\ingroup PkgShapeDetectionRGOnPoints
\brief Region type based on the quality of the least squares sphere
fit applied to 3D points.
This class fits a sphere to chunks of points in a 3D 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_sphere_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_3 = typename Traits::Point_3;
using Vector_3 = typename Traits::Vector_3;
using Plane_3 = typename Traits::Plane_3;
using Squared_length_3 = typename Traits::Compute_squared_length_3;
using Squared_distance_3 = typename Traits::Compute_squared_distance_3;
using Scalar_product_3 = typename Traits::Compute_scalar_product_3;
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 FT m_min_radius;
const FT m_max_radius;
const Point_map m_point_map;
const Normal_map m_normal_map;
const Squared_length_3 m_squared_length_3;
const Squared_distance_3 m_squared_distance_3;
const Scalar_product_3 m_scalar_product_3;
const Sqrt m_sqrt;
Point_3 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 3D points and \brief Region type based on the quality of the least squares sphere
corresponding 3D normal vectors fit applied to 3D points.
\param distance_threshold the maximum distance from a point to a This class fits a sphere to chunks of points in a 3D point set and
sphere. %Default is 1. 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 \tparam GeomTraits
between the normal of a point and the radius of a sphere. %Default a model of `Kernel`
is 25 degrees.
\param min_region_size the minimum number of 3D points a region \tparam InputRange
must have. %Default is 3. a model of `ConstRange` whose iterator type is `RandomAccessIterator`
\param minimum_radius the radius below which an estimated sphere \tparam PointMap
is considered as invalid and discarded. %Default is 0 (no limit). a model of `ReadablePropertyMap` whose key type is the value type of the input
range and value type is `Kernel::Point_3`
\param maximum_radius the radius above which an estimated sphere \tparam NormalMap
is considered as invalid and discarded. %Default is infinity (no a model of `ReadablePropertyMap` whose key type is the value type of the input
limit). range and value type is `Kernel::Vector_3`
\param point_map an instance of `PointMap` that maps an item from \cgalModels `RegionType`
`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_sphere_fit_region (const InputRange& input_range, template<
const FT distance_threshold = FT(1), typename GeomTraits,
const FT angle_threshold = FT(25), typename InputRange,
const std::size_t min_region_size = 3, typename PointMap,
const FT minimum_radius = FT(0), typename NormalMap>
const FT maximum_radius = std::numeric_limits<FT>::infinity(), class Least_squares_sphere_fit_region {
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_min_radius (minimum_radius)
, m_max_radius (maximum_radius)
, m_point_map(point_map)
, m_normal_map(normal_map)
, m_squared_length_3(traits.compute_squared_length_3_object())
, m_squared_distance_3(traits.compute_squared_distance_3_object())
, m_scalar_product_3(traits.compute_scalar_product_3_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);
}
/// @} public:
/// \name Types
/// @{
/// \name Access /// \cond SKIP_IN_MANUAL
/// @{ using Traits = GeomTraits;
using Input_range = InputRange;
using Point_map = PointMap;
using Normal_map = NormalMap;
/// \endcond
/*! /// Number type.
\brief implements `RegionType::is_part_of_region()`. typedef typename GeomTraits::FT FT;
This function controls if a point with the index `query_index` is /// @}
within the `distance_threshold` from the corresponding sphere and
if the angle between its normal and the sphere radius is within
the `angle_threshold`. If both conditions are satisfied, it
returns `true`, otherwise `false`.
\param query_index index of the query point private:
using Point_3 = typename Traits::Point_3;
using Vector_3 = typename Traits::Vector_3;
using Plane_3 = typename Traits::Plane_3;
\param indices indices of the inliers of the region using Squared_length_3 = typename Traits::Compute_squared_length_3;
using Squared_distance_3 = typename Traits::Compute_squared_distance_3;
using Scalar_product_3 = typename Traits::Compute_scalar_product_3;
The first parameter is not used in this implementation. using Get_sqrt = internal::Get_sqrt<Traits>;
using Sqrt = typename Get_sqrt::Sqrt;
\return Boolean `true` or `false` public:
/// \name Initialization
/// @{
\pre `query_index >= 0 && query_index < input_range.size()` /*!
*/ \brief initializes all internal data structures.
bool is_part_of_region (const std::size_t,
const std::size_t query_index, \tparam NamedParameters
const std::vector<std::size_t>& indices) const a sequence of \ref bgl_namedparameters "Named Parameters"
{
CGAL_precondition(query_index < m_input_range.size()); \param input_range
an instance of `InputRange` with 3D points and
corresponding 3D 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 sphere}
\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 sphere}
\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 3D points a region must have}
\cgalParamType{`std::size_t`}
\cgalParamDefault{3}
\cgalParamNEnd
\cgalParamNBegin{minimum_radius}
\cgalParamDescription{the radius below which an estimated sphere
is considered as invalid and discarded}
\cgalParamType{`GeomTraits::FT`}
\cgalParamDefault{0, no limit}
\cgalParamNEnd
\cgalParamNBegin{maximum_radius}
\cgalParamDescription{the radius above which an estimated sphere
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_3`}
\cgalParamDefault{`PointMap()`}
\cgalParamNEnd
\cgalParamNBegin{normal_map}
\cgalParamDescription{ an instance of `NormalMap` that maps an item from `input_range`
to `Kernel::Vector_3`}
\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<typename NamedParameters>
Least_squares_sphere_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_length_3(m_traits.compute_squared_length_3_object()),
m_squared_distance_3(m_traits.compute_squared_distance_3_object()),
m_scalar_product_3(m_traits.compute_scalar_product_3_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<FT>(std::cos(CGAL::to_double(
(max_angle * static_cast<FT>(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<double>::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 3D points and
corresponding 3D normal vectors
\param distance_threshold
the maximum distance from a point to a sphere. %Default is 1.
\param angle_threshold
the maximum accepted angle in degrees between the normal of a point
and the radius of a sphere. %Default is 25 degrees.
\param min_region_size
the minimum number of 3D points a region must have. %Default is 3.
\param minimum_radius
the radius below which an estimated sphere is considered as invalid
and discarded. %Default is 0 (no limit).
\param maximum_radius
the radius above which an estimated sphere 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_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`
\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_sphere_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<double>::max()),
const PointMap point_map = PointMap(),
const NormalMap normal_map = NormalMap(),
const GeomTraits traits = GeomTraits()) :
Least_squares_sphere_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 sphere and if the angle between
its normal and the sphere 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<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 sphere 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_3& query_point = get(m_point_map, key);
Vector_3 normal = get(m_normal_map, key);
const FT sq_dist = m_squared_distance_3(query_point, m_center);
if (std::isnan(sq_dist)) return false;
FT distance_to_center = m_sqrt (sq_dist);
FT distance_to_sphere = CGAL::abs (distance_to_center - m_radius);
if (distance_to_sphere > 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_3 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 sphere means something
if (indices.size() < 6)
return true; return true;
}
// TODO: Why do we get so many nan in this class? /*!
if (std::isnan(m_radius)) \brief implements `RegionType::is_valid_region()`.
return false;
// If radius is out of bound, nothing fits, early ending This function controls if the estimated radius is between `minimum_radius`
if (m_radius < m_min_radius || m_radius > m_max_radius) and `maximum_radius` and if the `region` contains at least `min_region_size` points.
return false;
const auto& key = *(m_input_range.begin() + query_index); \param region
const Point_3& query_point = get(m_point_map, key); indices of points included in the region
Vector_3 normal = get(m_normal_map, key);
const FT sq_dist = m_squared_distance_3(query_point, m_center); \return Boolean `true` or `false`
if (std::isnan(sq_dist)) return false; */
FT distance_to_center = m_sqrt (sq_dist); inline bool is_valid_region(const std::vector<std::size_t>& region) const {
FT distance_to_sphere = CGAL::abs (distance_to_center - m_radius); return (
(m_min_radius <= m_radius && m_radius <= m_max_radius) &&
(region.size() >= m_min_region_size)
);
}
if (distance_to_sphere > m_distance_threshold) /*!
return false; \brief implements `RegionType::update()`.
const FT sq_norm = normal * normal; This function fits the least squares sphere to all points from the `region`.
if (std::isnan(sq_norm)) return false;
normal = normal / m_sqrt (sq_norm);
Vector_3 ray (m_center, query_point); \param region
const FT sq_ray = ray * ray; indices of points included in the region
if (std::isnan(sq_ray)) return false;
ray = ray / m_sqrt (sq_ray);
if (CGAL::abs (normal * ray) < m_normal_threshold) \return Boolean `true` if the sphere fitting succeeded and `false` otherwise
return false;
return true; \pre `region.size() > 0`
} */
bool update(const std::vector<std::size_t>& region) {
/*! CGAL_precondition(region.size() > 0);
\brief implements `RegionType::is_valid_region()`. auto unary_function
= [&](const std::size_t& idx) -> const Point_3&
{
return get (m_point_map, *(m_input_range.begin() + idx));
};
This function controls if the estimated radius is between internal::sphere_fit
`minimum_radius` and `maximum_radius` and if the `region` contains (make_range (boost::make_transform_iterator
at least `min_region_size` points. (region.begin(), unary_function),
boost::make_transform_iterator
(region.end(), unary_function)),
m_sqrt, m_squared_distance_3, m_center, m_radius);
return true;
}
\param region /// @}
indices of points included in the region
\return Boolean `true` or `false` private:
*/ const Input_range& m_input_range;
inline bool is_valid_region(const std::vector<std::size_t>& region) const const Point_map m_point_map;
{ const Normal_map m_normal_map;
return ((m_min_radius <= m_radius && m_radius <= m_max_radius) const Traits m_traits;
&& (region.size() >= m_min_region_size));
}
/*! FT m_distance_threshold;
\brief implements `RegionType::update()`. 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 sphere to all points from the `region`. const Squared_length_3 m_squared_length_3;
const Squared_distance_3 m_squared_distance_3;
const Scalar_product_3 m_scalar_product_3;
const Sqrt m_sqrt;
\param region Point_3 m_center;
indices of points included in the region FT m_radius;
};
\pre `region.size() > 0`
*/
bool update(const std::vector<std::size_t>& region)
{
CGAL_precondition(region.size() > 0);
auto unary_function
= [&](const std::size_t& idx) -> const Point_3&
{
return get (m_point_map, *(m_input_range.begin() + idx));
};
internal::sphere_fit
(make_range (boost::make_transform_iterator
(region.begin(), unary_function),
boost::make_transform_iterator
(region.end(), unary_function)),
m_sqrt, m_squared_distance_3, m_center, m_radius);
return true;
}
/// @}
};
} // namespace Point_set } // namespace Point_set
} // namespace Shape_detection } // namespace Shape_detection

View File

@ -8,7 +8,7 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
// //
// //
// Author(s) : Simon Giraudot, Dmitry Anisimov // Author(s) : Simon Giraudot
// //
#ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_SPHERE_FIT_SORTING_H #ifndef CGAL_SHAPE_DETECTION_REGION_GROWING_POINT_SET_LEAST_SQUARES_SPHERE_FIT_SORTING_H
@ -16,202 +16,235 @@
#include <CGAL/license/Shape_detection.h> #include <CGAL/license/Shape_detection.h>
// STL includes.
#include <vector>
#include <algorithm>
// CGAL includes.
#include <CGAL/assertions.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
// Internal includes. // Internal includes.
#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
#include <CGAL/Shape_detection/Region_growing/internal/property_map.h> #include <CGAL/Shape_detection/Region_growing/internal/property_map.h>
namespace CGAL { namespace CGAL {
namespace Shape_detection { namespace Shape_detection {
namespace Point_set { namespace Point_set {
/*!
\ingroup PkgShapeDetectionRGOnPoints
\brief Sorting of 3D points with respect to the local sphere fit quality.
Indices of 3D input points are sorted with respect to the quality of the
least squares sphere fit applied to the neighboring points of each point.
\tparam GeomTraits
must be a model of `Kernel`.
\tparam InputRange
must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`.
\tparam NeighborQuery
must be 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_3`.
*/
template<typename GeomTraits,
typename InputRange,
typename NeighborQuery,
typename PointMap>
class Least_squares_sphere_fit_sorting
{
public:
/// \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
#ifdef DOXYGEN_RUNNING
/*! /*!
an `LvaluePropertyMap` whose key and value type is `std::size_t`. \ingroup PkgShapeDetectionRGOnPoints
This map provides an access to the ordered indices of input points.
\brief Sorting of 3D points with respect to the local sphere fit quality.
Indices of 3D input points are sorted with respect to the quality of the
least squares sphere fit applied to the neighboring points of each point.
\tparam GeomTraits
a model of `Kernel`
\tparam InputRange
a model of `ConstRange` whose iterator type is `RandomAccessIterator`
\tparam NeighborQuery
a model of `NeighborQuery`
\tparam PointMap
a model of `ReadablePropertyMap` whose key type is the value type of the input
range and value type is `Kernel::Point_3`
*/ */
typedef unspecified_type Seed_map; template<
#endif typename GeomTraits,
typename InputRange,
typename NeighborQuery,
typename PointMap>
class Least_squares_sphere_fit_sorting {
/// @} public:
/// \name Initialization /// \name Types
/// @{ /// @{
/*! /// \cond SKIP_IN_MANUAL
\brief initializes all internal data structures. using Traits = GeomTraits;
using Input_range = InputRange;
using Neighbor_query = NeighborQuery;
using Point_map = PointMap;
using Seed_map = internal::Seed_property_map;
/// \endcond
\param input_range #ifdef DOXYGEN_RUNNING
an instance of `InputRange` with 3D points /*!
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
\param neighbor_query /// @}
an instance of `NeighborQuery` that is used internally to
access point's neighbors
\param point_map private:
an instance of `PointMap` that maps an item from `input_range` // using FT = typename Traits::FT;
to `Kernel::Point_3` // using Compare_scores = internal::Compare_scores<FT>;
\pre `input_range.size() > 0` using Local_traits = Exact_predicates_inexact_constructions_kernel;
*/ using Local_FT = typename Local_traits::FT;
Least_squares_sphere_fit_sorting (const InputRange& input_range, using Local_point_3 = typename Local_traits::Point_3;
NeighborQuery& neighbor_query, using To_local_converter = Cartesian_converter<Traits, Local_traits>;
const PointMap point_map = PointMap()) using Compare_scores = internal::Compare_scores<Local_FT>;
: m_input_range(input_range)
, m_neighbor_query(neighbor_query)
, m_point_map(point_map)
, m_to_local_converter()
{
CGAL_precondition(input_range.size() > 0);
m_order.resize(m_input_range.size()); public:
for (std::size_t i = 0; i < m_input_range.size(); ++i) /// \name Initialization
m_order[i] = i; /// @{
m_scores.resize(m_input_range.size());
}
/// @} /*!
\brief initializes all internal data structures.
/// \name Sorting \tparam NamedParameters
/// @{ a sequence of \ref bgl_namedparameters "Named Parameters"
/*! \param input_range
\brief sorts indices of input points. an instance of `InputRange` with 3D points
*/
void sort()
{
compute_scores();
CGAL_postcondition(m_scores.size() > 0);
Compare_scores cmp(m_scores); \param neighbor_query
std::sort(m_order.begin(), m_order.end(), cmp); 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
/// \name Access \cgalNamedParamsBegin
/// @{ \cgalParamNBegin{point_map}
\cgalParamDescription{an instance of `PointMap` that maps an item from `input_range`
to `Kernel::Point_3`}
\cgalParamDefault{`PointMap()`}
\cgalParamNEnd
\cgalParamNBegin{geom_traits}
\cgalParamDescription{an instance of `GeomTraits`}
\cgalParamDefault{`GeomTraits()`}
\cgalParamNEnd
\cgalNamedParamsEnd
/*! \pre `input_range.size() > 0`
\brief returns an instance of `Seed_map` to access the ordered indices */
of input points. template<typename NamedParameters>
*/ Least_squares_sphere_fit_sorting(
Seed_map seed_map() const InputRange& input_range,
{ NeighborQuery& neighbor_query,
return Seed_map(m_order); const NamedParameters& np) :
} m_input_range(input_range),
m_neighbor_query(neighbor_query),
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());
private: std::iota(m_order.begin(), m_order.end(), 0);
m_scores.resize(m_input_range.size());
// Types.
using Local_traits = Exact_predicates_inexact_constructions_kernel;
using Local_FT = typename Local_traits::FT;
using Local_point_3 = typename Local_traits::Point_3;
using To_local_converter = Cartesian_converter<Traits, Local_traits>;
using Compare_scores = internal::Compare_scores<Local_FT>;
// Functions.
void compute_scores()
{
std::vector<std::size_t> neighbors;
std::vector<Local_point_3> points;
typename internal::Get_sqrt<Local_traits>::Sqrt sqrt;
typename Local_traits::Compute_squared_distance_3 squared_distance_3;
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_3 fitted_center;
Local_FT fitted_radius;
if (internal::sphere_fit (points, sqrt, squared_distance_3, fitted_center, fitted_radius))
{
// Score is min squared distance to sphere
m_scores[i] = Local_FT(0);
for (const Local_point_3& p : points)
m_scores[i] += abs (sqrt(squared_distance_3(p, fitted_center)) - fitted_radius);
}
else
m_scores[i] = Local_FT(std::numeric_limits<double>::infinity());
} }
}
// Fields. /*!
const Input_range& m_input_range; \brief initializes all internal data structures.
Neighbor_query& m_neighbor_query;
const Point_map m_point_map;
std::vector<std::size_t> m_order; \deprecated This constructor is deprecated since the version 5.5 of \cgal.
std::vector<Local_FT> m_scores;
const To_local_converter m_to_local_converter; \param input_range
}; an instance of `InputRange` with 3D 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_3`
\pre `input_range.size() > 0`
*/
CGAL_DEPRECATED_MSG("This constructor is deprecated since the version 5.5 of CGAL!")
Least_squares_sphere_fit_sorting(
const InputRange& input_range,
NeighborQuery& neighbor_query,
const PointMap point_map = PointMap()) :
Least_squares_sphere_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_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:
const Input_range& m_input_range;
Neighbor_query& m_neighbor_query;
const Point_map m_point_map;
const Traits m_traits;
std::vector<std::size_t> m_order;
std::vector<Local_FT> m_scores;
const To_local_converter m_to_local_converter;
void compute_scores() {
std::vector<std::size_t> neighbors;
std::vector<Local_point_3> points;
typename internal::Get_sqrt<Local_traits>::Sqrt sqrt;
typename Local_traits::Compute_squared_distance_3 squared_distance_3;
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_3 fitted_center;
Local_FT fitted_radius;
if (internal::sphere_fit (points, sqrt, squared_distance_3, fitted_center, fitted_radius))
{
// Score is min squared distance to sphere
m_scores[i] = Local_FT(0);
for (const Local_point_3& p : points)
m_scores[i] += abs (sqrt(squared_distance_3(p, fitted_center)) - fitted_radius);
}
else
m_scores[i] = Local_FT(std::numeric_limits<double>::infinity());
}
}
};
} // namespace Point_set } // namespace Point_set
} // namespace Shape_detection } // namespace Shape_detection