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

Merge branch 'CGAL-Solver-package-GF-old' into CGAL-Solver-package-GF

This commit is contained in:
Simon Giraudot 2015-09-23 10:29:14 +02:00
parent edbda537d4 3d3469524f
commit 36723443a4
105 changed files with 1781 additions and 1763 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Barycentric_coordinates_2/doc/Barycentric_coordinates_2/dependencies
View File

@ -8,4 +8,3 @@ Interpolation
STL_Extension
Triangulation_2
Algebraic_foundations
Surface_mesh_parameterization

74
Bounding_volumes/include/CGAL/Approximate_min_ellipsoid_d/Approximate_min_ellipsoid_d_impl.h
View File

@ -18,8 +18,8 @@
//
// Author(s) : Kaspar Fischer <fischerk@inf.ethz.ch>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <CGAL/Approximate_min_ellipsoid_d.h>
@ -150,34 +150,32 @@ namespace CGAL {
{
CGAL_APPEL_ASSERT(d==2);
typedef Simple_cartesian<double> K;
typedef Vector_2<K> Vector_2;
// write matrix M' as [ a, b; b, c ]:
const double matrix[3] = { E->matrix(0, 0), // a
E->matrix(0, 1), // b
E->matrix(1, 1) }; // c
std::pair<Vector_2, Vector_2> eigenvectors; // Note: not neces. normalized.
std::pair<double, double> eigenvalues; // Note: sorted descendent.
internal::eigen_symmetric_2<K>(matrix, eigenvectors, eigenvalues);
const CGAL::cpp11::array<double, 3> matrix = {{ E->matrix(0, 0), // a
E->matrix(0, 1), // b
E->matrix(1, 1) }}; // c
CGAL::cpp11::array<double, 4> eigenvectors; // Note: not neces. normalized.
CGAL::cpp11::array<double, 2> eigenvalues; // Note: sorted ascendent.
CGAL::Default_diagonalize_traits<double, 2>::diagonalize_selfadjoint_covariance_matrix
(matrix, eigenvalues, eigenvectors);
// normalize eigenvectors:
double l1=1.0/std::sqrt(eigenvectors.first.x()*eigenvectors.first.x()+
eigenvectors.first.y()*eigenvectors.first.y());
double l2=1.0/std::sqrt(eigenvectors.second.x()*eigenvectors.second.x()+
eigenvectors.second.y()*eigenvectors.second.y());
double l1=1.0/std::sqrt(eigenvectors[2]*eigenvectors[2]+
eigenvectors[3]*eigenvectors[3]);
double l2=1.0/std::sqrt(eigenvectors[0]*eigenvectors[0]+
eigenvectors[1]*eigenvectors[1]);
// store axes lengths:
lengths_.push_back(std::sqrt(factor/eigenvalues.first));
lengths_.push_back(std::sqrt(factor/eigenvalues.second));
lengths_.push_back(std::sqrt(factor/eigenvalues[1]));
lengths_.push_back(std::sqrt(factor/eigenvalues[0]));
// store directions:
directions_.resize(2);
directions_[0].push_back(eigenvectors.first.x()*l1);
directions_[0].push_back(eigenvectors.first.y()*l1);
directions_[1].push_back(eigenvectors.second.x()*l2);
directions_[1].push_back(eigenvectors.second.y()*l2);
directions_[0].push_back(eigenvectors[2]*l1);
directions_[0].push_back(eigenvectors[3]*l1);
directions_[1].push_back(eigenvectors[0]*l2);
directions_[1].push_back(eigenvectors[1]*l2);
}
template<class Traits>
@ -192,16 +190,18 @@ namespace CGAL {
// M' = [ b d e ]
// [ c e f ]
//
const double matrix[6] = { E->matrix(0, 0), // a
E->matrix(0, 1), // b
E->matrix(1, 1), // d
E->matrix(0, 2), // c
E->matrix(1, 2), // e
E->matrix(2, 2) }; // f
const CGAL::cpp11::array<double, 6> matrix = {{ E->matrix(0, 0), // a
E->matrix(0, 1), // b
E->matrix(0, 2), // c
E->matrix(1, 1), // d
E->matrix(1, 2), // e
E->matrix(2, 2) }}; // f
double eigenvectors[3 * 3]; // Note: not necessarily normalized.
double eigenvalues[3]; // Note: sorted descendent.
internal::eigen_symmetric<double>(matrix, 3, eigenvectors, eigenvalues);
CGAL::cpp11::array<double, 9> eigenvectors; // Note: not necessarily normalized.
CGAL::cpp11::array<double, 3> eigenvalues; // Note: sorted ascendent.
CGAL::Default_diagonalize_traits<double, 3>::diagonalize_selfadjoint_covariance_matrix
(matrix, eigenvalues, eigenvectors);
// normalize eigenvectors:
double l1 = 1.0/std::sqrt(eigenvectors[0] * eigenvectors[0]+ // x^2
@ -220,15 +220,15 @@ namespace CGAL {
// store directions:
directions_.resize(3);
directions_[0].push_back(eigenvectors[0]*l1);
directions_[0].push_back(eigenvectors[1]*l1);
directions_[0].push_back(eigenvectors[2]*l1);
directions_[0].push_back(eigenvectors[6]*l3);
directions_[0].push_back(eigenvectors[7]*l3);
directions_[0].push_back(eigenvectors[8]*l3);
directions_[1].push_back(eigenvectors[3]*l2);
directions_[1].push_back(eigenvectors[4]*l2);
directions_[1].push_back(eigenvectors[5]*l2);
directions_[2].push_back(eigenvectors[6]*l3);
directions_[2].push_back(eigenvectors[7]*l3);
directions_[2].push_back(eigenvectors[8]*l3);
directions_[2].push_back(eigenvectors[0]*l1);
directions_[2].push_back(eigenvectors[1]*l1);
directions_[2].push_back(eigenvectors[2]*l1);
}
template<class Traits>

1
Documentation/doc/Documentation/packages.txt
View File

@ -150,6 +150,7 @@ h1 {
\package_listing{STL_Extension}
\package_listing{BGL}
\package_listing{Solver_interface}
\package_listing{Point_set_processing_3/Property_map}
\package_listing{Circulator}
\package_listing{Generator}

1
Documentation/scripts/generate_how_to_cite.py
View File

@ -108,6 +108,7 @@ def gen_bib_entry(title, authors, bib, anchor):
def gen_txt_entry(title, authors, bib, anchor,k):
title_r=title.replace("Kernel","%Kernel").replace("Interval","%Interval").replace("Matrix","%Matrix").replace("Kinetic","%Kinetic").replace("CGAL","%CGAL").replace("Range","%Range")
authors=authors.replace("CGAL","%CGAL")
res="<tr valign=\"top\">\n\
<td align=\"right\" class=\"bibtexnumber\">\n\
[<a name=\""+bib+"-${CGAL_RELEASE_YEAR_ID}\">"+str(k)+"</a>]\n\

20
Installation/changes.html
View File

@ -135,9 +135,29 @@ and <code>src/</code> directories).
<!-- Geometry Processing -->
<!-- Spatial Searching and Sorting -->
<!-- Geometric Optimization -->
<h3>Principal Component Analysis</h3>
<ul>
<li>
Add a template parameter <code>DiagonalizeTraits</code> for
functions <code>CGAL::linear_least_squares_fitting_2()</code>
and <code>CGAL::linear_least_squares_fitting_3()</code>. This
allows to either choose the legacy internal diagonalization code
from CGAL or the Eigen implementation (or any class that is a
model of <code>DiagonalizeTraits</code>). Variants of the
function that automatically deduce the kernel also automatically
select the diagonalizer, so the API is mostly preserved.
</li>
</ul>
<!-- Interpolation -->
<!-- Kinetic Data Structures -->
<!-- Support Library -->
<h3>CGAL and Solvers</h3>
<ul>
<li>
This package now includes all CGAL concepts for solvers with
models using the third party Eigen library.
</li>
</ul>
<!-- Visualization -->
<!-- end of the div for 4.8 -->

2
Jet_fitting_3/doc/Jet_fitting_3/CGAL/Monge_via_jet_fitting.h
View File

@ -24,7 +24,7 @@ the geometric classes and tools required by local
computations.
\tparam SvdTraits features the linear
algebra algorithm required by the fitting method.
algebra algorithm required by the fitting method. The scalar type, `SvdTraits::FT`, must be the same as that of the `LocalKernel` concept : `LocalKernel::FT`.
\sa `Eigen_svd`
\sa `Monge_form`

6
Jet_fitting_3/doc/Jet_fitting_3/PackageDescription.txt
View File

@ -13,7 +13,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.3}
\cgalPkgDependsOn{Solvers from \ref thirdpartyEigen.}
\cgalPkgDependsOn{\ref PkgSolverSummary}
\cgalPkgBib{cgal:pc-eldp}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{Operations on Polyhedra,polyhedron_3.zip}
@ -25,12 +25,12 @@
## Concepts ##
- `DataKernel`
- `LocalKernel`
- `SvdTraits`
## Classes ##
- \link CGAL::Monge_via_jet_fitting::Monge_form `CGAL::Monge_via_jet_fitting< DataKernel, LocalKernel, SvdTraits>::Monge_form` \endlink
- `CGAL::Monge_via_jet_fitting<DataKernel, LocalKernel, SvdTraits>`
- `CGAL::Eigen_svd`
## Global Functions ##
The insert operator (operator<< ) is overloaded for the class \link CGAL::Monge_via_jet_fitting::Monge_form `Monge_form` \endlink.

1
Jet_fitting_3/doc/Jet_fitting_3/dependencies
View File

@ -5,3 +5,4 @@ Algebraic_foundations
Circulator
Stream_support
Number_types
Solver_interface

49
Jet_fitting_3/include/CGAL/Monge_via_jet_fitting.h
View File

@ -22,8 +22,8 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/circulator.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/eigen.h>
#include <CGAL/Cartesian_converter.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <math.h>
#include <utility>
#ifdef CGAL_EIGEN3_ENABLED
@ -362,27 +362,30 @@ compute_PCA(InputIterator begin, InputIterator end)
// assemble covariance matrix as a
// semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
FT covariance[6] = {xx,xy,yy,xz,yz,zz};
FT eigen_values[3];
FT eigen_vectors[9];
// 0 1 2
// 3 4
// 5
CGAL::cpp11::array<FT, 6> covariance = {{ xx,xy,xz,yy,yz,zz }};
CGAL::cpp11::array<FT, 3> eigen_values = {{ 0., 0., 0. }};
CGAL::cpp11::array<FT, 9> eigen_vectors = {{ 0., 0., 0. }};
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
CGAL::internal::eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
CGAL::Default_diagonalize_traits<FT,3>::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
//store in m_pca_basis
for (int i=0; i<3; i++)
{
m_pca_basis[i].first = eigen_values[i];
m_pca_basis[i].first = eigen_values[2-i];
}
Vector_3 v1(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
Vector_3 v1(eigen_vectors[6],eigen_vectors[7],eigen_vectors[8]);
m_pca_basis[0].second = v1;
Vector_3 v2(eigen_vectors[3],eigen_vectors[4],eigen_vectors[5]);
m_pca_basis[1].second = v2;
Vector_3 v3(eigen_vectors[6],eigen_vectors[7],eigen_vectors[8]);
Vector_3 v3(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
m_pca_basis[2].second = v3;
switch_to_direct_orientation(m_pca_basis[0].second,
m_pca_basis[1].second,
@ -524,15 +527,17 @@ compute_Monge_basis(const FT* A, Monge_form& monge_form)
//in the new orthonormal basis (Y,Z) of the tangent plane :
weingarten = inv *(1/det) * weingarten * change_XuXv2YZ;
//switch to eigen_symmetric algo for diagonalization of weingarten
FT W[3] = {weingarten(0,0), weingarten(1,0), weingarten(1,1)};
FT eval[2];
FT evec[4];
//eval in decreasing order
CGAL::internal::eigen_symmetric<FT>(W,2,evec,eval);
// diagonalization of weingarten
CGAL::cpp11::array<FT,3> W = {{ weingarten(0,0), weingarten(1,0), weingarten(1,1) }};
CGAL::cpp11::array<FT,2> eval = {{ 0., 0. }};
CGAL::cpp11::array<FT,4> evec = {{ 0., 0., 0., 0. }};
Vector_3 d_max = evec[0]*Y + evec[1]*Z,
d_min = evec[2]*Y + evec[3]*Z;
//eval in increasing order
CGAL::Default_diagonalize_traits<FT,2>::diagonalize_selfadjoint_covariance_matrix
(W, eval, evec);
Vector_3 d_max = evec[2]*Y + evec[3]*Z,
d_min = evec[0]*Y + evec[1]*Z;
switch_to_direct_orientation(d_max, d_min, normal);
Aff_transformation change_basis (d_max[0], d_max[1], d_max[2],
@ -548,8 +553,8 @@ compute_Monge_basis(const FT* A, Monge_form& monge_form)
monge_form.maximal_principal_direction() = L2D_converter(this->change_world2fitting.inverse()(d_max));
monge_form.minimal_principal_direction() = L2D_converter(this->change_world2fitting.inverse()(d_min));
monge_form.normal_direction() = L2D_converter(this->change_world2fitting.inverse()(normal));
monge_form.coefficients()[0] = L2D_NTconverter()(eval[0]);
monge_form.coefficients()[1] = L2D_NTconverter()(eval[1]);
monge_form.coefficients()[0] = L2D_NTconverter()(eval[1]);
monge_form.coefficients()[1] = L2D_NTconverter()(eval[0]);
}
//end else
}

34
Point_set_processing_3/doc/Point_set_processing_3/Concepts/VCMTraits.h
View File

@ -1,34 +0,0 @@
/*!
\ingroup PkgPointSetProcessingConcepts
\cgalConcept
Concept providing functions to extract eigenvectors and eigenvalue from
covariance matrices represented by an array `a` as follows:
<center>
\f$ \begin{bmatrix}
a[0] & a[1] & a[2] \\
a[1] & a[3] & a[4] \\
a[2] & a[4] & a[5] \\
\end{bmatrix}\f$
</center>
\cgalHasModel `CGAL::Eigen_vcm_traits`
*/
class VCMTraits
{
public:
/// fill `eigenvalues` with the eigenvalues of the covariance matrix represented by `cov`.
/// Eigenvalues are sorted by increasing order.
/// \return `true` if the operation was successful and `false` otherwise.
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<double,6>& cov,
cpp11::array<double, 3>& eigenvalues);
/// Extract the eigenvector associated to the largest eigenvalue
/// of the covariance matrix represented by `cov`.
/// \return `true` if the operation was successful and `false` otherwise.
static bool
extract_largest_eigenvector_of_covariance_matrix (
const cpp11::array<double,6>& cov,
cpp11::array<double,3> &normal);
};

2
Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt
View File

@ -13,7 +13,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.5}
\cgalPkgDependsOn{\ref thirdpartyEigen "Eigen"}
\cgalPkgDependsOn{\ref PkgSolverSummary}
\cgalPkgBib{cgal:ass-psp}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{See Polyhedral Surface,polyhedron_3.zip}

1
Point_set_processing_3/doc/Point_set_processing_3/dependencies
View File

@ -9,4 +9,5 @@ Property_map
Bounding_volumes
Principal_component_analysis
Jet_fitting_3
Solver_interface

115
Point_set_processing_3/include/CGAL/Eigen_vcm_traits.h
View File

@ -1,115 +0,0 @@
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
// Author(s) : Jocelyn Meyron and Quentin Mérigot
//
#ifndef CGAL_EIGEN_VCM_TRAITS_H
#define CGAL_EIGEN_VCM_TRAITS_H
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <CGAL/array.h>
namespace CGAL {
/// A model of the concept `VCMTraits` using \ref thirdpartyEigen.
/// \cgalModels `VCMTraits`
class Eigen_vcm_traits{
// Construct the covariance matrix
static Eigen::Matrix3f
construct_covariance_matrix (const cpp11::array<double,6>& cov) {
Eigen::Matrix3f m;
m(0,0) = cov[0]; m(0,1) = cov[1]; m(0,2) = cov[2];
m(1,1) = cov[3]; m(1,2) = cov[4];
m(2,2) = cov[5];
m(1, 0) = m(0,1); m(2, 0) = m(0, 2); m(2, 1) = m(1, 2);
return m;
}
// Diagonalize a selfadjoint matrix
static bool
diagonalize_selfadjoint_matrix (Eigen::Matrix3f &m,
Eigen::Matrix3f &eigenvectors,
Eigen::Vector3f &eigenvalues) {
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigensolver(m);
if (eigensolver.info() != Eigen::Success) {
return false;
}
eigenvalues = eigensolver.eigenvalues();
eigenvectors = eigensolver.eigenvectors();
return true;
}
public:
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<double,6>& cov,
cpp11::array<double, 3>& eigenvalues)
{
Eigen::Matrix3f m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
Eigen::Vector3f eigenvalues_;
Eigen::Matrix3f eigenvectors_;
bool res = diagonalize_selfadjoint_matrix(m, eigenvectors_, eigenvalues_);
if (res)
{
eigenvalues[0]=eigenvalues_[0];
eigenvalues[1]=eigenvalues_[1];
eigenvalues[2]=eigenvalues_[2];
}
return res;
}
// Extract the eigenvector associated to the largest eigenvalue
static bool
extract_largest_eigenvector_of_covariance_matrix (
const cpp11::array<double,6>& cov,
cpp11::array<double,3> &normal)
{
// Construct covariance matrix
Eigen::Matrix3f m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
Eigen::Vector3f eigenvalues;
Eigen::Matrix3f eigenvectors;
if (! diagonalize_selfadjoint_matrix(m, eigenvectors, eigenvalues)) {
return false;
}
// Eigenvalues are already sorted by increasing order
normal[0]=eigenvectors(0,0);
normal[1]=eigenvectors(1,0);
normal[2]=eigenvectors(2,0);
return true;
}
};
} // namespace CGAL
#endif // CGAL_EIGEN_VCM_TRAITS_H

12
Point_set_processing_3/include/CGAL/internal/Voronoi_covariance_3/voronoi_covariance_3.h
View File

@ -79,8 +79,8 @@ namespace CGAL {
template <class Point>
inline void operator () (const Point &a,
const Point &c,
const Point &b)
const Point &b,
const Point &c)
{
internal::covariance_matrix_tetrahedron (a[0], a[1], a[2],
b[0], b[1], b[2],
@ -109,8 +109,8 @@ namespace CGAL {
template <class Point>
inline void operator () (const Point &a,
const Point &c,
const Point &b)
const Point &b,
const Point &c)
{
const double vol = CGAL::volume(a, b, c, Point(CGAL::ORIGIN));
//std::cerr << "vol = " << vol << "\n";
@ -168,12 +168,10 @@ namespace CGAL {
typename Polyhedron::Halfedge_around_facet_circulator
h0 = it->facet_begin(), hf = h0--, hs = cpp11::next(hf);
while(1)
while(hs != h0)
{
f (h0->vertex()->point(), hf->vertex()->point(),
hs->vertex()->point());
if (hs == h0)
break;
++hs; ++hf;
}
}

21
Point_set_processing_3/include/CGAL/vcm_estimate_edges.h
View File

@ -39,18 +39,21 @@ namespace CGAL {
/// which however would result in selecting more points on sharper regions.
/// More details are provided in \cgalCite{cgal:mog-vbcfe-11}.
///
/// \tparam VCM_traits is a model of `VCMTraits`. If Eigen 3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined
/// then an overload using `Eigen_vcm_traits` is provided and this template parameter can be omitted.
/// \sa `CGAL::compute_vcm()`
/// \tparam VCMTraits is a model of `DiagonalizeTraits`. It can be
/// omitted: if Eigen 3 (or greater) is available and
/// `CGAL_EIGEN3_ENABLED` is defined then an overload using
/// `Eigen_diagonalize_traits` is provided. Otherwise, the internal
/// implementation `Diagonalize_traits` is used.
/// \sa CGAL::compute_vcm()`
///
template <class FT, class VCM_traits>
template <class FT, class VCMTraits>
bool
vcm_is_on_feature_edge (cpp11::array<FT,6> &cov,
double threshold,
VCM_traits)
VCMTraits)
{
cpp11::array<double,3> eigenvalues;
if (!VCM_traits::
if (!VCMTraits::
diagonalize_selfadjoint_covariance_matrix(cov, eigenvalues) )
{
return false;
@ -66,15 +69,15 @@ vcm_is_on_feature_edge (cpp11::array<FT,6> &cov,
#ifdef CGAL_EIGEN3_ENABLED
template <class FT>
bool
vcm_is_on_feature_edge (cpp11::array<FT,6> &cov,
double threshold)
{
return vcm_is_on_feature_edge(cov, threshold, Eigen_vcm_traits());
return vcm_is_on_feature_edge(cov, threshold,
CGAL::Default_diagonalize_traits<double, 3>());
}
#endif
} // namespace CGAL

59
Point_set_processing_3/include/CGAL/vcm_estimate_normals.h
View File

@ -32,9 +32,7 @@
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Fuzzy_sphere.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_vcm_traits.h>
#endif
#include <CGAL/Default_diagonalize_traits.h>
#include <iterator>
#include <vector>
@ -290,7 +288,7 @@ compute_vcm (ForwardIterator first,
/// this number of neighbors.
// This variant requires all of the parameters.
template < typename VCM_traits,
template < typename VCMTraits,
typename ForwardIterator,
typename PointPMap,
typename NormalPMap,
@ -345,7 +343,7 @@ vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input
int i = 0;
for (ForwardIterator it = first; it != beyond; ++it) {
cpp11::array<double, 3> enormal = {{ 0,0,0 }};
VCM_traits::extract_largest_eigenvector_of_covariance_matrix
VCMTraits::extract_largest_eigenvector_of_covariance_matrix
(cov[i], enormal);
typename Kernel::Vector_3 normal(enormal[0],
@ -368,15 +366,17 @@ vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input
/// @tparam ForwardIterator iterator over input points.
/// @tparam PointPMap is a model of `ReadablePropertyMap` with a value_type = `Kernel::Point_3`.
/// @tparam NormalPMap is a model of `WritablePropertyMap` with a value_type = `Kernel::Vector_3`.
/// \tparam VCM_traits is a model of `VCMTraits`. If Eigen 3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined
/// then an overload using `Eigen_vcm_traits` is provided and this template parameter can be omitted.
/// \tparam VCMTraits is a model of `DiagonalizeTraits`. It can be
/// omitted: if Eigen 3 (or greater) is available and
/// `CGAL_EIGEN3_ENABLED` is defined then an overload using
/// `Eigen_diagonalize_traits` is provided. Otherwise, the internal
/// implementation `Diagonalize_traits` is used.
// This variant deduces the kernel from the point property map
// and uses a radius for the convolution.
template < typename ForwardIterator,
typename PointPMap,
typename NormalPMap,
typename VCM_traits
typename VCMTraits
>
void
vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input point.
@ -385,16 +385,16 @@ vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input
NormalPMap normal_pmap, ///< property map: value_type of ForwardIterator -> Vector_3.
double offset_radius, ///< offset radius.
double convolution_radius, ///< convolution radius.
VCM_traits
VCMTraits
)
{
typedef typename boost::property_traits<PointPMap>::value_type Point;
typedef typename Kernel_traits<Point>::Kernel Kernel;
vcm_estimate_normals<VCM_traits>(first, beyond,
point_pmap, normal_pmap,
offset_radius, convolution_radius,
Kernel());
vcm_estimate_normals<VCMTraits>(first, beyond,
point_pmap, normal_pmap,
offset_radius, convolution_radius,
Kernel());
}
@ -409,15 +409,18 @@ vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input
/// @tparam ForwardIterator iterator over input points.
/// @tparam PointPMap is a model of `ReadablePropertyMap` with a value_type = `Kernel::Point_3`.
/// @tparam NormalPMap is a model of `WritablePropertyMap` with a value_type = `Kernel::Vector_3`.
/// \tparam VCM_traits is a model of `VCMTraits`. If Eigen 3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined
/// then an overload using `Eigen_vcm_traits` is provided and this template parameter can be omitted.
/// \tparam VCMTraits is a model of `DiagonalizeTraits`. It can be
/// omitted: if Eigen 3 (or greater) is available and
/// `CGAL_EIGEN3_ENABLED` is defined then an overload using
/// `Eigen_diagonalize_traits` is provided. Otherwise, the internal
/// implementation `Diagonalize_traits` is used.
// This variant deduces the kernel from the point property map
// and uses a number of neighbors for the convolution.
template < typename ForwardIterator,
typename PointPMap,
typename NormalPMap,
typename VCM_traits
typename VCMTraits
>
void
vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input point.
@ -426,21 +429,20 @@ vcm_estimate_normals (ForwardIterator first, ///< iterator over the first input
NormalPMap normal_pmap, ///< property map: value_type of ForwardIterator -> Vector_3.
double offset_radius, ///< offset radius.
unsigned int k, ///< number of neighbor points used for the convolution.
VCM_traits
VCMTraits
)
{
typedef typename boost::property_traits<PointPMap>::value_type Point;
typedef typename Kernel_traits<Point>::Kernel Kernel;
vcm_estimate_normals<VCM_traits>(first, beyond,
point_pmap, normal_pmap,
offset_radius, 0,
Kernel(),
k);
vcm_estimate_normals<VCMTraits>(first, beyond,
point_pmap, normal_pmap,
offset_radius, 0,
Kernel(),
k);
}
#ifdef CGAL_EIGEN3_ENABLED
template < typename ForwardIterator,
typename PointPMap,
typename NormalPMap
@ -453,7 +455,8 @@ vcm_estimate_normals (ForwardIterator first,
double offset_radius,
double convolution_radius)
{
vcm_estimate_normals(first, beyond, point_pmap, normal_pmap, offset_radius, convolution_radius, Eigen_vcm_traits());
vcm_estimate_normals(first, beyond, point_pmap, normal_pmap, offset_radius, convolution_radius,
CGAL::Default_diagonalize_traits<double, 3>());
}
template < typename ForwardIterator,
@ -468,9 +471,11 @@ vcm_estimate_normals (ForwardIterator first,
double offset_radius,
unsigned int nb_neighbors_convolve)
{
vcm_estimate_normals(first, beyond, point_pmap, normal_pmap, offset_radius, nb_neighbors_convolve, Eigen_vcm_traits());
vcm_estimate_normals(first, beyond, point_pmap, normal_pmap, offset_radius, nb_neighbors_convolve,
CGAL::Default_diagonalize_traits<double, 3>());
}
#endif
/// @cond SKIP_IN_MANUAL
// This variant creates a default point property map = Identity_property_map

4
Point_set_processing_3/test/Point_set_processing_3/vcm_plane_test.cpp
View File

@ -30,11 +30,11 @@ int main (void) {
std::cout << "Normal is " << points[0].second << std::endl;
// The normal at the origin should be (0, 0, -1)
// The normal at the origin should be (0, 0, 1)
double epsilon=2e-5;
assert(points[0].second.x() < epsilon && points[0].second.x() > -epsilon);
assert(points[0].second.y() < epsilon && points[0].second.y() > -epsilon);
assert(points[0].second.z() < -1+epsilon && points[0].second.z() > -1-epsilon);
assert(points[0].second.z() < 1+epsilon && points[0].second.z() > 1-epsilon);
return 0;
}

2
Poisson_surface_reconstruction_3/doc/Poisson_surface_reconstruction_3/PackageDescription.txt
View File

@ -10,7 +10,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.5}
\cgalPkgDependsOn{\ref thirdpartyEigen}
\cgalPkgDependsOn{\ref PkgSolverSummary}
\cgalPkgBib{cgal:asg-srps}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{See Polyhedral Surface,polyhedron_3.zip}

2
Poisson_surface_reconstruction_3/doc/Poisson_surface_reconstruction_3/dependencies
View File

@ -6,4 +6,4 @@ Circulator
Stream_support
Surface_mesher
Point_set_processing_3
Surface_mesh_parameterization
Solver_interface

2
Polygon_mesh_processing/doc/Polygon_mesh_processing/NamedParameters.txt
View File

@ -91,7 +91,7 @@ Here is the list of the named parameters available in this package:
\cgalNPBegin{sparse_linear_solver} \anchor PMP_sparse_linear_solver
is the solver used for fairing of polygon meshes.\n
\b Type: a class model of `SparseLinearAlgebraTraitsWithFactor_d`.\n
\b Type: a class model of `SparseLinearAlgebraWithFactorTraits_d`.\n
\b Default: if \ref thirdpartyEigen "Eigen" 3.2 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined, then the following overload of `Eigen_solver_traits` is provided as default value\n
\code
CGAL::Eigen_solver_traits<

2
Polygon_mesh_processing/doc/Polygon_mesh_processing/PackageDescription.txt
View File

@ -26,7 +26,7 @@ ranging from basic operations on simplices, to complex geometry processing algor
\cgalPkgShortInfoBegin
\cgalPkgSince{4.7}
\cgalPkgDependsOn{documented for each function; Sparse square solver such as those from \ref thirdpartyEigen}
\cgalPkgDependsOn{documented for each function; \ref PkgSolverSummary}
\cgalPkgBib{cgal:lty-pmp}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{Operations on Polyhedra,polyhedron_3.zip}

2
Polygon_mesh_processing/doc/Polygon_mesh_processing/dependencies
View File

@ -6,7 +6,7 @@ Circulator
Stream_support
Polyhedron
BGL
Surface_mesh_parameterization
Solver_interface
Surface_mesh
Surface_modeling
AABB_tree

7
Principal_component_analysis/doc/Principal_component_analysis/CGAL/linear_least_squares_fitting_2.h
View File

@ -13,6 +13,8 @@ The tag `tag` identifies the dimension to be considered from the objects. For po
The class `K` is the kernel in which the value type of the `InputIterator` is defined. It can be omitted and deduced automatically from the value type.
The class `DiagonalizeTraits_` is a model of `DiagonalizeTraits`. It can be omitted: if Eigen 3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined then an overload using `Eigen_diagonalize_traits` is provided. Otherwise, the internal implementation `Diagonalize_traits` is used.
\cgalHeading{Requirements}
<OL>
@ -28,14 +30,15 @@ omitted.
\pre first != beyond.
*/
template < typename InputIterator, typename K, typename Tag >
template < typename InputIterator, typename K, typename Tag, typename DiagonalizeTraits_ >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
typename K::Line_2 & line,
typename K::Point_2 & centroid,
const Tag & tag,
const K & k);
const K & k,
const DiagonalizeTraits_& diagonalize_traits);
} /* namespace CGAL */

23
Principal_component_analysis/doc/Principal_component_analysis/CGAL/linear_least_squares_fitting_3.h
View File

@ -64,6 +64,13 @@ The class `K` is the kernel in which the
value type of `InputIterator` is defined. It can be omitted and deduced
automatically from the value type.
The class `DiagonalizeTraits_` is a model of `DiagonalizeTraits`. It
can be omitted: if Eigen 3 (or greater) is available and
`CGAL_EIGEN3_ENABLED` is defined then an overload using
`Eigen_diagonalize_traits` is provided. Otherwise, the internal
implementation `Diagonalize_traits` is used.
\cgalHeading{Requirements}
<OL>
@ -77,14 +84,15 @@ automatically from the value type.
*/
template < typename InputIterator, typename K, typename Tag >
template < typename InputIterator, typename K, typename Tag, typename DiagonalizeTraits_ >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
typename K::Line_3& line,
typename K::Point_3& centroid,
const Tag& tag,
const K& k);
const K& k,
const DiagonalizeTraits_& diagonalize_traits);
/*!
\brief computes the best fitting 3D plane of a 3D object set in the
@ -100,6 +108,12 @@ of `InputIterator` is defined. It can be omitted and deduced
automatically from the value type. The tag `tag` identifies the
dimension to be considered from the objects (see above).
The class `DiagonalizeTraits_` is a model of `DiagonalizeTraits`. It
can be omitted: if Eigen 3 (or greater) is available and
`CGAL_EIGEN3_ENABLED` is defined then an overload using
`Eigen_diagonalize_traits` is provided. Otherwise, the internal
implementation `Diagonalize_traits` is used.
\cgalHeading{Requirements}
<OL>
@ -112,14 +126,15 @@ dimension to be considered from the objects (see above).
</OL>
*/
template < typename InputIterator, typename K, typename Tag >
template < typename InputIterator, typename K, typename Tag, typename DiagonalizeTraits_ >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
typename K::Plane_3& plane,
typename K::Point_3& centroid,
const Tag& tag,
const K& k);
const K& k,
const DiagonalizeTraits_& diagonalize_traits);
/// @}

1
Principal_component_analysis/doc/Principal_component_analysis/PackageDescription.txt
View File

@ -27,6 +27,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.2}
\cgalPkgDependsOn{\ref PkgSolverSummary}
\cgalPkgBib{cgal:ap-pcad}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{Principal Component Analysis,pca.zip,Operations on Polygons,polygon.zip,Operations on Polyhedra,polyhedron_3.zip}

1
Principal_component_analysis/doc/Principal_component_analysis/dependencies
View File

@ -4,3 +4,4 @@ STL_Extension
Algebraic_foundations
Circulator
Stream_support
Solver_interface

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

@ -20,7 +20,6 @@
#ifndef CGAL_LINEAR_LEAST_SQUARES_FITTING_UTIL_H
#define CGAL_LINEAR_LEAST_SQUARES_FITTING_UTIL_H
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/Dimension.h>
@ -52,7 +51,7 @@ template < typename InputIterator,
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K& , // kernel
const typename K::Point_3*, // used for indirection
@ -63,9 +62,9 @@ assemble_covariance_matrix_3(InputIterator first,
typedef typename K::Vector_3 Vector;
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
covariance[0] = covariance[1] = covariance[2] =
covariance[3] = covariance[4] = covariance[5] = (FT)0.0;
for(InputIterator it = first;
@ -76,8 +75,8 @@ assemble_covariance_matrix_3(InputIterator first,
Vector d = p - c;
covariance[0] += d.x() * d.x();
covariance[1] += d.x() * d.y();
covariance[2] += d.y() * d.y();
covariance[3] += d.x() * d.z();
covariance[2] += d.x() * d.z();
covariance[3] += d.y() * d.y();
covariance[4] += d.y() * d.z();
covariance[5] += d.z() * d.z();
}
@ -88,8 +87,8 @@ template < typename InputIterator,
typename K >
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
InputIterator beyond,
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K&, // kernel
const typename K::Triangle_3*,// used for indirection
@ -102,9 +101,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
//Final combined covariance matrix for all triangles and their combined mass
FT mass = 0.0;
@ -141,8 +140,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0];
covariance[1] += transformation[1][0];
covariance[2] += transformation[1][1];
covariance[3] += transformation[2][0];
covariance[2] += transformation[2][0];
covariance[3] += transformation[1][1];
covariance[4] += transformation[2][1];
covariance[5] += transformation[2][2];
@ -153,8 +152,8 @@ assemble_covariance_matrix_3(InputIterator first,
// 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[3] += mass * (-1.0 * c.z() * c.x());
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());
@ -165,8 +164,8 @@ template < typename InputIterator,
typename K >
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
InputIterator beyond,
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K& , // kernel
const typename K::Iso_cuboid_3*,// used for indirection
@ -179,9 +178,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
// final combined covariance matrix for all cuboids and their combined mass
FT mass = (FT)0.0;
@ -227,8 +226,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0] + volume * (2*x0*xav0 + x0*x0);
covariance[1] += transformation[1][0] + volume * (xav0*y0 + yav0*x0 + x0*y0);
covariance[2] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
covariance[3] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
covariance[2] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
covariance[3] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
covariance[4] += transformation[2][1] + volume * (yav0*z0 + y0*zav0 + z0*y0);
covariance[5] += transformation[2][2] + volume * (2*zav0*z0 + z0*z0);
@ -239,8 +238,8 @@ assemble_covariance_matrix_3(InputIterator first,
// the center of mass to get the covariance.
covariance[0] += mass * (- c.x() * c.x());
covariance[1] += mass * (- c.x() * c.y());
covariance[2] += mass * (- c.y() * c.y());
covariance[3] += mass * (- c.z() * c.x());
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());
}
@ -251,7 +250,7 @@ template < typename InputIterator,
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K& , // kernel
const typename K::Iso_cuboid_3*,// used for indirection
@ -264,9 +263,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
//Final combined covariance matrix for all cuboids and their combined mass
FT mass = (FT)0.0;
@ -319,8 +318,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0] + area * (2*x0*xav0 + x0*x0);
covariance[1] += transformation[1][0] + area * (xav0*y0 + yav0*x0 + x0*y0);
covariance[2] += transformation[1][1] + area * (2*y0*yav0 + y0*y0);
covariance[3] += transformation[2][0] + area * (x0*zav0 + xav0*z0 + x0*z0);
covariance[2] += transformation[2][0] + area * (x0*zav0 + xav0*z0 + x0*z0);
covariance[3] += transformation[1][1] + area * (2*y0*yav0 + y0*y0);
covariance[4] += transformation[2][1] + area * (yav0*z0 + y0*zav0 + z0*y0);
covariance[5] += transformation[2][2] + area * (2*zav0*z0 + z0*z0);
@ -331,8 +330,8 @@ assemble_covariance_matrix_3(InputIterator first,
// 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[3] += mass * (-1.0 * c.z() * c.x());
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());
@ -344,7 +343,7 @@ template < typename InputIterator,
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K&, // kernel
const typename K::Sphere_3*, // used for indirection
@ -357,9 +356,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
//Final combined covariance matrix for all spheres and their combined mass
FT mass = 0.0;
@ -402,8 +401,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0] + volume * x0*x0;
covariance[1] += transformation[1][0] + volume * x0*y0;
covariance[2] += transformation[1][1] + volume * y0*y0;
covariance[3] += transformation[2][0] + volume * x0*z0;
covariance[2] += transformation[2][0] + volume * x0*z0;
covariance[3] += transformation[1][1] + volume * y0*y0;
covariance[4] += transformation[2][1] + volume * z0*y0;
covariance[5] += transformation[2][2] + volume * z0*z0;
@ -414,8 +413,8 @@ assemble_covariance_matrix_3(InputIterator first,
// 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[3] += mass * (-1.0 * c.z() * c.x());
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());
@ -426,7 +425,7 @@ template < typename InputIterator,
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K&, // kernel
const typename K::Sphere_3*, // used for indirection
@ -439,9 +438,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
//Final combined covariance matrix for all spheres and their combined mass
FT mass = 0.0;
@ -485,8 +484,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0] + area * x0*x0;
covariance[1] += transformation[1][0] + area * x0*y0;
covariance[2] += transformation[1][1] + area * y0*y0;
covariance[3] += transformation[2][0] + area * x0*z0;
covariance[2] += transformation[2][0] + area * x0*z0;
covariance[3] += transformation[1][1] + area * y0*y0;
covariance[4] += transformation[2][1] + area * z0*y0;
covariance[5] += transformation[2][2] + area * z0*z0;
@ -497,8 +496,8 @@ assemble_covariance_matrix_3(InputIterator first,
// 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[3] += mass * (-1.0 * c.z() * c.x());
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());
@ -510,7 +509,7 @@ template < typename InputIterator,
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K& , // kernel
const typename K::Tetrahedron_3*,// used for indirection
@ -523,9 +522,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
//Final combined covariance matrix for all tetrahedrons and their combined mass
FT mass = 0.0;
@ -571,8 +570,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0] + volume * (2*x0*xav0 + x0*x0);
covariance[1] += transformation[1][0] + volume * (xav0*y0 + yav0*x0 + x0*y0);
covariance[2] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
covariance[3] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
covariance[2] += transformation[2][0] + volume * (x0*zav0 + xav0*z0 + x0*z0);
covariance[3] += transformation[1][1] + volume * (2*y0*yav0 + y0*y0);
covariance[4] += transformation[2][1] + volume * (yav0*z0 + y0*zav0 + z0*y0);
covariance[5] += transformation[2][2] + volume * (2*zav0*z0 + z0*z0);
@ -583,8 +582,8 @@ assemble_covariance_matrix_3(InputIterator first,
// 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[3] += mass * (-1.0 * c.z() * c.x());
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());
}
@ -595,7 +594,7 @@ template < typename InputIterator,
void
assemble_covariance_matrix_3(InputIterator first,
InputIterator beyond,
typename K::FT covariance[6], // covariance matrix
CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
const K& , // kernel
const typename K::Segment_3*,// used for indirection
@ -608,9 +607,9 @@ assemble_covariance_matrix_3(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 3 4 5
// 0 1 2
// 3 4
// 5
//Final combined covariance matrix for all segments and their combined mass
FT mass = 0.0;
@ -648,8 +647,8 @@ assemble_covariance_matrix_3(InputIterator first,
// and add to covariance matrix
covariance[0] += transformation[0][0];
covariance[1] += transformation[1][0];
covariance[2] += transformation[1][1];
covariance[3] += transformation[2][0];
covariance[2] += transformation[2][0];
covariance[3] += transformation[1][1];
covariance[4] += transformation[2][1];
covariance[5] += transformation[2][2];
@ -660,8 +659,8 @@ assemble_covariance_matrix_3(InputIterator first,
// 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[3] += mass * (-1.0 * c.z() * c.x());
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());
@ -671,23 +670,27 @@ assemble_covariance_matrix_3(InputIterator first,
// compute the eigen values and vectors of the covariance
// matrix and deduces the best linear fitting plane.
// returns fitting quality
template < typename K >
template < typename K, typename DiagonalizeTraits >
typename K::FT
fitting_plane_3(const typename K::FT covariance[6], // covariance matrix
fitting_plane_3(CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
typename K::Plane_3& plane, // best fit plane
const K& ) // kernel
const K&, // kernel
const DiagonalizeTraits& ) // Diagonalize traits
{
typedef typename K::FT FT;
typedef typename K::Plane_3 Plane;
typedef typename K::Vector_3 Vector;
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
FT eigen_values[3];
FT eigen_vectors[9];
eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
CGAL::cpp11::array<FT, 3> eigen_values = {{ 0. , 0., 0. }};
CGAL::cpp11::array<FT, 9> eigen_vectors = {{ 0., 0., 0.,
0., 0., 0.,
0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// degenerate case
if(eigen_values[0] == eigen_values[1] &&
@ -700,11 +703,11 @@ fitting_plane_3(const typename K::FT covariance[6], // covariance matrix
}
else // regular and line case
{
Vector normal(eigen_vectors[6],
eigen_vectors[7],
eigen_vectors[8]);
Vector normal(eigen_vectors[0],
eigen_vectors[1],
eigen_vectors[2]);
plane = Plane(c,normal);
return FT(1) - eigen_values[2] / eigen_values[1];
return FT(1) - eigen_values[0] / eigen_values[1];
} // end regular case
}
@ -712,23 +715,27 @@ fitting_plane_3(const typename K::FT covariance[6], // covariance matrix
// matrix and deduces the best linear fitting line
// (this is an internal function)
// returns fitting quality
template < typename K >
template < typename K, typename DiagonalizeTraits >
typename K::FT
fitting_line_3(const typename K::FT covariance[6], // covariance matrix
fitting_line_3(CGAL::cpp11::array<typename K::FT, 6>& covariance, // covariance matrix
const typename K::Point_3& c, // centroid
typename K::Line_3& line, // best fit line
const K&) // kernel
const K&, // kernel
const DiagonalizeTraits& ) // Diagonalize traits
{
typedef typename K::FT FT;
typedef typename K::Line_3 Line;
typedef typename K::Vector_3 Vector;
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
FT eigen_values[3];
FT eigen_vectors[9];
eigen_symmetric<FT>(covariance,3,eigen_vectors,eigen_values);
CGAL::cpp11::array<FT, 3> eigen_values = {{ 0. , 0., 0. }};
CGAL::cpp11::array<FT, 9> eigen_vectors = {{ 0., 0., 0.,
0., 0., 0.,
0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// isotropic case (infinite number of directions)
if(eigen_values[0] == eigen_values[1] &&
@ -742,9 +749,9 @@ fitting_line_3(const typename K::FT covariance[6], // covariance matrix
else
{
// regular case
Vector direction(eigen_vectors[0],eigen_vectors[1],eigen_vectors[2]);
Vector direction(eigen_vectors[6],eigen_vectors[7],eigen_vectors[8]);
line = Line(c,direction);
return (FT)1.0 - eigen_values[1] / eigen_values[0];
return (FT)1.0 - eigen_values[1] / eigen_values[2];
}
}

229
Principal_component_analysis/include/CGAL/eigen.h
View File

@ -1,229 +0,0 @@
// Copyright (c) 2005 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Bruno Levy, Pierre Alliez
#ifndef CGAL_EIGEN_H
#define CGAL_EIGEN_H
#include <cmath>
#include <CGAL/number_utils.h>
namespace CGAL {
namespace internal {
template <class FT>
void eigen_symmetric(const FT *mat,
const int n,
FT *eigen_vectors,
FT *eigen_values,
const int MAX_ITER = 100)
{
static const FT EPSILON = (FT)0.00001;
// number of entries in mat
int nn = (n*(n+1))/2;
// copy matrix
FT *a = new FT[nn];
int ij;
for(ij=0; ij<nn; ij++)
a[ij] = mat[ij];
// Fortran-porting
a--;
// init diagonalization matrix as the unit matrix
FT *v = new FT[n*n];
ij = 0;
int i;
for(i=0; i<n; i++)
for(int j=0; j<n; j++)
if(i==j)
v[ij++] = 1.0;
else
v[ij++] = 0.0;
// Fortran-porting
v--;
// compute weight of the non diagonal terms
ij = 1;
FT a_norm = 0.0;
for(i=1; i<=n; i++)
for(int j=1; j<=i; j++)
{
if( i!=j )
{
FT a_ij = a[ij];
a_norm += a_ij * a_ij;
}
ij++;
}
if(a_norm != 0.0)
{
FT a_normEPS = a_norm * EPSILON;
FT thr = a_norm;
// rotations
int nb_iter = 0;
while(thr > a_normEPS && nb_iter < MAX_ITER)
{
nb_iter++;
FT thr_nn = thr / nn;
for(int l=1; l< n; l++)
{
for(int m=l+1; m<=n; m++)
{
// compute sinx and cosx
int lq = (l*l-l)/2;
int mq = (m*m-m)/2;
int lm = l + mq;
FT a_lm = a[lm];
FT a_lm_2 = a_lm * a_lm;
if(a_lm_2 < thr_nn)
continue;
int ll = l + lq;
int mm = m + mq;
FT a_ll = a[ll];
FT a_mm = a[mm];
FT delta = a_ll - a_mm;
FT x;
if(delta == 0.0)
x = (FT) - CGAL_PI / 4;
else
x = (FT)(- std::atan( (a_lm+a_lm) / delta ) / 2.0);
FT sinx = std::sin(x);
FT cosx = std::cos(x);
FT sinx_2 = sinx * sinx;
FT cosx_2 = cosx * cosx;
FT sincos = sinx * cosx;
// rotate L and M columns
int ilv = n*(l-1);
int imv = n*(m-1);
int i;
for( i=1; i<=n;i++ )
{
if( (i!=l) && (i!=m) )
{
int iq = (i*i-i)/2;
int im;
if( i<m )
im = i + mq;
else
im = m + iq;
FT a_im = a[im];
int il;
if( i<l )
il = i + lq;
else
il = l + iq;
FT a_il = a[il];
a[il] = a_il * cosx - a_im * sinx;
a[im] = a_il * sinx + a_im * cosx;
}
ilv++;
imv++;
FT v_ilv = v[ilv];
FT v_imv = v[imv];
v[ilv] = cosx * v_ilv - sinx * v_imv;
v[imv] = sinx * v_ilv + cosx * v_imv;
}
x = a_lm * sincos;
x += x;
a[ll] = a_ll * cosx_2 + a_mm * sinx_2 - x;
a[mm] = a_ll * sinx_2 + a_mm * cosx_2 + x;
a[lm] = 0.0;
thr = CGAL::abs(thr - a_lm_2);
}
}
}
}
// convert indices and copy eigen values
a++;
for(i=0; i<n; i++)
{
int k = i + (i*(i+1))/2;
eigen_values[i] = a[k];
}
delete [] a;
// sort eigen values and vectors
int *index = new int[n];
for(i=0; i<n; i++)
index[i] = i;
for(i=0; i<(n-1); i++)
{
FT x = eigen_values[i];
int k = i;
for(int j=i+1; j<n; j++)
if(x < eigen_values[j])
{
k = j;
x = eigen_values[j];
}
eigen_values[k] = eigen_values[i];
eigen_values[i] = x;
int jj = index[k];
index[k] = index[i];
index[i] = jj;
}
// save eigen vectors
v++; // back to C++
ij = 0;
for(int k=0; k<n; k++ )
{
int ik = index[k]*n;
for(int i=0; i<n; i++)
eigen_vectors[ij++] = v[ik++];
}
delete [] v;
delete [] index;
}
} // end namespace internal
} //namespace CGAL
#endif // CGAL_EIGEN_H

121
Principal_component_analysis/include/CGAL/eigen_2.h
View File

@ -1,121 +0,0 @@
// Copyright (c) 2005 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
//
// Author(s) : Pierre Alliez
#ifndef CGAL_EIGEN_2_H
#define CGAL_EIGEN_2_H
#include <cmath>
#include <utility>
#include <CGAL/assertions.h>
namespace CGAL {
namespace internal {
// extract eigenvalues and eigenvectors from a 2x2 symmetric
// positive definite matrix.
// Note: computations involve a square root.
// Matrix numbering:
// a b
// b c
// Eigen values and vectors are sorted in descendent order.
template <typename K>
void eigen_symmetric_2(const typename K::FT *matrix, // a b c
std::pair<typename K::Vector_2,
typename K::Vector_2>& eigen_vectors,
std::pair<typename K::FT,
typename K::FT>& eigen_values)
{
// types
typedef typename K::FT FT;
typedef typename K::Vector_2 Vector;
// for better reading
FT a = matrix[0];
FT b = matrix[1];
FT c = matrix[2];
FT p = c*c - 2*a*c + 4*b*b + a*a;
CGAL_assertion(a >= 0.0 && c >= 0.0);
// degenerate or isotropic case
if(p == 0.0)
{
// unit eigen values by default
eigen_values.first = eigen_values.second = (FT)1.0;
// any vector is eigen vector
// the 2D canonical frame is output by default
eigen_vectors.first = Vector((FT)1.0,(FT)0.0);
eigen_vectors.second = Vector((FT)0.0,(FT)1.0);
}
else
{
if(b == 0.0)
{
if(a>=c)
{
eigen_values.first = a;
eigen_values.second = c;
eigen_vectors.first = Vector((FT)1.0, (FT)0.0);
eigen_vectors.second = Vector((FT)0.0, (FT)1.0);
}
else
{
eigen_values.first = c;
eigen_values.second = a;
eigen_vectors.first = Vector((FT)0.0, (FT)1.0);
eigen_vectors.second = Vector((FT)1.0, (FT)0.0);
}
}
else // generic case
{
FT l1 = (FT)(0.5 * ( -1*CGAL::sqrt(p) + c + a));
FT l2 = (FT)(0.5 * ( CGAL::sqrt(p) + c + a));
// all eigen values of a symmetric positive
// definite matrix must be real and positive
// we saturate the values if this is not the
// case for floating point computations.
l1 = (l1 < (FT)0.0) ? (FT)0.0 : l1;
l2 = (l2 < (FT)0.0) ? (FT)0.0 : l2;
// sort eigen values and vectors in descendent order.
if(l1 >= l2)
{
eigen_values.first = l1;
eigen_values.second = l2;
eigen_vectors.first = Vector((FT)1.0, (FT)(-(CGAL::sqrt(p)-c+a) / (2*b)));
eigen_vectors.second = Vector((FT)1.0, (FT)( (CGAL::sqrt(p)+c-a) / (2*b)));
}
else
{
eigen_values.first = l2;
eigen_values.second = l1;
eigen_vectors.first = Vector((FT)1.0, (FT)( (CGAL::sqrt(p)+c-a) / (2*b)));
eigen_vectors.second = Vector((FT)1.0, (FT)(-(CGAL::sqrt(p)-c+a) / (2*b)));
}
} // end generic case
} // end non-degenerate case
} // end eigen_symmetric_2
} // end namespace internal
} //namespace CGAL
#endif // CGAL_EIGEN_2_H

26
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_2.h
View File

@ -29,14 +29,15 @@
#include <CGAL/linear_least_squares_fitting_circles_2.h>
#include <CGAL/linear_least_squares_fitting_rectangles_2.h>
#include <CGAL/Dimension.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <iterator>
namespace CGAL {
template < typename InputIterator,
typename Kernel,
typename Tag>
typename Tag,
typename DiagonalizeTraits>
inline
typename Kernel::FT
linear_least_squares_fitting_2(InputIterator first,
@ -44,35 +45,39 @@ linear_least_squares_fitting_2(InputIterator first,
typename Kernel::Line_2& line,
typename Kernel::Point_2& centroid,
const Tag& tag,
const Kernel& kernel)
const Kernel& kernel,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
return internal::linear_least_squares_fitting_2(first, beyond, line,
centroid,(Value_type*)NULL,kernel,tag);
centroid,(Value_type*)NULL,kernel,tag,
diagonalize_traits);
}
// deduces the kernel from the points in container.
// Use default DiagonalizeTraits
template < typename InputIterator,
typename Line,
typename Point,
typename Tag>
typename Tag>
inline
typename Kernel_traits<Line>::Kernel::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
Line& line,
Point& centroid,
const Tag& tag)
const Tag& tag)
{
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
return CGAL::linear_least_squares_fitting_2(first,beyond,line,centroid,tag,Kernel());
return CGAL::linear_least_squares_fitting_2
(first,beyond,line,centroid,tag,Kernel(),
Default_diagonalize_traits<typename Kernel::FT, 2>());
}
template < typename InputIterator,
typename Line,
typename Tag >
typename Tag>
inline
typename Kernel_traits<Line>::Kernel::FT
linear_least_squares_fitting_2(InputIterator first,
@ -83,7 +88,8 @@ linear_least_squares_fitting_2(InputIterator first,
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
typename Kernel::Point_2 centroid; // unused
return CGAL::linear_least_squares_fitting_2(first,beyond,line,centroid,tag,Kernel());
return CGAL::linear_least_squares_fitting_2(first,beyond,line,centroid,tag,Kernel(),
Default_diagonalize_traits<typename Kernel::FT, 2>());
}

31
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_3.h
View File

@ -21,8 +21,6 @@
#define CGAL_LINEAR_LEAST_SQUARES_FITTING_3_H
#include <CGAL/basic.h>
//#include <CGAL/Algebraic_structure_traits.h>
//#include <CGAL/IO/io.h>
#include <CGAL/linear_least_squares_fitting_points_3.h>
#include <CGAL/linear_least_squares_fitting_segments_3.h>
@ -31,6 +29,8 @@
#include <CGAL/linear_least_squares_fitting_tetrahedra_3.h>
#include <CGAL/linear_least_squares_fitting_spheres_3.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <CGAL/Dimension.h>
#include <iterator>
@ -42,7 +42,8 @@ namespace CGAL {
template < typename InputIterator,
typename Object,
typename Kernel,
typename Tag >
typename Tag,
typename DiagonalizeTraits >
inline
typename Kernel::FT
linear_least_squares_fitting_3(InputIterator first,
@ -50,32 +51,38 @@ linear_least_squares_fitting_3(InputIterator first,
Object& object, // plane or line
typename Kernel::Point_3& centroid,
const Tag& tag, // dimension tag, ranges from 0 to 3
const Kernel& kernel)
const Kernel& kernel,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
return internal::linear_least_squares_fitting_3(first, beyond, object,
centroid, (Value_type*) NULL, kernel, tag);
centroid, (Value_type*) NULL, kernel, tag,
diagonalize_traits);
}
// deduces kernel from value type of input iterator
// use default DiagonalizeTraits
template < typename InputIterator,
typename Object,
typename Point,
typename Tag>
typename Point,
typename Tag >
inline
typename Kernel_traits<Object>::Kernel::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
Object& object, // plane or line
Point& centroid,
const Tag& tag) // dimension tag, ranges from 0 to 3
const Tag& tag) // dimension tag, ranges from 0 to 3
{
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,tag,Kernel());
return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,tag,Kernel(),
Default_diagonalize_traits<typename Kernel::FT, 3>());
}
// deduces kernel and does not write centroid
// use default DiagonalizeTraits
template < typename InputIterator,
typename Object,
typename Tag>
@ -84,12 +91,14 @@ typename Kernel_traits<Object>::Kernel::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
Object& object, // plane or line
const Tag& tag) // dimension tag, ranges from 0 to 3
const Tag& tag) // dimension tag, ranges from 0 to 3
{
typedef typename std::iterator_traits<InputIterator>::value_type Value_type;
typedef typename Kernel_traits<Value_type>::Kernel Kernel;
typename Kernel::Point_3 centroid; // not used by caller
return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,tag);
return CGAL::linear_least_squares_fitting_3(first,beyond,object,centroid,tag,Kernel(),
Default_diagonalize_traits<typename Kernel::FT, 3>());
}
} //namespace CGAL

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

@ -22,8 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/PCA_util.h>
@ -40,7 +38,7 @@ namespace internal {
// 0 is worst (isotropic case, returns a line with horizontal
// direction by default)
template < typename InputIterator, typename K >
template < typename InputIterator, typename K, typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -48,7 +46,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Circle_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits&)
{
// types
typedef typename K::FT FT;
@ -67,11 +66,11 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 0 1
// 2
//Final combined covariance matrix for all circles and their combined mass
FT mass = 0.0;
FT covariance[3] = {0.0,0.0,0.0};
CGAL::cpp11::array<FT, 3> covariance = {{ 0., 0., 0. }};
// assemble 2nd order moment about the origin.
FT temp[4] = {0.25, 0.0,
@ -120,22 +119,19 @@ linear_least_squares_fitting_2(InputIterator first,
covariance[2] += mass * (-1.0 * c.y() * c.y());
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
// internal::eigen_symmetric_2<K>(final_cov, eigen_vectors, eigen_values);
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(covariance,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
CGAL::cpp11::array<FT, 2> eigen_values = {{ 0. , 0. }};
CGAL::cpp11::array<FT, 4> eigen_vectors = {{ 0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
if(eigen_values[0] != eigen_values[1])
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
line = Line(c, Vector (eigen_vectors[2],eigen_vectors[3]));
return (FT)1.0 - eigen_values[0] / eigen_values[1];
}
else
{
@ -148,7 +144,7 @@ linear_least_squares_fitting_2(InputIterator first,
} // end linear_least_squares_fitting_2 for circle set with 2D tag
template < typename InputIterator, typename K >
template < typename InputIterator, typename K, typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -156,7 +152,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Circle_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& )
{
// types
typedef typename K::FT FT;
@ -174,11 +171,11 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
// 0 1
// 2
//Final combined covariance matrix for all circles and their combined mass
FT mass = 0.0;
FT covariance[3] = {0.0,0.0,0.0};
CGAL::cpp11::array<FT, 3> covariance = {{ 0., 0., 0. }};
// assemble 2nd order moment about the origin.
FT temp[4] = {1.0, 0.0,
@ -226,22 +223,19 @@ linear_least_squares_fitting_2(InputIterator first,
covariance[2] += mass * (-1.0 * c.y() * c.y());
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
// internal::eigen_symmetric_2<K>(final_cov, eigen_vectors, eigen_values);
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(covariance,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
CGAL::cpp11::array<FT, 2> eigen_values = {{ 0. , 0. }};
CGAL::cpp11::array<FT, 4> eigen_vectors = {{ 0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
if(eigen_values[1] != eigen_values[0])
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
line = Line(c, Vector(eigen_vectors[2],eigen_vectors[3]));
return (FT)1.0 - eigen_values[0] / eigen_values[1];
}
else
{

77
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
View File

@ -22,7 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen.h>
#include <CGAL/PCA_util.h>
#include <CGAL/linear_least_squares_fitting_points_3.h>
#include <CGAL/linear_least_squares_fitting_segments_3.h>
@ -36,7 +35,8 @@ namespace internal {
// fits a plane to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -44,7 +44,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<3>& tag)
const CGAL::Dimension_tag<3>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -56,18 +57,19 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,k,tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a plane to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -75,7 +77,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -87,18 +90,19 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,k,tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a plane to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -106,7 +110,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Segment_3 Segment;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -135,13 +140,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)NULL,k,tag);
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a plane to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -149,7 +156,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Point_3 Point;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -174,13 +182,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a line to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -188,7 +198,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<3>& tag)
const CGAL::Dimension_tag<3>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -200,17 +211,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,k,tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a line to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -218,7 +230,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -230,18 +243,19 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,k,tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Iso_cuboid*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a line to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -249,7 +263,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Segment_3 Segment;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -278,13 +293,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)NULL,k,tag);
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3
// fits a line to a 3D cuboid set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -292,7 +309,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Iso_cuboid_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Point_3 Point;
typedef typename K::Iso_cuboid_3 Iso_cuboid;
@ -317,7 +335,8 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_cuboids_3

31
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_2.h
View File

@ -22,8 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <iterator>
#include <cmath>
@ -37,7 +35,8 @@ namespace internal {
// 0 is worst (isotropic case, returns a line with horizontal
// direction by default).
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits>
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -45,7 +44,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Point_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits&)
{
// types
typedef typename K::FT FT;
@ -63,9 +63,11 @@ linear_least_squares_fitting_2(InputIterator first,
// assemble covariance matrix as a semi-definite matrix.
// Matrix numbering:
// 0
// 1 2
FT covariance[3] = {0,0,0};
// 0 1
// 2
CGAL::cpp11::array<FT, 3> covariance = {{ 0., 0., 0. }};
for(InputIterator it = first;
it != beyond;
it++)
@ -78,18 +80,19 @@ linear_least_squares_fitting_2(InputIterator first,
}
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
internal::eigen_symmetric_2<K>(covariance, eigen_vectors, eigen_values);
CGAL::cpp11::array<FT, 2> eigen_values = {{ 0. , 0. }};
CGAL::cpp11::array<FT, 4> eigen_vectors = {{ 0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
if(eigen_values[0] != eigen_values[1])
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
line = Line(c, Vector (eigen_vectors[2], eigen_vectors[3]));
return (FT)1.0 - eigen_values[0] / eigen_values[1];
}
else
{

21
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_points_3.h
View File

@ -22,7 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen.h>
#include <CGAL/PCA_util.h>
#include <iterator>
@ -35,7 +34,8 @@ namespace internal {
// 1 is best (zero variance orthogonally to the fitting line)
// 0 is worst (isotropic case, returns a plane with default direction)
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -43,7 +43,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Point_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Point_3 Point;
@ -55,11 +56,11 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6];
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Point*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end fit plane to point set
// fits a line to a 3D point set
@ -67,7 +68,8 @@ linear_least_squares_fitting_3(InputIterator first,
// 1 is best (zero variance orthogonally to the fitting line)
// 0 is worst (isotropic case, returns a line along x axis)
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -75,7 +77,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Point_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Point_3 Point;
@ -87,11 +90,11 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6];
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Point*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end fit line to point set
} // end namespace internal

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

@ -22,8 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/PCA_util.h>
@ -40,7 +38,7 @@ namespace internal {
// 0 is worst (isotropic case, returns a line with horizontal
// direction by default)
template < typename InputIterator, typename K >
template < typename InputIterator, typename K, typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -48,7 +46,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Iso_rectangle_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits&)
{
// types
typedef typename K::FT FT;
@ -70,7 +69,7 @@ linear_least_squares_fitting_2(InputIterator first,
// 1 2
//Final combined covariance matrix for all rectangles and their combined mass
FT mass = 0.0;
FT covariance[3] = {0.0,0.0,0.0};
CGAL::cpp11::array<FT, 3> covariance = {{ 0., 0., 0. }};
// assemble 2nd order moment about the origin.
FT temp[4] = {1/3.0, 0.25,
@ -123,22 +122,19 @@ linear_least_squares_fitting_2(InputIterator first,
covariance[2] += mass * (-1.0 * c.y() * c.y());
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(covariance,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
CGAL::cpp11::array<FT, 2> eigen_values = {{ 0. , 0. }};
CGAL::cpp11::array<FT, 4> eigen_vectors = {{ 0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
if(eigen_values[0] != eigen_values[1])
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
line = Line(c, Vector(eigen_vectors[2],eigen_vectors[3]));
return (FT)1.0 - eigen_values[0] / eigen_values[1];
}
else
{
@ -150,7 +146,7 @@ linear_least_squares_fitting_2(InputIterator first,
}
} // end linear_least_squares_fitting_2 for rectangle set with 2D tag
template < typename InputIterator, typename K >
template < typename InputIterator, typename K, typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -158,7 +154,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Iso_rectangle_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
// types
typedef typename K::Iso_rectangle_2 Iso_rectangle;
@ -179,13 +176,15 @@ linear_least_squares_fitting_2(InputIterator first,
segments.push_back(Segment_2(t[3],t[0]));
}
return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K(),tag);
return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,K(),tag,
diagonalize_traits);
} // end linear_least_squares_fitting_2 for rectangle set with 1D tag
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -193,7 +192,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Iso_rectangle_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
// types
typedef typename K::Iso_rectangle_2 Iso_rectangle;
@ -214,7 +214,8 @@ linear_least_squares_fitting_2(InputIterator first,
points.push_back(Point_2(t[3]));
}
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,K(),tag);
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,K(),tag,
diagonalize_traits);
} // end linear_least_squares_fitting_2 for rectangle set with 0D tag

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

@ -22,8 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/PCA_util.h>
@ -40,7 +38,7 @@ namespace internal {
// 0 is worst (isotropic case, returns a line with horizontal
// direction by default)
template < typename InputIterator, typename K >
template < typename InputIterator, typename K, typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -48,7 +46,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Segment_2*,// used for indirection
const K&, // kernel
const CGAL::Dimension_tag<1>& tag = CGAL::Dimension_tag<1>())
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits&)
{
// types
typedef typename K::FT FT;
@ -69,7 +68,7 @@ linear_least_squares_fitting_2(InputIterator first,
// 1 2
//Final combined covariance matrix for all segments and their combined mass
FT mass = 0.0;
FT covariance[3] = {0.0,0.0,0.0};
CGAL::cpp11::array<FT, 3> covariance = {{ 0., 0., 0. }};
// assemble 2nd order moment about the origin.
FT temp[4] = {1.0, 0.5, 0.5, 1.0};
@ -111,22 +110,19 @@ linear_least_squares_fitting_2(InputIterator first,
covariance[2] += mass * (-1.0 * c.y() * c.y());
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
// internal::eigen_symmetric_2<K>(covariance, eigen_vectors, eigen_values);
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(covariance,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
CGAL::cpp11::array<FT, 2> eigen_values = {{ 0. , 0. }};
CGAL::cpp11::array<FT, 4> eigen_vectors = {{ 0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
if(eigen_values[0] != eigen_values[1])
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
line = Line(c, Vector(eigen_vectors[2],eigen_vectors[3]));
return (FT)1.0 - eigen_values[0] / eigen_values[1];
}
else
{
@ -138,7 +134,7 @@ linear_least_squares_fitting_2(InputIterator first,
}
} // end linear_least_squares_fitting_2 for segment set with 1D tag
template < typename InputIterator, typename K >
template < typename InputIterator, typename K, typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -146,7 +142,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename K::Point_2& c, // centroid
const typename K::Segment_2*,// used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
// types
typedef typename K::Point_2 Point;
@ -164,7 +161,8 @@ linear_least_squares_fitting_2(InputIterator first,
points.push_back(s[0]);
points.push_back(s[1]);
}
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,k,(Point*)NULL,tag);
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,k,(Point*)NULL,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_2 for segment set with 1D tag

39
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
View File

@ -22,7 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen.h>
#include <CGAL/PCA_util.h>
#include <CGAL/linear_least_squares_fitting_points_3.h>
@ -35,7 +34,8 @@ namespace internal {
// fits a plane to a 3D segment set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -43,7 +43,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Segment_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Segment_3 Segment;
@ -55,17 +56,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,k,tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Segment*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_segments_3
// fits a plane to a 3D segment set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -73,7 +75,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Segment_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Segment_3 Segment;
typedef typename K::Point_3 Point;
@ -92,13 +95,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_segments_3
// fits a line to a 3D segment set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -106,7 +111,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Segment_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Segment_3 Segment;
@ -118,17 +124,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,k,tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Segment*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_segments_3
// fits a plane to a 3D segment set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -136,7 +143,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Segment_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Segment_3 Segment;
typedef typename K::Point_3 Point;
@ -155,7 +163,8 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_segments_3

41
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_spheres_3.h
View File

@ -22,7 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen.h>
#include <CGAL/PCA_util.h>
#include <iterator>
@ -33,7 +32,8 @@ namespace internal {
// fits a plane to a set of 3D balls (3D)
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -41,7 +41,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Sphere_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<3>& tag)
const CGAL::Dimension_tag<3>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Sphere_3 Sphere;
@ -53,17 +54,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_spheres_3
// fits a plane to a 3D sphere set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -71,7 +73,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Sphere_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Sphere_3 Sphere;
@ -83,18 +86,19 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_spheres_3
// fits a line to a 3D ball set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -102,7 +106,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Sphere_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<3>& tag)
const CGAL::Dimension_tag<3>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Sphere_3 Sphere;
@ -114,18 +119,19 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_spheres_3
// fits a line to a 3D sphere set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -133,7 +139,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Sphere_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Sphere_3 Sphere;
@ -145,11 +152,11 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Sphere*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_spheres_3

75
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
View File

@ -22,7 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen.h>
#include <CGAL/PCA_util.h>
#include <CGAL/linear_least_squares_fitting_points_3.h>
@ -37,7 +36,8 @@ namespace internal {
// fits a plane to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -45,7 +45,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<3>& tag)
const CGAL::Dimension_tag<3>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Tetrahedron_3 Tetrahedron;
@ -57,16 +58,17 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Tetrahedron*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a plane to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -74,7 +76,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Tetrahedron_3 Tetrahedron;
typedef typename K::Triangle_3 Triangle;
@ -95,13 +98,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,c,(Triangle*)NULL,k,tag);
return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,c,(Triangle*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a plane to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -109,7 +114,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Tetrahedron_3 Tetrahedron;
typedef typename K::Segment_3 Segment;
@ -133,13 +139,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)NULL,k,tag);
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a plane to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -147,7 +155,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Tetrahedron_3 Tetrahedron;
typedef typename K::Point_3 Point;
@ -168,13 +177,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a line to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -182,7 +193,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<3>& tag)
const CGAL::Dimension_tag<3>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Tetrahedron_3 Tetrahedron;
@ -194,17 +206,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Tetrahedron*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a line to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -212,7 +225,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Tetrahedron_3 Tetrahedron;
typedef typename K::Triangle_3 Triangle;
@ -233,13 +247,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,c,(Triangle*)NULL,k,tag);
return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,c,(Triangle*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a line to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -247,7 +263,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Tetrahedron_3 Tetrahedron;
typedef typename K::Segment_3 Segment;
@ -270,13 +287,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)NULL,k,tag);
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedrons_3
// fits a line to a 3D tetrahedron set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -284,7 +303,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Tetrahedron_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Tetrahedron_3 Tetrahedron;
typedef typename K::Point_3 Point;
@ -305,7 +325,8 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_tetrahedra_3

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

@ -22,8 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen_2.h>
#include <CGAL/eigen.h>
#include <CGAL/Linear_algebraCd.h>
#include <CGAL/PCA_util.h>
@ -41,7 +39,7 @@ namespace internal {
// direction by default)
template < typename InputIterator,
typename Kernel >
typename Kernel, typename DiagonalizeTraits >
typename Kernel::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -49,7 +47,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename Kernel::Point_2& c, // centroid
const typename Kernel::Triangle_2*,// used for indirection
const Kernel&, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits&)
{
// types
typedef typename Kernel::FT FT;
@ -71,7 +70,7 @@ linear_least_squares_fitting_2(InputIterator first,
// 1 2
//Final combined covariance matrix for all triangles and their combined mass
FT mass = 0.0;
FT covariance[3] = {0.0, 0.0, 0.0};
CGAL::cpp11::array<FT, 3> covariance = {{ 0., 0., 0. }};
// assemble the 2nd order moment about the origin.
FT temp[4] = {1/12.0, 1/24.0,
@ -123,21 +122,19 @@ linear_least_squares_fitting_2(InputIterator first,
// std::cout<<"cov: "<<covariance[0]*covariance[2]<<" =? "<<covariance[1]*covariance[1]<<std::endl;
// solve for eigenvalues and eigenvectors.
// eigen values are sorted in descending order,
// eigen values are sorted in ascending order,
// eigen vectors are sorted in accordance.
std::pair<FT,FT> eigen_values;
std::pair<Vector,Vector> eigen_vectors;
FT eigen_vectors1[4];
FT eigen_values1[2];
eigen_symmetric<FT>(covariance,2, eigen_vectors1, eigen_values1);
eigen_values = std::make_pair(eigen_values1[0],eigen_values1[1]);
eigen_vectors = std::make_pair(Vector(eigen_vectors1[0],eigen_vectors1[1]),Vector(eigen_vectors1[2],eigen_vectors1[3]));
CGAL::cpp11::array<FT, 2> eigen_values = {{ 0. , 0. }};
CGAL::cpp11::array<FT, 4> eigen_vectors = {{ 0., 0., 0. }};
DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix
(covariance, eigen_values, eigen_vectors);
// check unicity and build fitting line accordingly
if(eigen_values.first != eigen_values.second)
if(eigen_values[0] != eigen_values[1])
{
// regular case
line = Line(c, eigen_vectors.first);
return (FT)1.0 - eigen_values.second / eigen_values.first;
line = Line(c, Vector(eigen_vectors[2],eigen_vectors[3]));
return (FT)1.0 - eigen_values[0] / eigen_values[1];
}
else
{
@ -150,7 +147,8 @@ linear_least_squares_fitting_2(InputIterator first,
} // end linear_least_squares_fitting_2 for triangle set with 2D tag
template < typename InputIterator,
typename Kernel >
typename Kernel,
typename DiagonalizeTraits >
typename Kernel::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -158,7 +156,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename Kernel::Point_2& c, // centroid
const typename Kernel::Triangle_2*,// used for indirection
const Kernel&, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
// types
typedef typename Kernel::Triangle_2 Triangle;
@ -178,12 +177,14 @@ linear_least_squares_fitting_2(InputIterator first,
segments.push_back(Segment(t[2],t[0]));
}
return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,tag,Kernel());
return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,tag,Kernel(),
diagonalize_traits);
} // end linear_least_squares_fitting_2 for triangle set with 1D tag
template < typename InputIterator,
typename Kernel >
typename Kernel,
typename DiagonalizeTraits >
typename Kernel::FT
linear_least_squares_fitting_2(InputIterator first,
InputIterator beyond,
@ -191,7 +192,8 @@ linear_least_squares_fitting_2(InputIterator first,
typename Kernel::Point_2& c, // centroid
const typename Kernel::Triangle_2*,// used for indirection
const Kernel&, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
// types
@ -212,7 +214,8 @@ linear_least_squares_fitting_2(InputIterator first,
points.push_back(Point(t[2]));
}
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,tag,Kernel());
return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,tag,Kernel(),
diagonalize_traits);
} // end linear_least_squares_fitting_2 for triangle set with 0D tag

57
Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
View File

@ -22,7 +22,6 @@
#include <CGAL/basic.h>
#include <CGAL/centroid.h>
#include <CGAL/eigen.h>
#include <CGAL/PCA_util.h>
#include <list>
@ -34,7 +33,8 @@ namespace internal {
// fits a plane to a 3D triangle set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -42,7 +42,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Triangle_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Triangle_3 Triangle;
@ -54,17 +55,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Triangle*) NULL,tag);
// compute fitting plane
return fitting_plane_3(covariance,c,plane,k);
return fitting_plane_3(covariance,c,plane,k,diagonalize_traits);
} // end linear_least_squares_fitting_triangles_3
// fits a plane to a 3D triangle set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -72,7 +74,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Triangle_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Triangle_3 Triangle;
typedef typename K::Segment_3 Segment;
@ -92,13 +95,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)NULL,k,tag);
return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_triangles_3
// fits a plane to a 3D triangle set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -106,7 +111,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Triangle_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Triangle_3 Triangle;
typedef typename K::Point_3 Point;
@ -125,13 +131,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting plane
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_triangles_3
// fits a line to a 3D triangle set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -139,7 +147,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Triangle_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<2>& tag)
const CGAL::Dimension_tag<2>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::FT FT;
typedef typename K::Triangle_3 Triangle;
@ -151,17 +160,18 @@ linear_least_squares_fitting_3(InputIterator first,
c = centroid(first,beyond,K(),tag);
// assemble covariance matrix
FT covariance[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
assemble_covariance_matrix_3(first,beyond,covariance,c,k,(Triangle*) NULL,tag);
// compute fitting line
return fitting_line_3(covariance,c,line,k);
return fitting_line_3(covariance,c,line,k,diagonalize_traits);
} // end linear_least_squares_fitting_triangles_3
// fits a line to a 3D triangle set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -169,7 +179,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Triangle_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<1>& tag)
const CGAL::Dimension_tag<1>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Triangle_3 Triangle;
typedef typename K::Segment_3 Segment;
@ -189,13 +200,15 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)NULL,k,tag);
return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_triangles_3
// fits a line to a 3D triangle set
template < typename InputIterator,
typename K >
typename K,
typename DiagonalizeTraits >
typename K::FT
linear_least_squares_fitting_3(InputIterator first,
InputIterator beyond,
@ -203,7 +216,8 @@ linear_least_squares_fitting_3(InputIterator first,
typename K::Point_3& c, // centroid
const typename K::Triangle_3*, // used for indirection
const K& k, // kernel
const CGAL::Dimension_tag<0>& tag)
const CGAL::Dimension_tag<0>& tag,
const DiagonalizeTraits& diagonalize_traits)
{
typedef typename K::Triangle_3 Triangle;
typedef typename K::Point_3 Point;
@ -222,7 +236,8 @@ linear_least_squares_fitting_3(InputIterator first,
}
// compute fitting line
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag);
return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)NULL,k,tag,
diagonalize_traits);
} // end linear_least_squares_fitting_triangles_3

8
Principal_component_analysis/test/Principal_component_analysis/test_linear_least_squares_fitting_points_2.cpp
View File

@ -3,6 +3,7 @@
#include <CGAL/algorithm.h>
#include <CGAL/linear_least_squares_fitting_2.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <vector>
#include <cassert>
@ -31,7 +32,9 @@ void test_2D()
FT quality;
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),k);
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),k,
CGAL::Default_diagonalize_traits<FT,2>());
std::cout << "done (quality: " << quality << ")" << std::endl;
if(!line.is_horizontal())
@ -71,7 +74,8 @@ void test_2D_point_set(const unsigned int nb_points)
FT quality;
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),k);
quality = linear_least_squares_fitting_2(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),k,
CGAL::Default_diagonalize_traits<FT,2>());
std::cout << "done (quality: " << quality << ")" << std::endl;

10
Principal_component_analysis/test/Principal_component_analysis/test_linear_least_squares_fitting_points_3.cpp
View File

@ -3,6 +3,8 @@
#include <CGAL/Cartesian.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <cassert>
#include <cstdlib>
@ -43,13 +45,17 @@ void fit_point_set(std::list<Point>& points,
std::cout << "fit 3D line...";
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),kernel);
quality = linear_least_squares_fitting_3(points.begin(),points.end(),line,centroid,CGAL::Dimension_tag<0>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
std::cout << "done (quality: " << quality << ")" << std::endl;
std::cout << "fit 3D plane...";
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,centroid,CGAL::Dimension_tag<0>());
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,centroid,CGAL::Dimension_tag<0>(),kernel);
quality = linear_least_squares_fitting_3(points.begin(),points.end(),plane,centroid,CGAL::Dimension_tag<0>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
std::cout << "done (quality: " << quality << ")" << std::endl;
}

9
Principal_component_analysis/test/Principal_component_analysis/test_linear_least_squares_fitting_segments_3.cpp
View File

@ -5,6 +5,7 @@
#include <CGAL/algorithm.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/point_generators_3.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <vector>
#include <cassert>
@ -34,13 +35,15 @@ FT fit_set(std::list<Segment>& segments,
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,centroid,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,centroid,CGAL::Dimension_tag<1>(),kernel);
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),line,centroid,CGAL::Dimension_tag<1>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,centroid,CGAL::Dimension_tag<1>());
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,centroid,CGAL::Dimension_tag<1>(),kernel);
quality = linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,centroid,CGAL::Dimension_tag<1>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
return quality;
return quality;
}

15
Principal_component_analysis/test/Principal_component_analysis/test_linear_least_squares_fitting_spheres_3.cpp
View File

@ -2,6 +2,7 @@
#include <CGAL/Cartesian.h>
#include <CGAL/linear_least_squares_fitting_3.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <list>
@ -34,19 +35,23 @@ int main(void)
Kernel kernel;
Point centroid;
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<3>(),kernel);
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<2>(),kernel);
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<3>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),line,centroid,CGAL::Dimension_tag<2>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<2>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<3>(),kernel);
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<2>(),kernel);
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<3>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
linear_least_squares_fitting_3(spheres.begin(),spheres.end(),plane,centroid,CGAL::Dimension_tag<2>(),kernel,
CGAL::Default_diagonalize_traits<FT,3>());
return 0;
}

2
Ridges_3/doc/Ridges_3/PackageDescription.txt
View File

@ -15,7 +15,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.3}
\cgalPkgDependsOn{ \ref PkgJet_fitting_3Summary and solvers from \ref thirdpartyEigen.}
\cgalPkgDependsOn{\ref PkgSolverSummary}
\cgalPkgBib{cgal:cp-arutsm}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgShortInfoEnd

136
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Concepts/WeightedApproximation_3.h
View File

@ -1,136 +0,0 @@
//Orthogonal projection of a point onto a weighted PCA plane.
//Copyright (C) 2014 INRIA - Sophia Antipolis
//
//This program is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//This program is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
//GNU General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
/// A concept for computing an approximation of a weighted point set.
/** \ingroup PkgScaleSpaceReconstruction3Concepts
* \cgalConcept
* This approximation can be used to fit other points to the point set.
*
* \cgalHasModel `CGAL::Weighted_PCA_approximation_3`
*/
class WeightedApproximation_3 {
public:
/// \name Types
/// \{
typedef unspecified_type FT; ///< defines the field number type.
typedef unspecified_type Point; ///< defines the point type.
/// \}
public:
/// \name Constructors
/// \{
/// constructs an approximation of a undefined point set.
/** The point set holds a fixed number of points with undefined coordinates.
*
* \param size is the size of the points set.
*
* \note this does not compute the approximation.
*/
WeightedApproximation_3( unsigned int size );
/// \}
// computes the approximation of a point set.
/* Similar to constructing a approximation and calling
* <code>[set_points(points_begin, points_end, weights_begin)](\ref WeightedApproximation_3::set_points )</code>
*
* \tparam PointIterator is an input iterator over the point collection.
* The value type of the iterator must be a `Point`.
* \tparam WeightIterator is an input iterator over the collection of
* weights. The value type of the iterator must be an `FT`.
*
* \param points_begin is an iterator to the first point.
* \param points_end is a past-the-end oterator for the points.
* \param weights_begin is an iterator to the weight of the first point.
*
* \pre The number of weights must be at least equal to the number of
* points.
*/
template < typename PointIterator, typename WeightIterator >
WeightedApproximation_3( PointIterator points_begin, PointIterator points_end, WeightIterator weights_begin );
public:
/// \name Point Set.
/// \{
/// changes a weighted point in the set.
/** This invalidates the approximation. `compute()` should be called
* after all points have been set.
* \pre i must be smaller than the total size of the point set.
*/
void set_point( unsigned int i, const Point& p, const FT& w );
/// gives the size of the weighted point set.
std::size_t size() const;
/// \}
// changes the weighted point set.
/* After these points are set, the approximation is immediately computed.
*
* \tparam PointIterator is an input iterator over the point collection.
* The value type of the iterator must be a `Point`.
* \tparam WeightIterator is an input iterator over the collection of
* weights. The value type of the iterator must be an `FT`.
*
* \param points_begin is an iterator to the first point.
* \param points_end is a past-the-end iterator for the points.
* \param weights_begin is an iterator to the weight of the first point.
*
* \return whether the approximation converges. If the approximation does
* not converge this may indicate that the point set is too small, or the
* affine hull of the points cannot contain the approximation.
*
* \pre The points must fit in the approximation.
* \pre The number of weights must be at least equal to the number of
* points.
*/
template < typename PointIterator, typename WeightIterator >
bool set_points( PointIterator points_begin, PointIterator points_end,
WeightIterator weights_begin );
public:
/// \name Approximation
/// \{
/// computes the approximation.
/** \return whether the approximation converges. If the approximation does
* not converge this may indicate that the point set is too small, or the
* affine hull of the points cannot contain the approximation.
*/
bool compute();
/// checks whether the approximation has been computed successfully.
bool is_computed() const;
/// \}
public:
/// \name Fitting
/// \{
/// fits a point to the approximation.
/** \param p is the point to fit.
*
* \return the point on the approximation closest to `p`.
*
* \pre The approximation must have been computed.
*/
Point fit( const Point& p );
/// \}
}; // class WeightedApproximation_3

7
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/PackageDescription.txt
View File

@ -20,20 +20,15 @@
\cgalPkgSince{4.6}
\cgalPkgBib{cgal:ssr3}
\cgalPkgLicense{\ref licensesGPL "GPL" }
\cgalPkgDependsOn{\ref PkgAlphaShapes3Summary, \ref PkgSpatialSearchingDSummary, eigenvector solver such as \ref thirdpartyEigen 3.1.2 (or greater)}
\cgalPkgDependsOn{\ref PkgAlphaShapes3Summary, \ref PkgSpatialSearchingDSummary, \ref PkgSolverSummary}
\cgalPkgShortInfoEnd
\cgalPkgDescriptionEnd
\cgalClassifedRefPages
## Concepts ##
- `WeightedApproximation_3`
## Classes ##
- `CGAL::Scale_space_surface_reconstruction_3<Gt,FS,Sh,wA,Ct>`
- `CGAL::Weighted_PCA_approximation_3<Kernel>`
*/

2
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/Scale_space_reconstruction_3.txt
View File

@ -90,7 +90,7 @@ The main class `Scale_space_surface_reconstruction_3` contains all the functiona
The neighborhood size is estimated using `Orthogonal_k_neighbor_search`. The point set is generally stored in a `Orthogonal_k_neighbor_search::Tree`. When the neighborhood size is estimated, this tree is searched for nearest neighbors.
The scale-space is constructed at the original scale of the points. An iteration of increasing the scale is computed using a weighted PCA procedure. As described by Digne <i>et al.</i> \cgalCite{cgal:dmsl-ssmrp-11}, unlike similar methods this procedure does not lead to an <em>undesirable clustering effect</em>. By default the efficient \ref thirdpartyEigen library is used for this procedure. It is also possible to provide your own model for the `WeightedApproximation_3` concept. The weighted PCA procedure is performed locally per point, so it can be performed with parallel computing (linking with Intel TBB and passing the `Parallel_tag` to the reconstruction class is required).
The scale-space is constructed at the original scale of the points. An iteration of increasing the scale is computed using a weighted PCA procedure. As described by Digne <i>et al.</i> \cgalCite{cgal:dmsl-ssmrp-11}, unlike similar methods this procedure does not lead to an <em>undesirable clustering effect</em>. By default the efficient \ref thirdpartyEigen library is used for this procedure if available. Otherwise, the internal fallback `Diagonalize_traits` is used. It is also possible to provide your own model for the `DiagonalizeTraits` concept. The weighted PCA procedure is performed locally per point, so it can be performed with parallel computing (linking with Intel TBB and passing the `Parallel_tag` to the reconstruction class is required).
The mesh reconstruction is performed by filtering a \ref Chapter_3D_Alpha_Shapes "3D alpha shape" of the point set at a fixed scale. This filtering constructs a triangle for each regular facet; each singular facet results in two triangles facing opposite directions.

1
Scale_space_reconstruction_3/doc/Scale_space_reconstruction_3/dependencies
View File

@ -4,3 +4,4 @@ STL_Extension
Triangulation_3
Alpha_shapes_3
Spatial_searching
Solver_interface

53
Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Scale_space_surface_reconstruction_3_impl.h
View File

@ -162,21 +162,58 @@ public:
return;
Static_search search(_tree, _pts[i], _nn[i]);
Approximation pca( _nn[i] );
Point barycenter (0., 0., 0.);
FT weight_sum = 0.;
unsigned int column = 0;
// Compute total weight
for( typename Static_search::iterator nit = search.begin();
nit != search.end() && column < _nn[i];
++nit, ++column ) {
pca.set_point( column, boost::get<0>(nit->first), 1.0 / _nn[boost::get<1>(nit->first)] );
}
++nit, ++column )
weight_sum += (1.0 / _nn[boost::get<1>(nit->first)]);
column = 0;
// Compute barycenter
for( typename Static_search::iterator nit = search.begin();
nit != search.end() && column < _nn[i];
++nit, ++column )
{
Vector v (CGAL::ORIGIN, boost::get<0>(nit->first));
barycenter = barycenter + ((1.0 / _nn[boost::get<1>(nit->first)]) / weight_sum) * v;
}
CGAL::cpp11::array<FT, 6> covariance = {{ 0., 0., 0., 0., 0., 0. }};
column = 0;
// Compute covariance matrix of Weighted PCA
for( typename Static_search::iterator nit = search.begin();
nit != search.end() && column < _nn[i];
++nit, ++column )
{
Vector v (barycenter, boost::get<0>(nit->first));
FT w = (1.0 / _nn[boost::get<1>(nit->first)]);
v = w*v;
covariance[0] += w * v.x () * v.x ();
covariance[1] += w * v.x () * v.y ();
covariance[2] += w * v.x () * v.z ();
covariance[3] += w * v.y () * v.y ();
covariance[4] += w * v.y () * v.z ();
covariance[5] += w * v.z () * v.z ();
}
// Compute the weighted least-squares planar approximation of the point set.
if( !pca.compute() )
return;
CGAL::cpp11::array<FT, 9> eigenvectors = {{ 0., 0., 0.,
0., 0., 0.,
0., 0., 0. }};
CGAL::cpp11::array<FT, 3> eigenvalues = {{ 0., 0., 0. }};
wA::diagonalize_selfadjoint_covariance_matrix
(covariance, eigenvalues, eigenvectors);
// The vertex is moved by projecting it onto the plane
// through the barycenter and orthogonal to the Eigen vector with smallest Eigen value.
_pts[i] = pca.fit( _pts[i] );
Vector norm (eigenvectors[0], eigenvectors[1], eigenvectors[2]);
Vector b2p (barycenter, _pts[i]);
_pts[i] = barycenter + b2p - ((norm * b2p) * norm);
}
}; // class AdvanceSS

177
Scale_space_reconstruction_3/include/CGAL/Scale_space_reconstruction_3/Weighted_PCA_approximation_3.h
View File

@ -1,177 +0,0 @@
//Copyright (C) 2014 INRIA - Sophia Antipolis
//
//This program is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//This program is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
//GNU General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author(s): Thijs van Lankveld
#ifndef CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_PROJECTION_3_H
#define CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_PROJECTION_3_H
#include <Eigen/Dense>
namespace CGAL {
/// approximates a point set using a weighted least-squares plane.
/** \ingroup PkgScaleSpaceReconstruction3Classes
* The weighted least-squares planar approximation contains the barycenter
* of the points and is orthogonal to the eigenvector corresponding to the
* smallest eigenvalue.
*
* Point are fitted to this approximation by projecting them orthogonally onto
* the weighted least-squares plane.
*
* This class requires the eigenvector solvers of \ref thirdpartyEigen version
* 3.1.2 (or greater).
*
* \tparam Kernel the geometric traits class of the input and output. It
* must have a `RealEmbeddable` field number type.
*
* \note Irrespective of the geometric traits class, the approximation is
* always estimated up to double precision.
*
* \cgalModels `WeightedApproximation_3`
*/
template < class Kernel >
class Weighted_PCA_approximation_3 {
public:
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
private:
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Matrix3D; // 3-by-dynamic double-value matrix.
typedef Eigen::Array<double, 1, Eigen::Dynamic> Array1D; // 1-by-dynamic double-value array.
typedef Eigen::Matrix3d Matrix3; // 3-by-3 double-value matrix.
typedef Eigen::Vector3d Vector3; // 3-by-1 double-value matrix.
typedef Eigen::SelfAdjointEigenSolver< Matrix3 > EigenSolver; // Eigen value decomposition solver.
private:
bool _comp; // Whether the eigenvalues are computed.
Matrix3D _pts; // points.
Array1D _wts; // weights.
Vector3 _bary; // barycenter.
Vector3 _norm; // normal.
public:
// constructs an default approximation to hold the points.
/* \param size is the number of points that will be added.
*/
Weighted_PCA_approximation_3( unsigned int size ): _comp(false), _pts(3,size), _wts(1,size) {}
// computes the weighted least-squares planar approximation of a point set.
/* Similar to constructing an empty projection and calling
* <code>[set_points(points_begin, points_end, weights_begin)](\ref WeightedApproximation_3::set_points )</code>
*
* \tparam PointIterator is an input iterator over the point collection.
* The value type of the iterator must be a `Point`.
* \tparam WeightIterator is an input iterator over the collection of
* weights. The value type of the iterator must be an `FT`.
*
* \param points_begin is an iterator to the first point.
* \param points_end is a past-the-end oterator for the points.
* \param weights_begin is an iterator to the weight of the first point.
*/
template < typename PointIterator, typename WeightIterator >
Weighted_PCA_approximation_3( PointIterator points_begin, PointIterator points_end, WeightIterator weights_begin )
: _comp(false) {
std::size_t size = std::distance( points_begin, points_end );
_pts = Matrix3D(3,size);
_wts = Array1D(1,size);
set_points( points_begin, points_end, weights_begin );
}
public:
// sets a weighted point in the collection.
void set_point( unsigned int i, const Point& p, const FT& w ) {
_pts( 0, i ) = CGAL::to_double( p[0] );
_pts( 1, i ) = CGAL::to_double( p[1] );
_pts( 2, i ) = CGAL::to_double( p[2] );
_wts( i ) = CGAL::to_double( w );
_comp = false;
}
// sets the weighted points and compute PCA.
template < typename PointIterator, typename WeightIterator >
bool set_points( PointIterator points_begin, PointIterator points_end, WeightIterator weights_begin ) {
for( unsigned int column = 0; points_begin != points_end; ++column, ++points_begin, ++weights_begin ) {
_pts( 0, column ) = CGAL::to_double( (*points_begin)[0] );
_pts( 1, column ) = CGAL::to_double( (*points_begin)[1] );
_pts( 2, column ) = CGAL::to_double( (*points_begin)[2] );
_wts( column ) = CGAL::to_double( *weights_begin );
}
return compute();
}
// gives the size of the weighted point set.
std::size_t size() const { return _pts.cols(); }
// compute weighted PCA.
bool compute() {
// Construct the barycenter.
_bary = ( _pts.array().rowwise() * _wts ).rowwise().sum() / _wts.sum();
// Replace the points by their difference with the barycenter.
_pts = ( _pts.colwise() - _bary ).array().rowwise() * _wts;
// Construct the weighted covariance matrix.
Matrix3 covar = Matrix3::Zero();
for( int column = 0; column < _pts.cols(); ++column )
covar += _wts.matrix()(column) * _pts.col(column) * _pts.col(column).transpose();
// Construct the Eigen system.
EigenSolver solver( covar );
// If the Eigen system does not converge, the vertex is not moved.
if( solver.info() != Eigen::Success )
return false;
// Find the Eigen vector with the smallest Eigen value.
std::ptrdiff_t index;
solver.eigenvalues().minCoeff( &index );
if( solver.eigenvectors().col( index ).imag() != Vector3::Zero() ) {
// This should actually never happen,
// because the covariance matrix is symmetric!
CGAL_assertion( false );
return false;
}
// The normal is the Eigen vector with smallest Eigen value.
_norm = solver.eigenvectors().col( index ).real();
_comp = true;
return true;
}
// checks whether the weighted least-squares approximating plane has been computed.
bool is_computed() const { return _comp; }
public:
// projects a point onto the weighted PCA plane.
Point fit( const Point& p ) {
CGAL_assertion( _comp );
// The point is moved by projecting it onto the plane through the
// barycenter and orthogonal to the normal.
Vector3 to_p = Vector3( to_double( p[0] ),
to_double( p[1] ),
to_double( p[2] ) ) - _bary;
Vector3 proj = _bary + to_p - ( _norm.dot(to_p) * _norm );
return Point( proj(0), proj(1), proj(2) );
}
}; // class Weighted_PCA_approximation_3
} // namespace CGAL
#endif // CGAL_SCALE_SPACE_RECONSTRUCTION_3_WEIGHTED_PCA_PROJECTION_3_H

26
Scale_space_reconstruction_3/include/CGAL/Scale_space_surface_reconstruction_3.h
View File

@ -36,13 +36,10 @@
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Random.h>
#include <CGAL/Default_diagonalize_traits.h>
#include <CGAL/Scale_space_reconstruction_3/Shape_construction_3.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Scale_space_reconstruction_3/Weighted_PCA_approximation_3.h>
#endif // CGAL_EIGEN3_ENABLED
#include <boost/mpl/and.hpp>
namespace CGAL {
@ -99,25 +96,22 @@ namespace CGAL {
* only has an impact on the run-time.
* \tparam Sh determines whether to collect the surface per shell. It
* must be a `Boolean_tag` type. The default value is `Tag_true`.
* \tparam wA must be a model of `WeightedPCAProjection_3` and determines how
* to approximate a weighted point set. If \ref thirdpartyEigen 3.1.2 (or
* greater) is available and CGAL_EIGEN3_ENABLED is defined, then
* `Weighted_PCA_approximation_3<DelaunayTriangulationTraits_3>` is used by default.
* \tparam wA must be a model of `DiagonalizeTraits` and determines
* how to diagonalize a weighted covariance matrix to approximate a
* weighted point set. It can be omitted: if Eigen 3 (or greater) is
* available and `CGAL_EIGEN3_ENABLED` is defined then an overload
* using `Eigen_diagonalize_traits` is provided. Otherwise, the
* internal implementation `Diagonalize_traits` is used.
* \tparam Ct indicates whether to use concurrent processing. It must be
* either `Sequential_tag` or `Parallel_tag` (the default value).
*/
template < class Gt, class FS = Tag_true, class Sh = Tag_true, class wA = Default, class Ct = Parallel_tag >
template < class Gt, class FS = Tag_true, class Sh = Tag_true,
class wA = Default_diagonalize_traits<typename Gt::FT, 3>, class Ct = Parallel_tag >
class Scale_space_surface_reconstruction_3 {
typedef typename Default::Get< wA,
#ifdef CGAL_EIGEN3_ENABLED
Weighted_PCA_approximation_3<Gt>
#else // CGAL_EIGEN3_ENABLED
void
#endif // CGAL_EIGEN3_ENABLED
>::type Approximation;
public:
typedef typename Gt::Point_3 Point; ///< defines the point type.
typedef typename Gt::Vector_3 Vector; ///< defines the vector type.
typedef boost::tuple<Point,int> Point_and_int;
private:

23
Solver_interface/doc/Solver_interface/CGAL/Diagonalize_traits.h Normal file
View File

@ -0,0 +1,23 @@
namespace CGAL {
/*!
\ingroup PkgSolver
The class `Diagonalize_traits` provides an internal
implementation for the diagonalization of Variance-Covariance
Matrices.
\cgalModels `DiagonalizeTraits`
*/
template <typename FT, unsigned int dim = 3>
class Diagonalize_traits{
};
} // namespace CGAL

31
Solver_interface/doc/Solver_interface/CGAL/Eigen_diagonalize_traits.h Normal file
View File

@ -0,0 +1,31 @@
namespace CGAL {
/*!
\ingroup PkgSolver
The class `Eigen_diagonalize_traits` provides an interface to the
diagonalization of covariance matrices of \ref thirdpartyEigen
"Eigen".
The version 3.1 (or greater) of \ref thirdpartyEigen "Eigen" must be
available on the system.
\cgalModels `DiagonalizeTraits`
\sa http://eigen.tuxfamily.org
*/
template <typename FT, unsigned int dim = 3>
class Eigen_diagonalize_traits{
};
} // namespace CGAL

49
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/include/CGAL/Eigen_matrix.h → Solver_interface/doc/Solver_interface/CGAL/Eigen_matrix.h
View File

@ -2,24 +2,18 @@
namespace CGAL {
/*!
\ingroup PkgSurfaceParameterizationAlgebra
\ingroup PkgSolver
The class `Eigen_sparse_matrix` is a C++ wrapper around \ref thirdpartyEigen "Eigen" matrix type `Eigen::SparseMatrix`
The class `Eigen_sparse_matrix` is a wrapper around \ref thirdpartyEigen "Eigen" matrix type <a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::SparseMatrix` </a>
that represents general matrices, be they symmetric or not.
The version 3.1 (or greater) of \ref thirdpartyEigen "Eigen" must be available on the system.
\cgalModels `SparseLinearAlgebraTraits_d::Matrix`
Parameters
--------------
`T`: Number type.
\tparam T Number type.
\sa `CGAL::Eigen_solver_traits<T>`
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
\sa `CGAL::Eigen_vector<T>`
\sa http://eigen.tuxfamily.org
*/
template< typename T >
class Eigen_sparse_matrix {
@ -41,23 +35,19 @@ typedef unspecified_type EigenType;
namespace CGAL {
/*!
\ingroup PkgSurfaceParameterizationAlgebra
\ingroup PkgSolver
The class `Eigen_sparse_symmetric_matrix` is a C++ wrapper around \ref thirdpartyEigen "Eigen" matrix type `Eigen::SparseMatrix`.
The class `Eigen_sparse_symmetric_matrix` is a wrapper around \ref thirdpartyEigen "Eigen" matrix type <a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::SparseMatrix` </a>
As the matrix is symmetric only the lower triangle part is stored.
\cgalModels `SparseLinearAlgebraTraits_d::Matrix`
Parameters
--------------
`T`: Number type.
\tparam T Number type.
\sa `CGAL::Eigen_solver_traits<T>`
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
\sa `CGAL::Eigen_sparse_matrix<T>`
\sa `CGAL::Eigen_vector<T>`
\sa http://eigen.tuxfamily.org
*/
template< typename T >
@ -75,4 +65,29 @@ typedef unspecified_type EigenType;
/// @}
}; /* end Eigen_sparse_symmetric_matrix */
/*!
\ingroup PkgSolver
The class `Eigen_matrix` is a wrapper around \ref thirdpartyEigen "Eigen"
matrix type
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html">`Eigen::Matrix`</a>.
\cgalModels `SvdTraits::Matrix`
\tparam T Number type.
\sa `CGAL::Eigen_svd`
\sa `CGAL::Eigen_vector<T>`
*/
template< typename T >
class Eigen_matrix {
public:
};
} /* end namespace CGAL */

9
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/include/CGAL/Eigen_solver_traits.h → Solver_interface/doc/Solver_interface/CGAL/Eigen_solver_traits.h
View File

@ -2,18 +2,15 @@
namespace CGAL {
/*!
\ingroup PkgSurfaceParameterizationAlgebra
\ingroup PkgSolver
The class `Eigen_solver_traits` provides an interface to the sparse solvers of \ref thirdpartyEigen "Eigen".
The version 3.1 (or greater) of \ref thirdpartyEigen "Eigen" must be available on the system.
\cgalModels `#SparseLinearAlgebraTraitsWithFactor_d` and `#NormalEquationSparseLinearAlgebraTraits_d`
\cgalModels `SparseLinearAlgebraWithFactorTraits_d` and `NormalEquationSparseLinearAlgebraTraits_d`
Parameters
--------------
`T`: a sparse solver of \ref thirdpartyEigen "Eigen". The default solver is the iterative bi-congugate gradient stabilized solver
`Eigen::BiCGSTAB` for `double`.
\tparam T A sparse solver of \ref thirdpartyEigen "Eigen". The default solver is the iterative bi-congugate gradient stabilized solver `Eigen::BiCGSTAB` for `double`.
\sa `CGAL::Eigen_sparse_matrix<T>`
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`

10
Jet_fitting_3/doc/Jet_fitting_3/CGAL/Eigen_svd.h → Solver_interface/doc/Solver_interface/CGAL/Eigen_svd.h
View File

@ -2,11 +2,11 @@
namespace CGAL {
/*!
\ingroup PkgJet_fitting_3
\ingroup PkgSolver
The class `Eigen_svd` provides an algorithm to solve in the least
square sense a linear system with a singular value decomposition using
\ref thirdpartyEigen. The field type is `double`.
\ref thirdpartyEigen.
\cgalModels `SvdTraits`
@ -15,7 +15,11 @@ square sense a linear system with a singular value decomposition using
class Eigen_svd {
public:
/// @}
typedef double FT;
typedef Eigen_vector<FT> Vector;
typedef Eigen_matrix<FT> Matrix;
}; /* end Eigen_svd */
} /* end namespace CGAL */

13
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/include/CGAL/Eigen_vector.h → Solver_interface/doc/Solver_interface/CGAL/Eigen_vector.h
View File

@ -2,22 +2,21 @@
namespace CGAL {
/*!
\ingroup PkgSurfaceParameterizationAlgebra
\ingroup PkgSolver
The class `Eigen_vector` is a C++ wrapper around \ref thirdpartyEigen "Eigen" vector, which is a simple array of numbers.
The version 3.1 (or greater) of \ref thirdpartyEigen "Eigen" must be available on the system.
The class `Eigen_vector` is a wrapper around \ref thirdpartyEigen "Eigen" vector
type <a href="http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html"> </a>,
which is a simple array of numbers.
\cgalModels `SvdTraits::Vector`
\cgalModels `SparseLinearAlgebraTraits_d::Vector`.
Parameters
--------------
`T`: Number type.
\tparam T Number type.
\sa `CGAL::Eigen_solver_traits<T>`
\sa `CGAL::Eigen_sparse_matrix<T>`
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
\sa http://eigen.tuxfamily.org
*/
template< typename T >

57
Solver_interface/doc/Solver_interface/Concepts/DiagonalizeTraits.h Normal file
View File

@ -0,0 +1,57 @@
/*!
\ingroup PkgSolverConcepts
\cgalConcept
Concept providing functions to extract eigenvectors and eigenvalues
from covariance matrices represented by an array `a`, using symmetric
diagonalization. For example, a matrix of dimension 3 is defined as
follows:
<center>
\f$ \begin{bmatrix}
a[0] & a[1] & a[2] \\
a[1] & a[3] & a[4] \\
a[2] & a[4] & a[5] \\
\end{bmatrix}\f$
</center>
\tparam FT Number type
\tparam dim Dimension of the matrices and vectors
\cgalHasModel `CGAL::Eigen_diagonalize_traits`
\cgalHasModel `CGAL::Diagonalize_traits`
*/
template <typename FT, unsigned int dim = 3>
class DiagonalizeTraits
{
public:
/// fill `eigenvalues` with the eigenvalues of the covariance matrix represented by `cov`.
/// Eigenvalues are sorted by increasing order.
/// \return `true` if the operation was successful and `false` otherwise.
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues);
/// fill `eigenvalues` with the eigenvalues and `eigenvectors` with
/// the eigenvectors of the covariance matrix represented by `cov`.
/// Eigenvalues are sorted by increasing order.
/// \return `true` if the operation was successful and `false`
/// otherwise.
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues,
cpp11::array<FT, dim * dim>& eigenvectors);
/// Extract the eigenvector associated to the largest eigenvalue
/// of the covariance matrix represented by `cov`.
/// \return `true` if the operation was successful and `false` otherwise.
static bool
extract_largest_eigenvector_of_covariance_matrix (
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT,dim> &normal);
};

5
Surface_mesh_skeletonization/doc/Surface_mesh_skeletonization/Concepts/NormalEquationSparseLinearAlgebraTraits_d.h → Solver_interface/doc/Solver_interface/Concepts/NormalEquationSparseLinearAlgebraTraits_d.h
View File

@ -1,6 +1,7 @@
/*!
\ingroup PkgMeanCurvatureSkeleton3Concepts
\ingroup PkgSolverConcepts
\cgalConcept
@ -72,5 +73,5 @@ bool normal_equation_solver(const Matrix& A, const Vector& B, Vector& X);
/// @}
}; /* end SparseLinearAlgebraTraitsWithFactor_d */
}; /* end NormalEquationSparseLinearAlgebraTraits_d */

28
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/Concepts/SparseLinearAlgebraTraits_d.h → Solver_interface/doc/Solver_interface/Concepts/SparseLinearAlgebraTraits_d.h
View File

@ -1,18 +1,14 @@
/*!
\ingroup PkgSurfaceParameterizationConcepts
\ingroup PkgSolverConcepts
\cgalConcept
The concept `SparseLinearAlgebraTraits_d` is used to solve sparse linear systems <I>A\f$ \times \f$ X = B</I>.
\cgalRefines `LinearAlgebraTraits_d`
\cgalHasModel `CGAL::Eigen_solver_traits<T>`
\cgalHasModel `OpenNL::DefaultLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` in OpenNL package
\cgalHasModel `OpenNL::SymmetricLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` in OpenNL package
\sa `SparseLinearAlgebraTraits_d::Matrix`
\sa `SparseLinearAlgebraTraits_d::Vector`
*/
@ -66,23 +62,23 @@ bool linear_solver(const Matrix& A, const Vector& B, Vector& X, NT& D);
/// @}
}; /* end SparseLinearAlgebraTraits_d */
/*!
\ingroup PkgSurfaceParameterizationConcepts
\cgalConcept
`SparseLinearAlgebraTraits_d::Vector` is a concept of a vector that can be multiplied by a sparse matrix.
\cgalRefines `LinearAlgebraTraits_d::Vector`
\cgalHasModel `CGAL::Eigen_vector<T>`
\cgalHasModel `OpenNL::FullVector<T>` in `OpenNL` package
\sa `SparseLinearAlgebraTraits_d`
\sa `SparseLinearAlgebraTraits_d::Matrix`
*/
class Vector {
class SparseLinearAlgebraTraits_d::Vector {
public:
/// \name Types
@ -142,23 +138,20 @@ NT& operator[](int row);
}; /* end Vector */
/*!
\ingroup PkgSurfaceParameterizationConcepts
\cgalConcept
`SparseLinearAlgebraTraits_d::Matrix` is a concept of a sparse matrix class.
\cgalRefines `LinearAlgebraTraits_d::Matrix`
\cgalHasModel `CGAL::Eigen_sparse_matrix<T>`
\cgalHasModel `CGAL::Eigen_sparse_symmetric_matrix<T>`
\cgalHasModel `OpenNL::SparseMatrix<T>` in `OpenNL` package
\sa `SparseLinearAlgebraTraits_d`
\sa `SparseLinearAlgebraTraits_d::Vector`
*/
class Matrix {
class SparseLinearAlgebraTraits_d::Matrix {
public:
/// \name Types
@ -243,6 +236,3 @@ void set_coef(int row, int column, NT value, bool new_coef = false);
}; /* end Matrix */
}; /* end SparseLinearAlgebraTraits_d */

18
Surface_modeling/doc/Surface_modeling/Concepts/SparseLinearAlgebraTraitsWithFactor_d.h → Solver_interface/doc/Solver_interface/Concepts/SparseLinearAlgebraWithFactorTraits_d.h
View File

@ -1,6 +1,6 @@
/*!
\ingroup PkgSurfaceModelingConcepts
\ingroup PkgSolverConcepts
\cgalConcept
@brief Concept describing the set of requirements for a direct sparse linear system solver with factorization.
@ -12,7 +12,7 @@ method to solve the system for different right-hand vectors.
\cgalHasModel `CGAL::Eigen_solver_traits<T>`
*/
class SparseLinearAlgebraTraitsWithFactor_d {
class SparseLinearAlgebraWithFactorTraits_d {
public:
/// \name Creation
@ -23,7 +23,7 @@ public:
Default constructor.
*/
SparseLinearAlgebraTraitsWithFactor_d();
SparseLinearAlgebraWithFactorTraits_d();
/// @}
@ -31,20 +31,20 @@ SparseLinearAlgebraTraitsWithFactor_d();
/// @{
/*!
Factorize the sparse matrix A.
This factorization is used in SparseLinearAlgebraTraitsWithFactor_d::linear_solver to solve the system for different right-hand side vectors.
Factorize the sparse matrix `A`.
This factorization is used in `SparseLinearAlgebraWithFactorTraits_d::linear_solver()` to solve the system for different right-hand side vectors.
See `::SparseLinearAlgebraTraits_d::linear_solver()` for the description of `D`.
@return true if the factorization is successful
@return `true` if the factorization is successful
*/
bool factor(const Matrix& A, NT& D);
/*!
Solve the sparse linear system \f$ A \times X = B\f$, with \f$ A \f$ being the matrix provided in SparseLinearAlgebraTraitsWithFactor_d::factor.
@return true if the solver is successful
Solve the sparse linear system \f$ A \times X = B\f$, with \f$ A \f$ being the matrix provided in `SparseLinearAlgebraWithFactorTraits_d::factor()`.
@return `true` if the solver is successful
*/
bool linear_solver(const Vector& B, Vector& X);
/// @}
}; /* end SparseLinearAlgebraTraitsWithFactor_d */
}; /* end SparseLinearAlgebraWithFactorTraits_d */

40
Jet_fitting_3/doc/Jet_fitting_3/Concepts/SvdTraits.h → Solver_interface/doc/Solver_interface/Concepts/SvdTraits.h
View File

@ -1,24 +1,12 @@
/*!
\ingroup PkgJet_fitting_3Concepts
\ingroup PkgSolverConcepts
\cgalConcept
The concept `SvdTraits` describes the set of requirements to be
fulfilled by any class used to instantiate the third template
parameter of the class
`CGAL::Monge_via_jet_fitting<DataKernel,LocalKernel,SvdTraits>`.
It describes the linear algebra types and algorithms needed by the
class `CGAL::Monge_via_jet_fitting`.
\cgalHeading{Requirements}
The scalar type, `SvdTraits::FT`, must be the same as that of
the `LocalKernel` concept : `LocalKernel::FT`.
The concept `SvdTraits` describes the linear algebra types and algorithms needed
to solve in the least square sense a linear system with a singular value decomposition
\cgalHasModel `CGAL::Eigen_svd`
\sa `LocalKernel`
*/
class SvdTraits {
public:
@ -39,7 +27,7 @@ public:
/*!
The matrix type, model of the concept `SvdTraits::Matrix`.
*/
typedef unspecified_type matrix;
typedef unspecified_type Matrix;
/// @}
@ -60,9 +48,10 @@ public:
}; /* end SvdTraits */
/*!
\ingroup PkgJet_fitting_3Concepts
\cgalConcept
Concept of vector type used by the concept SvdTraits.
Concept of vector type used by the concept `SvdTraits`.
\cgalHasModel `CGAL::Eigen_vector<T>`
*/
class SvdTraits::Vector {
public:
@ -76,12 +65,12 @@ public:
size_t size();
/*!
return the \f$ i^{th}\f$ entry, \f$ i\f$ from \f$ 0\f$ to \f$ size()-1\f$.
return the `i`th entry, `i` from `0` to `size()-1`.
*/
FT operator()(size_t i);
/*!
set the \f$ i^{th}\f$ entry to `value`.
set the `i`'th entry to `value`.
*/
void set(size_t i, const FT value);
@ -93,9 +82,10 @@ public:
/*!
\ingroup PkgJet_fitting_3Concepts
\cgalConcept
Concept of matrix type used by the concept SvdTraits.
Concept of matrix type used by the concept `SvdTraits`.
\cgalHasModel `CGAL::Eigen_matrix<T>`
*/
class SvdTraits::Matrix {
public:
@ -115,13 +105,13 @@ public:
size_t number_of_columns();
/*!
return the entry at row \f$ i\f$ and column \f$ j\f$, \f$ i\f$ from \f$ 0\f$ to `number_of_rows - 1`,
\f$ j\f$ from \f$ 0\f$ to `number_of_columns - 1`.
return the entry at row `i` and column `j`, `i` from `0` to `number_of_rows - 1`,
`j` from `0` to `number_of_columns - 1`.
*/
FT operator()(size_t i, size_t j);
/*!
set the entry at row \f$ i\f$ and column \f$ j\f$ to \f$ value\f$.
set the entry at row `i` and column `j` to `value`.
*/
void set(size_t i, size_t j, const FT value);
};

7
Solver_interface/doc/Solver_interface/Doxyfile.in Normal file
View File

@ -0,0 +1,7 @@
@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
PROJECT_NAME = "CGAL ${CGAL_CREATED_VERSION_NUM} - CGAL and Solvers"
INPUT = ${CMAKE_SOURCE_DIR}/Solver_interface/doc/Solver_interface/
EXAMPLE_PATH = ${CMAKE_SOURCE_DIR}/Solver_interface/examples/Solver_interface/

50
Solver_interface/doc/Solver_interface/PackageDescription.txt Normal file
View File

@ -0,0 +1,50 @@
/// \defgroup PkgSolver CGAL and Solvers Reference
/// \defgroup PkgSolverConcepts Concepts
/// \ingroup PkgSolver
///
/*!
\addtogroup PkgSolver
\cgalPkgDescriptionBegin{CGAL and Solvers,PkgSolverSummary}
\cgalPkgPicture{solver.png}
\cgalPkgSummaryBegin
\cgalPkgAuthors{Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, and Laurent Saboret}
\cgalPkgDesc{This package provides concepts and models for solving linear systems with dense or sparse matrices.}
\cgalPkgManuals{Chapter_CGAL_and_Solvers,PkgSolver}
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{4.8}
\cgalPkgBib{cgal:eb-solver}
\cgalPkgLicense{\ref licensesLGPL "LGPL"}
\cgalPkgShortInfoEnd
\cgalPkgDescriptionEnd
\cgalClassifedRefPages
## Concepts ##
- `DiagonalizeTraits`
- `SparseLinearAlgebraTraits_d`
- `SparseLinearAlgebraWithFactorTraits_d`
- `SvdTraits`
## Classes ##
- `CGAL::Eigen_solver_traits`
- `CGAL::Eigen_svd`
- `CGAL::Eigen_diagonalize_traits`
- `CGAL::Eigen_matrix`
- `CGAL::Eigen_vector`
- `CGAL::Eigen_sparse_matrix`
- `CGAL::Eigen_sparse_symmetric_matrix`
- `CGAL::Diagonalize_traits`
*/

115
Solver_interface/doc/Solver_interface/Solver_interface.txt Normal file
View File

@ -0,0 +1,115 @@
namespace CGAL {
/*!
\mainpage User Manual
\anchor Chapter_CGAL_and_Solvers
\anchor chapterSolvers
\cgalAutoToc
\authors Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, and Laurent Saboret
Several \cgal packages have to solve linear systems with dense or
sparse matrices. This package provides concepts and models for that
purpose.
We generally provide models using the \ref thirdpartyEigen
library. Wrappers for the Eigen classes `Eigen_matrix` and
`Eigen_vector` are also provided when needed.
It is straightforward to develop equivalent models for
other solvers, for example those found in the <a
href="https://software.intel.com/en-us/intel-mkl">Intel Math Kernel
Library (MKL)</a>.
\section SectionSolverDiagonalize Matrix Diagonalization
The concept `DiagonalizeTraits<T,dim>` defines an interface for the
diagonalization and computation of eigenvectors and eigenvalues of a
symmetric matrix. `T` is the number type and `dim` is the dimension of
the matrices and vector (set to 3 by default). We provide two models
for this concept:
- `Eigen_diagonalize_traits<T,dim>` uses the \ref thirdpartyEigen library.
- `Diagonalize_traits<T,dim>` is an internal implementation that does not
depend on another library.
Although both models achieve the same computation,
`Eigen_diagonalize_traits<T,dim>` is faster and should thus be used if the
\ref thirdpartyEigen library is available. The eigenvalues are stored
in ascending order and eigenvectors are stored in accordance.
This is an example of an eigen decomposition of a matrix using this
class:
\cgalExample{Solver_interface/diagonalize_matrix.cpp}
\section SectionSolverSVD Singular Value Decomposition
The concept `SvdTraits` defines an interface for solving in the least
square sense a linear system with a singular value decomposition. The
field type is `double`. We provide the model `Eigen_svd` that uses the
\ref thirdpartyEigen library.
Here is a simple example that shows how to handle matrices, vectors
and this solver:
\cgalExample{Solver_interface/singular_value_decomposition.cpp}
\section SectionSolverSparse Sparse Solvers
We define 3 concepts for sparse linear algebra:
- `SparseLinearAlgebraTraits_d`
- `SparseLinearAlgebraWithFactorTraits_d`
- `NormalEquationSparseLinearAlgebraTraits_d`
An interface to the sparse solvers from the \ref thirdpartyEigen
library is provided as a model for these 3 concepts through the class
`Eigen_solver_traits<T>`. This solver traits class can be used for an
iterative or a direct, symmetric or general sparse solvers. The
specific solver to be used must be given as template parameter.
Each \cgal package using a sparse solver specifies which type of matrix
and solver is required:
\code{.cpp}
typedef CGAL::Eigen_sparse_matrix<double>::EigenType EigenMatrix;
//iterative general solver
typedef CGAL::Eigen_solver_traits< Eigen::BiCGSTAB<EigenMatrix> > Iterative_general_solver;
//iterative symmetric solver
typedef CGAL::Eigen_solver_traits< Eigen::ConjugateGradient<EigenMatrix> > Iterative_symmetric_solver;
//direct symmetric solver
typedef CGAL::Eigen_solver_traits< Eigen::SimplicialCholesky<EigenMatrix> > Direct_symmetric_solver;
\endcode
Here is an example that shows how to fill the sparse matrix and call
the solver:
\cgalExample{Solver_interface/sparse_solvers.cpp}
\section SolversHistory Implementation History
This package is the result of the increasing needs for linear solvers in \cgal.
The first packages that introduced the solver concepts were \ref PkgSurfaceParameterization,
\ref PkgSurfaceReconstructionFromPointSets and \ref PkgJet_fitting_3.
At that time, these packages were relying on \sc{Taucs}, \sc{LAPACK}, \sc{BLAS} and \sc{OpenNL}. Gaël Guennebaud
then introduced new models using the \ref thirdpartyEigen library that became the
only supported models by \cgal.
Later on the packages \ref PkgMeanCurvatureSkeleton3Summary and \ref PkgSurfaceModelingSummary
extended the existing concepts.
Simon Giraudot was responsible for gathering all concepts and classes, and also wrote this user manual
with the help of Andreas Fabri.
*/
} /* namespace CGAL */

10
Solver_interface/doc/Solver_interface/dependencies Normal file
View File

@ -0,0 +1,10 @@
Manual
Kernel_23
STL_Extension
Algebraic_foundations
Miscellany
Surface_mesh_parameterization
Surface_mesh_skeletonization
Surface_modeling
Jet_fitting_3
Surface_reconstruction_points_3

5
Solver_interface/doc/Solver_interface/examples.txt Normal file
View File

@ -0,0 +1,5 @@
/*!
\example Solver_interface/diagonalize_matrix.cpp
\example Solver_interface/singular_value_decomposition.cpp
\example Solver_interface/sparse_solvers.cpp
*/

BIN
Solver_interface/doc/Solver_interface/fig/solver.png Executable file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 909 B

35
Solver_interface/examples/Solver_interface/CMakeLists.txt Normal file
View File

@ -0,0 +1,35 @@
# Created by the script cgal_create_cmake_script
# This is the CMake script for compiling a CGAL application.
project( examples_ )
cmake_minimum_required(VERSION 2.8.10)
find_package(CGAL QUIET COMPONENTS Core )
if ( CGAL_FOUND )
include( ${CGAL_USE_FILE} )
include( CGAL_CreateSingleSourceCGALProgram )
include_directories (BEFORE "../include")
# Use Eigen
find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater)
if (EIGEN3_FOUND)
include( ${EIGEN3_USE_FILE} )
create_single_source_cgal_program( "singular_value_decomposition.cpp" )
create_single_source_cgal_program( "sparse_solvers.cpp" )
endif()
create_single_source_cgal_program( "diagonalize_matrix.cpp" )
else()
message(STATUS "This program requires the CGAL library, and will not be compiled.")
endif()

51
Solver_interface/examples/Solver_interface/diagonalize_matrix.cpp Normal file
View File

@ -0,0 +1,51 @@
#include <iostream>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#else
#include <CGAL/Diagonalize_traits.h>
#endif
typedef double FT;
typedef CGAL::cpp11::array<FT, 6> Eigen_matrix;
typedef CGAL::cpp11::array<FT, 3> Eigen_vector;
typedef CGAL::cpp11::array<FT, 9> Eigen_three_vectors;
// If Eigen is enabled, use it, otherwise fallback to the internal model
#ifdef CGAL_EIGEN3_ENABLED
typedef CGAL::Eigen_diagonalize_traits<FT, 3> Diagonalize_traits;
#else
typedef CGAL::Diagonalize_traits<FT, 3> Diagonalize_traits;
#endif
int main(void)
{
Eigen_matrix covariance = {{ 0., 0., 0., 0., 0., 0. }};
// Fill matrix with random numbers
for (std::size_t i = 0; i < 6; ++ i)
covariance[i] = rand ();
Eigen_vector eigenvalues;
Eigen_three_vectors eigenvectors;
if (!(Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix (covariance,
eigenvalues,
eigenvectors)))
{
std::cerr << "Error: cannot diagonalize matrix" << std::endl;
return -1;
}
// Print result
for (std::size_t i = 0; i < 3; ++ i)
{
std::cout << "Eigenvalue " << i+1 << " = " << eigenvalues[i] << std::endl
<< " with eigenvector [ ";
for (std::size_t j = 0; j < 3; ++ j)
std::cout << eigenvectors[3*i + j] << " ";
std::cout << "]" << std::endl;
}
return 0;
}

37
Solver_interface/examples/Solver_interface/singular_value_decomposition.cpp Normal file
View File

@ -0,0 +1,37 @@
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_matrix.h>
#include <CGAL/Eigen_vector.h>
#include <CGAL/Eigen_svd.h>
typedef CGAL::Eigen_svd Svd;
#endif
typedef Svd::FT FT;
typedef Svd::Vector Eigen_vector;
typedef Svd::Matrix Eigen_matrix;
int main(void)
{
std::size_t degree = 3;
Eigen_vector B (degree);
Eigen_matrix M (degree, degree);
// Fill B and M with random numbers
for (std::size_t i = 0; i < degree; ++ i)
{
B.set (i, rand());
for (std::size_t j = 0; j < degree; ++ j)
M.set (i, j, rand());
}
// Solve MX=B
std::cout << Svd::solve(M, B) << std::endl;
// Print result
std::cout << "Solution of SVD = [ ";
for (std::size_t i = 0; i < degree; ++ i)
std::cout << B.vector()[i] << " ";
std::cout << "]" << std::endl;
return 0;
}

49
Solver_interface/examples/Solver_interface/sparse_solvers.cpp Normal file
View File

@ -0,0 +1,49 @@
#include <CGAL/Eigen_solver_traits.h>
#include <CGAL/Eigen_matrix.h>
typedef CGAL::Eigen_solver_traits<> Eigen_solver;
typedef Eigen_solver::NT FT;
typedef Eigen_solver::Matrix Eigen_matrix;
typedef Eigen_solver::Vector Eigen_vector;
int main(void)
{
srand (static_cast<unsigned int>(time (NULL)));
std::size_t degree = 3000;
std::size_t nb_nonzero_coef = 100;
Eigen_vector B (degree); // Zero vector
Eigen_matrix A (degree);
Eigen_vector diag (degree);
// Randomly make some coefficients of the matrix non-zero
for (std::size_t i = 0; i < nb_nonzero_coef; ++ i)
{
int x = rand () % degree;
int y = rand () % degree;
FT value = rand () / (FT)RAND_MAX;
A.add_coef (x, y, value);
A.add_coef (y, x, value);
}
// Create sparse matrix
A.assemble_matrix();
Eigen_vector X (degree);
FT d;
Eigen_solver solver;
if (!(solver.linear_solver (A, B, X, d)))
{
std::cerr << "Error: linear solver failed" << std::endl;
return -1;
}
std::cerr << "Linear solve succeeded" << std::endl;
return 0;
}

58
Solver_interface/include/CGAL/Default_diagonalize_traits.h Normal file
View File

@ -0,0 +1,58 @@
#ifndef CGAL_DEFAULT_DIAGONALIZE_TRAITS_H
#define CGAL_DEFAULT_DIAGONALIZE_TRAITS_H
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_diagonalize_traits.h>
#else
#include <CGAL/Diagonalize_traits.h>
#endif
namespace CGAL {
// Wrapper class designed to automatically use
// Eigen_diagonalize_traits if Eigen is available and otherwise use
// the fallback Diagonalize_traits class.
template <typename FT, unsigned int dim = 3>
class Default_diagonalize_traits{
#ifdef CGAL_EIGEN3_ENABLED
typedef Eigen_diagonalize_traits<FT, dim> Base;
#else
typedef Diagonalize_traits<FT, dim> Base;
#endif
public:
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues)
{
return Base::diagonalize_selfadjoint_covariance_matrix (cov, eigenvalues);
}
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues,
cpp11::array<FT, dim * dim>& eigenvectors)
{
return Base::diagonalize_selfadjoint_covariance_matrix (cov, eigenvalues, eigenvectors);
}
// Extract the eigenvector associated to the largest eigenvalue
static bool
extract_largest_eigenvector_of_covariance_matrix (
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT,dim> &normal)
{
return Base::extract_largest_eigenvector_of_covariance_matrix (cov, normal);
}
};
} // namespace CGAL
#endif // CGAL_DEFAULT_DIAGONALIZE_TRAITS_H

247
Solver_interface/include/CGAL/Diagonalize_traits.h Normal file
View File

@ -0,0 +1,247 @@
#ifndef CGAL_DIAGONALIZE_TRAITS_H
#define CGAL_DIAGONALIZE_TRAITS_H
#include <CGAL/array.h>
namespace CGAL {
/// A model of the concept `DiagonalizeTraits`
/// \cgalModels `DiagonalizeTraits`
template <typename FT, unsigned int dim = 3>
class Diagonalize_traits{
public:
static bool
diagonalize_selfadjoint_covariance_matrix
(const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues)
{
cpp11::array<FT, dim * dim> eigenvectors;
return diagonalize_selfadjoint_covariance_matrix (cov, eigenvalues, eigenvectors);
}
// Extract the eigenvector associated to the largest eigenvalue
static bool
extract_largest_eigenvector_of_covariance_matrix
(const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT,dim> &normal)
{
cpp11::array<FT, dim> eigenvalues;
cpp11::array<FT, dim * dim> eigenvectors;
diagonalize_selfadjoint_covariance_matrix (cov, eigenvalues, eigenvectors);
for (std::size_t i = 0; i < dim; ++ i)
normal[i] = static_cast<FT> (eigenvectors[(dim*(dim-1))+i]);
return true;
}
static bool diagonalize_selfadjoint_covariance_matrix
(const cpp11::array<FT, (dim * (dim+1) / 2)>& mat,
cpp11::array<FT, dim>& eigen_values,
cpp11::array<FT, dim * dim>& eigen_vectors)
{
const int n = dim;
const int MAX_ITER = 100;
static const FT EPSILON = (FT)0.00001;
// number of entries in mat
int nn = (n*(n+1))/2;
// copy matrix
FT *a = new FT[nn];
int ij;
// This function requires lower triangular, so we switch
for (int i = 0; i < n; ++ i)
for (int j = i; j < n; ++ j)
a[(n * i) + j - ((i * (i+1)) / 2)]
= mat[i + (j * (j+1) / 2)];
// Fortran-porting
a--;
// init diagonalization matrix as the unit matrix
FT *v = new FT[n*n];
ij = 0;
int i;
for(i=0; i<n; i++)
for(int j=0; j<n; j++)
if(i==j)
v[ij++] = 1.0;
else
v[ij++] = 0.0;
// Fortran-porting
v--;
// compute weight of the non diagonal terms
ij = 1;
FT a_norm = 0.0;
for(i=1; i<=n; i++)
for(int j=1; j<=i; j++)
{
if( i!=j )
{
FT a_ij = a[ij];
a_norm += a_ij * a_ij;
}
ij++;
}
if(a_norm != 0.0)
{
FT a_normEPS = a_norm * EPSILON;
FT thr = a_norm;
// rotations
int nb_iter = 0;
while(thr > a_normEPS && nb_iter < MAX_ITER)
{
nb_iter++;
FT thr_nn = thr / nn;
for(int l=1; l< n; l++)
{
for(int m=l+1; m<=n; m++)
{
// compute sinx and cosx
int lq = (l*l-l)/2;
int mq = (m*m-m)/2;
int lm = l + mq;
FT a_lm = a[lm];
FT a_lm_2 = a_lm * a_lm;
if(a_lm_2 < thr_nn)
continue;
int ll = l + lq;
int mm = m + mq;
FT a_ll = a[ll];
FT a_mm = a[mm];
FT delta = a_ll - a_mm;
FT x;
if(delta == 0.0)
x = (FT) - CGAL_PI / 4;
else
x = (FT)(- std::atan( (a_lm+a_lm) / delta ) / 2.0);
FT sinx = std::sin(x);
FT cosx = std::cos(x);
FT sinx_2 = sinx * sinx;
FT cosx_2 = cosx * cosx;
FT sincos = sinx * cosx;
// rotate L and M columns
int ilv = n*(l-1);
int imv = n*(m-1);
int i;
for( i=1; i<=n;i++ )
{
if( (i!=l) && (i!=m) )
{
int iq = (i*i-i)/2;
int im;
if( i<m )
im = i + mq;
else
im = m + iq;
FT a_im = a[im];
int il;
if( i<l )
il = i + lq;
else
il = l + iq;
FT a_il = a[il];
a[il] = a_il * cosx - a_im * sinx;
a[im] = a_il * sinx + a_im * cosx;
}
ilv++;
imv++;
FT v_ilv = v[ilv];
FT v_imv = v[imv];
v[ilv] = cosx * v_ilv - sinx * v_imv;
v[imv] = sinx * v_ilv + cosx * v_imv;
}
x = a_lm * sincos;
x += x;
a[ll] = a_ll * cosx_2 + a_mm * sinx_2 - x;
a[mm] = a_ll * sinx_2 + a_mm * cosx_2 + x;
a[lm] = 0.0;
thr = CGAL::abs(thr - a_lm_2);
}
}
}
}
// convert indices and copy eigen values
a++;
for(i=0; i<n; i++)
{
int k = i + (i*(i+1))/2;
eigen_values[i] = a[k];
}
delete [] a;
// sort eigen values and vectors
int *index = new int[n];
for(i=0; i<n; i++)
index[i] = i;
for(i=0; i<(n-1); i++)
{
FT x = eigen_values[i];
int k = i;
for(int j=i+1; j<n; j++)
if(x > eigen_values[j])
{
k = j;
x = eigen_values[j];
}
eigen_values[k] = eigen_values[i];
eigen_values[i] = x;
int jj = index[k];
index[k] = index[i];
index[i] = jj;
}
// save eigen vectors
v++; // back to C++
ij = 0;
for(int k=0; k<n; k++ )
{
int ik = index[k]*n;
for(int i=0; i<n; i++)
eigen_vectors[ij++] = v[ik++];
}
delete [] v;
delete [] index;
return true;
}
};
} // namespace CGAL
#endif // CGAL_DIAGONALIZE_TRAITS_H

149
Solver_interface/include/CGAL/Eigen_diagonalize_traits.h Normal file
View File

@ -0,0 +1,149 @@
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL$
// $Id$
//
// Author(s) : Jocelyn Meyron and Quentin Mérigot
//
#ifndef CGAL_EIGEN_DIAGONALIZE_TRAITS_H
#define CGAL_EIGEN_DIAGONALIZE_TRAITS_H
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <CGAL/array.h>
namespace CGAL {
/// A model of the concept `DiagonalizeTraits` using \ref thirdpartyEigen.
/// \cgalModels `DiagonalizeTraits`
template <typename FT, unsigned int dim = 3>
class Eigen_diagonalize_traits{
typedef Eigen::Matrix<FT, dim, dim> Matrix;
typedef Eigen::Matrix<FT, dim, 1> Vector;
// Construct the covariance matrix
static Matrix
construct_covariance_matrix
(const cpp11::array<FT, (dim * (dim+1) / 2)>& cov) {
Matrix m;
for (std::size_t i = 0; i < dim; ++ i)
for (std::size_t j = i; j < dim; ++ j)
{
m(i,j) = static_cast<float>(cov[(dim * i) + j - ((i * (i+1)) / 2)]);
if (i != j)
m(j,i) = m(i,j);
}
return m;
}
// Diagonalize a selfadjoint matrix
static bool
diagonalize_selfadjoint_matrix (Matrix& m,
Matrix& eigenvectors,
Vector& eigenvalues) {
Eigen::SelfAdjointEigenSolver<Matrix> eigensolver(m);
if (eigensolver.info() != Eigen::Success) {
return false;
}
eigenvalues = eigensolver.eigenvalues();
eigenvectors = eigensolver.eigenvectors();
return true;
}
public:
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues)
{
Matrix m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
Vector eigenvalues_;
Matrix eigenvectors_;
bool res = diagonalize_selfadjoint_matrix(m, eigenvectors_, eigenvalues_);
if (res)
{
for (std::size_t i = 0; i < dim; ++ i)
eigenvalues[i] = static_cast<FT>(eigenvalues_[i]);
}
return res;
}
static bool
diagonalize_selfadjoint_covariance_matrix(
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT, dim>& eigenvalues,
cpp11::array<FT, dim * dim>& eigenvectors)
{
Matrix m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
Vector eigenvalues_;
Matrix eigenvectors_;
bool res = diagonalize_selfadjoint_matrix(m, eigenvectors_, eigenvalues_);
if (res)
{
for (std::size_t i = 0; i < dim; ++ i)
{
eigenvalues[i] = static_cast<FT>(eigenvalues_[i]);
for (std::size_t j = 0; j < dim; ++ j)
eigenvectors[dim*i + j]=static_cast<FT>(eigenvectors_(j,i));
}
}
return res;
}
// Extract the eigenvector associated to the largest eigenvalue
static bool
extract_largest_eigenvector_of_covariance_matrix (
const cpp11::array<FT, (dim * (dim+1) / 2)>& cov,
cpp11::array<FT,dim> &normal)
{
// Construct covariance matrix
Matrix m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
Vector eigenvalues;
Matrix eigenvectors;
if (! diagonalize_selfadjoint_matrix(m, eigenvectors, eigenvalues)) {
return false;
}
// Eigenvalues are sorted by increasing order
for (unsigned int i = 0; i < dim; ++ i)
normal[i] = static_cast<FT> (eigenvectors(i,dim-1));
return true;
}
};
} // namespace CGAL
#endif // CGAL_EIGEN_DIAGONALIZE_TRAITS_H

8
Solver_interface/include/CGAL/Eigen_solver_traits.h
View File

@ -69,9 +69,9 @@ namespace internal {
/// is a generic traits class for solving asymmetric or symmetric positive definite (SPD)
/// sparse linear systems using one of the Eigen solvers.
/// The default solver is the iterative bi-congugate gradient stabilized solver
/// Eigen::BiCGSTAB for double.
/// <a href="http://eigen.tuxfamily.org/dox/classEigen_1_1BiCGSTAB.html">Eigen::BiCGSTAB</a> for double.
///
/// \cgalModels `SparseLinearAlgebraTraitsWithFactor_d`.
/// \cgalModels `SparseLinearAlgebraWithFactorTraits_d`.
template<class EigenSolverT = Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType> >
class Eigen_solver_traits
@ -79,6 +79,7 @@ class Eigen_solver_traits
typedef typename EigenSolverT::Scalar Scalar;
// Public types
public:
typedef EigenSolverT Solver;
typedef Scalar NT;
typedef typename internal::Get_eigen_matrix<EigenSolverT,NT>::type Matrix;
typedef Eigen_vector<Scalar> Vector;
@ -159,7 +160,7 @@ protected:
};
//specilization of the solver for BiCGSTAB as for surface parameterization, the
//specialization of the solver for BiCGSTAB as for surface parameterization, the
//intializer should be a vector of one's (this was the case in 3.1-alpha but not in the official 3.1).
template<>
class Eigen_solver_traits< Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType> >
@ -168,6 +169,7 @@ class Eigen_solver_traits< Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenTyp
typedef EigenSolverT::Scalar Scalar;
// Public types
public:
typedef EigenSolverT Solver;
typedef Scalar NT;
typedef internal::Get_eigen_matrix<EigenSolverT,NT>::type Matrix;
typedef Eigen_vector<Scalar> Vector;

0
Jet_fitting_3/include/CGAL/Lapack/Linear_algebra_lapack.h → Solver_interface/include/CGAL/Lapack_svd.h
View File

27
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/PackageDescription.txt
View File

@ -12,7 +12,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{3.2}
\cgalPkgDependsOn{Solvers as \ref thirdpartyEigen.}
\cgalPkgDependsOn{\ref PkgSolverSummary}
\cgalPkgBib{cgal:sal-pptsm2}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{Operations on Polyhedra,polyhedron_3.zip}
@ -31,7 +31,6 @@
- `BorderParameterizer_3`
- `ParameterizationMesh_3`
- `ParameterizationPatchableMesh_3`
- `SparseLinearAlgebraTraits_d`
- `PolyhedronTraitsWithKernel_3`
## Surface Parameterization Methods ##
@ -55,7 +54,7 @@ This \cgal package implements several parameterization methods:
- `CGAL::Barycentric_mapping_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
- `CGAL::Discrete_authalic_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
- `CGAL::Discrete_conformal_map_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
- `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
- `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3>`
- `CGAL::Mean_value_coordinates_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
## Border Parameterization Methods ##
@ -122,16 +121,6 @@ The package provides an interface with `CGAL::Polyhedron_3<Traits>`:
A `(u,v)` pair is computed for each inner vertex (i.e.\ its halfedges share the same `(u,v)` pair), while a `(u,v)` pair is computed for each border halfedge. The user must iterate over the mesh halfedges to get the result.
## Sparse Linear Algebra ##
Since parameterizing meshes requires
efficient representation of sparse matrices and efficient iterative or
direct linear solvers, we provide an interface to several
sparse linear solvers:
- <span class="textsc">Eigen</span> 3.1 (or greater) is the library recommended by %CGAL solving sparse systems.
- OpenNL (authored by Bruno Lévy) is shipped with %CGAL is the default solver.
- `OpenNL::DefaultLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` in OpenNL package
- `OpenNL::SymmetricLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` in OpenNL package
## Helper Classes ##
@ -250,18 +239,6 @@ The package provides an interface with `CGAL::Polyhedron_3<Traits>`:
/*! \defgroup PkgSurfaceParameterizationAlgebra Sparse Linear Algebra
\ingroup PkgSurfaceParameterization
Since parameterizing meshes requires
efficient representation of sparse matrices and efficient iterative or
direct linear solvers, we provide an interface to several
sparse linear solvers:
- <span class="textsc">Eigen</span> 3.1 (or greater) is the library recommended by %CGAL solving sparse systems.
- OpenNL (authored by Bruno Lévy) is shipped with %CGAL is the default solver.
- `OpenNL::DefaultLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` in OpenNL package
- `OpenNL::SymmetricLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` in OpenNL package
*/
/*!
\defgroup PkgSurfaceParameterizationHelper Helper Class

153
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/Surface_mesh_parameterization.txt
View File

@ -40,9 +40,8 @@ data structure.
Since parameterizing meshes require efficient representation of sparse
matrices and efficient iterative or direct linear solvers, we provide
a unified interface to a linear solver library (\ref thirdpartyEigen "Eigen"),
and propose a separate package devoted to OpenNL sparse
linear solver.
a unified interface to linear solvers as described in Chapter
\ref PkgSolverSummary.
Note that linear solvers commonly use double precision floating point
numbers. Therefore, this package is intended to be used with a \cgal %Cartesian kernel with doubles.
@ -69,8 +68,8 @@ Parameterizer_traits_3<ParameterizationMesh_3>::Error_code parameterize (Paramet
The function `parameterize()` applies a default surface parameterization
method, namely Floater Mean Value Coordinates \cgalCite{cgal:f-mvc-03}, with an
arc-length circular border parameterization, and using OpenNL sparse
linear solver \cgalCite{cgal:l-nmdgp-05}. The `ParameterizationMesh_3` concept defines the input meshes handled by `parameterize()`. See Section \ref secInputMeshforparameterize. The result is stored into the `(u,v)` fields of the mesh halfedges.
arc-length circular border parameterization, and using a sparse
linear solver from the \ref thirdpartyEigen library. The `ParameterizationMesh_3` concept defines the input meshes handled by `parameterize()`. See Section \ref secInputMeshforparameterize. The result is stored into the `(u,v)` fields of the mesh halfedges.
The mesh must be a triangle mesh surface with one connected component.
Note: `Parameterizer_traits_3<ParameterizationMesh_3>` is the (pure virtual) superclass of all surface parameterizations and defines the error codes.
@ -155,7 +154,7 @@ and \ref secBorderParameterizationsforFreeMethods.
\subsection Surface_mesh_parameterizationTheSparseLinearAlgebraTraitsd The SparseLinearAlgebraTraits_d Concept
This package solves sparse linear systems using solvers which are models
of `SparseLinearAlgebraTraits_d`. See Section \ref secSparseLinearAlgebra.
of `SparseLinearAlgebraTraits_d`. See Chapter \ref PkgSolverSummary.
\subsection Surface_mesh_parameterizationTheParameterizationMesh3 The ParameterizationMesh_3 and ParameterizationPatchableMesh_3 Concepts
@ -419,65 +418,7 @@ Parameterizer::Error_code err = CGAL::parameterize(mesh_adaptor, Parameterizer()
\endcode
\section secSparseLinearAlgebra Sparse Linear Algebra
Parameterizing triangle meshes requires both efficient representation
of sparse matrices and efficient iterative or direct linear
solvers. We provide links to \ref thirdpartyEigen "Eigen" library
and include a separate package devoted to OpenNL sparse linear solver.
\subsection Surface_mesh_parameterizationListofSolvers List of Solvers
We provide an interface to several sparse linear solvers, as models
of the `SparseLinearAlgebraTraits_d` concept:
- An interface to sparse solvers from the `OpenNL` library \cgalCite{cgal:l-nmdgp-05} is provided through classes
`OpenNL::DefaultLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` and
`OpenNL::SymmetricLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>`. The OpenNL library version shipped with \cgal is a lightweight default sparse linear solver. It does not support large systems, but it is portable and
supports exact number types.
- An interface to all sparse solvers from the \ref thirdpartyEigen "Eigen" library is provided through the class
`Eigen_solver_traits<T>`. This solver traits class can be used for an iterative or a direct,
symmetric or general sparse solvers. The \ref thirdpartyEigen "Eigen" solver to be used must be given as template parameter.
\subsection Surface_mesh_parameterizationEigenSolver Eigen Solver Example
The example \ref Surface_mesh_parameterization/Eigen_parameterization.cpp "Eigen_parameterization.cpp" computes the
default parameterization method (Floater mean value coordinates with a circular border),
but specifically instantiates an \ref thirdpartyEigen "Eigen" solver. Specifying a specific solver
instead of the default one (OpenNL) means using the third parameter of
`Mean_value_coordinates_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`. The differences with the first
example \ref Surface_mesh_parameterization/Simple_parameterization.cpp "Simple_parameterization.cpp" are:
\code{.cpp}
#include <CGAL/Eigen_solver_traits.h>
...
//***************************************
// Floater Mean Value Coordinates parameterization
// (circular border) with Eigen solver
//***************************************
// Circular border parameterizer (the default)
typedef CGAL::Circular_border_arc_length_parameterizer_3<Parameterization_polyhedron_adaptor>
Border_parameterizer;
// Eigen solver
typedef CGAL::Eigen_solver_traits<> Solver;
// Floater Mean Value Coordinates parameterization
// (circular border) with Eigen solver
typedef CGAL::Mean_value_coordinates_parameterizer_3<Parameterization_polyhedron_adaptor,
Border_parameterizer,
Solver>
Parameterizer;
Parameterizer::Error_code err = CGAL::parameterize(mesh_adaptor, Parameterizer());
...
\endcode
\section secCuttingaMesh Cutting a Mesh
@ -605,63 +546,7 @@ and is therefore non-singular (Gram theorem).
</UL>
\subsection Surface_mesh_parameterizationPrecision Precision
Two algorithms of this package construct the sparse linear system(s)
using trigonometric functions, and are this incompatible with exact arithmetic:
<UL>
<LI>Floater mean value coordinates
<LI>Circular border parameterization
</UL>
On the other hand, linear solvers commonly use double precision floating point
numbers.
OpenNL's BICGSTAB solver (accessible through the
`OpenNL::DefaultLinearSolverTraits<COEFFTYPE, MATRIX, VECTOR, SOLVER>` interface)
is the only solver supported by this package which
computes exact results, when used with an exact arithmetic. This package is
intended to be used mainly with a \cgal %Cartesian kernel with doubles.
\subsection Surface_mesh_parameterizationOpenNLsBICGSTAB OpenNL's BICGSTAB Solver with an Exact Arithmetic
The BICGSTAB conjugate gradient is in disguise a direct solver.
In a nutshell, it computes a vector basis
orthogonal with respect to the matrix, and the coordinates of the solution in this vector basis.
Each iteration computes one component of the basis and one coordinate, therefore the algorithm
converges to the solution in \f$ n\f$ iterations, where \f$ n\f$ is the dimension of the matrix. More precisely, it is shown to converge in \f$ k\f$ iteration, where \f$ k\f$ is the number of distinct eigenvalues of the matrix.
\subsection Surface_mesh_parameterizationSolverswith Solvers with a Floating Point Arithmetic
<I>OpenNL's BICGSTAB example:</I>
When inexact numerical types are used (e.g. doubles), accumulated errors slow down convergence
(in practice, it requires approximately \f$ 5k\f$ iterations to converge).
The required number of iterations depends on the eigenvalues of the matrix, and these eigenvalues depend
on the shape of the triangles. The optimum is when the triangles are equilateral (then the solver converges
in less than 10 iterations). The worst case is obtained when the mesh has a large number of skinny triangles (near-singular Jacobian matrix of the triangle). In this case, the spectrum of the matrix
is wide (many different eigenvalues), and the solver requires nearly \f$ 5n\f$ iterations to converge.
\subsection Surface_mesh_parameterizationAlgorithmic Algorithmic Complexity
In this package, we focus on piecewise linear mappings onto a planar
domain. All surface parameterization methods are based on solving one (or two)
sparse linear system(s).
The algorithmic complexity is dominated by the resolution of the sparse linear system(s).
<I>OpenNL's BICGSTAB example:</I>
At each iteration, the operation of highest complexity is the product between the sparse-matrix and a vector.
The sparse matrix has a fixed number of non-zero coefficients per row,
therefore the matrix / vector product has \f$ O(n)\f$ complexity.
Since convergence is reached after \f$ k\f$ iterations, the complexity is \f$ O(k.n)\f$
(where \f$ k\f$ is the number of distinct eigenvalues of the matrix).
Therefore, best case complexity is \f$ O(n)\f$ (equilateral triangles),
and worst case complexity is \f$ O(n^2)\f$ (skinny triangles).
\section Surface_mesh_parameterizationSoftware Software Design
@ -801,7 +686,7 @@ The purpose of such a model is to:
<LI>Support several kind of meshes.
<LI>Hide the implementation of extra fields specific to the parameterization domain
(`index`, `u`, `v`, `is_parameterized`).
<LI>Handle in the mesh type the complexity of <I>virtually</I> cutting a mesh
<LI>%Handle in the mesh type the complexity of <I>virtually</I> cutting a mesh
to make it homeomorphic to a disk (instead of duplicating this
code in each parameterization method).
</OL>
@ -834,15 +719,6 @@ This mainly means that:
can be set <I>per corner</I> (which is a more general way of saying <I>per half-edge</I>).
</UL>
\subsection Surface_mesh_parameterizationSparseLinearAlgebraTraitsd SparseLinearAlgebraTraits_d Concept
This package solves sparse linear systems using solvers which are models
of `SparseLinearAlgebraTraits_d`.
`SparseLinearAlgebraTraits_d` is a sub-concept of the `LinearAlgebraTraits_d` concept
in `Kernel_d`.
The goal is to adapt easily code written for dense matrices to sparse ones,
and vice-versa.
\subsection Surface_mesh_parameterizationCuttingaMesh_1 Cutting a Mesh
@ -855,23 +731,7 @@ not intend to cover this topic at the moment.
\section Surface_mesh_parameterizationExtendingthe Extending the Package and Reusing Code
\subsection Surface_mesh_parameterizationReusingMesh Reusing Mesh Adaptors
`ParameterizationMesh_3` defines a concept to access to a
general polyhedral mesh.
It is optimized for the `Surface_mesh_parameterization` package
only in the sense that it
defines the accessors to fields specific to the parameterization domain
(`index`, `u`, `v`, `is_parameterized`).
It may be easily generalized.
\subsection Surface_mesh_parameterizationReusingSparse Reusing Sparse Linear Algebra
The `SparseLinearAlgebraTraits_d` concept and the traits classes
for \ref thirdpartyEigen "Eigen" and OpenNL are independent of the rest of the
`Surface_mesh_parameterization` package, and may be reused by
\cgal developers for other purposes.
\subsection Surface_mesh_parameterizationAddingNewParameterization Adding New Parameterization Methods
@ -900,6 +760,7 @@ Square, circular and two-points border parameterizations are good starting point
Obviously, this package would benefit of having robust algorithms
which transform arbitrary meshes into topological disks.
*/
} /* namespace CGAL */

5
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/dependencies
View File

@ -6,5 +6,6 @@ Circulator
Stream_support
Polyhedron
Kernel_d
Surface_modeling
Surface_mesh_skeletonization
Solver_interface
Convex_hull_2

11
Surface_mesh_parameterization/doc/Surface_mesh_parameterization/examples.txt
View File

@ -1,5 +1,6 @@
/// \example Surface_mesh_parameterization/Simple_parameterization.cpp
/// \example Surface_mesh_parameterization/Mesh_cutting_parameterization.cpp
/// \example Surface_mesh_parameterization/Complete_parameterization_example.cpp
/// \example Surface_mesh_parameterization/Eigen_parameterization.cpp
/// \example Surface_mesh_parameterization/polyhedron_ex_parameterization.cpp
/*!
\example Surface_mesh_parameterization/Simple_parameterization.cpp
\example Surface_mesh_parameterization/Mesh_cutting_parameterization.cpp
\example Surface_mesh_parameterization/Complete_parameterization_example.cpp
\example Surface_mesh_parameterization/polyhedron_ex_parameterization.cpp
*/

3
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Authalic_parameterization.cpp
View File

@ -27,7 +27,6 @@ int main(int argc, char * argv[])
std::cerr << "PARAMETERIZATION" << std::endl;
std::cerr << " Discrete Authalic Parameterization" << std::endl;
std::cerr << " circle border" << std::endl;
std::cerr << " OpenNL solver" << std::endl;
//***************************************
// decode parameters
@ -68,7 +67,7 @@ int main(int argc, char * argv[])
//***************************************
// Discrete Authalic Parameterization
// (defaults are circular border and OpenNL solver)
// (defaults are circular border and Eigen solver)
//***************************************
typedef CGAL::Discrete_authalic_parameterizer_3<Parameterization_polyhedron_adaptor>

18
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/CMakeLists.txt
View File

@ -62,25 +62,25 @@ if ( CGAL_FOUND )
list(APPEND CGAL_3RD_PARTY_LIBRARIES ${Boost_PROGRAM_OPTIONS_LIBRARY})
endif()
# Executables that do *not* require Eigen
find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater)
if (EIGEN3_FOUND)
# Executables that require Eigen 3.1
include( ${EIGEN3_USE_FILE} )
create_single_source_cgal_program( "Authalic_parameterization.cpp" )
create_single_source_cgal_program( "Mesh_cutting_parameterization.cpp" )
create_single_source_cgal_program( "Simple_parameterization.cpp" )
create_single_source_cgal_program( "Square_border_parameterization.cpp" )
find_package(Eigen3 3.1.0) #(requires 3.1.0 or greater)
if (EIGEN3_FOUND)
# Executables that require Eigen 3.1
include( ${EIGEN3_USE_FILE} )
create_single_source_cgal_program( "Complete_parameterization_example.cpp" )
create_single_source_cgal_program( "Eigen_parameterization.cpp" )
create_single_source_cgal_program( "polyhedron_ex_parameterization.cpp" )
else(EIGEN3_FOUND)
message(STATUS "NOTICE: Some examples require Eigen 3.1 (or greater) and will not be compiled.")
endif(EIGEN3_FOUND)
create_single_source_cgal_program( "polyhedron_ex_parameterization.cpp" )
else()
message(STATUS "NOTICE: This program requires the CGAL library, and will not be compiled.")

6
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Complete_parameterization_example.cpp
View File

@ -248,14 +248,10 @@ int main(int argc, char * argv[])
// Border parameterizer
typedef CGAL::Square_border_arc_length_parameterizer_3<Mesh_patch_polyhedron>
Border_parameterizer;
// Eigen solver
typedef CGAL::Eigen_solver_traits<> Solver;
// Discrete Authalic Parameterization (square border)
// with Eigen solver
typedef CGAL::Discrete_authalic_parameterizer_3<Mesh_patch_polyhedron,
Border_parameterizer,
Solver> Parameterizer;
Border_parameterizer> Parameterizer;
Parameterizer::Error_code err = CGAL::parameterize(mesh_patch, Parameterizer());
switch(err) {

1
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cmd
View File

@ -1 +0,0 @@
data/mannequin-devil.off

124
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Eigen_parameterization.cpp
View File

@ -1,124 +0,0 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Parameterization_polyhedron_adaptor_3.h>
#include <CGAL/parameterize.h>
#include <CGAL/Eigen_solver_traits.h>
#include <iostream>
#include <fstream>
#include <cstdlib>
// ----------------------------------------------------------------------------
// Private types
// ----------------------------------------------------------------------------
typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
// ----------------------------------------------------------------------------
// main()
// ----------------------------------------------------------------------------
int main(int argc, char * argv[])
{
std::cerr << "PARAMETERIZATION" << std::endl;
std::cerr << " Floater parameterization" << std::endl;
std::cerr << " Circle border" << std::endl;
std::cerr << " Eigen solver" << std::endl;
//***************************************
// decode parameters
//***************************************
if (argc-1 != 1)
{
std::cerr << "Usage: " << argv[0] << " input_file.off" << std::endl;
return(EXIT_FAILURE);
}
// File name is:
const char* input_filename = argv[1];
//***************************************
// Read the mesh
//***************************************
// Read the mesh
std::ifstream stream(input_filename);
Polyhedron mesh;
stream >> mesh;
if(!stream || !mesh.is_valid() || mesh.empty())
{
std::cerr << "Error: cannot read OFF file " << input_filename << std::endl;
return EXIT_FAILURE;
}
//***************************************
// Create Polyhedron adaptor
// Note: no cutting => we support only
// meshes that are topological disks
//***************************************
typedef CGAL::Parameterization_polyhedron_adaptor_3<Polyhedron>
Parameterization_polyhedron_adaptor;
Parameterization_polyhedron_adaptor mesh_adaptor(mesh);
//***************************************
// Floater Mean Value Coordinates parameterization
// (circular border) with Eigen solver
//***************************************
// Circular border parameterizer (the default)
typedef CGAL::Circular_border_arc_length_parameterizer_3<Parameterization_polyhedron_adaptor>
Border_parameterizer;
// Eigen solver
typedef CGAL::Eigen_solver_traits<> Solver;
// Floater Mean Value Coordinates parameterization
// (circular border) with Eigen solver
typedef CGAL::Mean_value_coordinates_parameterizer_3<Parameterization_polyhedron_adaptor,
Border_parameterizer,
Solver>
Parameterizer;
Parameterizer::Error_code err = CGAL::parameterize(mesh_adaptor, Parameterizer());
switch(err) {
case Parameterizer::OK: // Success
break;
case Parameterizer::ERROR_EMPTY_MESH: // Input mesh not supported
case Parameterizer::ERROR_NON_TRIANGULAR_MESH:
case Parameterizer::ERROR_NO_TOPOLOGICAL_DISC:
case Parameterizer::ERROR_BORDER_TOO_SHORT:
std::cerr << "Input mesh not supported: " << Parameterizer::get_error_message(err) << std::endl;
return EXIT_FAILURE;
break;
default: // Error
std::cerr << "Error: " << Parameterizer::get_error_message(err) << std::endl;
return EXIT_FAILURE;
break;
};
//***************************************
// Output
//***************************************
// Raw output: dump (u,v) pairs
Polyhedron::Vertex_const_iterator pVertex;
for (pVertex = mesh.vertices_begin();
pVertex != mesh.vertices_end();
pVertex++)
{
// (u,v) pair is stored in any halfedge
double u = mesh_adaptor.info(pVertex->halfedge())->uv().x();
double v = mesh_adaptor.info(pVertex->halfedge())->uv().y();
std::cout << "(u,v) = (" << u << "," << v << ")" << std::endl;
}
return EXIT_SUCCESS;
}

1
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Mesh_cutting_parameterization.cpp
View File

@ -93,7 +93,6 @@ int main(int argc, char * argv[])
std::cerr << "PARAMETERIZATION" << std::endl;
std::cerr << " Floater parameterization" << std::endl;
std::cerr << " Circle border" << std::endl;
std::cerr << " OpenNL solver" << std::endl;
std::cerr << " Very simple cut if model is not a topological disk" << std::endl;
//***************************************

3
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Simple_parameterization.cpp
View File

@ -25,7 +25,6 @@ int main(int argc, char * argv[])
std::cerr << "PARAMETERIZATION" << std::endl;
std::cerr << " Floater parameterization" << std::endl;
std::cerr << " Circle border" << std::endl;
std::cerr << " OpenNL solver" << std::endl;
//***************************************
// decode parameters
@ -66,7 +65,7 @@ int main(int argc, char * argv[])
//***************************************
// Floater Mean Value Coordinates parameterization
// (defaults are circular border and OpenNL solver)
// (defaults are circular border and Eigen solver)
//***************************************
typedef CGAL::Parameterizer_traits_3<Parameterization_polyhedron_adaptor>

1
Surface_mesh_parameterization/examples/Surface_mesh_parameterization/Square_border_parameterization.cpp
View File

@ -26,7 +26,6 @@ int main(int argc, char * argv[])
std::cerr << "PARAMETERIZATION" << std::endl;
std::cerr << " Floater parameterization" << std::endl;
std::cerr << " square border" << std::endl;
std::cerr << " OpenNL solver" << std::endl;
//***************************************
// decode parameters

5
Surface_mesh_parameterization/include/CGAL/Barycentric_mapping_parameterizer_3.h
View File

@ -24,6 +24,7 @@
#include <CGAL/Fixed_border_parameterizer_3.h>
#include <CGAL/surface_mesh_parameterization_assertions.h>
#include <CGAL/Eigen_solver_traits.h>
/// \file Barycentric_mapping_parameterizer_3.h
@ -58,7 +59,7 @@ namespace CGAL {
\sa `CGAL::Fixed_border_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
\sa `CGAL::Discrete_authalic_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
\sa `CGAL::Discrete_conformal_map_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
\sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
\sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3>`
\sa `CGAL::Mean_value_coordinates_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
*/
@ -68,7 +69,7 @@ template
class BorderParameterizer_3
= Circular_border_arc_length_parameterizer_3<ParameterizationMesh_3>,
class SparseLinearAlgebraTraits_d
= OpenNL::DefaultLinearSolverTraits<typename ParameterizationMesh_3::NT>
= Eigen_solver_traits<Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT< double > > >
>
class Barycentric_mapping_parameterizer_3
: public Fixed_border_parameterizer_3<ParameterizationMesh_3,

5
Surface_mesh_parameterization/include/CGAL/Discrete_authalic_parameterizer_3.h
View File

@ -24,6 +24,7 @@
#include <CGAL/Fixed_border_parameterizer_3.h>
#include <CGAL/surface_mesh_parameterization_assertions.h>
#include <CGAL/Eigen_solver_traits.h>
/// \file Discrete_authalic_parameterizer_3.h
@ -55,7 +56,7 @@ namespace CGAL {
/// \sa `CGAL::Fixed_border_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::Barycentric_mapping_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::Discrete_conformal_map_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3>`
/// \sa `CGAL::Mean_value_coordinates_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
template
@ -64,7 +65,7 @@ template
class BorderParameterizer_3 ///< Strategy to parameterize the surface border
= Circular_border_arc_length_parameterizer_3<ParameterizationMesh_3>,
class SparseLinearAlgebraTraits_d ///< Traits class to solve a sparse linear system
= OpenNL::DefaultLinearSolverTraits<typename ParameterizationMesh_3::NT>
= Eigen_solver_traits<Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT< double > > >
>
class Discrete_authalic_parameterizer_3
: public Fixed_border_parameterizer_3<ParameterizationMesh_3,

5
Surface_mesh_parameterization/include/CGAL/Discrete_conformal_map_parameterizer_3.h
View File

@ -24,6 +24,7 @@
#include <CGAL/Fixed_border_parameterizer_3.h>
#include <CGAL/surface_mesh_parameterization_assertions.h>
#include <CGAL/Eigen_solver_traits.h>
/// \file Discrete_conformal_map_parameterizer_3.h
@ -61,7 +62,7 @@ namespace CGAL {
/// \sa `CGAL::Fixed_border_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::Barycentric_mapping_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::Discrete_authalic_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3>`
/// \sa `CGAL::Mean_value_coordinates_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
template
@ -70,7 +71,7 @@ template
class BorderParameterizer_3
= Circular_border_arc_length_parameterizer_3<ParameterizationMesh_3>,
class SparseLinearAlgebraTraits_d
= OpenNL::DefaultLinearSolverTraits<typename ParameterizationMesh_3::NT>
= Eigen_solver_traits<Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT< double > > >
>
class Discrete_conformal_map_parameterizer_3
: public Fixed_border_parameterizer_3<ParameterizationMesh_3,

7
Surface_mesh_parameterization/include/CGAL/Fixed_border_parameterizer_3.h
View File

@ -24,7 +24,8 @@
#include <CGAL/circulator.h>
#include <CGAL/Timer.h>
#include <CGAL/OpenNL/linear_solver.h>
#include <CGAL/Eigen_solver_traits.h>
#include <CGAL/Parameterizer_traits_3.h>
#include <CGAL/Circular_border_parameterizer_3.h>
@ -80,8 +81,8 @@ template
class BorderParameterizer_3 ///< Strategy to parameterize the surface border
= Circular_border_arc_length_parameterizer_3<ParameterizationMesh_3>,
class SparseLinearAlgebraTraits_d ///< Traits class to solve a sparse linear system
= OpenNL::DefaultLinearSolverTraits<typename ParameterizationMesh_3::NT>
>
= Eigen_solver_traits<Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT< double > > > >
class Fixed_border_parameterizer_3
: public Parameterizer_traits_3<ParameterizationMesh_3>
{

25
Surface_mesh_parameterization/include/CGAL/LSCM_parameterizer_3.h
View File

@ -43,15 +43,17 @@ namespace CGAL {
/// \ingroup PkgSurfaceParameterizationMethods
///
/// The class LSCM_parameterizer_3 implements the
/// The class `LSCM_parameterizer_3` implements the
/// *Least Squares Conformal Maps (LSCM)* parameterization \cgalCite{cgal:lprm-lscm-02}.
///
/// This is a conformal parameterization, i.e. it attempts to preserve angles.
///
/// This is a free border parameterization. No need to map the surface's border
/// This is a free border parameterization. No need to map the border of the surface
/// onto a convex polygon (only two pinned vertices are needed to ensure a
/// unique solution), but one-to-one mapping is *not* guaranteed.
///
/// Note that his parametrization method has no template parameter for changing the sparse solver.
///
/// \cgalModels `ParameterizerTraits_3`
///
///
@ -66,14 +68,17 @@ template
<
class ParameterizationMesh_3, ///< 3D surface mesh.
class BorderParameterizer_3
= Two_vertices_parameterizer_3<ParameterizationMesh_3>,
= Two_vertices_parameterizer_3<ParameterizationMesh_3>
#ifndef DOXYGEN_RUNNING
,
///< Strategy to parameterize the surface border.
///< The minimum is to parameterize two vertices.
class SparseLinearAlgebraTraits_d
class SparseLinearAlgebraTraits_d
= OpenNL::SymmetricLinearSolverTraits<typename ParameterizationMesh_3::NT>
///< Traits class to solve a sparse linear system.
///< We may use a symmetric definite positive solver because LSCM
///< solves the system in the least squares sense.
//< Traits class to solve a sparse linear system.
//< We may use a symmetric definite positive solver because LSCM
//< solves the system in the least squares sense.
#endif
>
class LSCM_parameterizer_3
: public Parameterizer_traits_3<ParameterizationMesh_3>
@ -94,11 +99,11 @@ public:
/// Export BorderParameterizer_3 template parameter.
typedef BorderParameterizer_3 Border_param;
/// Export SparseLinearAlgebraTraits_d template parameter.
// Private types
private:
typedef SparseLinearAlgebraTraits_d Sparse_LA;
// Private types
private:
// Mesh_Adaptor_3 subtypes:
typedef typename Adaptor::NT NT;
typedef typename Adaptor::Point_2 Point_2;

6
Surface_mesh_parameterization/include/CGAL/Mean_value_coordinates_parameterizer_3.h
View File

@ -24,6 +24,7 @@
#include <CGAL/Fixed_border_parameterizer_3.h>
#include <CGAL/surface_mesh_parameterization_assertions.h>
#include <CGAL/Eigen_solver_traits.h>
/// \file Mean_value_coordinates_parameterizer_3.h
@ -56,7 +57,7 @@ namespace CGAL {
/// \sa `CGAL::Barycentric_mapping_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::Discrete_authalic_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::Discrete_conformal_map_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3, SparseLinearAlgebraTraits_d>`
/// \sa `CGAL::LSCM_parameterizer_3<ParameterizationMesh_3, BorderParameterizer_3>`
template
<
@ -64,8 +65,9 @@ template
class BorderParameterizer_3 ///< Strategy to parameterize the surface border
= Circular_border_arc_length_parameterizer_3<ParameterizationMesh_3>,
class SparseLinearAlgebraTraits_d ///< Traits class to solve a sparse linear system
= OpenNL::DefaultLinearSolverTraits<typename ParameterizationMesh_3::NT>
= Eigen_solver_traits<Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType, Eigen::IncompleteLUT< double > > >
>
class Mean_value_coordinates_parameterizer_3
: public Fixed_border_parameterizer_3<ParameterizationMesh_3,
BorderParameterizer_3,

7
Surface_mesh_parameterization/test/Surface_mesh_parameterization/extensive_parameterization_test.cpp
View File

@ -407,10 +407,9 @@ int main(int argc, char * argv[])
err = CGAL::parameterize(
mesh_patch,
CGAL::Discrete_authalic_parameterizer_3<
Mesh_patch_polyhedron,
CGAL::Square_border_arc_length_parameterizer_3<Mesh_patch_polyhedron>,
CGAL::Eigen_solver_traits<>
>());
Mesh_patch_polyhedron,
CGAL::Square_border_arc_length_parameterizer_3<Mesh_patch_polyhedron> >()
);
switch(err) {
case Parameterizer::OK: // Success
break;

2
Surface_mesh_skeletonization/doc/Surface_mesh_skeletonization/PackageDescription.txt
View File

@ -18,7 +18,7 @@
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{4.7}
\cgalPkgDependsOn{Sparse symmetric solver such as those from \ref thirdpartyEigen}
\cgalPkgDependsOn{ \ref PkgSolverSummary}
\cgalPkgBib{cgal:glt-tsms}
\cgalPkgLicense{\ref licensesGPL "GPL"}
\cgalPkgDemo{See Polyhedral Surface,polyhedron_3.zip}

Some files were not shown because too many files have changed in this diff Show More