/usr/include/CGAL/Eigen_diagonalize_traits.h is in libcgal-dev 4.11-2build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 | // 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 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) : Jocelyn Meyron and Quentin Mérigot
//
#ifndef CGAL_EIGEN_DIAGONALIZE_TRAITS_H
#define CGAL_EIGEN_DIAGONALIZE_TRAITS_H
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
// If the matrix to diagonalize is of dimension 2x2 or 3x3, Eigen
// provides a faster implementation using a closed-form
// algorithm. However, it offers less precision. See:
// https://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html
// This is usually acceptable for CGAL algorithms but one might want
// to use the slower but more accurate version. In that case, just
// uncomment the following line:
//#define DO_NOT_USE_EIGEN_COMPUTEDIRECT_FOR_DIAGONALIZATION
#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{
public:
typedef cpp11::array<FT, dim> Vector;
typedef cpp11::array<FT, dim*dim> Matrix;
typedef cpp11::array<FT, (dim * (dim+1) / 2)> Covariance_matrix;
private:
typedef Eigen::Matrix<FT, dim, dim> EigenMatrix;
typedef Eigen::Matrix<FT, dim, 1> EigenVector;
// Construct the covariance matrix
static EigenMatrix
construct_covariance_matrix
(const Covariance_matrix& cov) {
EigenMatrix 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 (EigenMatrix& m,
EigenMatrix& eigenvectors,
EigenVector& eigenvalues) {
Eigen::SelfAdjointEigenSolver<EigenMatrix> eigensolver;
#ifndef DO_NOT_USE_EIGEN_COMPUTEDIRECT_FOR_DIAGONALIZATION
if (dim == 2 || dim == 3)
eigensolver.computeDirect(m);
else
#endif
eigensolver.compute(m);
if (eigensolver.info() != Eigen::Success) {
return false;
}
eigenvalues = eigensolver.eigenvalues();
eigenvectors = eigensolver.eigenvectors();
return true;
}
public:
static bool
diagonalize_selfadjoint_covariance_matrix(
const Covariance_matrix& cov,
Vector& eigenvalues)
{
EigenMatrix m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
EigenVector eigenvalues_;
EigenMatrix 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 Covariance_matrix& cov,
Vector& eigenvalues,
Matrix& eigenvectors)
{
EigenMatrix m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
EigenVector eigenvalues_;
EigenMatrix 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 Covariance_matrix& cov,
Vector& normal)
{
// Construct covariance matrix
EigenMatrix m = construct_covariance_matrix(cov);
// Diagonalizing the matrix
EigenVector eigenvalues;
EigenMatrix 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
|