songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Allow user to assemble covariance matrix on triangles using Matrix::FT != Kernel::FT

This commit is contained in:
Simon Giraudot 2018-07-06 10:10:19 +02:00
parent 7467130cfd
commit f3baef9c42
5 changed files with 47 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
Principal_component_analysis/include/CGAL/PCA_util.h
View File

@ -32,13 +32,13 @@ namespace CGAL {
namespace internal { namespace internal {
// Initialize a matrix in n dimension by an array or numbers // Initialize a matrix in n dimension by an array or numbers
template <typename K> template <typename FT>
typename CGAL::Linear_algebraCd<typename K::FT>::Matrix typename CGAL::Linear_algebraCd<FT>::Matrix
init_matrix(const int n, init_matrix(const int n,
typename K::FT entries[]) FT entries[])
{ {
CGAL_assertion(n > 1); // dimension > 1 CGAL_assertion(n > 1); // dimension > 1
typedef typename CGAL::Linear_algebraCd<typename K::FT>::Matrix Matrix; typedef typename CGAL::Linear_algebraCd<FT>::Matrix Matrix;
Matrix m(n); Matrix m(n);
int i,j; int i,j;
@ -102,24 +102,28 @@ assemble_covariance_matrix_3(InputIterator first,
const CGAL::Dimension_tag<2>&, const CGAL::Dimension_tag<2>&,
const DiagonalizeTraits&) const DiagonalizeTraits&)
{ {
typedef typename K::FT FT; // Use FT from the DiagonalizeTraits to avoid warnings if the Kernel
// uses double but the user wants to compute diagonalization using floats
typedef typename DiagonalizeTraits::Vector::value_type FT;
typedef typename K::Triangle_3 Triangle; typedef typename K::Triangle_3 Triangle;
typedef typename CGAL::Linear_algebraCd<FT> LA; typedef typename CGAL::Linear_algebraCd<FT> LA;
typedef typename LA::Matrix Matrix; typedef typename LA::Matrix Matrix;
// assemble covariance matrix as a semi-definite matrix. // assemble covariance matrix as a semi-definite matrix.
// Matrix numbering: // Matrix numbering:
// 0 1 2 // 0 1 2
// 3 4 // 3 4
// 5 // 5
//Final combined covariance matrix for all triangles and their combined mass //Final combined covariance matrix for all triangles and their combined mass
FT mass = 0.0; FT mass = FT(0.0);
// assemble 2nd order moment about the origin. // assemble 2nd order moment about the origin.
FT temp[9] = {1.0/12.0, 1.0/24.0, 1.0/24.0, FT temp[9] = {FT(1.0/12.0), FT(1.0/24.0), FT(1.0/24.0),
1.0/24.0, 1.0/12.0, 1.0/24.0, FT(1.0/24.0), FT(1.0/12.0), FT(1.0/24.0),
1.0/24.0, 1.0/24.0, 1.0/12.0}; FT(1.0/24.0), FT(1.0/24.0), FT(1.0/12.0)};
Matrix moment = init_matrix<K>(3,temp); Matrix moment = init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -130,11 +134,11 @@ assemble_covariance_matrix_3(InputIterator first,
const Triangle& t = *it; const Triangle& t = *it;
// defined for convenience. // defined for convenience.
FT delta[9] = {t[0].x(), t[1].x(), t[2].x(), FT delta[9] = {FT(t[0].x()), FT(t[1].x()), FT(t[2].x()),
t[0].y(), t[1].y(), t[2].y(), FT(t[0].y()), FT(t[1].y()), FT(t[2].y()),
t[0].z(), t[1].z(), t[2].z()}; FT(t[0].z()), FT(t[1].z()), FT(t[2].z())};
Matrix transformation = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT area = std::sqrt(t.squared_area()); FT area = FT(std::sqrt(t.squared_area()));
// skip zero measure primitives // skip zero measure primitives
if(area == (FT)0.0) if(area == (FT)0.0)
@ -158,12 +162,12 @@ assemble_covariance_matrix_3(InputIterator first,
// Translate the 2nd order moment calculated about the origin to // Translate the 2nd order moment calculated about the origin to
// the center of mass to get the covariance. // the center of mass to get the covariance.
covariance[0] += mass * (-1.0 * c.x() * c.x()); covariance[0] += mass * FT(-1.0 * c.x() * c.x());
covariance[1] += mass * (-1.0 * c.x() * c.y()); covariance[1] += mass * FT(-1.0 * c.x() * c.y());
covariance[2] += mass * (-1.0 * c.z() * c.x()); covariance[2] += mass * FT(-1.0 * c.z() * c.x());
covariance[3] += mass * (-1.0 * c.y() * c.y()); covariance[3] += mass * FT(-1.0 * c.y() * c.y());
covariance[4] += mass * (-1.0 * c.z() * c.y()); covariance[4] += mass * FT(-1.0 * c.z() * c.y());
covariance[5] += mass * (-1.0 * c.z() * c.z()); covariance[5] += mass * FT(-1.0 * c.z() * c.z());
} }
@ -198,7 +202,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT temp[9] = {(FT)(1.0/3.0), (FT)(1.0/4.0), (FT)(1.0/4.0), FT temp[9] = {(FT)(1.0/3.0), (FT)(1.0/4.0), (FT)(1.0/4.0),
(FT)(1.0/4.0), (FT)(1.0/3.0), (FT)(1.0/4.0), (FT)(1.0/4.0), (FT)(1.0/3.0), (FT)(1.0/4.0),
(FT)(1.0/4.0), (FT)(1.0/4.0), (FT)(1.0/3.0)}; (FT)(1.0/4.0), (FT)(1.0/4.0), (FT)(1.0/3.0)};
Matrix moment = init_matrix<K>(3,temp); Matrix moment = init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -216,7 +220,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0, 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].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 = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT volume = t.volume(); FT volume = t.volume();
// skip zero measure primitives // skip zero measure primitives
@ -285,7 +289,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT temp[9] = {(FT)(7.0/3.0), (FT)1.5, (FT)1.5, FT temp[9] = {(FT)(7.0/3.0), (FT)1.5, (FT)1.5,
(FT)1.5, (FT)(7.0/3.0), (FT)1.5, (FT)1.5, (FT)(7.0/3.0), (FT)1.5,
(FT)1.5, (FT)1.5, (FT)(7.0/3.0)}; (FT)1.5, (FT)1.5, (FT)(7.0/3.0)};
Matrix moment = init_matrix<K>(3,temp); Matrix moment = init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -302,7 +306,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0, 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].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 = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT area = std::pow(delta[0]*delta[0] + delta[3]*delta[3] + 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[6]*delta[6],1/3.0)*std::pow(delta[1]*delta[1] +
delta[4]*delta[4] + delta[7]*delta[7],1/3.0)*2 + delta[4]*delta[4] + delta[7]*delta[7],1/3.0)*2 +
@ -380,7 +384,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT temp[9] = {4.0/15.0, 0.0, 0.0, FT temp[9] = {4.0/15.0, 0.0, 0.0,
0.0, 4.0/15.0, 0.0, 0.0, 4.0/15.0, 0.0,
0.0, 0.0, 4.0/15.0}; 0.0, 0.0, 4.0/15.0};
Matrix moment = init_matrix<K>(3,temp); Matrix moment = init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -395,7 +399,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT delta[9] = {radius, 0.0, 0.0, FT delta[9] = {radius, 0.0, 0.0,
0.0, radius, 0.0, 0.0, radius, 0.0,
0.0, 0.0, radius}; 0.0, 0.0, radius};
Matrix transformation = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT volume = (FT)(4.0/3.0) * radius * t.squared_radius(); FT volume = (FT)(4.0/3.0) * radius * t.squared_radius();
// skip zero measure primitives // skip zero measure primitives
@ -464,7 +468,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT temp[9] = {4.0/3.0, 0.0, 0.0, FT temp[9] = {4.0/3.0, 0.0, 0.0,
0.0, 4.0/3.0, 0.0, 0.0, 4.0/3.0, 0.0,
0.0, 0.0, 4.0/3.0}; 0.0, 0.0, 4.0/3.0};
Matrix moment = init_matrix<K>(3,temp); Matrix moment = init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -480,7 +484,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT delta[9] = {radius, 0.0, 0.0, FT delta[9] = {radius, 0.0, 0.0,
0.0, radius, 0.0, 0.0, radius, 0.0,
0.0, 0.0, radius}; 0.0, 0.0, radius};
Matrix transformation = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT area = (FT)4.0 * t.squared_radius(); FT area = (FT)4.0 * t.squared_radius();
// skip zero measure primitives // skip zero measure primitives
@ -550,7 +554,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT temp[9] = {1.0/60.0, 1.0/120.0, 1.0/120.0, FT temp[9] = {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/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/60.0};
Matrix moment = init_matrix<K>(3,temp); Matrix moment = init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -568,7 +572,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT delta[9] = {t[1].x()-x0, t[2].x()-x0, t[3].x()-x0, FT delta[9] = {t[1].x()-x0, t[2].x()-x0, t[3].x()-x0,
t[1].y()-y0, t[2].y()-y0, t[3].y()-y0, t[1].y()-y0, t[2].y()-y0, t[3].y()-y0,
t[1].z()-z0, t[2].z()-z0, t[3].z()-z0}; t[1].z()-z0, t[2].z()-z0, t[3].z()-z0};
Matrix transformation = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT volume = t.volume(); FT volume = t.volume();
// skip zero measure primitives // skip zero measure primitives
@ -637,7 +641,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT temp[9] = {1.0, 0.5, 0.0, FT temp[9] = {1.0, 0.5, 0.0,
0.5, 1.0, 0.0, 0.5, 1.0, 0.0,
0.0, 0.0, 0.0}; 0.0, 0.0, 0.0};
Matrix moment = (FT)(1.0/3.0) * init_matrix<K>(3,temp); Matrix moment = (FT)(1.0/3.0) * init_matrix<FT>(3,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -652,7 +656,7 @@ assemble_covariance_matrix_3(InputIterator first,
FT delta[9] = {t[0].x(), t[1].x(), 0.0, FT delta[9] = {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};
Matrix transformation = init_matrix<K>(3,delta); Matrix transformation = init_matrix<FT>(3,delta);
FT length = std::sqrt(t.squared_length()); FT length = std::sqrt(t.squared_length());
// skip zero measure primitives // skip zero measure primitives

8
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
View File

@ -79,7 +79,7 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble 2nd order moment about the origin. // assemble 2nd order moment about the origin.
FT temp[4] = {0.25, 0.0, FT temp[4] = {0.25, 0.0,
0.0, 0.25}; 0.0, 0.25};
Matrix moment = init_matrix<K>(2,temp); Matrix moment = init_matrix<FT>(2,temp);
// Matrix moment = Matrix(2,true,PI); // Matrix moment = Matrix(2,true,PI);
for(InputIterator it = first; for(InputIterator it = first;
@ -95,7 +95,7 @@ linear_least_squares_fitting_2(InputIterator first,
FT radius = std::sqrt(t.squared_radius()); FT radius = std::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<K>(2,delta); Matrix transformation = init_matrix<FT>(2,delta);
FT area = t.squared_radius(); FT area = t.squared_radius();
CGAL_assertion(area != 0.0); CGAL_assertion(area != 0.0);
@ -184,7 +184,7 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble 2nd order moment about the origin. // assemble 2nd order moment about the origin.
FT temp[4] = {1.0, 0.0, FT temp[4] = {1.0, 0.0,
0.0, 1.0}; 0.0, 1.0};
Matrix moment = init_matrix<K>(2,temp); Matrix moment = init_matrix<FT>(2,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -199,7 +199,7 @@ linear_least_squares_fitting_2(InputIterator first,
FT radius = std::sqrt(t.squared_radius()); FT radius = std::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<K>(2,delta); Matrix transformation = init_matrix<FT>(2,delta);
FT length = 2 * radius; FT length = 2 * radius;
CGAL_assertion(length != 0.0); CGAL_assertion(length != 0.0);

4
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
View File

@ -78,7 +78,7 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble 2nd order moment about the origin. // assemble 2nd order moment about the origin.
FT temp[4] = {1/3.0, 0.25, FT temp[4] = {1/3.0, 0.25,
0.25, 1/3.0}; 0.25, 1/3.0};
Matrix moment = init_matrix<K>(2,temp); Matrix moment = init_matrix<FT>(2,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -98,7 +98,7 @@ linear_least_squares_fitting_2(InputIterator first,
FT delta[4] = {x1-x0, 0.0, FT delta[4] = {x1-x0, 0.0,
0.0, y2-y0}; 0.0, y2-y0};
Matrix transformation = init_matrix<K>(2,delta); Matrix transformation = init_matrix<FT>(2,delta);
FT area = (x1-x0)*(y2-y0); FT area = (x1-x0)*(y2-y0);
CGAL_assertion(area != 0.0); CGAL_assertion(area != 0.0);

4
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
View File

@ -76,7 +76,7 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble 2nd order moment about the origin. // assemble 2nd order moment about the origin.
FT temp[4] = {1.0, 0.5, 0.5, 1.0}; FT temp[4] = {1.0, 0.5, 0.5, 1.0};
Matrix moment = (1.0/3.0) * init_matrix<K>(2,temp); Matrix moment = (1.0/3.0) * init_matrix<FT>(2,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -90,7 +90,7 @@ linear_least_squares_fitting_2(InputIterator first,
// FT example = CGAL::to_double(t[0].x()); // FT example = CGAL::to_double(t[0].x());
FT delta[4] = {t[0].x(), t[1].x(), FT delta[4] = {t[0].x(), t[1].x(),
t[0].y(), t[1].y()}; t[0].y(), t[1].y()};
Matrix transformation = init_matrix<K>(2,delta); Matrix transformation = init_matrix<FT>(2,delta);
FT length = std::sqrt(t.squared_length()); FT length = std::sqrt(t.squared_length());
CGAL_assertion(length != 0.0); CGAL_assertion(length != 0.0);

4
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
View File

@ -80,7 +80,7 @@ linear_least_squares_fitting_2(InputIterator first,
FT temp[4] = {1/12.0, 1/24.0, FT temp[4] = {1/12.0, 1/24.0,
1/24.0, 1/12.0}; 1/24.0, 1/12.0};
Matrix moment = init_matrix<Kernel>(2,temp); Matrix moment = init_matrix<FT>(2,temp);
for(InputIterator it = first; for(InputIterator it = first;
it != beyond; it != beyond;
@ -96,7 +96,7 @@ linear_least_squares_fitting_2(InputIterator first,
FT delta[4] = {t[1].x() - x0, t[2].x() - x0, FT delta[4] = {t[1].x() - x0, t[2].x() - x0,
t[1].y() - y0, t[2].y() - y0}; t[1].y() - y0, t[2].y() - y0};
Matrix transformation = init_matrix<Kernel>(2,delta); Matrix transformation = init_matrix<FT>(2,delta);
FT area = 0.5 * std::abs(LA::determinant(transformation)); FT area = 0.5 * std::abs(LA::determinant(transformation));
CGAL_assertion(area!=0); CGAL_assertion(area!=0);