mirror of https://github.com/CGAL/cgal
Make PCA compatible with Epeck
This commit is contained in:
parent
e58561aeca
commit
35a096011e
|
|
@ -22,6 +22,12 @@ namespace CGAL {
|
||||||
|
|
||||||
namespace internal {
|
namespace internal {
|
||||||
|
|
||||||
|
template <typename FT>
|
||||||
|
FT approximate_cbrt (const FT& x)
|
||||||
|
{
|
||||||
|
return static_cast<FT>(std::cbrt (CGAL::to_double(x)));
|
||||||
|
}
|
||||||
|
|
||||||
// assemble covariance matrix from a triangle set
|
// assemble covariance matrix from a triangle set
|
||||||
template < typename InputIterator,
|
template < typename InputIterator,
|
||||||
typename K >
|
typename K >
|
||||||
|
|
@ -67,7 +73,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
||||||
t[0].y(), t[1].y(), t[2].y(),
|
t[0].y(), t[1].y(), t[2].y(),
|
||||||
t[0].z(), t[1].z(), t[2].z();
|
t[0].z(), t[1].z(), t[2].z();
|
||||||
|
|
||||||
FT area = std::sqrt(t.squared_area());
|
FT area = CGAL::approximate_sqrt(t.squared_area());
|
||||||
|
|
||||||
// skip zero measure primitives
|
// skip zero measure primitives
|
||||||
if(area == (FT)0.0)
|
if(area == (FT)0.0)
|
||||||
|
|
@ -236,15 +242,15 @@ assemble_covariance_matrix_3(InputIterator first,
|
||||||
t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
|
t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
|
||||||
t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
|
t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
|
||||||
Matrix transformation (delta);
|
Matrix transformation (delta);
|
||||||
FT area = std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
|
FT area = approximate_cbrt(delta[0]*delta[0] + delta[3]*delta[3] +
|
||||||
delta[6]*delta[6],1/3.0)*std::pow(delta[1]*delta[1] +
|
delta[6]*delta[6])*approximate_cbrt(delta[1]*delta[1] +
|
||||||
delta[4]*delta[4] + delta[7]*delta[7],1/3.0)*2 +
|
delta[4]*delta[4] + delta[7]*delta[7])*2 +
|
||||||
std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
|
approximate_cbrt(delta[0]*delta[0] + delta[3]*delta[3] +
|
||||||
delta[6]*delta[6],1/3.0)*std::pow(delta[2]*delta[2] +
|
delta[6]*delta[6])*approximate_cbrt(delta[2]*delta[2] +
|
||||||
delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2 +
|
delta[5]*delta[5] + delta[8]*delta[8])*2 +
|
||||||
std::pow(delta[1]*delta[1] + delta[4]*delta[4] +
|
approximate_cbrt(delta[1]*delta[1] + delta[4]*delta[4] +
|
||||||
delta[7]*delta[7],1/3.0)*std::pow(delta[2]*delta[2] +
|
delta[7]*delta[7])*approximate_cbrt(delta[2]*delta[2] +
|
||||||
delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2;
|
delta[5]*delta[5] + delta[8]*delta[8])*2;
|
||||||
|
|
||||||
// skip zero measure primitives
|
// skip zero measure primitives
|
||||||
if(area == (FT)0.0)
|
if(area == (FT)0.0)
|
||||||
|
|
@ -324,7 +330,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
||||||
const Sphere& t = *it;
|
const Sphere& t = *it;
|
||||||
|
|
||||||
// defined for convenience.
|
// defined for convenience.
|
||||||
FT radius = std::sqrt(t.squared_radius());
|
FT radius = CGAL::approximate_sqrt(t.squared_radius());
|
||||||
Matrix transformation;
|
Matrix transformation;
|
||||||
transformation << radius, 0.0, 0.0,
|
transformation << radius, 0.0, 0.0,
|
||||||
0.0, radius, 0.0,
|
0.0, radius, 0.0,
|
||||||
|
|
@ -409,7 +415,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
||||||
|
|
||||||
// defined for convenience.
|
// defined for convenience.
|
||||||
// FT example = CGAL::to_double(t[0].x());
|
// FT example = CGAL::to_double(t[0].x());
|
||||||
FT radius = std::sqrt(t.squared_radius());
|
FT radius = CGAL::approximate_sqrt(t.squared_radius());
|
||||||
Matrix transformation;
|
Matrix transformation;
|
||||||
transformation << radius, 0.0, 0.0,
|
transformation << radius, 0.0, 0.0,
|
||||||
0.0, radius, 0.0,
|
0.0, radius, 0.0,
|
||||||
|
|
@ -576,8 +582,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
||||||
transformation << t[0].x(), t[1].x(), 0.0,
|
transformation << t[0].x(), t[1].x(), 0.0,
|
||||||
t[0].y(), t[1].y(), 0.0,
|
t[0].y(), t[1].y(), 0.0,
|
||||||
t[0].z(), t[1].z(), 1.0;
|
t[0].z(), t[1].z(), 1.0;
|
||||||
using std::sqrt;
|
FT length = CGAL::approximate_sqrt(t.squared_length());
|
||||||
FT length = sqrt(t.squared_length());
|
|
||||||
|
|
||||||
// skip zero measure primitives
|
// skip zero measure primitives
|
||||||
if(length == (FT)0.0)
|
if(length == (FT)0.0)
|
||||||
|
|
|
||||||
|
|
@ -83,7 +83,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
||||||
|
|
||||||
// defined for convenience.
|
// defined for convenience.
|
||||||
// FT example = CGAL::to_double(t[0].x());
|
// FT example = CGAL::to_double(t[0].x());
|
||||||
FT radius = std::sqrt(t.squared_radius());
|
FT radius = CGAL::approximate_sqrt(t.squared_radius());
|
||||||
FT delta[4] = {radius, 0.0,
|
FT delta[4] = {radius, 0.0,
|
||||||
0.0, radius};
|
0.0, radius};
|
||||||
Matrix transformation = init_matrix<FT>(2,delta);
|
Matrix transformation = init_matrix<FT>(2,delta);
|
||||||
|
|
@ -189,7 +189,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
||||||
|
|
||||||
// defined for convenience.
|
// defined for convenience.
|
||||||
// FT example = CGAL::to_double(t[0].x());
|
// FT example = CGAL::to_double(t[0].x());
|
||||||
FT radius = std::sqrt(t.squared_radius());
|
FT radius = CGAL::approximate_sqrt(t.squared_radius());
|
||||||
FT delta[4] = {radius, 0.0,
|
FT delta[4] = {radius, 0.0,
|
||||||
0.0, radius};
|
0.0, radius};
|
||||||
Matrix transformation = init_matrix<FT>(2,delta);
|
Matrix transformation = init_matrix<FT>(2,delta);
|
||||||
|
|
@ -199,7 +199,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
||||||
// Find the 2nd order moment for the circle wrt to the origin by an affine transformation.
|
// Find the 2nd order moment for the circle wrt to the origin by an affine transformation.
|
||||||
|
|
||||||
// Transform the standard 2nd order moment using the transformation matrix
|
// Transform the standard 2nd order moment using the transformation matrix
|
||||||
transformation = 0.5 * length * transformation * moment * LA::transpose(transformation);
|
transformation = FT(0.5) * length * transformation * moment * LA::transpose(transformation);
|
||||||
|
|
||||||
// Translate the 2nd order moment to the center of the circle.
|
// Translate the 2nd order moment to the center of the circle.
|
||||||
FT x0 = t.center().x();
|
FT x0 = t.center().x();
|
||||||
|
|
|
||||||
|
|
@ -83,7 +83,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
||||||
t[0].y(), t[1].y()};
|
t[0].y(), t[1].y()};
|
||||||
Matrix transformation = init_matrix<FT>(2,delta);
|
Matrix transformation = init_matrix<FT>(2,delta);
|
||||||
using std::sqrt;
|
using std::sqrt;
|
||||||
FT length = sqrt(t.squared_length());
|
FT length = CGAL::approximate_sqrt(t.squared_length());
|
||||||
CGAL_assertion(length != 0.0);
|
CGAL_assertion(length != 0.0);
|
||||||
|
|
||||||
// Find the 2nd order moment for the segment wrt to the origin by an affine transformation.
|
// Find the 2nd order moment for the segment wrt to the origin by an affine transformation.
|
||||||
|
|
|
||||||
|
|
@ -88,7 +88,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
||||||
t[1].y() - y0, t[2].y() - y0};
|
t[1].y() - y0, t[2].y() - y0};
|
||||||
|
|
||||||
Matrix transformation = init_matrix<FT>(2,delta);
|
Matrix transformation = init_matrix<FT>(2,delta);
|
||||||
FT area = 0.5 * std::abs(LA::determinant(transformation));
|
FT area = 0.5 * CGAL::abs(LA::determinant(transformation));
|
||||||
CGAL_assertion(area!=0);
|
CGAL_assertion(area!=0);
|
||||||
|
|
||||||
// Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
|
// Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
|
||||||
|
|
|
||||||
|
|
@ -77,7 +77,7 @@ linear_least_squares_fitting_3(InputIterator first,
|
||||||
// precondition: at least one element in the container.
|
// precondition: at least one element in the container.
|
||||||
CGAL_precondition(first != beyond);
|
CGAL_precondition(first != beyond);
|
||||||
|
|
||||||
std::list<Segment> segments;
|
std::list<Segment> segments;
|
||||||
for(InputIterator it = first;
|
for(InputIterator it = first;
|
||||||
it != beyond;
|
it != beyond;
|
||||||
it++)
|
it++)
|
||||||
|
|
@ -88,6 +88,7 @@ linear_least_squares_fitting_3(InputIterator first,
|
||||||
segments.push_back(Segment(t[2],t[0]));
|
segments.push_back(Segment(t[2],t[0]));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
// compute fitting plane
|
// compute fitting plane
|
||||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
|
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
|
||||||
diagonalize_traits);
|
diagonalize_traits);
|
||||||
|
|
|
||||||
|
|
@ -120,9 +120,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Segment& s = *it;
|
const Segment& s = *it;
|
||||||
using std::abs;
|
FT length = CGAL::approximate_sqrt(CGAL::abs(s.squared_length()));
|
||||||
using std::sqrt;
|
|
||||||
FT length = sqrt(abs(s.squared_length()));
|
|
||||||
Point c = K().construct_midpoint_2_object()(s[0],s[1]);
|
Point c = K().construct_midpoint_2_object()(s[0],s[1]);
|
||||||
v = v + length * (c - ORIGIN);
|
v = v + length * (c - ORIGIN);
|
||||||
sum_lengths += length;
|
sum_lengths += length;
|
||||||
|
|
@ -216,7 +214,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Triangle& triangle = *it;
|
const Triangle& triangle = *it;
|
||||||
FT unsigned_area = std::abs(triangle.area());
|
FT unsigned_area = CGAL::abs(triangle.area());
|
||||||
Point c = K().construct_centroid_2_object()(triangle[0],triangle[1],triangle[2]);
|
Point c = K().construct_centroid_2_object()(triangle[0],triangle[1],triangle[2]);
|
||||||
v = v + unsigned_area * (c - ORIGIN);
|
v = v + unsigned_area * (c - ORIGIN);
|
||||||
sum_areas += unsigned_area;
|
sum_areas += unsigned_area;
|
||||||
|
|
@ -252,7 +250,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Circle& s = *it;
|
const Circle& s = *it;
|
||||||
FT radius = std::sqrt(s.squared_radius());
|
FT radius = CGAL::approximate_sqrt(s.squared_radius());
|
||||||
Point c = s.center();
|
Point c = s.center();
|
||||||
v = v + radius * (c - ORIGIN);
|
v = v + radius * (c - ORIGIN);
|
||||||
sum_lengths += radius;
|
sum_lengths += radius;
|
||||||
|
|
@ -381,7 +379,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Iso_rectangle& r = *it;
|
const Iso_rectangle& r = *it;
|
||||||
FT unsigned_area = std::abs(r.area());
|
FT unsigned_area = CGAL::abs(r.area());
|
||||||
Point c = K().construct_centroid_2_object()(r[0],r[1],r[2],r[3]);
|
Point c = K().construct_centroid_2_object()(r[0],r[1],r[2],r[3]);
|
||||||
v = v + unsigned_area * (c - ORIGIN);
|
v = v + unsigned_area * (c - ORIGIN);
|
||||||
sum_areas += unsigned_area;
|
sum_areas += unsigned_area;
|
||||||
|
|
@ -442,8 +440,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Segment& s = *it;
|
const Segment& s = *it;
|
||||||
using std::sqrt;
|
FT length = CGAL::approximate_sqrt(s.squared_length());
|
||||||
FT length = sqrt(s.squared_length());
|
|
||||||
Point c = CGAL::midpoint(s.source(),s.target());
|
Point c = CGAL::midpoint(s.source(),s.target());
|
||||||
// Point c = K().construct_midpoint_3_object()(s[0],s[1]);
|
// Point c = K().construct_midpoint_3_object()(s[0],s[1]);
|
||||||
//Point c = Point((s[0][0] + s[1][0])/2.0, (s[0][1] + s[1][1])/2.0, (s[0][2] + s[1][2])/2.0);
|
//Point c = Point((s[0][0] + s[1][0])/2.0, (s[0][1] + s[1][1])/2.0, (s[0][2] + s[1][2])/2.0);
|
||||||
|
|
@ -539,7 +536,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Triangle& triangle = *it;
|
const Triangle& triangle = *it;
|
||||||
FT unsigned_area = std::sqrt(triangle.squared_area());
|
FT unsigned_area = CGAL::approximate_sqrt(triangle.squared_area());
|
||||||
Point c = K().construct_centroid_3_object()(triangle[0],triangle[1],triangle[2]);
|
Point c = K().construct_centroid_3_object()(triangle[0],triangle[1],triangle[2]);
|
||||||
v = v + unsigned_area * (c - ORIGIN);
|
v = v + unsigned_area * (c - ORIGIN);
|
||||||
sum_areas += unsigned_area;
|
sum_areas += unsigned_area;
|
||||||
|
|
@ -608,7 +605,7 @@ centroid(InputIterator begin,
|
||||||
it++)
|
it++)
|
||||||
{
|
{
|
||||||
const Sphere& sphere = *it;
|
const Sphere& sphere = *it;
|
||||||
FT unsigned_volume = sphere.squared_radius() * std::sqrt(sphere.squared_radius());
|
FT unsigned_volume = sphere.squared_radius() * CGAL::approximate_sqrt(sphere.squared_radius());
|
||||||
Point c = sphere.center();
|
Point c = sphere.center();
|
||||||
v = v + unsigned_volume * (c - ORIGIN);
|
v = v + unsigned_volume * (c - ORIGIN);
|
||||||
sum_volumes += unsigned_volume;
|
sum_volumes += unsigned_volume;
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue