mirror of https://github.com/CGAL/cgal
removed fitting.h and using np in circle classes
This commit is contained in:
Dmitry Anisimov
parent
b28f32269b
commit
0fbc3ea1bb
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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<double>::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<unsigned char>
|
||||
|
|
|
|||
|
|
@ -16,294 +16,401 @@
|
|||
|
||||
#include <CGAL/license/Shape_detection.h>
|
||||
|
||||
#include <cmath>
|
||||
#include <vector>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
#include <CGAL/number_utils.h>
|
||||
#include <CGAL/Cartesian_converter.h>
|
||||
|
||||
// Internal includes.
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/fitting.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Shape_detection {
|
||||
namespace Point_set {
|
||||
|
||||
/*!
|
||||
\ingroup PkgShapeDetectionRGOnPoints
|
||||
|
||||
\brief Region type based on the quality of the least squares circle
|
||||
fit applied to 2D points.
|
||||
|
||||
This class fits a circle to chunks of points in a 2D point set and
|
||||
controls the quality of this fit. If all quality conditions are
|
||||
satisfied, the chunk is accepted as a valid region, otherwise
|
||||
rejected.
|
||||
|
||||
\tparam GeomTraits
|
||||
must be a model of `Kernel`.
|
||||
|
||||
\tparam InputRange
|
||||
must be a model of `ConstRange` whose iterator type is `RandomAccessIterator`.
|
||||
|
||||
\tparam PointMap
|
||||
must be an `LvaluePropertyMap` whose key type is the value type of the input
|
||||
range and value type is `Kernel::Point_3`.
|
||||
|
||||
\tparam NormalMap
|
||||
must be an `LvaluePropertyMap` whose key type is the value type of the input
|
||||
range and value type is `Kernel::Vector_3`.
|
||||
|
||||
\cgalModels `RegionType`
|
||||
*/
|
||||
template<typename GeomTraits,
|
||||
typename InputRange,
|
||||
typename PointMap,
|
||||
typename NormalMap>
|
||||
class Least_squares_circle_fit_region
|
||||
{
|
||||
|
||||
public:
|
||||
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
using Traits = GeomTraits;
|
||||
using Input_range = InputRange;
|
||||
using Point_map = PointMap;
|
||||
using Normal_map = NormalMap;
|
||||
/// \endcond
|
||||
|
||||
/// Number type.
|
||||
typedef typename GeomTraits::FT FT;
|
||||
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
using Point_2 = typename Traits::Point_2;
|
||||
using Vector_2 = typename Traits::Vector_2;
|
||||
|
||||
using Squared_distance_2 = typename Traits::Compute_squared_distance_2;
|
||||
|
||||
using Get_sqrt = internal::Get_sqrt<Traits>;
|
||||
using Sqrt = typename Get_sqrt::Sqrt;
|
||||
|
||||
private:
|
||||
|
||||
const Input_range& m_input_range;
|
||||
|
||||
const FT m_distance_threshold;
|
||||
const FT m_normal_threshold;
|
||||
const std::size_t m_min_region_size;
|
||||
const 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<FT>::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<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_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<Traits>;
|
||||
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<std::size_t>& 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<typename NamedParameters>
|
||||
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<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 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<double>::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<std::size_t>& indices) const {
|
||||
|
||||
CGAL_precondition(query_index < m_input_range.size());
|
||||
|
||||
// First, we need to integrate at least 6 points so that the
|
||||
// computed circle means something
|
||||
if (indices.size() < 6)
|
||||
return true;
|
||||
|
||||
// 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<std::size_t>& 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<std::size_t>& 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<std::size_t>& 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<std::size_t>& region)
|
||||
{
|
||||
CGAL_precondition(region.size() > 0);
|
||||
|
||||
auto unary_function
|
||||
= [&](const std::size_t& idx) -> const Point_2&
|
||||
{
|
||||
return get (m_point_map, *(m_input_range.begin() + idx));
|
||||
};
|
||||
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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 <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.
|
||||
#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/fitting.h>
|
||||
|
||||
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<FT>;
|
||||
|
||||
/// \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<Traits, Local_traits>;
|
||||
using Compare_scores = internal::Compare_scores<Local_FT>;
|
||||
|
||||
/*!
|
||||
\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<typename NamedParameters>
|
||||
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<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_2> points;
|
||||
|
||||
typename internal::Get_sqrt<Local_traits>::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<double>::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<std::size_t> m_order;
|
||||
std::vector<Local_FT> 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<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_2> points;
|
||||
|
||||
typename internal::Get_sqrt<Local_traits>::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<double>::infinity());
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Point_set
|
||||
} // namespace Shape_detection
|
||||
|
|
|
|||
|
|
@ -24,7 +24,6 @@
|
|||
#include <CGAL/Cartesian_converter.h>
|
||||
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/fitting.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Shape_detection {
|
||||
|
|
|
|||
|
|
@ -28,7 +28,6 @@
|
|||
// 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/fitting.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Shape_detection {
|
||||
|
|
|
|||
|
|
@ -24,7 +24,6 @@
|
|||
#include <CGAL/Cartesian_converter.h>
|
||||
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/fitting.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Shape_detection {
|
||||
|
|
|
|||
|
|
@ -28,7 +28,6 @@
|
|||
// 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/fitting.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Shape_detection {
|
||||
|
|
|
|||
|
|
@ -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 <CGAL/license/Shape_detection.h>
|
||||
|
||||
#include <CGAL/Eigen_diagonalize_traits.h>
|
||||
#include <CGAL/property_map.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Shape_detection {
|
||||
namespace internal {
|
||||
|
||||
template <typename PointRange,
|
||||
typename Sqrt, typename Squared_distance_2,
|
||||
typename Point_2, typename FT>
|
||||
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<FT, 4>;
|
||||
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<FT>(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 <typename PointRange,
|
||||
typename Sqrt, typename Squared_distance_3,
|
||||
typename Point_3, typename FT>
|
||||
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<FT, 5>;
|
||||
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<FT>(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 <typename PointRange, typename PointMap, typename NormalMap,
|
||||
typename Sqrt, typename Squared_distance_3,
|
||||
typename Line_3, typename FT>
|
||||
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<PointMap>::value_type;
|
||||
using Vector_3 = typename boost::property_traits<NormalMap>::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<FT>(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
|
||||
|
|
@ -16,15 +16,6 @@
|
|||
|
||||
#include <CGAL/license/Shape_detection.h>
|
||||
|
||||
// STL includes.
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
#include <type_traits>
|
||||
|
||||
// CGAL includes.
|
||||
#include <CGAL/assertions.h>
|
||||
#include <CGAL/number_utils.h>
|
||||
|
||||
// Internal includes.
|
||||
#include <CGAL/Shape_detection/Region_growing/internal/utils.h>
|
||||
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@
|
|||
// STL includes.
|
||||
#include <set>
|
||||
#include <map>
|
||||
#include <cmath>
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
#include <type_traits>
|
||||
|
|
@ -345,6 +346,240 @@ namespace internal {
|
|||
return std::make_pair(plane, static_cast<FT>(score));
|
||||
}
|
||||
|
||||
template <typename PointRange,
|
||||
typename Sqrt, typename Squared_distance_2,
|
||||
typename Point_2, typename FT>
|
||||
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<FT, 4>;
|
||||
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<FT>(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 <typename PointRange,
|
||||
typename Sqrt, typename Squared_distance_3,
|
||||
typename Point_3, typename FT>
|
||||
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<FT, 5>;
|
||||
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<FT>(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 <typename PointRange, typename PointMap, typename NormalMap,
|
||||
typename Sqrt, typename Squared_distance_3,
|
||||
typename Line_3, typename FT>
|
||||
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<PointMap>::value_type;
|
||||
using Vector_3 = typename boost::property_traits<NormalMap>::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<FT>(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
|
||||
|
|
|
|||
Loading…
Reference in New Issue