mirror of https://github.com/CGAL/cgal
Merge pull request #5130 from sgiraudot/PCA-Several_enhancements-GF
[Principal Component Analysis] Kernel compatibility, diagonalization precision and overload efficiency
This commit is contained in:
Sebastien Loriot
GitHub
commit
db383e15b5
|
|
@ -17,11 +17,18 @@
|
|||
#include <CGAL/Eigen_diagonalize_traits.h>
|
||||
#include <CGAL/Default_diagonalize_traits.h>
|
||||
#include <CGAL/Dimension.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
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
|
||||
template < typename InputIterator,
|
||||
typename K >
|
||||
|
|
@ -49,9 +56,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
|
||||
// assemble 2nd order moment about the origin.
|
||||
Matrix moment;
|
||||
moment << 1.0/12.0, 1.0/24.0, 1.0/24.0,
|
||||
1.0/24.0, 1.0/12.0, 1.0/24.0,
|
||||
1.0/24.0, 1.0/24.0, 1.0/12.0;
|
||||
moment << FT(1.0/12.0), FT(1.0/24.0), FT(1.0/24.0),
|
||||
FT(1.0/24.0), FT(1.0/12.0), FT(1.0/24.0),
|
||||
FT(1.0/24.0), FT(1.0/24.0), FT(1.0/12.0);
|
||||
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
|
|
@ -67,7 +74,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
t[0].y(), t[1].y(), t[2].y(),
|
||||
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
|
||||
if(area == (FT)0.0)
|
||||
|
|
@ -89,14 +96,16 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
mass += area;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.z() * c.x());
|
||||
covariance[3] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[4] += mass * (-1.0 * c.z() * c.y());
|
||||
covariance[5] += mass * (-1.0 * c.z() * c.z());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.z() * c.x());
|
||||
covariance[3] += -mass * (c.y() * c.y());
|
||||
covariance[4] += -mass * (c.z() * c.y());
|
||||
covariance[5] += -mass * (c.z() * c.z());
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -141,18 +150,25 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
|
||||
// defined for convenience.
|
||||
// FT example = CGAL::to_double(t[0].x());
|
||||
FT x0 = t[0].x();
|
||||
FT y0 = t[0].y();
|
||||
FT z0 = t[0].z();
|
||||
FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0,
|
||||
t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
|
||||
t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
|
||||
Matrix transformation (delta);
|
||||
FT volume = t.volume();
|
||||
FT x0 = t.xmin();
|
||||
FT y0 = t.ymin();
|
||||
FT z0 = t.zmin();
|
||||
|
||||
FT x1 = t.xmax();
|
||||
FT y1 = t.ymax();
|
||||
FT z1 = t.zmax();
|
||||
|
||||
Matrix transformation;
|
||||
transformation << x1 - x0, 0 , 0 ,
|
||||
0 , y1 - y0, 0 ,
|
||||
0 , 0 , z1 - z0;
|
||||
|
||||
FT volume = (x1-x0) * (y1-y0) * (z1-z0);
|
||||
|
||||
// skip zero measure primitives
|
||||
if(volume == (FT)0.0)
|
||||
continue;
|
||||
CGAL_assertion(volume > 0.0);
|
||||
|
||||
// Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
|
||||
|
||||
|
|
@ -160,9 +176,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
transformation = volume * transformation * moment * transformation.transpose();
|
||||
|
||||
// Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
|
||||
FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
|
||||
FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
|
||||
FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
|
||||
FT xav0 = (x1 - x0) / FT(2.0);
|
||||
FT yav0 = (y1 - y0) / FT(2.0);
|
||||
FT zav0 = (z1 - z0) / FT(2.0);
|
||||
|
||||
// and add to covariance matrix
|
||||
covariance[0] += transformation(0,0) + volume * (2*x0*xav0 + x0*x0);
|
||||
|
|
@ -175,6 +191,8 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
mass += volume;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (- c.x() * c.x());
|
||||
|
|
@ -198,6 +216,32 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
const CGAL::Dimension_tag<2>&,
|
||||
const Eigen_diagonalize_traits<typename K::FT, 3>&)
|
||||
{
|
||||
#if 1
|
||||
typedef typename K::FT FT;
|
||||
typedef typename K::Iso_cuboid_3 Iso_cuboid;
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
auto converter = [](const Iso_cuboid& c, int idx) -> Triangle
|
||||
{
|
||||
// Decomposition of 6 faces of the cuboid into 12 triangles
|
||||
static constexpr std::array<std::array<int, 3>, 12 > indices
|
||||
= {{ { 0, 1, 2 }, { 0, 2, 3 }, { 2, 3, 4 }, { 2, 4, 7 },
|
||||
{ 3, 4, 5 }, { 3, 5, 0 }, { 4, 5, 6 }, { 4, 6, 7 },
|
||||
{ 5, 6, 1 }, { 5, 1, 0 }, { 6, 7, 2 }, { 6, 2, 1 } }};
|
||||
return Triangle (c[indices[idx][0]], c[indices[idx][1]], c[indices[idx][2]]);
|
||||
};
|
||||
|
||||
assemble_covariance_matrix_3
|
||||
(make_subiterator<Triangle, 12> (first, converter),
|
||||
make_subiterator<Triangle, 12> (beyond),
|
||||
covariance, c, K(), (Triangle*)nullptr, CGAL::Dimension_tag<2>(),
|
||||
Eigen_diagonalize_traits<FT, 3>());
|
||||
|
||||
#else
|
||||
// This variant uses the standard formulas but seems to be broken
|
||||
// (line/plane estimated appear to be wrong). In the absence of a
|
||||
// reliable fix so far, the above workaround applying PCA to a
|
||||
// decomposition of the cuboid into triangles is used.
|
||||
|
||||
typedef typename K::FT FT;
|
||||
typedef typename K::Iso_cuboid_3 Iso_cuboid;
|
||||
typedef typename Eigen::Matrix<FT, 3, 3> Matrix;
|
||||
|
|
@ -225,26 +269,25 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
const Iso_cuboid& t = *it;
|
||||
|
||||
// defined for convenience.
|
||||
FT x0 = t[0].x();
|
||||
FT y0 = t[0].y();
|
||||
FT z0 = t[0].z();
|
||||
FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0,
|
||||
t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
|
||||
t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
|
||||
Matrix transformation (delta);
|
||||
FT area = std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
|
||||
delta[6]*delta[6],1/3.0)*std::pow(delta[1]*delta[1] +
|
||||
delta[4]*delta[4] + delta[7]*delta[7],1/3.0)*2 +
|
||||
std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
|
||||
delta[6]*delta[6],1/3.0)*std::pow(delta[2]*delta[2] +
|
||||
delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2 +
|
||||
std::pow(delta[1]*delta[1] + delta[4]*delta[4] +
|
||||
delta[7]*delta[7],1/3.0)*std::pow(delta[2]*delta[2] +
|
||||
delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2;
|
||||
FT x0 = t.xmin();
|
||||
FT y0 = t.ymin();
|
||||
FT z0 = t.zmin();
|
||||
|
||||
FT x1 = t.xmax();
|
||||
FT y1 = t.ymax();
|
||||
FT z1 = t.zmax();
|
||||
|
||||
Matrix transformation;
|
||||
transformation << x1 - x0, 0 , 0 ,
|
||||
0 , y1 - y0, 0 ,
|
||||
0 , 0 , z1 - z0;
|
||||
|
||||
FT area = FT(2) * ((x1-x0)*(y1-y0) + (x1-x0)*(z1-z0) + (y1-y0)*(z1-z0));
|
||||
|
||||
// skip zero measure primitives
|
||||
if(area == (FT)0.0)
|
||||
continue;
|
||||
CGAL_assertion(area > 0.0);
|
||||
|
||||
// Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
|
||||
|
||||
|
|
@ -252,9 +295,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
transformation = area * transformation * moment * transformation.transpose();
|
||||
|
||||
// Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
|
||||
FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
|
||||
FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
|
||||
FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
|
||||
FT xav0 = (x1 - x0) / (2.0);
|
||||
FT yav0 = (y1 - y0) / (2.0);
|
||||
FT zav0 = (z1 - z0) / (2.0);
|
||||
|
||||
// and add to covariance matrix
|
||||
covariance[0] += transformation(0,0) + area * (2*x0*xav0 + x0*x0);
|
||||
|
|
@ -267,15 +310,17 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
mass += area;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.z() * c.x());
|
||||
covariance[3] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[4] += mass * (-1.0 * c.z() * c.y());
|
||||
covariance[5] += mass * (-1.0 * c.z() * c.z());
|
||||
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.z() * c.x());
|
||||
covariance[3] += -mass * (c.y() * c.y());
|
||||
covariance[4] += -mass * (c.z() * c.y());
|
||||
covariance[5] += -mass * (c.z() * c.z());
|
||||
#endif
|
||||
}
|
||||
|
||||
// assemble covariance matrix from a sphere set
|
||||
|
|
@ -305,9 +350,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
|
||||
// assemble 2nd order moment about the origin.
|
||||
Matrix moment;
|
||||
moment << 4.0/15.0, 0.0, 0.0,
|
||||
0.0, 4.0/15.0, 0.0,
|
||||
0.0, 0.0, 4.0/15.0;
|
||||
moment << FT(4.0/15.0), FT(0.0), FT(0.0),
|
||||
FT(0.0), FT(4.0/15.0), FT(0.0),
|
||||
FT(0.0), FT(0.0), FT(4.0/15.0);
|
||||
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
|
|
@ -318,7 +363,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
const Sphere& t = *it;
|
||||
|
||||
// defined for convenience.
|
||||
FT radius = std::sqrt(t.squared_radius());
|
||||
FT radius = CGAL::approximate_sqrt(t.squared_radius());
|
||||
Matrix transformation;
|
||||
transformation << radius, 0.0, 0.0,
|
||||
0.0, radius, 0.0,
|
||||
|
|
@ -350,14 +395,16 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
mass += volume;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.z() * c.x());
|
||||
covariance[3] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[4] += mass * (-1.0 * c.z() * c.y());
|
||||
covariance[5] += mass * (-1.0 * c.z() * c.z());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.z() * c.x());
|
||||
covariance[3] += -mass * (c.y() * c.y());
|
||||
covariance[4] += -mass * (c.z() * c.y());
|
||||
covariance[5] += -mass * (c.z() * c.z());
|
||||
|
||||
}
|
||||
// assemble covariance matrix from a sphere set
|
||||
|
|
@ -387,9 +434,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
|
||||
// assemble 2nd order moment about the origin.
|
||||
Matrix moment;
|
||||
moment << 4.0/3.0, 0.0, 0.0,
|
||||
0.0, 4.0/3.0, 0.0,
|
||||
0.0, 0.0, 4.0/3.0;
|
||||
moment << FT(4.0/3.0), FT(0.0), FT(0.0),
|
||||
FT(0.0), FT(4.0/3.0), FT(0.0),
|
||||
FT(0.0), FT(0.0), FT(4.0/3.0);
|
||||
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
|
|
@ -401,7 +448,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
|
||||
// defined for convenience.
|
||||
// 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;
|
||||
transformation << radius, 0.0, 0.0,
|
||||
0.0, radius, 0.0,
|
||||
|
|
@ -433,14 +480,16 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
mass += area;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.z() * c.x());
|
||||
covariance[3] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[4] += mass * (-1.0 * c.z() * c.y());
|
||||
covariance[5] += mass * (-1.0 * c.z() * c.z());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.z() * c.x());
|
||||
covariance[3] += -mass * (c.y() * c.y());
|
||||
covariance[4] += -mass * (c.z() * c.y());
|
||||
covariance[5] += -mass * (c.z() * c.z());
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -471,9 +520,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
// 5
|
||||
// assemble 2nd order moment about the origin.
|
||||
Matrix moment;
|
||||
moment << 1.0/60.0, 1.0/120.0, 1.0/120.0,
|
||||
1.0/120.0, 1.0/60.0, 1.0/120.0,
|
||||
1.0/120.0, 1.0/120.0, 1.0/60.0;
|
||||
moment << FT(1.0/60.0), FT(1.0/120.0), FT(1.0/120.0),
|
||||
FT(1.0/120.0), FT(1.0/60.0), FT(1.0/120.0),
|
||||
FT(1.0/120.0), FT(1.0/120.0), FT(1.0/60.0);
|
||||
|
||||
Matrix accum; // zero by default
|
||||
accum << 0, 0, 0, 0, 0, 0, 0, 0, 0;
|
||||
|
|
@ -548,9 +597,9 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
|
||||
// assemble 2nd order moment about the origin.
|
||||
Matrix moment;
|
||||
moment << 1.0/3.0, 0.5/3.0, 0.0,
|
||||
0.5/3.0, 1.0/3.0, 0.0,
|
||||
0.0, 0.0, 0.0;
|
||||
moment << FT(1.0/3.0), FT(0.5/3.0), FT(0.0),
|
||||
FT(0.5/3.0), FT(1.0/3.0), FT(0.0),
|
||||
FT(0.0), FT(0.0), FT(0.0);
|
||||
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
|
|
@ -566,8 +615,7 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
transformation << t[0].x(), t[1].x(), 0.0,
|
||||
t[0].y(), t[1].y(), 0.0,
|
||||
t[0].z(), t[1].z(), 1.0;
|
||||
using std::sqrt;
|
||||
FT length = sqrt(t.squared_length());
|
||||
FT length = CGAL::approximate_sqrt(t.squared_length());
|
||||
|
||||
// skip zero measure primitives
|
||||
if(length == (FT)0.0)
|
||||
|
|
@ -589,14 +637,16 @@ assemble_covariance_matrix_3(InputIterator first,
|
|||
mass += length;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.z() * c.x());
|
||||
covariance[3] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[4] += mass * (-1.0 * c.z() * c.y());
|
||||
covariance[5] += mass * (-1.0 * c.z() * c.z());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.z() * c.x());
|
||||
covariance[3] += -mass * (c.y() * c.y());
|
||||
covariance[4] += -mass * (c.z() * c.y());
|
||||
covariance[5] += -mass * (c.z() * c.z());
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,99 @@
|
|||
// Copyright (c) 2020 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org)
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Simon Giraudot
|
||||
//
|
||||
|
||||
#ifndef CGAL_PCA_SUBITERATOR_H
|
||||
#define CGAL_PCA_SUBITERATOR_H
|
||||
|
||||
#include <boost/iterator/iterator_facade.hpp>
|
||||
|
||||
namespace CGAL
|
||||
{
|
||||
|
||||
template <typename InputIterator, typename ValueType, int Size>
|
||||
class Subiterator
|
||||
: public boost::iterator_facade<Subiterator<InputIterator, ValueType, Size>,
|
||||
ValueType,
|
||||
std::input_iterator_tag>
|
||||
{
|
||||
public:
|
||||
using Self = Subiterator<InputIterator, ValueType, Size>;
|
||||
using Facade = boost::iterator_facade<Self, ValueType, std::input_iterator_tag>;
|
||||
|
||||
using Input_type = typename std::iterator_traits<InputIterator>::value_type;
|
||||
using Output_type = ValueType;
|
||||
|
||||
using Converter = std::function<Output_type(const Input_type&, int)>;
|
||||
|
||||
private:
|
||||
|
||||
Converter m_converter;
|
||||
int m_next_index;
|
||||
InputIterator m_base;
|
||||
mutable Output_type m_current;
|
||||
|
||||
public:
|
||||
|
||||
Subiterator() { }
|
||||
|
||||
Subiterator(InputIterator begin, const Converter& converter)
|
||||
: m_converter(converter), m_next_index(0), m_base(begin)
|
||||
{ }
|
||||
|
||||
Subiterator(InputIterator end)
|
||||
: m_next_index(0), m_base(end)
|
||||
{ }
|
||||
|
||||
private:
|
||||
|
||||
friend class boost::iterator_core_access;
|
||||
|
||||
void increment()
|
||||
{
|
||||
++ m_next_index;
|
||||
if (m_next_index == Size)
|
||||
{
|
||||
++ m_base;
|
||||
m_next_index = 0;
|
||||
}
|
||||
}
|
||||
|
||||
bool equal(const Self& other) const
|
||||
{
|
||||
return this->m_base == other.m_base && this->m_next_index == other.m_next_index;
|
||||
}
|
||||
|
||||
Output_type& dereference() const
|
||||
{
|
||||
m_current = m_converter (*m_base, m_next_index);
|
||||
return const_cast<Output_type&>(m_current);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename ValueType, int Size, typename InputIterator>
|
||||
Subiterator<InputIterator, ValueType, Size>
|
||||
make_subiterator (InputIterator begin,
|
||||
const typename Subiterator<InputIterator, ValueType, Size>::Converter& converter)
|
||||
{
|
||||
return Subiterator<InputIterator, ValueType, Size>(begin, converter);
|
||||
}
|
||||
|
||||
template <typename ValueType, int Size, typename InputIterator>
|
||||
Subiterator<InputIterator, ValueType, Size>
|
||||
make_subiterator (InputIterator end)
|
||||
{
|
||||
return Subiterator<InputIterator, ValueType, Size>(end);
|
||||
}
|
||||
|
||||
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_PCA_SUBITERATOR_H
|
||||
|
|
@ -83,7 +83,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
|
||||
// defined for convenience.
|
||||
// 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,
|
||||
0.0, radius};
|
||||
Matrix transformation = init_matrix<FT>(2,delta);
|
||||
|
|
@ -107,11 +107,13 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
mass += area;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.y() * c.y());
|
||||
|
||||
// solve for eigenvalues and eigenvectors.
|
||||
// eigen values are sorted in ascending order,
|
||||
|
|
@ -133,7 +135,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// isotropic case (infinite number of directions)
|
||||
// by default: assemble a line that goes through
|
||||
// the centroid and with a default horizontal vector.
|
||||
line = Line(c, Vector(1.0, 0.0));
|
||||
line = Line(c, Vector(FT(1), FT(0)));
|
||||
return (FT)0.0;
|
||||
}
|
||||
|
||||
|
|
@ -187,7 +189,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
|
||||
// defined for convenience.
|
||||
// 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,
|
||||
0.0, radius};
|
||||
Matrix transformation = init_matrix<FT>(2,delta);
|
||||
|
|
@ -197,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.
|
||||
|
||||
// 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.
|
||||
FT x0 = t.center().x();
|
||||
|
|
@ -211,11 +213,13 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
mass += length;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.y() * c.y());
|
||||
|
||||
// solve for eigenvalues and eigenvectors.
|
||||
// eigen values are sorted in ascending order,
|
||||
|
|
@ -237,7 +241,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// isotropic case (infinite number of directions)
|
||||
// by default: assemble a line that goes through
|
||||
// the centroid and with a default horizontal vector.
|
||||
line = Line(c, Vector(1.0, 0.0));
|
||||
line = Line(c, Vector(FT(1), FT(0)));
|
||||
return (FT)0.0;
|
||||
}
|
||||
} // end linear_least_squares_fitting_2 for circle set with 1D tag
|
||||
|
|
|
|||
|
|
@ -20,6 +20,7 @@
|
|||
#include <CGAL/PCA_util.h>
|
||||
#include <CGAL/linear_least_squares_fitting_points_3.h>
|
||||
#include <CGAL/linear_least_squares_fitting_segments_3.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
#include <list>
|
||||
#include <iterator>
|
||||
|
|
@ -112,29 +113,22 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Iso_cuboid& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[2],t[3]));
|
||||
segments.push_back(Segment(t[3],t[0]));
|
||||
segments.push_back(Segment(t[4],t[5]));
|
||||
segments.push_back(Segment(t[5],t[6]));
|
||||
segments.push_back(Segment(t[6],t[7]));
|
||||
segments.push_back(Segment(t[7],t[4]));
|
||||
segments.push_back(Segment(t[5],t[0]));
|
||||
segments.push_back(Segment(t[1],t[6]));
|
||||
segments.push_back(Segment(t[2],t[7]));
|
||||
segments.push_back(Segment(t[3],t[4]));
|
||||
}
|
||||
auto converter = [](const Iso_cuboid& c, int idx) -> Segment
|
||||
{
|
||||
if (idx < 7)
|
||||
return Segment (c[idx], c[idx+1]);
|
||||
if (idx < 10)
|
||||
return Segment (c[idx - 6], c[(idx-1)%8]);
|
||||
if (idx == 10)
|
||||
return Segment (c[5], c[0]);
|
||||
CGAL_assertion (idx == 11);
|
||||
return Segment (c[7], c[4]);
|
||||
};
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Segment, 12> (first, converter),
|
||||
make_subiterator<Segment, 12> (beyond),
|
||||
plane,c,(Segment*)nullptr,k,tag,diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_cuboids_3
|
||||
|
||||
|
|
@ -154,29 +148,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Point_3 Point;
|
||||
typedef typename K::Iso_cuboid_3 Iso_cuboid;
|
||||
auto converter = [](const Iso_cuboid& c, int idx) -> Point { return c[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Iso_cuboid& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
points.push_back(t[2]);
|
||||
points.push_back(t[3]);
|
||||
points.push_back(t[4]);
|
||||
points.push_back(t[5]);
|
||||
points.push_back(t[6]);
|
||||
points.push_back(t[7]);
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 8> (first, converter),
|
||||
make_subiterator<Point, 8> (beyond),
|
||||
plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_cuboids_3
|
||||
|
||||
|
|
@ -259,33 +240,26 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Segment_3 Segment;
|
||||
typedef typename K::Iso_cuboid_3 Iso_cuboid;
|
||||
auto converter = [](const Iso_cuboid& c, int idx) -> Segment
|
||||
{
|
||||
if (idx < 7)
|
||||
return Segment (c[idx], c[idx+1]);
|
||||
if (idx < 10)
|
||||
return Segment (c[idx - 6], c[(idx-1)%8]);
|
||||
if (idx == 10)
|
||||
return Segment (c[5], c[0]);
|
||||
CGAL_assertion (idx == 11);
|
||||
return Segment (c[7], c[4]);
|
||||
};
|
||||
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Iso_cuboid& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[2],t[3]));
|
||||
segments.push_back(Segment(t[3],t[0]));
|
||||
segments.push_back(Segment(t[4],t[5]));
|
||||
segments.push_back(Segment(t[5],t[6]));
|
||||
segments.push_back(Segment(t[6],t[7]));
|
||||
segments.push_back(Segment(t[7],t[4]));
|
||||
segments.push_back(Segment(t[5],t[0]));
|
||||
segments.push_back(Segment(t[1],t[6]));
|
||||
segments.push_back(Segment(t[2],t[7]));
|
||||
segments.push_back(Segment(t[3],t[4]));
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Segment, 12> (first, converter),
|
||||
make_subiterator<Segment, 12> (beyond),
|
||||
line,c,(Segment*)nullptr,k,tag,diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_cuboids_3
|
||||
|
||||
|
|
@ -305,29 +279,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Point_3 Point;
|
||||
typedef typename K::Iso_cuboid_3 Iso_cuboid;
|
||||
auto converter = [](const Iso_cuboid& c, int idx) -> Point { return c[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Iso_cuboid& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
points.push_back(t[2]);
|
||||
points.push_back(t[3]);
|
||||
points.push_back(t[4]);
|
||||
points.push_back(t[5]);
|
||||
points.push_back(t[6]);
|
||||
points.push_back(t[7]);
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 8> (first, converter),
|
||||
make_subiterator<Point, 8> (beyond),
|
||||
line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_cuboids_3
|
||||
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@
|
|||
#include <CGAL/centroid.h>
|
||||
#include <CGAL/Linear_algebraCd.h>
|
||||
#include <CGAL/PCA_util.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
#include <iterator>
|
||||
#include <list>
|
||||
|
|
@ -67,8 +68,8 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0. }};
|
||||
|
||||
// assemble 2nd order moment about the origin.
|
||||
FT temp[4] = {1/3.0, 0.25,
|
||||
0.25, 1/3.0};
|
||||
FT temp[4] = {FT(1/3.0), FT(0.25),
|
||||
FT(0.25), FT(1/3.0)};
|
||||
Matrix moment = init_matrix<FT>(2,temp);
|
||||
|
||||
for(InputIterator it = first;
|
||||
|
|
@ -100,8 +101,8 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
transformation = area * transformation * moment * LA::transpose(transformation);
|
||||
|
||||
// Translate the 2nd order moment to the center of the rectangle.
|
||||
FT xav0 = (x1-x0)/2.0;
|
||||
FT yav0 = (y2-y0)/2.0;
|
||||
FT xav0 = (x1-x0)/FT(2);
|
||||
FT yav0 = (y2-y0)/FT(2);
|
||||
// and add to covariance matrix
|
||||
covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);
|
||||
covariance[1] += transformation[0][1] + area * (x0*yav0 + xav0*y0 + x0*y0);
|
||||
|
|
@ -110,11 +111,13 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
mass += area;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.y() * c.y());
|
||||
|
||||
// solve for eigenvalues and eigenvectors.
|
||||
// eigen values are sorted in ascending order,
|
||||
|
|
@ -136,7 +139,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// isotropic case (infinite number of directions)
|
||||
// by default: assemble a line that goes through
|
||||
// the centroid and with a default horizontal vector.
|
||||
line = Line(c, Vector(1.0, 0.0));
|
||||
line = Line(c, Vector(FT(1), FT(0)));
|
||||
return (FT)0.0;
|
||||
}
|
||||
} // end linear_least_squares_fitting_2 for rectangle set with 2D tag
|
||||
|
|
@ -155,25 +158,16 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// types
|
||||
typedef typename K::Iso_rectangle_2 Iso_rectangle;
|
||||
typedef typename K::Segment_2 Segment_2;
|
||||
auto converter = [](const Iso_rectangle& r, int idx) -> Segment_2 { return Segment_2(r[idx], r[(idx+1)%4]); };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment_2> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Iso_rectangle& t = *it;
|
||||
segments.push_back(Segment_2(t[0],t[1]));
|
||||
segments.push_back(Segment_2(t[1],t[2]));
|
||||
segments.push_back(Segment_2(t[2],t[3]));
|
||||
segments.push_back(Segment_2(t[3],t[0]));
|
||||
}
|
||||
|
||||
return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,
|
||||
(Segment_2*)nullptr, K(),tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_2
|
||||
(make_subiterator<Segment_2, 4> (first, converter),
|
||||
make_subiterator<Segment_2, 4> (beyond),
|
||||
line,c,(Segment_2*)nullptr,K(),tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_2 for rectangle set with 1D tag
|
||||
|
||||
|
|
@ -194,25 +188,16 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// types
|
||||
typedef typename K::Iso_rectangle_2 Iso_rectangle;
|
||||
typedef typename K::Point_2 Point_2;
|
||||
auto converter = [](const Iso_rectangle& r, int idx) -> Point_2 { return r[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point_2> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Iso_rectangle& t = *it;
|
||||
points.push_back(Point_2(t[0]));
|
||||
points.push_back(Point_2(t[1]));
|
||||
points.push_back(Point_2(t[2]));
|
||||
points.push_back(Point_2(t[3]));
|
||||
}
|
||||
|
||||
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,
|
||||
(Point_2*)nullptr, K(),tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_2
|
||||
(make_subiterator<Point_2, 4> (first, converter),
|
||||
make_subiterator<Point_2, 4> (beyond),
|
||||
line,c,(Point_2*)nullptr,K(),tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_2 for rectangle set with 0D tag
|
||||
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@
|
|||
#include <CGAL/centroid.h>
|
||||
#include <CGAL/Linear_algebraCd.h>
|
||||
#include <CGAL/PCA_util.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
#include <iterator>
|
||||
#include <list>
|
||||
|
|
@ -67,7 +68,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
|
||||
// assemble 2nd order moment about the origin.
|
||||
FT temp[4] = {1.0, 0.5, 0.5, 1.0};
|
||||
Matrix moment = (1.0/3.0) * init_matrix<FT>(2,temp);
|
||||
Matrix moment = FT(1.0/3.0) * init_matrix<FT>(2,temp);
|
||||
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
|
|
@ -83,7 +84,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
t[0].y(), t[1].y()};
|
||||
Matrix transformation = init_matrix<FT>(2,delta);
|
||||
using std::sqrt;
|
||||
FT length = sqrt(t.squared_length());
|
||||
FT length = CGAL::approximate_sqrt(t.squared_length());
|
||||
CGAL_assertion(length != 0.0);
|
||||
|
||||
// Find the 2nd order moment for the segment wrt to the origin by an affine transformation.
|
||||
|
|
@ -99,11 +100,13 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
mass += length;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[0] += -mass * ( c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.y() * c.y());
|
||||
|
||||
// solve for eigenvalues and eigenvectors.
|
||||
// eigen values are sorted in ascending order,
|
||||
|
|
@ -125,7 +128,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// isotropic case (infinite number of directions)
|
||||
// by default: assemble a line that goes through
|
||||
// the centroid and with a default horizontal vector.
|
||||
line = Line(c, Vector(1.0, 0.0));
|
||||
line = Line(c, Vector(FT(1), FT(0)));
|
||||
return (FT)0.0;
|
||||
}
|
||||
} // end linear_least_squares_fitting_2 for segment set with 1D tag
|
||||
|
|
@ -144,21 +147,16 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// types
|
||||
typedef typename K::Point_2 Point;
|
||||
typedef typename K::Segment_2 Segment;
|
||||
auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Segment& s = *it;
|
||||
points.push_back(s[0]);
|
||||
points.push_back(s[1]);
|
||||
}
|
||||
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_2
|
||||
(make_subiterator<Point, 2> (first, converter),
|
||||
make_subiterator<Point, 2> (beyond),
|
||||
line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_2 for segment set with 1D tag
|
||||
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@
|
|||
#include <CGAL/centroid.h>
|
||||
#include <CGAL/PCA_util.h>
|
||||
#include <CGAL/linear_least_squares_fitting_points_3.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
#include <list>
|
||||
#include <iterator>
|
||||
|
|
@ -74,24 +75,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Segment_3 Segment;
|
||||
typedef typename K::Point_3 Point;
|
||||
auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Segment& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 2> (first, converter),
|
||||
make_subiterator<Point, 2> (beyond),
|
||||
plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_segments_3
|
||||
|
||||
// fits a line to a 3D segment set
|
||||
|
|
@ -141,24 +134,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Segment_3 Segment;
|
||||
typedef typename K::Point_3 Point;
|
||||
auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Segment& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 2> (first, converter),
|
||||
make_subiterator<Point, 2> (beyond),
|
||||
line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_segments_3
|
||||
|
||||
} // end namespace internal
|
||||
|
|
|
|||
|
|
@ -22,6 +22,7 @@
|
|||
#include <CGAL/linear_least_squares_fitting_points_3.h>
|
||||
#include <CGAL/linear_least_squares_fitting_segments_3.h>
|
||||
#include <CGAL/linear_least_squares_fitting_triangles_3.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
#include <iterator>
|
||||
|
||||
|
|
@ -75,26 +76,17 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
auto converter = [](const Tetrahedron& t, int idx) -> Triangle
|
||||
{ return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Triangle> triangles;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Tetrahedron& t = *it;
|
||||
triangles.push_back(Triangle(t[0],t[1],t[2]));
|
||||
triangles.push_back(Triangle(t[0],t[2],t[3]));
|
||||
triangles.push_back(Triangle(t[0],t[3],t[1]));
|
||||
triangles.push_back(Triangle(t[3],t[1],t[2]));
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,c,(Triangle*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Triangle, 4> (first, converter),
|
||||
make_subiterator<Triangle, 4> (beyond),
|
||||
plane,c,(Triangle*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_tetrahedrons_3
|
||||
|
||||
// fits a plane to a 3D tetrahedron set
|
||||
|
|
@ -113,28 +105,21 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename K::Segment_3 Segment;
|
||||
auto converter = [](const Tetrahedron& t, int idx) -> Segment
|
||||
{
|
||||
if (idx < 4)
|
||||
return Segment (t[idx], t[(idx+1)%4]);
|
||||
return Segment (t[idx-4], t[idx-2]);
|
||||
};
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Tetrahedron& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[1],t[3]));
|
||||
segments.push_back(Segment(t[2],t[3]));
|
||||
segments.push_back(Segment(t[0],t[2]));
|
||||
segments.push_back(Segment(t[0],t[3]));
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Segment, 6> (first, converter),
|
||||
make_subiterator<Segment, 6> (beyond),
|
||||
plane,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_tetrahedrons_3
|
||||
|
||||
|
|
@ -154,26 +139,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename K::Point_3 Point;
|
||||
auto converter = [](const Tetrahedron& t, int idx) -> Point { return t[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Tetrahedron& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
points.push_back(t[2]);
|
||||
points.push_back(t[3]);
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 4> (first, converter),
|
||||
make_subiterator<Point, 4> (beyond),
|
||||
plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_tetrahedrons_3
|
||||
|
||||
// fits a line to a 3D tetrahedron set
|
||||
|
|
@ -223,26 +198,17 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
auto converter = [](const Tetrahedron& t, int idx) -> Triangle
|
||||
{ return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Triangle> triangles;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Tetrahedron& t = *it;
|
||||
triangles.push_back(Triangle(t[0],t[1],t[2]));
|
||||
triangles.push_back(Triangle(t[0],t[2],t[3]));
|
||||
triangles.push_back(Triangle(t[0],t[3],t[1]));
|
||||
triangles.push_back(Triangle(t[3],t[1],t[2]));
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,c,(Triangle*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Triangle, 4> (first, converter),
|
||||
make_subiterator<Triangle, 4> (beyond),
|
||||
line,c,(Triangle*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_tetrahedrons_3
|
||||
|
||||
// fits a line to a 3D tetrahedron set
|
||||
|
|
@ -261,28 +227,21 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename K::Segment_3 Segment;
|
||||
auto converter = [](const Tetrahedron& t, int idx) -> Segment
|
||||
{
|
||||
if (idx < 4)
|
||||
return Segment (t[idx], t[(idx+1)%4]);
|
||||
return Segment (t[idx-4], t[idx-2]);
|
||||
};
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Tetrahedron& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[1],t[3]));
|
||||
segments.push_back(Segment(t[2],t[3]));
|
||||
segments.push_back(Segment(t[0],t[2]));
|
||||
segments.push_back(Segment(t[0],t[3]));
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Segment, 6> (first, converter),
|
||||
make_subiterator<Segment, 6> (beyond),
|
||||
line,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_tetrahedrons_3
|
||||
|
||||
// fits a line to a 3D tetrahedron set
|
||||
|
|
@ -301,26 +260,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Tetrahedron_3 Tetrahedron;
|
||||
typedef typename K::Point_3 Point;
|
||||
auto converter = [](const Tetrahedron& t, int idx) -> Point { return t[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Tetrahedron& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
points.push_back(t[2]);
|
||||
points.push_back(t[3]);
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 4> (first, converter),
|
||||
make_subiterator<Point, 4> (beyond),
|
||||
line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_tetrahedra_3
|
||||
|
||||
} // end namespace internal
|
||||
|
|
|
|||
|
|
@ -68,8 +68,8 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0. }};
|
||||
|
||||
// assemble the 2nd order moment about the origin.
|
||||
FT temp[4] = {1/12.0, 1/24.0,
|
||||
1/24.0, 1/12.0};
|
||||
FT temp[4] = {FT(1/12.0), FT(1/24.0),
|
||||
FT(1/24.0), FT(1/12.0)};
|
||||
|
||||
Matrix moment = init_matrix<FT>(2,temp);
|
||||
|
||||
|
|
@ -88,7 +88,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
t[1].y() - y0, t[2].y() - y0};
|
||||
|
||||
Matrix transformation = init_matrix<FT>(2,delta);
|
||||
FT area = 0.5 * std::abs(LA::determinant(transformation));
|
||||
FT area = FT(0.5) * CGAL::abs(LA::determinant(transformation));
|
||||
CGAL_assertion(area!=0);
|
||||
|
||||
// Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
|
||||
|
|
@ -97,8 +97,8 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
transformation = 2 * area * transformation * moment * LA::transpose(transformation);
|
||||
|
||||
// Translate the 2nd order moment to (x0,y0).
|
||||
FT xav0 = (delta[0]+delta[1])/3.0;
|
||||
FT yav0 = (delta[2]+delta[3])/3.0;
|
||||
FT xav0 = (delta[0]+delta[1])/FT(3);
|
||||
FT yav0 = (delta[2]+delta[3])/FT(3);
|
||||
|
||||
// and add to the covariance matrix
|
||||
covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);
|
||||
|
|
@ -108,11 +108,13 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
mass += area;
|
||||
}
|
||||
|
||||
CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
|
||||
|
||||
// Translate the 2nd order moment calculated about the origin to
|
||||
// the center of mass to get the covariance.
|
||||
covariance[0] += mass * (-1.0 * c.x() * c.x());
|
||||
covariance[1] += mass * (-1.0 * c.x() * c.y());
|
||||
covariance[2] += mass * (-1.0 * c.y() * c.y());
|
||||
covariance[0] += -mass * (c.x() * c.x());
|
||||
covariance[1] += -mass * (c.x() * c.y());
|
||||
covariance[2] += -mass * (c.y() * c.y());
|
||||
|
||||
// std::cout<<"cov: "<<covariance[0]*covariance[2]<<" =? "<<covariance[1]*covariance[1]<<std::endl;
|
||||
|
||||
|
|
@ -136,7 +138,7 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// isotropic case (infinite number of directions)
|
||||
// by default: assemble a line that goes through
|
||||
// the centroid and with a default horizontal vector.
|
||||
line = Line(c, Vector(1.0, 0.0));
|
||||
line = Line(c, Vector(FT(1), FT(0)));
|
||||
return (FT)0.0;
|
||||
}
|
||||
} // end linear_least_squares_fitting_2 for triangle set with 2D tag
|
||||
|
|
@ -157,23 +159,16 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
// types
|
||||
typedef typename Kernel::Triangle_2 Triangle;
|
||||
typedef typename Kernel::Segment_2 Segment;
|
||||
auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Triangle& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[2],t[0]));
|
||||
}
|
||||
|
||||
return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,tag,Kernel(),
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_2
|
||||
(make_subiterator<Segment, 3> (first, converter),
|
||||
make_subiterator<Segment, 3> (beyond),
|
||||
line,c,(Segment*)nullptr,Kernel(),tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_2 for triangle set with 1D tag
|
||||
|
||||
|
|
@ -194,24 +189,16 @@ linear_least_squares_fitting_2(InputIterator first,
|
|||
|
||||
typedef typename Kernel::Triangle_2 Triangle;
|
||||
typedef typename Kernel::Point_2 Point;
|
||||
auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Triangle& t = *it;
|
||||
points.push_back(Point(t[0]));
|
||||
points.push_back(Point(t[1]));
|
||||
points.push_back(Point(t[2]));
|
||||
}
|
||||
|
||||
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,tag,Kernel(),
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_2
|
||||
(make_subiterator<Point, 3> (first, converter),
|
||||
make_subiterator<Point, 3> (beyond),
|
||||
line,c,(Point*)nullptr,Kernel(),tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_2 for triangle set with 0D tag
|
||||
|
||||
} // end namespace internal
|
||||
|
|
|
|||
|
|
@ -18,6 +18,7 @@
|
|||
#include <CGAL/basic.h>
|
||||
#include <CGAL/centroid.h>
|
||||
#include <CGAL/PCA_util.h>
|
||||
#include <CGAL/Subiterator.h>
|
||||
|
||||
#include <list>
|
||||
#include <iterator>
|
||||
|
|
@ -73,24 +74,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
typedef typename K::Segment_3 Segment;
|
||||
auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Triangle& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[2],t[0]));
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Segment, 3> (first, converter),
|
||||
make_subiterator<Segment, 3> (beyond),
|
||||
plane,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_triangles_3
|
||||
|
||||
|
|
@ -110,23 +103,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
typedef typename K::Point_3 Point;
|
||||
auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Triangle& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
points.push_back(t[2]);
|
||||
}
|
||||
|
||||
// compute fitting plane
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 3> (first, converter),
|
||||
make_subiterator<Point, 3> (beyond),
|
||||
plane,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
} // end linear_least_squares_fitting_triangles_3
|
||||
|
||||
|
|
@ -177,25 +163,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
typedef typename K::Segment_3 Segment;
|
||||
auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
|
||||
std::list<Segment> segments;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Triangle& t = *it;
|
||||
segments.push_back(Segment(t[0],t[1]));
|
||||
segments.push_back(Segment(t[1],t[2]));
|
||||
segments.push_back(Segment(t[2],t[0]));
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Segment, 3> (first, converter),
|
||||
make_subiterator<Segment, 3> (beyond),
|
||||
line,c,(Segment*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_triangles_3
|
||||
|
||||
// fits a line to a 3D triangle set
|
||||
|
|
@ -214,24 +191,16 @@ linear_least_squares_fitting_3(InputIterator first,
|
|||
{
|
||||
typedef typename K::Triangle_3 Triangle;
|
||||
typedef typename K::Point_3 Point;
|
||||
auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
|
||||
|
||||
// precondition: at least one element in the container.
|
||||
CGAL_precondition(first != beyond);
|
||||
std::list<Point> points;
|
||||
for(InputIterator it = first;
|
||||
it != beyond;
|
||||
it++)
|
||||
{
|
||||
const Triangle& t = *it;
|
||||
points.push_back(t[0]);
|
||||
points.push_back(t[1]);
|
||||
points.push_back(t[2]);
|
||||
}
|
||||
|
||||
// compute fitting line
|
||||
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
|
||||
return linear_least_squares_fitting_3
|
||||
(make_subiterator<Point, 3> (first, converter),
|
||||
make_subiterator<Point, 3> (beyond),
|
||||
line,c,(Point*)nullptr,k,tag,
|
||||
diagonalize_traits);
|
||||
} // end linear_least_squares_fitting_triangles_3
|
||||
|
||||
} // end namespace internal
|
||||
|
|
|
|||
|
|
@ -0,0 +1,107 @@
|
|||
#include <CGAL/Simple_cartesian.h>
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/linear_least_squares_fitting_2.h>
|
||||
#include <CGAL/linear_least_squares_fitting_3.h>
|
||||
|
||||
template <typename K>
|
||||
typename K::Point_2 point_2 (int x, int y)
|
||||
{ return typename K::Point_2 (typename K::FT(x), typename K::FT(y)); }
|
||||
template <typename K>
|
||||
typename K::Point_3 point_3 (int x, int y, int z)
|
||||
{ return typename K::Point_3 (typename K::FT(x), typename K::FT(y), typename K::FT(z)); }
|
||||
|
||||
// Generate default objects so that the test compiles *AND* runs fine
|
||||
template <typename K, typename O> O default_object(const O&) { return O(); }
|
||||
template <typename K> typename K::Point_2 default_object(const typename K::Point_2&)
|
||||
{ return point_2<K>(0,0); }
|
||||
template <typename K> typename K::Segment_2 default_object(const typename K::Segment_2&)
|
||||
{ return typename K::Segment_2(point_2<K>(0,0), point_2<K>(0,1)); }
|
||||
template <typename K> typename K::Circle_2 default_object(const typename K::Circle_2&)
|
||||
{ return typename K::Circle_2(point_2<K>(0,0), typename K::FT(1.0)); }
|
||||
template <typename K> typename K::Triangle_2 default_object(const typename K::Triangle_2&)
|
||||
{ return typename K::Triangle_2(point_2<K>(0,0), point_2<K>(0,1), point_2<K>(1,0)); }
|
||||
template <typename K> typename K::Iso_rectangle_2 default_object(const typename K::Iso_rectangle_2&)
|
||||
{ return typename K::Iso_rectangle_2(point_2<K>(0,0), point_2<K>(1,1)); }
|
||||
template <typename K> typename K::Point_3 default_object(const typename K::Point_3&)
|
||||
{ return point_3<K>(0,0,0); }
|
||||
template <typename K> typename K::Segment_3 default_object(const typename K::Segment_3&)
|
||||
{ return typename K::Segment_3(point_3<K>(0,0,0), point_3<K>(0,0,1)); }
|
||||
template <typename K> typename K::Sphere_3 default_object(const typename K::Sphere_3&)
|
||||
{ return typename K::Sphere_3(point_3<K>(0,0,0), typename K::FT(1.0)); }
|
||||
template <typename K> typename K::Triangle_3 default_object(const typename K::Triangle_3&)
|
||||
{ return typename K::Triangle_3(point_3<K>(0,0,0), point_3<K>(0,0,1), point_3<K>(0,1,0)); }
|
||||
template <typename K> typename K::Tetrahedron_3 default_object(const typename K::Tetrahedron_3&)
|
||||
{ return typename K::Tetrahedron_3(point_3<K>(0,0,0), point_3<K>(0,0,1),
|
||||
point_3<K>(0,1,0), point_3<K>(1,0,0)); }
|
||||
template <typename K> typename K::Iso_cuboid_3 default_object(const typename K::Iso_cuboid_3&)
|
||||
{ return typename K::Iso_cuboid_3(point_3<K>(0,0,0), point_3<K>(1,1,1)); }
|
||||
|
||||
template <typename Kernel, typename Object, int dim>
|
||||
void test_2d()
|
||||
{
|
||||
std::array<Object, 1> dummy = { default_object<Kernel>(Object()) };
|
||||
typename Kernel::Line_2 line;
|
||||
typename Kernel::Point_2 centroid;
|
||||
CGAL::linear_least_squares_fitting_2 (dummy.begin(), dummy.end(), line, centroid,
|
||||
CGAL::Dimension_tag<dim>(), Kernel(),
|
||||
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 2>());
|
||||
}
|
||||
|
||||
template <typename Kernel, typename Object, int dim>
|
||||
void test_3d()
|
||||
{
|
||||
std::array<Object, 1> dummy = { default_object<Kernel>(Object()) };
|
||||
typename Kernel::Line_3 line;
|
||||
typename Kernel::Plane_3 plane;
|
||||
typename Kernel::Point_3 centroid;
|
||||
CGAL::linear_least_squares_fitting_3 (dummy.begin(), dummy.end(), line, centroid,
|
||||
CGAL::Dimension_tag<dim>(), Kernel(),
|
||||
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>());
|
||||
CGAL::linear_least_squares_fitting_3 (dummy.begin(), dummy.end(), plane, centroid,
|
||||
CGAL::Dimension_tag<dim>(), Kernel(),
|
||||
CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>());
|
||||
}
|
||||
|
||||
template <typename Kernel>
|
||||
void test_kernel()
|
||||
{
|
||||
test_2d<Kernel, typename Kernel::Point_2, 0>();
|
||||
test_2d<Kernel, typename Kernel::Segment_2, 1>();
|
||||
test_2d<Kernel, typename Kernel::Segment_2, 0>();
|
||||
test_2d<Kernel, typename Kernel::Circle_2, 2>();
|
||||
test_2d<Kernel, typename Kernel::Circle_2, 1>();
|
||||
test_2d<Kernel, typename Kernel::Triangle_2, 2>();
|
||||
test_2d<Kernel, typename Kernel::Triangle_2, 1>();
|
||||
test_2d<Kernel, typename Kernel::Triangle_2, 0>();
|
||||
test_2d<Kernel, typename Kernel::Iso_rectangle_2, 2>();
|
||||
test_2d<Kernel, typename Kernel::Iso_rectangle_2, 1>();
|
||||
test_2d<Kernel, typename Kernel::Iso_rectangle_2, 0>();
|
||||
test_3d<Kernel, typename Kernel::Point_3, 0>();
|
||||
test_3d<Kernel, typename Kernel::Segment_3, 1>();
|
||||
test_3d<Kernel, typename Kernel::Segment_3, 0>();
|
||||
test_3d<Kernel, typename Kernel::Sphere_3, 3>();
|
||||
test_3d<Kernel, typename Kernel::Sphere_3, 2>();
|
||||
test_3d<Kernel, typename Kernel::Triangle_3, 2>();
|
||||
test_3d<Kernel, typename Kernel::Triangle_3, 1>();
|
||||
test_3d<Kernel, typename Kernel::Triangle_3, 0>();
|
||||
test_3d<Kernel, typename Kernel::Tetrahedron_3, 3>();
|
||||
test_3d<Kernel, typename Kernel::Tetrahedron_3, 2>();
|
||||
test_3d<Kernel, typename Kernel::Tetrahedron_3, 1>();
|
||||
test_3d<Kernel, typename Kernel::Tetrahedron_3, 0>();
|
||||
test_3d<Kernel, typename Kernel::Iso_cuboid_3, 3>();
|
||||
test_3d<Kernel, typename Kernel::Iso_cuboid_3, 2>();
|
||||
test_3d<Kernel, typename Kernel::Iso_cuboid_3, 1>();
|
||||
test_3d<Kernel, typename Kernel::Iso_cuboid_3, 0>();
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test_kernel<CGAL::Simple_cartesian<float> >();
|
||||
test_kernel<CGAL::Simple_cartesian<double> >();
|
||||
test_kernel<CGAL::Exact_predicates_inexact_constructions_kernel>();
|
||||
test_kernel<CGAL::Exact_predicates_exact_constructions_kernel>();
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -225,7 +225,7 @@ void assert_quality (const std::vector<Point_3>& points, const Fitted& fitted)
|
|||
std::cerr << "mean distance = " << mean_dist << std::endl;
|
||||
|
||||
CGAL_assertion_code
|
||||
(double limit = 1e-3 * std::sqrt (CGAL::squared_distance (points.front(), points.back())));
|
||||
(double limit = 1e-5 * std::sqrt (CGAL::squared_distance (points.front(), points.back())));
|
||||
CGAL_assertion (mean_dist < limit);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -120,9 +120,7 @@ centroid(InputIterator begin,
|
|||
it++)
|
||||
{
|
||||
const Segment& s = *it;
|
||||
using std::abs;
|
||||
using std::sqrt;
|
||||
FT length = sqrt(abs(s.squared_length()));
|
||||
FT length = CGAL::approximate_sqrt(CGAL::abs(s.squared_length()));
|
||||
Point c = K().construct_midpoint_2_object()(s[0],s[1]);
|
||||
v = v + length * (c - ORIGIN);
|
||||
sum_lengths += length;
|
||||
|
|
@ -216,7 +214,7 @@ centroid(InputIterator begin,
|
|||
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]);
|
||||
v = v + unsigned_area * (c - ORIGIN);
|
||||
sum_areas += unsigned_area;
|
||||
|
|
@ -252,7 +250,7 @@ centroid(InputIterator begin,
|
|||
it++)
|
||||
{
|
||||
const Circle& s = *it;
|
||||
FT radius = std::sqrt(s.squared_radius());
|
||||
FT radius = CGAL::approximate_sqrt(s.squared_radius());
|
||||
Point c = s.center();
|
||||
v = v + radius * (c - ORIGIN);
|
||||
sum_lengths += radius;
|
||||
|
|
@ -381,7 +379,7 @@ centroid(InputIterator begin,
|
|||
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]);
|
||||
v = v + unsigned_area * (c - ORIGIN);
|
||||
sum_areas += unsigned_area;
|
||||
|
|
@ -442,8 +440,7 @@ centroid(InputIterator begin,
|
|||
it++)
|
||||
{
|
||||
const Segment& s = *it;
|
||||
using std::sqrt;
|
||||
FT length = sqrt(s.squared_length());
|
||||
FT length = CGAL::approximate_sqrt(s.squared_length());
|
||||
Point c = CGAL::midpoint(s.source(),s.target());
|
||||
// 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);
|
||||
|
|
@ -539,7 +536,7 @@ centroid(InputIterator begin,
|
|||
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]);
|
||||
v = v + unsigned_area * (c - ORIGIN);
|
||||
sum_areas += unsigned_area;
|
||||
|
|
@ -608,7 +605,7 @@ centroid(InputIterator begin,
|
|||
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();
|
||||
v = v + unsigned_volume * (c - ORIGIN);
|
||||
sum_volumes += unsigned_volume;
|
||||
|
|
@ -717,7 +714,7 @@ centroid(InputIterator begin,
|
|||
{
|
||||
const Iso_cuboid& cuboid = *it;
|
||||
FT unsigned_area = 2 * ((cuboid.xmax()-cuboid.xmin())*(cuboid.ymax()-cuboid.ymin()) + (cuboid.xmax()-cuboid.xmin())*(cuboid.zmax()-cuboid.zmin()) + (cuboid.ymax()-cuboid.ymin())*(cuboid.zmax()-cuboid.zmin()));
|
||||
Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
|
||||
Point c = K().construct_midpoint_3_object()(cuboid[0],cuboid[7]);
|
||||
v = v + unsigned_area * (c - ORIGIN);
|
||||
sum_areas += unsigned_area;
|
||||
}
|
||||
|
|
@ -750,7 +747,7 @@ centroid(InputIterator begin,
|
|||
{
|
||||
const Iso_cuboid& cuboid = *it;
|
||||
FT unsigned_volume = cuboid.volume();
|
||||
Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
|
||||
Point c = K().construct_midpoint_3_object()(cuboid[0],cuboid[7]);
|
||||
v = v + unsigned_volume * (c - ORIGIN);
|
||||
sum_volumes += unsigned_volume;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -26,10 +26,20 @@
|
|||
|
||||
//#define DO_NOT_USE_EIGEN_COMPUTEDIRECT_FOR_DIAGONALIZATION
|
||||
|
||||
#include <CGAL/array.h>
|
||||
#include <CGAL/number_utils.h>
|
||||
#include <CGAL/double.h>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace internal {
|
||||
|
||||
template <typename FT>
|
||||
struct Restricted_FT { typedef double type; };
|
||||
template <>
|
||||
struct Restricted_FT<float> { typedef float type; };
|
||||
|
||||
}
|
||||
|
||||
/// \ingroup PkgSolverInterfaceRef
|
||||
///
|
||||
/// The class `Eigen_diagonalize_traits` provides an interface to the
|
||||
|
|
@ -47,14 +57,15 @@ namespace CGAL {
|
|||
template <typename FT, unsigned int dim = 3>
|
||||
class Eigen_diagonalize_traits
|
||||
{
|
||||
typedef typename internal::Restricted_FT<FT>::type RFT;
|
||||
public:
|
||||
typedef std::array<FT, dim> Vector;
|
||||
typedef std::array<FT, dim*dim> Matrix;
|
||||
typedef std::array<FT, (dim * (dim+1) / 2)> Covariance_matrix;
|
||||
|
||||
private:
|
||||
typedef Eigen::Matrix<FT, dim, dim> EigenMatrix;
|
||||
typedef Eigen::Matrix<FT, dim, 1> EigenVector;
|
||||
typedef Eigen::Matrix<RFT, dim, dim> EigenMatrix;
|
||||
typedef Eigen::Matrix<RFT, dim, 1> EigenVector;
|
||||
|
||||
/// Construct the covariance matrix
|
||||
static EigenMatrix construct_covariance_matrix(const Covariance_matrix& cov)
|
||||
|
|
@ -65,7 +76,7 @@ private:
|
|||
{
|
||||
for(std::size_t j=i; j<dim; ++j)
|
||||
{
|
||||
m(i,j) = static_cast<float>(cov[(dim * i) + j - ((i * (i+1)) / 2)]);
|
||||
m(i,j) = static_cast<RFT>(CGAL::to_double(cov[(dim * i) + j - ((i * (i+1)) / 2)]));
|
||||
|
||||
if(i != j)
|
||||
m(j,i) = m(i,j);
|
||||
|
|
|
|||
Loading…
Reference in New Issue