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

Rewrite most of OBB code to fix design issues and bugs

This commit is contained in:
Mael Rouxel-Labbé 2019-12-09 12:25:38 +01:00
parent cd4123bf2c
commit 51e965e1dc
11 changed files with 655 additions and 1139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Optimal_bounding_box/benchmarks/Optimal_bounding_box/bench_obb.cpp
View File

@ -32,15 +32,6 @@ void gather_mesh_points(SurfaceMesh& mesh, std::vector<Point>& points)
points.push_back(get(pmap, v));
}
template <typename Point>
double calculate_volume(std::vector<Point> points)
{
CGAL::Bbox_3 bbox;
bbox = bbox_3(points.begin(), points.end());
K::Iso_cuboid_3 ic(bbox);
return ic.volume();
}
void bench_finding_obb(std::string fname)
{
std::ifstream input(fname);

228
Optimal_bounding_box/include/CGAL/Optimal_bounding_box/evolution.h
View File

@ -1,23 +1,16 @@
// Copyright (c) 2018 GeometryFactory (France).
// Copyright (c) 2018-2019 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Konstantinos Katrioplas
// Mael Rouxel-Labbé
//
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_EVOLUTION_H
#define CGAL_OPTIMAL_BOUNDING_BOX_EVOLUTION_H
@ -29,141 +22,124 @@
#include <CGAL/number_utils.h>
#include <algorithm>
#include <array>
#include <iostream>
#include <vector>
namespace CGAL {
namespace Optimal_bounding_box {
template <typename Linear_algebra_traits>
template <typename PointRange, typename Traits>
class Evolution
{
typedef typename Linear_algebra_traits::MatrixXd MatrixXd;
typedef typename Linear_algebra_traits::Matrix3d Matrix3d;
typedef typename Linear_algebra_traits::Vector3d Vector3d;
public:
Evolution(Population<Linear_algebra_traits>& pop, MatrixXd& points)
: population(pop), point_data(points)
{}
typedef typename Traits::FT FT;
typedef typename Traits::Matrix Matrix;
void genetic_algorithm()
typedef Optimal_bounding_box::Population<Traits> Population;
typedef typename Population::Simplex Simplex;
typedef Optimal_bounding_box::Fitness_map<Population, PointRange, Traits> Fitness_map;
Evolution(const PointRange& points,
CGAL::Random& rng,
const Traits& traits)
:
m_population(traits),
m_rng(rng), // @todo just a parameter of genetic_algorithm() ?
m_points(points),
m_traits(traits)
{ }
void genetic_algorithm(const std::size_t population_size = 50)
{
// random permutations
const std::size_t m = population.size();
m_population.initialize(population_size, m_rng);
//groups 1,2 : size m/2 groups 3,4 : size (m - m/2). m/2 is floored
const std::size_t size_first_group = m/2;
const std::size_t size_second_group = m - size_first_group;
const std::size_t m = population_size;
const std::size_t first_group_size = m / 2;
const std::size_t second_group_size = m - first_group_size;
std::vector<std::size_t> ids1(size_first_group), ids2(size_first_group);
std::vector<std::size_t> ids3(size_second_group), ids4(size_second_group);
std::vector<std::size_t> group1(first_group_size), group2(first_group_size);
std::vector<std::size_t> group3(second_group_size), group4(second_group_size);
CGAL::Random rng;
int im = static_cast<int>(m);
std::generate(ids1.begin(), ids1.end(), [&rng, &im] () { return rng.get_int(0, im); });
std::generate(ids2.begin(), ids2.end(), [&rng, &im] () { return rng.get_int(0, im); });
std::generate(ids3.begin(), ids3.end(), [&rng, &im] () { return rng.get_int(0, im); });
std::generate(ids4.begin(), ids4.end(), [&rng, &im] () { return rng.get_int(0, im); });
Population<Linear_algebra_traits> group1(size_first_group), group2(size_first_group);
Population<Linear_algebra_traits> group3(size_second_group), group4(size_second_group);
for(std::size_t i=0; i<size_first_group; ++i)
group1[i] = population[ids1[i]];
for(std::size_t i=0; i<size_first_group; ++i)
group2[i] = population[ids2[i]];
for(std::size_t i=0; i<size_second_group; ++i)
group3[i] = population[ids3[i]];
for(std::size_t i=0; i<size_second_group; ++i)
group4[i] = population[ids4[i]];
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
check_det(group1);
check_det(group2);
check_det(group3);
check_det(group4);
#endif
std::generate(group1.begin(), group1.end(), [&]{ return m_rng.get_int(0, im); });
std::generate(group2.begin(), group2.end(), [&]{ return m_rng.get_int(0, im); });
std::generate(group3.begin(), group3.end(), [&]{ return m_rng.get_int(0, im); });
std::generate(group4.begin(), group4.end(), [&]{ return m_rng.get_int(0, im); });
// crossover I
Population<Linear_algebra_traits> offspringsA(size_first_group);
double bias = 0.1;
for(std::size_t i=0; i<size_first_group; ++i)
std::vector<Simplex> new_simplices(m);
for(std::size_t i=0; i<first_group_size; ++i)
{
std::vector<Matrix3d> offspring(4);
std::array<Matrix, 4> offspring;
for(int j=0; j<4; ++j)
{
double r = rng.get_double();
double fitnessA = compute_fitness<Linear_algebra_traits>(group1[i][j], point_data);
double fitnessB = compute_fitness<Linear_algebra_traits>(group2[i][j], point_data);
double threshold;
const double r = m_rng.get_double();
const double fitnessA = compute_fitness<Traits>(m_population[group1[i]][j], m_points);
const double fitnessB = compute_fitness<Traits>(m_population[group2[i]][j], m_points);
const double threshold = (fitnessA < fitnessB) ? (0.5 + bias) : (0.5 - bias);
if(fitnessA < fitnessB)
threshold = 0.5 + bias;
if(r < threshold)
offspring[j] = m_population[group1[i]][j];
else
threshold = 0.5 - bias;
if(r < threshold) // choose A
offspring[j] = group1[i][j];
else // choose B
offspring[j] = group2[i][j];
offspring[j] = m_population[group2[i]][j];
}
offspringsA[i] = offspring;
}
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
std::cout << "offspringsA: \n" ;
check_det(offspringsA);
#endif
new_simplices[i] = std::move(offspring);
}
// crossover II
Population<Linear_algebra_traits> offspringsB(size_second_group);
bias = 0.1;
bias = 0.1; // @fixme should the bias change? What should be the initial value?
for(std::size_t i=0; i<size_second_group; ++i)
for(std::size_t i=0; i<second_group_size; ++i)
{
std::vector<Matrix3d> offspring(4);
std::array<Matrix, 4> offspring;
for(int j=0; j<4; ++j)
{
double fitnessA = compute_fitness<Linear_algebra_traits>(group3[i][j], point_data);
double fitnessB = compute_fitness<Linear_algebra_traits>(group4[i][j], point_data);
double lambda;
if(fitnessA < fitnessB)
lambda = 0.5 + bias;
else
lambda = 0.5 - bias;
const double fitnessA = compute_fitness<Traits>(m_population[group3[i]][j], m_points);
const double fitnessB = compute_fitness<Traits>(m_population[group4[i]][j], m_points);
const double lambda = (fitnessA < fitnessB) ? (0.5 + bias) : (0.5 - bias);
const double rambda = 1 - lambda; // the 'l' in 'lambda' stands for left
// combine information from A and B
offspring[j] = lambda * group3[i][j] + lambda * group4[i][j];
Matrix new_vertex(3, 3);
const Matrix& lm = m_population[group3[i]][j];
const Matrix& rm = m_population[group4[i]][j];
// just avoiding having to add matrix sums and scalar multiplications to the concept
new_vertex.set(0, 0, lambda * lm(0, 0) + rambda * rm(0, 0));
new_vertex.set(0, 1, lambda * lm(0, 1) + rambda * rm(0, 1));
new_vertex.set(0, 2, lambda * lm(0, 2) + rambda * rm(0, 2));
new_vertex.set(1, 0, lambda * lm(1, 0) + rambda * rm(1, 0));
new_vertex.set(1, 1, lambda * lm(1, 1) + rambda * rm(1, 1));
new_vertex.set(1, 2, lambda * lm(1, 2) + rambda * rm(1, 2));
new_vertex.set(2, 0, lambda * lm(2, 0) + rambda * rm(2, 0));
new_vertex.set(2, 1, lambda * lm(2, 1) + rambda * rm(2, 1));
new_vertex.set(2, 2, lambda * lm(2, 2) + rambda * rm(2, 2));
offspring[j] = m_traits.qr_factorization(new_vertex);
}
// qr factorization of the offspring
Linear_algebra_traits::qr_factorization(offspring);
offspringsB[i] = offspring;
new_simplices[first_group_size + i] = std::move(offspring);
}
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
std::cout << "offspringsB: \n" ;
check_det(offspringsB);
#endif
CGAL_assertion(offspringsA.size() == size_first_group);
CGAL_assertion(offspringsB.size() == size_second_group);
CGAL_assertion(offspringsA.size() + offspringsB.size() == population.size());
// next generatrion
for(std::size_t i=0; i<size_first_group; ++i)
population[i] = offspringsA[i];
for(std::size_t i=0; i<size_second_group; ++i)
population[size_first_group + i] = offspringsB[i];
m_population.simplices() = std::move(new_simplices);
}
void evolve(std::size_t generations)
void evolve(const std::size_t generations)
{
// hardcoded population size
const std::size_t population_size = 50;
// hardcoded nelder_mead_iterations
const std::size_t nelder_mead_iterations = 20;
@ -176,10 +152,10 @@ public:
for(std::size_t t=0; t<generations; ++t)
{
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
std::cout << "generation= " << t << "\n";
std::cout << "generation = " << t << "\n";
#endif
genetic_algorithm();
genetic_algorithm(population_size);
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
//std::cout << "pop after genetic" << std::endl;
@ -187,22 +163,22 @@ public:
//std::cout << std::endl;
#endif
for(std::size_t s = 0; s < population.size(); ++s)
nelder_mead<Linear_algebra_traits>(population[s], point_data, nelder_mead_iterations);
for(std::size_t s=0; s<population_size; ++s)
nelder_mead(m_population[s], m_points, nelder_mead_iterations, m_traits);
// stopping criteria
Fitness_map fitness_map(m_population, m_points);
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
//std::cout << "pop after nelder mead: " << std::endl;
//pop.show_population();
//std::cout << std::endl;
Fitness_map<Linear_algebra_traits, Matrix3d, MatrixXd> fitness_map_debug(population, point_data);
Matrix3d R_now = fitness_map_debug.get_best();
std::cout << "det = " << Linear_algebra_traits::determinant(R_now) << std::endl;
const Matrix& R_now = fitness_map.get_best();
std::cout << "det = " << m_traits.determinant(R_now) << std::endl;
#endif
// stopping criteria
Fitness_map<Linear_algebra_traits, Matrix3d, MatrixXd> fitness_map(population, point_data);
new_fit_value = fitness_map.get_best_fitness_value();
double difference = new_fit_value - prev_fit_value;
const double difference = new_fit_value - prev_fit_value;
if(CGAL::abs(difference) < tolerance * new_fit_value)
++stale;
@ -214,29 +190,19 @@ public:
}
}
const Matrix3d get_best()
const Matrix& get_best()
{
Fitness_map<Linear_algebra_traits, Matrix3d, MatrixXd> fitness_map(population, point_data);
Fitness_map fitness_map(m_population, m_points);
return fitness_map.get_best();
}
private:
Population<Linear_algebra_traits> population;
MatrixXd point_data;
Population m_population;
CGAL::Random m_rng;
const PointRange& m_points;
const Traits& m_traits;
};
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
template <typename Simplex>
void check_det(Population<Simplex>& pop)
{
for(std::size_t i=0, ps=pop.size(); i<ps; ++i)
{
for(std::size_t j=0; j<4; ++j)
std::cout << Linear_algebra_traits::determinant(pop[i][j]) << std::endl;
}
}
#endif
} // end namespace Optimal_bounding_box
} // end namespace CGAL

114
Optimal_bounding_box/include/CGAL/Optimal_bounding_box/fitness_function.h
View File

@ -1,96 +1,85 @@
// Copyright (c) 2018 GeometryFactory (France).
// Copyright (c) 2018-2019 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Konstantinos Katrioplas
// Mael Rouxel-Labbé
//
#ifndef CGAL_OPTIMAL_BOUNDING_FITNESS_FUNCTION_H
#define CGAL_OPTIMAL_BOUNDING_FITNESS_FUNCTION_H
#include <CGAL/Optimal_bounding_box/population.h>
#include <CGAL/assertions.h>
#include <algorithm>
#include <limits>
namespace CGAL {
namespace Optimal_bounding_box {
template <typename Linear_algebra_traits, typename Vertex, typename Matrix>
double compute_fitness(const Vertex& R, const Matrix& data)
template <typename Traits, typename PointRange>
typename Traits::FT
compute_fitness(const typename Traits::Matrix& R, // rotation matrix
const PointRange& points)
{
// R: rotation matrix
CGAL_assertion(R.cols() == 3);
CGAL_assertion(R.rows() == 3);
// data: points
CGAL_assertion(data.cols() == 3);
CGAL_assertion(data.rows() >= 3);
typedef typename Traits::FT FT;
typedef typename Traits::Point_3 Point;
typedef typename Linear_algebra_traits::Vector3d Vector3d;
typedef typename Linear_algebra_traits::Index Index;
CGAL_assertion(R.number_of_rows() == 3 && R.number_of_columns() == 3);
CGAL_assertion(points.size() >= 3);
double xmin, xmax, ymin, ymax, zmin, zmax;
for(Index i=0; i < static_cast<Index>(data.rows()); ++i)
FT xmin, ymin, zmin, xmax, ymax, zmax;
xmin = ymin = zmin = std::numeric_limits<double>::max();
xmax = ymax = zmax = std::numeric_limits<double>::lowest();
for(const Point& pt : points)
{
Vector3d vec = Linear_algebra_traits::row3(data, i);
vec = R * vec;
const FT x = pt.x(), y = pt.y(), z = pt.z();
if(i == 0)
{
xmin = xmax = vec.coeff(0);
ymin = ymax = vec.coeff(1);
zmin = zmax = vec.coeff(2);
}
else
{
if(vec.coeff(0) < xmin) xmin = vec.coeff(0);
if(vec.coeff(1) < ymin) ymin = vec.coeff(1);
if(vec.coeff(2) < zmin) zmin = vec.coeff(2);
if(vec.coeff(0) > xmax) xmax = vec.coeff(0);
if(vec.coeff(1) > ymax) ymax = vec.coeff(1);
if(vec.coeff(2) > zmax) zmax = vec.coeff(2);
}
const FT rx = x*R(0, 0) + y*R(0, 1) + z*R(0, 2);
const FT ry = x*R(1, 0) + y*R(1, 1) + z*R(1, 2);
const FT rz = x*R(2, 0) + y*R(2, 1) + z*R(2, 2);
xmin = (std::min)(xmin, rx);
ymin = (std::min)(ymin, ry);
zmin = (std::min)(zmin, rz);
xmax = (std::max)(xmax, rx);
ymax = (std::max)(ymax, ry);
zmax = (std::max)(zmax, rz);
}
CGAL_assertion(xmax > xmin);
CGAL_assertion(ymax > ymin);
CGAL_assertion(zmax > zmin);
// volume
return ((xmax - xmin) * (ymax - ymin) * (zmax - zmin));
}
template <typename Linear_algebra_traits, typename Vertex, typename Matrix>
template <typename Population, typename PointRange, typename Traits>
struct Fitness_map
{
Fitness_map(Population<Linear_algebra_traits>& p, Matrix& points)
: pop(p), points(points)
{}
typedef typename Traits::FT FT;
typedef typename Population::Vertex Vertex;
const Vertex get_best()
Fitness_map(const Population& population,
const PointRange& points)
:
m_pop(population),
m_points(points)
{ }
const Vertex& get_best() const // @todo any point caching this?
{
std::size_t simplex_id, vertex_id;
double best_fitness = std::numeric_limits<int>::max();
for(std::size_t i=0; i<pop.size(); ++i)
FT best_fitness = std::numeric_limits<double>::max();
for(std::size_t i=0, ps=m_pop.size(); i<ps; ++i)
{
for(std::size_t j=0; j<4; ++j)
{
const Vertex vertex = pop[i][j];
const double fitness = compute_fitness<Linear_algebra_traits>(vertex, points);
const Vertex& vertex = m_pop[i][j];
const FT fitness = compute_fitness<Traits>(vertex, m_points);
if(fitness < best_fitness)
{
simplex_id = i;
@ -100,17 +89,18 @@ struct Fitness_map
}
}
return pop[simplex_id][vertex_id];
return m_pop[simplex_id][vertex_id];
}
double get_best_fitness_value()
FT get_best_fitness_value() const
{
const Vertex best_mat = get_best();
return compute_fitness<Linear_algebra_traits>(best_mat, points);
const Vertex& best_mat = get_best();
return compute_fitness<Traits>(best_mat, m_points);
}
Population<Linear_algebra_traits>& pop;
const Matrix& points;
private:
const Population& m_pop;
const PointRange& m_points;
};
} // end namespace Optimal_bounding_box

118
Optimal_bounding_box/include/CGAL/Optimal_bounding_box/helper.h
View File

@ -1,24 +1,16 @@
// Copyright (c) 2018 GeometryFactory (France).
// Copyright (c) 2018-2019 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Konstantinos Katrioplas
// Mael Rouxel-Labbé
//
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_HELPER_H
#define CGAL_OPTIMAL_BOUNDING_BOX_HELPER_H
@ -30,6 +22,11 @@
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_matrix.h>
#include <Eigen/QR>
#endif
#include <fstream>
#include <string>
#include <vector>
@ -37,74 +34,45 @@
namespace CGAL {
namespace Optimal_bounding_box {
template <typename Matrix, typename Point>
void fill_matrix(const std::vector<Point>& v_points, Matrix& points_mat)
#ifdef CGAL_EIGEN3_ENABLED
// @tmp move the header include too
template <typename K>
class Optimal_bounding_box_traits
: public K
{
points_mat.resize(v_points.size(), 3);
for(std::size_t i = 0; i < v_points.size(); ++i)
public:
typedef typename K::FT FT;
typedef CGAL::Eigen_matrix<FT, 3, 3> Matrix;
private:
typedef typename Matrix::EigenType EigenType;
public:
explicit Optimal_bounding_box_traits(const K& k = K()) : K(k) { }
/// Get the transpose of a `Matrix<NT>` matrix
Matrix transpose(const Matrix& mat) const
{
Point p = v_points[i];
points_mat.set_coef(i, 0, CGAL::to_double(p.x()));
points_mat.set_coef(i, 1, CGAL::to_double(p.y()));
points_mat.set_coef(i, 2, CGAL::to_double(p.z()));
}
}
template <typename SurfaceMesh, typename Matrix>
void sm_to_matrix(SurfaceMesh& sm, Matrix& mat)
{
typedef typename boost::property_map<SurfaceMesh, boost::vertex_point_t>::const_type Vpm;
typedef typename boost::property_traits<Vpm>::reference Point_ref;
typedef typename boost::graph_traits<SurfaceMesh>::vertex_descriptor vertex_descriptor;
Vpm vpm = get(boost::vertex_point, sm);
mat.resize(vertices(sm).size(), 3);
std::size_t i = 0;
for(vertex_descriptor v : vertices(sm))
{
Point_ref p = get(vpm, v);
mat.set_coef(i, 0, CGAL::to_double(p.x()));
mat.set_coef(i, 1, CGAL::to_double(p.y()));
mat.set_coef(i, 2, CGAL::to_double(p.z()));
++i;
}
}
template <typename Point>
double calculate_volume(const std::vector<Point>& points)
{
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
CGAL::Bbox_3 bbox = bbox_3(points.begin(), points.end());
K::Iso_cuboid_3 ic(bbox);
return ic.volume();
}
// it is called after post processing in debug only
template <typename Matrix>
void matrix_to_mesh_and_draw(Matrix& data_points, std::string filename)
{
typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Mesh;
// Simplex -> std::vector
std::vector<Point> points;
for(int i = 0; i < data_points.rows(); ++i)
{
Point p(data_points(i, 0), data_points(i, 1), data_points(i, 2));
points.push_back(p);
return Matrix(mat.eigen_object().transpose());
}
Mesh mesh;
CGAL::make_hexahedron(points[0], points[1], points[2], points[3],
points[4], points[5], points[6], points[7], mesh);
/// Get the determinant of a `Matrix<NT>` matrix
FT determinant(const Matrix& matrix) const
{
return matrix.eigen_object().determinant();
}
std::ofstream out(filename);
out << mesh;
out.close();
}
/// Performs QR decomposition of matrix A to a unitary matrix and an upper triagonal
/// and returns the unitary matrix.
Matrix qr_factorization(const Matrix& A) const
{
Eigen::HouseholderQR<EigenType> qr(A.eigen_object());
CGAL_assertion(CGAL::abs(determinant(Matrix(EigenType(qr.householderQ()))) - 1.) < 0.000001);
return Matrix(EigenType(qr.householderQ()));
}
};
#endif
} // end namespace Optimal_bounding_box
} // end namespace CGAL

186
Optimal_bounding_box/include/CGAL/Optimal_bounding_box/nelder_mead_functions.h
View File

@ -1,188 +1,160 @@
// Copyright (c) 2018 GeometryFactory (France).
// Copyright (c) 2018-2019 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Konstantinos Katrioplas
// Mael Rouxel-Labbé
//
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_NEALDER_MEAD_FUNCTIONS_H
#define CGAL_OPTIMAL_BOUNDING_BOX_NEALDER_MEAD_FUNCTIONS_H
#include <CGAL/assertions.h>
#include <CGAL/Optimal_bounding_box/fitness_function.h>
#include <boost/iterator/counting_iterator.hpp>
#include <CGAL/assertions.h>
#include <vector>
#include <algorithm>
#include <array>
namespace CGAL {
namespace Optimal_bounding_box {
template <typename Linear_algebra_traits, typename Matrix>
const Matrix reflection(const Matrix& S_centroid,
const Matrix& S_worst)
template <typename Matrix, typename Traits>
Matrix reflection(const Matrix& S_centroid,
const Matrix& S_worst,
const Traits& traits)
{
CGAL_assertion(S_centroid.rows() == 3);
CGAL_assertion(S_centroid.rows() == 3);
CGAL_assertion(S_worst.cols() == 3);
CGAL_assertion(S_worst.cols() == 3);
CGAL_assertion(S_centroid.number_of_rows() == 3 && S_centroid.number_of_columns() == 3);
CGAL_assertion(S_worst.number_of_rows() == 3 && S_worst.number_of_columns() == 3);
return S_centroid * Linear_algebra_traits::transpose(S_worst) * S_centroid;
return S_centroid * traits.transpose(S_worst) * S_centroid;
}
template <typename Linear_algebra_traits, typename Matrix>
const Matrix expansion(const Matrix& S_centroid,
const Matrix& S_worst,
const Matrix& S_reflection)
template <typename Matrix, typename Traits>
Matrix expansion(const Matrix& S_centroid,
const Matrix& S_worst,
const Matrix& S_reflection,
const Traits& traits)
{
CGAL_assertion(S_centroid.rows() == 3);
CGAL_assertion(S_centroid.rows() == 3);
CGAL_assertion(S_worst.cols() == 3);
CGAL_assertion(S_worst.cols() == 3);
CGAL_assertion(S_reflection.cols() == 3);
CGAL_assertion(S_reflection.cols() == 3);
CGAL_assertion(S_centroid.number_of_rows() == 3 && S_centroid.number_of_columns() == 3);
CGAL_assertion(S_worst.number_of_rows() == 3 && S_worst.number_of_columns() == 3);
CGAL_assertion(S_reflection.number_of_rows() == 3 && S_reflection.number_of_columns() == 3);
return S_centroid * Linear_algebra_traits::transpose(S_worst) * S_reflection;
return S_centroid * traits.transpose(S_worst) * S_reflection;
}
template <typename Linear_algebra_traits, typename Matrix>
template <typename Matrix, typename Traits>
Matrix mean(const Matrix& m1,
const Matrix& m2)
const Matrix& m2,
const Traits& traits)
{
// same API for reduction
CGAL_assertion(m1.rows() == 3);
CGAL_assertion(m1.rows() == 3);
CGAL_assertion(m2.cols() == 3);
CGAL_assertion(m2.cols() == 3);
CGAL_assertion(m1.number_of_rows() == 3 && m1.number_of_columns() == 3);
CGAL_assertion(m2.number_of_rows() == 3 && m2.number_of_columns() == 3);
Matrix reduction = 0.5 * m1 + 0.5 * m2;
Matrix Q = Linear_algebra_traits::qr_factorization(reduction);
double det = Linear_algebra_traits::determinant(Q);
return Q / det;
const Matrix reduction = 0.5 * m1 + 0.5 * m2;
const Matrix Q = traits.qr_factorization(reduction);
const typename Traits::FT det = traits.determinant(Q);
return (1. / det) * Q;
}
template <typename Linear_algebra_traits, typename Matrix>
template <typename Matrix, typename Traits>
const Matrix nm_centroid(const Matrix& S1,
const Matrix& S2,
const Matrix& S3)
const Matrix& S3,
const Traits& traits)
{
Matrix mean = (S1 + S2 + S3) / 3.0;
Matrix Q = Linear_algebra_traits::qr_factorization(mean);
double det = Linear_algebra_traits::determinant(Q);
return Q / det;
const Matrix mean = (1./3.) * (S1 + S2 + S3);
const Matrix Q = traits.qr_factorization(mean);
const typename Traits::FT det = traits.determinant(Q);
return (1. / det) * Q;
}
// needed in nelder mead algorithm
struct Comparator
// It's a 3D simplex with 4 rotation matrices as vertices
template <typename Simplex, typename PointRange, typename Traits>
void nelder_mead(Simplex& simplex,
const PointRange& points,
const std::size_t nelder_mead_iterations,
const Traits& traits)
{
Comparator(const std::vector<double>& in) : fitness(in) { }
typedef typename Traits::FT FT;
typedef typename Traits::Matrix Matrix;
inline bool operator() (std::size_t& i, std::size_t& j) {
return fitness[i] < fitness[j];
}
std::array<FT, 4> fitness;
std::array<std::size_t, 4> indices = {{ 0, 1, 2, 3 }};
const std::vector<double>& fitness;
};
// simplex: 4 rotation matrices are its vertices
template <typename Linear_algebra_traits>
void nelder_mead(std::vector<typename Linear_algebra_traits::Matrix3d>& simplex,
const typename Linear_algebra_traits::MatrixXd& point_data,
std::size_t nelder_mead_iterations)
{
CGAL_assertion(simplex.size() == 4); // tetrahedron
typedef typename Linear_algebra_traits::Matrix3d Matrix3d;
std::vector<double> fitness(4);
std::vector<std::size_t> indices(boost::counting_iterator<std::size_t>(0),
boost::counting_iterator<std::size_t>(simplex.size()));
for(std::size_t t = 0; t < nelder_mead_iterations; ++t)
for(std::size_t t=0; t<nelder_mead_iterations; ++t)
{
for(std::size_t i=0; i<4; ++i)
{
fitness[i] = compute_fitness<Linear_algebra_traits>(simplex[i], point_data);
}
CGAL_assertion(fitness.size() == 4);
CGAL_assertion(indices.size() == 4);
fitness[i] = compute_fitness<Traits>(simplex[i], points);
// get indices of sorted sequence
Comparator compare_indices(fitness); // @todo lambda this stuff
std::sort(indices.begin(), indices.end(), compare_indices);
std::sort(indices.begin(), indices.end(),
[&fitness](const std::size_t i, const std::size_t j) -> bool
{ return fitness[i] < fitness[j]; });
// new sorted simplex & fitness
std::vector<Matrix3d> s_simplex(4);
std::vector<double> s_fitness(4);
for(int i = 0; i < 4; ++i)
Simplex s_simplex;
std::array<FT, 4> s_fitness;
for(int i=0; i<4; ++i)
{
s_simplex[i] = simplex[indices[i]];
s_fitness[i] = fitness[indices[i]];
}
simplex = s_simplex;
fitness = s_fitness;
simplex = std::move(s_simplex);
fitness = std::move(s_fitness);
// centroid
const Matrix3d v_centroid = nm_centroid<Linear_algebra_traits>(simplex[0], simplex[1], simplex[2]);
const Matrix v_centroid = nm_centroid(simplex[0], simplex[1], simplex[2], traits);
// find worst's vertex reflection
const Matrix3d v_worst = simplex[3];
const Matrix3d v_refl = reflection<Linear_algebra_traits>(v_centroid, v_worst);
const double f_refl = compute_fitness<Linear_algebra_traits>(v_refl, point_data);
const Matrix& v_worst = simplex[3];
const Matrix v_refl = reflection(v_centroid, v_worst, traits);
const FT f_refl = compute_fitness<Traits>(v_refl, points);
if(f_refl < fitness[2])
{
if(f_refl >= fitness[0]) // if reflected point is not better than the best
{
// do reflection
simplex[3] = v_refl;
// reflection
simplex[3] = std::move(v_refl);
}
else
{
// expansion
const Matrix3d v_expand = expansion<Linear_algebra_traits>(v_centroid, v_worst, v_refl);
const double f_expand = compute_fitness<Linear_algebra_traits>(v_expand, point_data);
const Matrix v_expand = expansion(v_centroid, v_worst, v_refl, traits);
const FT f_expand = compute_fitness<Traits>(v_expand, points);
if(f_expand < f_refl)
simplex[3] = v_expand;
simplex[3] = std::move(v_expand);
else
simplex[3] = v_refl;
simplex[3] = std::move(v_refl);
}
}
else // if reflected vertex is not better
{
const Matrix3d v_mean = mean<Linear_algebra_traits>(v_centroid, v_worst);
const double f_mean = compute_fitness<Linear_algebra_traits>(v_mean, point_data);
const Matrix v_mean = mean(v_centroid, v_worst, traits);
const FT f_mean = compute_fitness<Traits>(v_mean, points);
if(f_mean <= fitness[3])
{
// contraction of worst
simplex[3] = v_mean;
simplex[3] = std::move(v_mean);
}
else
{
// reduction: move all vertices towards the best
for(std::size_t i=1; i<4; ++i)
{
simplex[i] = mean<Linear_algebra_traits>(simplex[i], simplex[0]);
}
simplex[i] = mean(simplex[i], simplex[0], traits);
}
}
CGAL_assertion(simplex.size() == 4); // tetrahedron
} // iterations
} // nelder mead iterations
}
} // end namespace Optimal_bounding_box

312
Optimal_bounding_box/include/CGAL/Optimal_bounding_box/optimal_bounding_box.h
View File

@ -1,24 +1,16 @@
// Copyright (c) 2018 GeometryFactory (France).
// Copyright (c) 2018-2019 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Konstantinos Katrioplas
// Mael Rouxel-Labbé
//
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_OBB_H
#define CGAL_OPTIMAL_BOUNDING_BOX_OBB_H
@ -26,122 +18,88 @@
#include <CGAL/Optimal_bounding_box/evolution.h>
#include <CGAL/Optimal_bounding_box/helper.h>
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h> // @tmp
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <CGAL/assertions.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/boost/graph/helpers.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Iso_cuboid_3.h>
#include <CGAL/Random.h>
#include <CGAL/Simple_cartesian.h>
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
#include <CGAL/Timer.h>
#endif
#if defined(CGAL_EIGEN3_ENABLED)
#include <CGAL/Eigen_linear_algebra_traits.h>
#include <CGAL/Real_timer.h>
#endif
#include <array>
#include <iostream>
#include <iterator>
#include <type_traits>
#include <vector>
namespace CGAL {
namespace Optimal_bounding_box {
namespace internal {
// works on matrices only
/// \cond SKIP_IN_MANUAL
template <typename Linear_algebra_traits, typename Vertex, typename Matrix>
void post_processing(const Matrix& points,
Vertex& R,
Matrix& obb)
template <typename PointRange, typename Traits>
void post_processing(const typename Traits::Matrix& R,
std::array<typename Traits::Point_3, 8>& obb,
const PointRange& points,
const Traits& traits)
{
CGAL_assertion(points.cols() == 3);
CGAL_assertion(R.rows() == 3);
CGAL_assertion(R.cols() == 3);
CGAL_assertion(obb.rows() == 8);
CGAL_assertion(obb.cols() == 3);
typedef typename Traits::FT FT;
typedef typename Traits::Point_3 Point;
typedef typename Traits::Matrix Matrix;
// 1) rotate points with R
Matrix rotated_points(points.rows(), points.cols());
rotated_points = points * Linear_algebra_traits::transpose(R);
CGAL_assertion(R.number_of_rows() == 3 && R.number_of_columns() == 3);
// 2) get AABB from rotated points
typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_3 Point;
typedef typename Linear_algebra_traits::Index index;
const Matrix Rt = traits.transpose(R);
// Simplex -> std::vector
std::vector<Point> v_points;
for(index i=0; i<static_cast<index>(rotated_points.rows()); ++i)
{
Point p(rotated_points(i, 0), rotated_points(i, 1), rotated_points(i, 2));
v_points.push_back(p);
}
CGAL::Bbox_3 bbox;
bbox = bbox_3(v_points.begin(), v_points.end());
K::Iso_cuboid_3 ic(bbox);
Matrix aabb(8, 3);
for(std::size_t i = 0; i<8; ++i)
for(const Point& pt : points)
{
aabb.set_coef(i, 0, ic[i].x());
aabb.set_coef(i, 1, ic[i].y());
aabb.set_coef(i, 2, ic[i].z());
// @fixme should it be R here Rt at the other one... ?
const FT x = pt.x(), y = pt.y(), z = pt.z();
Point rotated_pt(x*R(0, 0) + y*R(0, 1) + z*R(0, 2),
x*R(1, 0) + y*R(1, 1) + z*R(1, 2),
x*R(2, 0) + y*R(2, 1) + z*R(2, 2));
bbox += traits.construct_bbox_3_object()(rotated_pt);
}
// @todo could avoid building a cuboid
typename Traits::Iso_cuboid_3 ic(bbox);
// 3) apply inverse rotation to rotated AABB
obb = aabb * R;
for(std::size_t i = 0; i<8; ++i)
{
const FT x = ic[i].x(), y = ic[i].y(), z = ic[i].z();
obb[i] = Point(x*Rt(0, 0) + y*Rt(0, 1) + z*Rt(0, 2),
x*Rt(1, 0) + y*Rt(1, 1) + z*Rt(1, 2),
x*Rt(2, 0) + y*Rt(2, 1) + z*Rt(2, 2));
}
}
/// \endcond
/// \ingroup OBB_grp
/// calculates the optimal bounding box.
///
/// @tparam Point the point type
/// @tparam LinearAlgebraTraits a model of `LinearAlgebraTraits`. If no instance of `LinearAlgebraTraits`
/// is provided, then `CGAL::Eigen_linear_algebra_traits` is used.
///
/// @param points the input points that are included in the optimal bounding box.
/// @param obb_points the eight points of the optimal bounding box to be calculated.
/// @param use_ch a bool flag to indicating whether to use the convex hull of the input points
/// as an optimization step.
template <typename Point, typename LinearAlgebraTraits>
void compute_optimal_bounding_box(const std::vector<Point>& points,
std::vector<Point>& obb_points,
LinearAlgebraTraits&,
bool use_ch)
template <typename PointRange, typename Traits>
void construct_optimal_bounding_box(std::array<typename Traits::Point_3, 8>& obb_points,
const PointRange& points,
CGAL::Random& rng,
const Traits& traits)
{
CGAL_assertion(points.size() >= 3);
typedef typename Traits::Matrix Matrix;
if(obb_points.size() != 8)
obb_points.resize(8);
// eigen linear algebra traits
typedef typename LinearAlgebraTraits::MatrixXd MatrixXd;
typedef typename LinearAlgebraTraits::Matrix3d Matrix3d;
MatrixXd points_mat;
if(use_ch) // get the ch3
{
std::vector<Point> ch_points;
CGAL::extreme_points_3(points, std::back_inserter(ch_points));
CGAL::Optimal_bounding_box::fill_matrix(ch_points, points_mat);
}
else
{
CGAL::Optimal_bounding_box::fill_matrix(points, points_mat);
}
std::size_t max_generations = 100;
Population<LinearAlgebraTraits> pop(50);
const std::size_t max_generations = 100;
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
CGAL::Timer timer;
CGAL::Real_timer timer;
timer.start();
#endif
CGAL::Optimal_bounding_box::Evolution<LinearAlgebraTraits> search_solution(pop, points_mat);
Evolution<PointRange, Traits> search_solution(points, rng, traits);
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
std::cout << "constructor: " << timer.time() << std::endl;
@ -155,105 +113,163 @@ void compute_optimal_bounding_box(const std::vector<Point>& points,
timer.reset();
#endif
Matrix3d rotation = search_solution.get_best();
const Matrix& rotation = search_solution.get_best();
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
std::cout << "get best: " << timer.time() << std::endl;
#endif
MatrixXd obb; // may be preallocated at compile time
obb.resize(8, 3);
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
timer.reset();
#endif
post_processing<LinearAlgebraTraits>(points_mat, rotation, obb);
post_processing(rotation, obb_points, points, traits);
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
std::cout << "post porcessing: " << timer.time() << std::endl;
#endif
// matrix -> vector
for(std::size_t i=0; i<8; ++i)
obb_points.emplace(obb_points.begin() + i, obb(i, 0), obb(i, 1), obb(i, 2));
}
template <typename Point>
void compute_optimal_bounding_box(const std::vector<Point>& points,
std::vector<Point>& obb_points,
bool use_ch)
template <typename PointRange, typename Traits>
void construct_optimal_bounding_box(std::array<typename Traits::Point_3, 8>& obb_points,
const bool use_ch,
const PointRange& points,
CGAL::Random& rng,
const Traits& traits)
{
#if defined(CGAL_EIGEN3_ENABLED)
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
#else
#pragma message("Error: You must either provide linear traits or link CGAL with the Eigen library")
Linear_algebra_traits; // no parameter provided, and Eigen is not enabled --> don't compile!
#endif
typedef typename Traits::Matrix Matrix;
typedef typename Traits::Point_3 Point;
Linear_algebra_traits la_traits;
compute_optimal_bounding_box(points, obb_points, la_traits, use_ch);
CGAL_static_assertion((std::is_same<typename boost::range_value<PointRange>::type, Point>::value));
if(use_ch) // construct the convex hull to reduce the number of points
{
std::vector<Point> ch_points;
extreme_points_3(points, std::back_inserter(ch_points));
return construct_optimal_bounding_box(obb_points, ch_points, rng, traits);
}
else
{
return construct_optimal_bounding_box(obb_points, points, rng, traits);
}
}
} // namespace Optimal_bounding_box
} // namespace internal
/// \ingroup OBB_grp
/// calculates the optimal bounding box.
///
/// @tparam PolygonMesh a model of `FaceListGraph`
/// @tparam LinearAlgebraTraits a model of `LinearAlgebraTraits`. If no instance of `LinearAlgebraTraits`
/// is provided, then `CGAL::Eigen_linear_algebra_traits` is used.
/// \tparam PointRange a model of `Range` with value type `Point`
/// \tparam NamedParameters
///
/// @param pmesh the input mesh.
/// @param obbmesh the hexaedron mesh to be built out of the optimal bounding box.
/// @param la_traits an instance of the linear algebra traits.
/// @param use_ch a bool flag to indicating whether to use the convex hull of the input points
/// as an optimization step.
template <typename PolygonMesh, typename LinearAlgebraTraits>
void compute_optimal_bounding_box(const PolygonMesh& pmesh,
PolygonMesh& obbmesh,
LinearAlgebraTraits& la_traits,
bool use_ch)
/// \param points the input points that are included in the optimal bounding box.
/// \param obb_points the eight points of the optimal bounding box to be calculated.
///
template <typename PointRange,
typename Point,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
void optimal_bounding_box(const PointRange& points,
std::array<Point, 8>& obb_points,
const CGAL_PMP_NP_CLASS& np)
{
using CGAL::parameters::choose_parameter;
using CGAL::parameters::get_parameter;
#if 0 // @fixme
typedef typename GetSvdTraits<CGAL_PMP_NP_CLASS>::type SVD_traits;
CGAL_static_assertion_msg(
!(std::is_same<SVD_traits, typename GetSvdTraits<CGAL_PMP_NP_CLASS>::NoTraits>::value),
"Error: must provide traits or enable Eigen");
const SVD_traits traits = choose_parameter(get_parameter(np, svd_traits), SVD_traits());
#else
typedef typename CGAL::Kernel_traits<Point>::type K;
typedef Optimal_bounding_box::Optimal_bounding_box_traits<K> Traits;
Traits traits;
#endif
// @fixme
const bool use_ch = true; //choose_parameter(get_parameter(np, internal_np::use_convex_hull), false);
unsigned int seed = choose_parameter(get_parameter(np, internal_np::random_seed), 0);
CGAL::Random rng(seed);
return Optimal_bounding_box::internal::construct_optimal_bounding_box(obb_points, use_ch, points, rng, traits);
}
template <typename PolygonMesh, typename Point>
void optimal_bounding_box(const PolygonMesh& pmesh,
std::array<Point, 8>& obb_mesh)
{
return optimal_bounding_box(pmesh, obb_mesh, CGAL::parameters::all_default());
}
///////////////////////////////////////////////////////////////////////////////////////////////////
/// With polygon meshes
/////////////////////////////////////////////////////////////////////////////////////////////////
/// \ingroup OBB_grp
/// calculates the optimal bounding box.
///
/// \tparam PolygonMesh a model of `FaceListGraph`
/// \tparam NamedParameters
///
/// \param pmesh the input mesh.
/// \param obb_mesh the hexaedron mesh to be built out of the optimal bounding box.
///
template <typename PolygonMesh,
typename CGAL_PMP_NP_TEMPLATE_PARAMETERS>
void optimal_bounding_box(const PolygonMesh& pmesh,
PolygonMesh& obb_mesh,
const CGAL_PMP_NP_CLASS& np)
{
namespace PMP = CGAL::Polygon_mesh_processing;
using CGAL::parameters::choose_parameter;
using CGAL::parameters::get_parameter;
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename PMP::GetVertexPointMap<PolygonMesh, CGAL_PMP_NP_CLASS>::const_type VPM;
typedef typename boost::property_traits<VPM>::value_type Point;
CGAL_assertion(vertices(pmesh).size() >= 3);
if(vertices(pmesh).size() <= 3)
{
std::cerr << "Not enough points in the mesh!\n";
// @fixme is that even true
std::cerr << "The optimal bounding box cannot be computed for a mesh with fewer than 4 vertices!\n";
return;
}
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type Vpm;
typedef typename boost::property_traits<Vpm>::value_type Point;
VPM vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
get_const_property_map(vertex_point, pmesh));
std::vector<Point> points;
Vpm pmap = get(boost::vertex_point, pmesh);
for(vertex_descriptor v : vertices(pmesh))
points.push_back(get(pmap, v));
points.reserve(num_vertices(pmesh));
std::vector<Point> obb_points;
compute_optimal_bounding_box(points, obb_points, la_traits, use_ch);
for(vertex_descriptor v : vertices(pmesh))
points.push_back(get(vpm, v));
std::array<Point, 8> obb_points;
optimal_bounding_box(points, obb_points, np);
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
obb_points[4], obb_points[5], obb_points[6], obb_points[7], obbmesh);
obb_points[4], obb_points[5], obb_points[6], obb_points[7], obb_mesh);
}
///////////////////////////////////////////////////////////////////////////////////////////////////
/// Convenience overloads for polygon meshes
/////////////////////////////////////////////////////////////////////////////////////////////////
template <typename PolygonMesh>
void compute_optimal_bounding_box(const PolygonMesh& pmesh,
PolygonMesh& obbmesh,
bool use_ch)
void optimal_bounding_box(const PolygonMesh& pmesh,
PolygonMesh& obb_mesh)
{
#if defined(CGAL_EIGEN3_ENABLED)
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
#else
#pragma message("Error: You must either provide linear traits or link CGAL with the Eigen library")
Linear_algebra_traits; // no parameter provided, and Eigen is not enabled --> don't compile!
#endif
Linear_algebra_traits la_traits;
compute_optimal_bounding_box(pmesh, obbmesh, la_traits, use_ch);
return optimal_bounding_box(pmesh, obb_mesh, CGAL::parameters::all_default());
}
} // end namespace Optimal_bounding_box
} // end namespace CGAL
#endif // CGAL_OPTIMAL_BOUNDING_BOX_OBB_H

144
Optimal_bounding_box/include/CGAL/Optimal_bounding_box/population.h
View File

@ -1,113 +1,87 @@
// Copyright (c) 2018 GeometryFactory (France).
// Copyright (c) 2018-2019 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Konstantinos Katrioplas
// Mael Rouxel-Labbé
//
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_POPULATION_H
#define CGAL_OPTIMAL_BOUNDING_BOX_POPULATION_H
#include <CGAL/assertions.h>
#include <CGAL/Random.h>
#include <utility>
#include <vector>
namespace CGAL {
namespace Optimal_bounding_box {
template<typename Linear_algebra_traits>
template<typename Traits>
class Population
{
typedef typename Linear_algebra_traits::Matrix3d Matrix;
typedef std::vector<Matrix> Simplex;
public:
Population(std::size_t size)
:
n(size),
random_generator(CGAL::Random())
{
// reserve pop space
pop.reserve(n);
typedef typename Traits::FT FT;
typedef typename Traits::Matrix Matrix;
typedef Matrix Vertex;
// create simplices
for(std::size_t i=0 ; i<n; ++i)
{
Simplex simplex(4);
create_simplex(simplex);
CGAL_assertion(simplex.size() == 4);
pop.push_back(simplex);
}
typedef std::array<Matrix, 4> Simplex;
typedef std::vector<Simplex> Simplex_container;
private:
Matrix create_vertex(CGAL::Random& rng) const
{
Matrix R;
for(std::size_t i=0; i<3; ++i)
for(std::size_t j=0; j<3; ++j)
R.set(i, j, FT(rng.get_double()));
return R;
}
// create random population
Simplex create_simplex(CGAL::Random& rng) const
{
Simplex simplex;
for(std::size_t i=0; i<4; ++i)
simplex[i] = m_traits.qr_factorization(create_vertex(rng));
return simplex;
}
public:
Population(const Traits& traits) : m_traits(traits) { }
void initialize(std::size_t population_size,
CGAL::Random& rng)
{
m_pop.clear();
m_pop.reserve(population_size);
for(std::size_t i=0; i<population_size; ++i)
m_pop.emplace_back(create_simplex(rng));
}
// Access
std::size_t size() const { return m_pop.size(); }
Simplex& operator[](const std::size_t i) { CGAL_assertion(i < m_pop.size()); return m_pop[i]; }
const Simplex& operator[](const std::size_t i) const { CGAL_assertion(i < m_pop.size()); return m_pop[i]; }
Simplex_container& simplices() { return m_pop; }
// Debug
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
void show_population();
#endif
std::size_t size() const { return n; }
// access simplex
Simplex& operator[](const std::size_t i)
{
CGAL_assertion(i < n);
return pop[i];
}
const Simplex& operator[](const std::size_t i) const
{
CGAL_assertion(i < n);
return pop[i];
}
private:
// create random population
void create_simplex(Simplex& simplex)
{
CGAL_assertion(simplex.size() == 4);
for(std::size_t i=0; i<4; ++i)
{
Matrix R;
if(R.cols() == 0 || R.rows() == 0)
R.resize(3, 3);
create_vertex(R);
Matrix Q = Linear_algebra_traits::qr_factorization(R);
simplex[i] = Q;
}
CGAL_assertion(simplex.size() == 4);
}
void create_vertex(Matrix& R)
{
CGAL_assertion(R.rows() == 3);
CGAL_assertion(R.cols() == 3);
for(std::size_t i=0; i<3; ++i)
{
for(std::size_t j=0; j<3; ++j)
R.set_coef(i, j, random_generator.get_double());
}
}
std::size_t n;
CGAL::Random random_generator;
std::vector<Simplex> pop;
std::vector<Simplex> m_pop;
const Traits& m_traits;
};
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
@ -115,15 +89,11 @@ template <typename Matrix>
void Population<Matrix>::show_population()
{
std::size_t id = 0;
for(const Simplex& i : pop)
for(const Simplex& simplex : m_pop)
{
CGAL_assertion(i.size() == 4);
std:: cout << "Simplex: "<< id++ << std::endl;
for(const Matrix& R : i)
{
std::cout << R; // eigen out
std::cout << "\n\n";
}
std::cout << "Simplex: " << id++ << std::endl;
for(const Matrix& R : simplex)
std::cout << R << "\n\n";
std:: cout << std:: endl;
}
}

260
Optimal_bounding_box/test/Optimal_bounding_box/test_optimization_algorithms.cpp
View File

@ -21,13 +21,12 @@ bool assert_doubles(double d1, double d2, double epsilon)
void test_nelder_mead()
{
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
typedef Linear_algebra_traits::Matrix3d Matrix3d;
typedef Linear_algebra_traits::MatrixXd MatrixXd;
typedef Linear_algebra_traits::Matrix Matrix;
MatrixXd data_points(4,3);
Matrix data_points(4,3);
data_points(0,0) = 0.866802;
data_points(0,1) = 0.740808,
data_points(0,2) = 0.895304,
data_points(0,1) = 0.740808;
data_points(0,2) = 0.895304;
data_points(1,0) = 0.912651;
data_points(1,1) = 0.761565;
data_points(1,2) = 0.160330;
@ -35,19 +34,19 @@ void test_nelder_mead()
data_points(2,1) = 0.892578;
data_points(2,2) = 0.737412;
data_points(3,0) = 0.166461;
data_points(3,1) = 0.149912,
data_points(3,1) = 0.149912;
data_points(3,2) = 0.364944;
// one simplex
std::vector<Matrix3d> simplex(4);
Matrix3d v0(3,3);
Matrix3d v1(3,3);
Matrix3d v2(3,3);
Matrix3d v3(3,3);
std::vector<Matrix> simplex(4);
Matrix v0(3,3);
Matrix v1(3,3);
Matrix v2(3,3);
Matrix v3(3,3);
v0(0,0) = -0.2192721;
v0(0,1) = 0.2792986,
v0(0,2) = -0.9348326,
v0(0,1) = 0.2792986;
v0(0,2) = -0.9348326;
v0(1,0) = -0.7772152;
v0(1,1) = -0.6292092;
v0(1,2) = -0.0056861;
@ -93,58 +92,58 @@ void test_nelder_mead()
std::size_t nm_iterations = 19;
CGAL::Optimal_bounding_box::nelder_mead<Linear_algebra_traits>(simplex, data_points, nm_iterations);
CGAL_assertion_code(double epsilon = 1e-5);
CGAL_assertion_code(Matrix3d v0_new = simplex[0]);
CGAL_assertion(assert_doubles(v0_new(0,0), -0.288975, epsilon));
CGAL_assertion(assert_doubles(v0_new(0,1), 0.7897657, epsilon));
CGAL_assertion(assert_doubles(v0_new(0,2), -0.541076, epsilon));
CGAL_assertion(assert_doubles(v0_new(1,0), -0.9407046, epsilon));
CGAL_assertion(assert_doubles(v0_new(1,1), -0.3391466, epsilon));
CGAL_assertion(assert_doubles(v0_new(1,2), 0.0073817, epsilon));
CGAL_assertion(assert_doubles(v0_new(2,0), -0.1776743, epsilon));
CGAL_assertion(assert_doubles(v0_new(2,1), 0.5111260, epsilon));
CGAL_assertion(assert_doubles(v0_new(2,2), 0.84094, epsilon));
double epsilon = 1e-5;
const Matrix& v0_new = simplex[0];
assert(assert_doubles(v0_new(0,0), -0.288975, epsilon));
assert(assert_doubles(v0_new(0,1), 0.7897657, epsilon));
assert(assert_doubles(v0_new(0,2), -0.541076, epsilon));
assert(assert_doubles(v0_new(1,0), -0.9407046, epsilon));
assert(assert_doubles(v0_new(1,1), -0.3391466, epsilon));
assert(assert_doubles(v0_new(1,2), 0.0073817, epsilon));
assert(assert_doubles(v0_new(2,0), -0.1776743, epsilon));
assert(assert_doubles(v0_new(2,1), 0.5111260, epsilon));
assert(assert_doubles(v0_new(2,2), 0.84094, epsilon));
CGAL_assertion_code(Matrix3d v1_new = simplex[1]);
CGAL_assertion(assert_doubles(v1_new(0,0), -0.458749, epsilon));
CGAL_assertion(assert_doubles(v1_new(0,1), 0.823283, epsilon));
CGAL_assertion(assert_doubles(v1_new(0,2), -0.334296, epsilon));
CGAL_assertion(assert_doubles(v1_new(1,0), -0.885235, epsilon));
CGAL_assertion(assert_doubles(v1_new(1,1), -0.455997, epsilon));
CGAL_assertion(assert_doubles(v1_new(1,2), 0.091794, epsilon));
CGAL_assertion(assert_doubles(v1_new(2,0), -0.076866, epsilon));
CGAL_assertion(assert_doubles(v1_new(2,1), 0.338040, epsilon));
CGAL_assertion(assert_doubles(v1_new(2,2), 0.937987, epsilon));
const Matrix& v1_new = simplex[1];
assert(assert_doubles(v1_new(0,0), -0.458749, epsilon));
assert(assert_doubles(v1_new(0,1), 0.823283, epsilon));
assert(assert_doubles(v1_new(0,2), -0.334296, epsilon));
assert(assert_doubles(v1_new(1,0), -0.885235, epsilon));
assert(assert_doubles(v1_new(1,1), -0.455997, epsilon));
assert(assert_doubles(v1_new(1,2), 0.091794, epsilon));
assert(assert_doubles(v1_new(2,0), -0.076866, epsilon));
assert(assert_doubles(v1_new(2,1), 0.338040, epsilon));
assert(assert_doubles(v1_new(2,2), 0.937987, epsilon));
CGAL_assertion_code(Matrix3d v2_new = simplex[2]);
CGAL_assertion(assert_doubles(v2_new(0,0), -0.346582, epsilon));
CGAL_assertion(assert_doubles(v2_new(0,1), 0.878534, epsilon));
CGAL_assertion(assert_doubles(v2_new(0,2), -0.328724, epsilon));
CGAL_assertion(assert_doubles(v2_new(1,0), -0.936885, epsilon));
CGAL_assertion(assert_doubles(v2_new(1,1), -0.341445, epsilon));
CGAL_assertion(assert_doubles(v2_new(1,2), 0.075251, epsilon));
CGAL_assertion(assert_doubles(v2_new(2,0), -0.046131, epsilon));
CGAL_assertion(assert_doubles(v2_new(2,1), 0.334057, epsilon));
CGAL_assertion(assert_doubles(v2_new(2,2), 0.941423, epsilon));
const Matrix& v2_new = simplex[2];
assert(assert_doubles(v2_new(0,0), -0.346582, epsilon));
assert(assert_doubles(v2_new(0,1), 0.878534, epsilon));
assert(assert_doubles(v2_new(0,2), -0.328724, epsilon));
assert(assert_doubles(v2_new(1,0), -0.936885, epsilon));
assert(assert_doubles(v2_new(1,1), -0.341445, epsilon));
assert(assert_doubles(v2_new(1,2), 0.075251, epsilon));
assert(assert_doubles(v2_new(2,0), -0.046131, epsilon));
assert(assert_doubles(v2_new(2,1), 0.334057, epsilon));
assert(assert_doubles(v2_new(2,2), 0.941423, epsilon));
CGAL_assertion_code(Matrix3d v3_new = simplex[3]);
CGAL_assertion(assert_doubles(v3_new(0,0), -0.394713, epsilon));
CGAL_assertion(assert_doubles(v3_new(0,1), 0.791782, epsilon));
CGAL_assertion(assert_doubles(v3_new(0,2), -0.466136, epsilon));
CGAL_assertion(assert_doubles(v3_new(1,0), -0.912112, epsilon));
CGAL_assertion(assert_doubles(v3_new(1,1), -0.398788, epsilon));
CGAL_assertion(assert_doubles(v3_new(1,2), 0.094972, epsilon));
CGAL_assertion(assert_doubles(v3_new(2,0), -0.110692, epsilon));
CGAL_assertion(assert_doubles(v3_new(2,1), 0.462655, epsilon));
CGAL_assertion(assert_doubles(v3_new(2,2), 0.879601, epsilon));
const Matrix& v3_new = simplex[3];
assert(assert_doubles(v3_new(0,0), -0.394713, epsilon));
assert(assert_doubles(v3_new(0,1), 0.791782, epsilon));
assert(assert_doubles(v3_new(0,2), -0.466136, epsilon));
assert(assert_doubles(v3_new(1,0), -0.912112, epsilon));
assert(assert_doubles(v3_new(1,1), -0.398788, epsilon));
assert(assert_doubles(v3_new(1,2), 0.094972, epsilon));
assert(assert_doubles(v3_new(2,0), -0.110692, epsilon));
assert(assert_doubles(v3_new(2,1), 0.462655, epsilon));
assert(assert_doubles(v3_new(2,2), 0.879601, epsilon));
}
void test_genetic_algorithm()
{
CGAL::Eigen_dense_matrix<double, -1, -1> data_points(4, 3); // -1 = dynamic size at run time
data_points(0,0) = 0.866802;
data_points(0,1) = 0.740808,
data_points(0,2) = 0.895304,
data_points(0,1) = 0.740808;
data_points(0,2) = 0.895304;
data_points(1,0) = 0.912651;
data_points(1,1) = 0.761565;
data_points(1,2) = 0.160330;
@ -152,14 +151,14 @@ void test_genetic_algorithm()
data_points(2,1) = 0.892578;
data_points(2,2) = 0.737412;
data_points(3,0) = 0.166461;
data_points(3,1) = 0.149912,
data_points(3,1) = 0.149912;
data_points(3,2) = 0.364944;
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
CGAL::Optimal_bounding_box::Population<Linear_algebra_traits> pop(5);
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> evolution(pop, data_points);
evolution.genetic_algorithm();
CGAL_assertion(pop.size() == 5);
assert(pop.size() == 5);
}
void test_random_unit_tetra()
@ -169,8 +168,8 @@ void test_random_unit_tetra()
// points are on their convex hull
data_points(0,0) = 0.866802;
data_points(0,1) = 0.740808,
data_points(0,2) = 0.895304,
data_points(0,1) = 0.740808;
data_points(0,2) = 0.895304;
data_points(1,0) = 0.912651;
data_points(1,1) = 0.761565;
data_points(1,2) = 0.160330;
@ -178,7 +177,7 @@ void test_random_unit_tetra()
data_points(2,1) = 0.892578;
data_points(2,2) = 0.737412;
data_points(3,0) = 0.166461;
data_points(3,1) = 0.149912,
data_points(3,1) = 0.149912;
data_points(3,2) = 0.364944;
typedef CGAL::Simple_cartesian<double> K;
@ -199,31 +198,31 @@ void test_random_unit_tetra()
out.close();
#endif
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
CGAL_assertion_code(typedef Linear_algebra_traits::Matrix3d Matrix3d);
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
typedef Linear_algebra_traits::Matrix Matrix;
std::size_t generations = 10;
CGAL::Optimal_bounding_box::Population<Linear_algebra_traits> pop(50);
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> evolution(pop, data_points);
evolution.evolve(generations);
CGAL_assertion_code(Matrix3d R = evolution.get_best());
CGAL_assertion_code(double epsilon = 1e-3);
CGAL_assertion(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
CGAL_assertion(assert_doubles(R(0,0), -0.25791, epsilon));
CGAL_assertion(assert_doubles(R(0,1), 0.796512, epsilon));
CGAL_assertion(assert_doubles(R(0,2), -0.546855, epsilon));
CGAL_assertion(assert_doubles(R(1,0), -0.947128, epsilon));
CGAL_assertion(assert_doubles(R(1,1), -0.320242, epsilon));
CGAL_assertion(assert_doubles(R(1,2), -0.0197553, epsilon));
CGAL_assertion(assert_doubles(R(2,0), -0.190861, epsilon));
CGAL_assertion(assert_doubles(R(2,1), 0.512847, epsilon));
CGAL_assertion(assert_doubles(R(2,2), 0.836992, epsilon));
Matrix R = evolution.get_best();
const double epsilon = 1e-3;
assert(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
assert(assert_doubles(R(0,0), -0.25791, epsilon));
assert(assert_doubles(R(0,1), 0.796512, epsilon));
assert(assert_doubles(R(0,2), -0.546855, epsilon));
assert(assert_doubles(R(1,0), -0.947128, epsilon));
assert(assert_doubles(R(1,1), -0.320242, epsilon));
assert(assert_doubles(R(1,2), -0.0197553, epsilon));
assert(assert_doubles(R(2,0), -0.190861, epsilon));
assert(assert_doubles(R(2,1), 0.512847, epsilon));
assert(assert_doubles(R(2,2), 0.836992, epsilon));
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
// postprocessing
CGAL::Eigen_dense_matrix<double> obb(8, 3);
CGAL::Optimal_bounding_box::post_processing(data_points, R, obb);
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, "data/random_unit_tetra_result.off");
#endif
}
@ -236,12 +235,11 @@ void test_reference_tetrahedron(const char* fname)
std::exit(1);
}
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
typedef Linear_algebra_traits::MatrixXd MatrixXd;
CGAL_assertion_code(typedef Linear_algebra_traits::Matrix3d Matrix3d);
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
typedef Linear_algebra_traits::Matrix Matrix;
// points in a matrix
MatrixXd points;
Matrix points;
CGAL::Optimal_bounding_box::sm_to_matrix(mesh, points);
std::size_t generations = 10;
@ -249,33 +247,32 @@ void test_reference_tetrahedron(const char* fname)
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> experiment(pop, points);
experiment.evolve(generations);
CGAL_assertion_code(Matrix3d R = experiment.get_best());
CGAL_assertion_code(double epsilon = 1e-5);
CGAL_assertion(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
Matrix R = experiment.get_best();
double epsilon = 1e-5;
assert(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
// postprocessing
MatrixXd obb(8, 3);
Matrix obb(8, 3);
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, "data/OBB.off");
#endif
#endif
}
void test_long_tetrahedron(std::string fname)
void test_long_tetrahedron(const std::string fname)
{
std::ifstream input(fname);
CGAL::Surface_mesh<K::Point_3> mesh;
if (!input || !(input >> mesh) || mesh.is_empty()) {
if (!input || !(input >> mesh) || mesh.is_empty())
{
std::cerr << fname << " is not a valid off file.\n";
std::exit(1);
}
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
typedef Linear_algebra_traits::MatrixXd MatrixXd;
CGAL_assertion_code(typedef Linear_algebra_traits::Matrix3d Matrix3d);
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
typedef Linear_algebra_traits::Matrix Matrix;
// points in a matrix
MatrixXd points;
Matrix points;
CGAL::Optimal_bounding_box::sm_to_matrix(mesh, points);
std::size_t max_generations = 10;
@ -283,42 +280,45 @@ void test_long_tetrahedron(std::string fname)
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> experiment(pop, points);
experiment.evolve(max_generations);
CGAL_assertion_code(Matrix3d R = experiment.get_best());
CGAL_assertion_code(double epsilon = 1e-3);
CGAL_assertion(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
CGAL_assertion(assert_doubles(R(0,0), -1, epsilon));
CGAL_assertion(assert_doubles(R(0,1), 0, epsilon));
CGAL_assertion(assert_doubles(R(0,2), 0, epsilon));
CGAL_assertion(assert_doubles(R(1,0), 0, epsilon));
CGAL_assertion(assert_doubles(R(1,1), -0.707107, epsilon));
CGAL_assertion(assert_doubles(R(1,2), 0.707106, epsilon) ||
assert_doubles(R(1,2), -0.707106, epsilon));
CGAL_assertion(assert_doubles(R(2,0), 0, epsilon));
CGAL_assertion(assert_doubles(R(2,1), 0.707106, epsilon) ||
assert_doubles(R(1,2), -0.707106, epsilon));
CGAL_assertion(assert_doubles(R(2,2), 0.707107, epsilon));
Matrix R = experiment.get_best();
double epsilon = 1e-3;
assert(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
assert(assert_doubles(R(0,0), -1, epsilon));
assert(assert_doubles(R(0,1), 0, epsilon));
assert(assert_doubles(R(0,2), 0, epsilon));
assert(assert_doubles(R(1,0), 0, epsilon));
assert(assert_doubles(R(1,1), -0.707107, epsilon));
assert(assert_doubles(R(1,2), 0.707106, epsilon) ||
assert_doubles(R(1,2), -0.707106, epsilon));
assert(assert_doubles(R(2,0), 0, epsilon));
assert(assert_doubles(R(2,1), 0.707106, epsilon) ||
assert_doubles(R(1,2), -0.707106, epsilon));
assert(assert_doubles(R(2,2), 0.707107, epsilon));
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
// postprocessing
MatrixXd obb(8, 3);
Matrix obb(8, 3);
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, fname + "result.off");
#endif
#endif
}
void test_compute_obb_evolution(std::string fname)
void test_compute_obb_evolution(const std::string fname)
{
std::ifstream input(fname);
typedef CGAL::Surface_mesh<K::Point_3> SMesh;
SMesh mesh;
if (!input || !(input >> mesh) || mesh.is_empty()) {
if (!input || !(input >> mesh) || mesh.is_empty())
{
std::cerr << fname << " is not a valid off file.\n";
std::exit(1);
}
// get mesh points
std::vector<K::Point_3> sm_points;
typedef typename boost::graph_traits<SMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<SMesh, boost::vertex_point_t>::const_type PointPMap;
typedef typename boost::graph_traits<SMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<SMesh, boost::vertex_point_t>::const_type PointPMap;
PointPMap pmap = get(boost::vertex_point, mesh);
for(vertex_descriptor v : vertices(mesh))
sm_points.push_back(get(pmap, v));
@ -328,26 +328,26 @@ void test_compute_obb_evolution(std::string fname)
std::vector<K::Point_3> obb_points;
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points, la_traits, true);
CGAL_assertion_code(double epsilon = 1e-3);
CGAL_assertion_code(double vol = CGAL::Optimal_bounding_box::calculate_volume(obb_points));
CGAL_assertion(assert_doubles(vol, 0.883371, epsilon));
double epsilon = 1e-3;
double vol = CGAL::Optimal_bounding_box::calculate_volume(obb_points);
assert(assert_doubles(vol, 0.883371, epsilon));
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
/*
for(int i = 0; i < 8; ++i)
std::cout << obb_points[i].x() << " " << obb_points[i].y() << " " << obb_points[i].z() << "\n" ;
*/
CGAL::Surface_mesh<K::Point_3> result_mesh;
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
obb_points[4], obb_points[5], obb_points[6], obb_points[7], result_mesh);
obb_points[4], obb_points[5], obb_points[6], obb_points[7], result_mesh);
std::ofstream out("data/obb_result.off");
out << result_mesh;
out.close();
#endif
#endif
}
void test_compute_obb_mesh(std::string fname)
void test_compute_obb_mesh(const std::string fname)
{
std::ifstream input(fname);
CGAL::Surface_mesh<K::Point_3> mesh;
@ -361,14 +361,15 @@ void test_compute_obb_mesh(std::string fname)
CGAL::Surface_mesh< K::Point_3> obbmesh;
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(mesh, obbmesh, la_traits, true);
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
std::ofstream out("/tmp/result_elephant.off");
out << obbmesh;
out.close();
#endif
#endif
}
void test_function_defaults_traits(std::string fname1, std::string fname2)
void test_function_defaults_traits(const std::string fname1,
const std::string fname2)
{
std::ifstream input1(fname1);
CGAL::Surface_mesh<K::Point_3> mesh1;
@ -389,6 +390,7 @@ void test_function_defaults_traits(std::string fname1, std::string fname2)
// test one
std::vector<K::Point_3> sm_points;
typedef CGAL::Surface_mesh<K::Point_3> SMesh;
typedef typename boost::graph_traits<SMesh>::vertex_descriptor vertex_descriptor;
typedef typename boost::property_map<SMesh, boost::vertex_point_t>::const_type PointPMap;
PointPMap pmap = get(boost::vertex_point, mesh1);
@ -399,12 +401,12 @@ void test_function_defaults_traits(std::string fname1, std::string fname2)
std::vector<K::Point_3> obb_points;
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points, true);
CGAL_assertion_code(double epsilon = 1e-3);
CGAL_assertion_code(double vol = CGAL::Optimal_bounding_box::calculate_volume(obb_points));
CGAL_assertion(assert_doubles(vol, 0.883371, epsilon));
const double epsilon = 1e-3;
const double vol = CGAL::Optimal_bounding_box::calculate_volume(obb_points);
assert(assert_doubles(vol, 0.883371, epsilon));
// test two
CGAL::Surface_mesh< K::Point_3> obbmesh;
CGAL::Surface_mesh<K::Point_3> obbmesh;
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(mesh2, obbmesh, true);
}

2
Polyhedron/demo/Polyhedron/Plugins/PCA/CMakeLists.txt
View File

@ -18,7 +18,7 @@ polyhedron_demo_plugin(create_bbox_mesh_plugin Create_bbox_mesh_plugin)
target_link_libraries(create_bbox_mesh_plugin PUBLIC scene_surface_mesh_item)
polyhedron_demo_plugin(create_obb_mesh_plugin Create_obb_mesh_plugin)
target_link_libraries(create_obb_mesh_plugin PUBLIC scene_surface_mesh_item scene_polyhedron_item scene_polyhedron_selection_item scene_points_with_normal_item)
target_link_libraries(create_obb_mesh_plugin PUBLIC scene_surface_mesh_item scene_selection_item scene_points_with_normal_item)
qt5_wrap_ui( volumesUI_FILES Basic_generator_widget.ui)
polyhedron_demo_plugin(basic_generator_plugin Basic_generator_plugin ${volumesUI_FILES} KEYWORDS PolygonMesh PointSetProcessing)

99
Polyhedron/demo/Polyhedron/Plugins/PCA/Create_obb_mesh_plugin.cpp
View File

@ -1,29 +1,25 @@
#include <QtCore/qglobal.h>
#include <CGAL/Three/Scene_item.h>
#include <CGAL/Three/Scene_interface.h>
#include <QAction>
#include <QMainWindow>
#include <QMessageBox>
#include <QApplication>
#include "Scene_surface_mesh_item.h"
#include "Polyhedron_type.h"
#include "Scene_polyhedron_item.h"
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
#include "Scene_polyhedron_selection_item.h"
#include "Scene_points_with_normal_item.h"
#include <boost/graph/graph_traits.hpp>
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
#include <CGAL/Optimal_bounding_box/optimal_bounding_box.h>
#include <CGAL/Eigen_linear_algebra_traits.h>
//typedef Scene_surface_mesh_item Scene_facegraph_item;
//typedef Scene_facegraph_item Scene_facegraph_item;
typedef Scene_facegraph_item::Face_graph FaceGraph;
typedef Polyhedron::Point_3 Point_3;
using namespace CGAL::Three;
typedef Scene_surface_mesh_item Scene_facegraph_item;
class Create_obb_mesh_plugin :
public QObject,
public CGAL::Three::Polyhedron_demo_plugin_interface
@ -36,36 +32,21 @@ public:
void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface*);
QList<QAction*> actions() const;
bool applicable(QAction*) const {
/*
if (scene->selectionIndices().size() == 1)
{
return qobject_cast<Scene_facegraph_item*>(scene->item(scene->mainSelectionIndex()))
|| qobject_cast<Scene_polyhedron_selection_item*>(scene->item(scene->mainSelectionIndex()));
}
Q_FOREACH(int index, scene->selectionIndices())
{
if (qobject_cast<Scene_facegraph_item*>(scene->item(index)))
return true;
}
return false;
*/
bool applicable(QAction*) const
{
if(scene->mainSelectionIndex() != -1
&& scene->item(scene->mainSelectionIndex())->isFinite())
return true;
return false;
}
return false;
}
protected:
void gather_mesh_points(std::vector<Point_3>& points);
void obb();
public Q_SLOTS:
void createObb() {
void createObb()
{
QApplication::setOverrideCursor(Qt::WaitCursor);
obb();
QApplication::restoreOverrideCursor();
@ -75,11 +56,8 @@ private:
Scene_interface* scene;
QMainWindow* mw;
QAction* actionObb;
}; // end Create_obb_mesh_plugin class
void Create_obb_mesh_plugin::init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface*)
{
scene = scene_interface;
@ -89,17 +67,16 @@ void Create_obb_mesh_plugin::init(QMainWindow* mainWindow, Scene_interface* scen
connect(actionObb, SIGNAL(triggered()), this, SLOT(createObb()));
}
QList<QAction*> Create_obb_mesh_plugin::actions() const {
QList<QAction*> Create_obb_mesh_plugin::actions() const
{
return QList<QAction*>() << actionObb;
}
void Create_obb_mesh_plugin::gather_mesh_points(std::vector<Point_3>& points)
{
const Scene_interface::Item_id index = scene->mainSelectionIndex();
Scene_facegraph_item* poly_item =
qobject_cast<Scene_facegraph_item*>(scene->item(index));
Scene_facegraph_item* item = qobject_cast<Scene_facegraph_item*>(scene->item(index));
Scene_polyhedron_selection_item* selection_item =
qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index));
@ -107,7 +84,7 @@ void Create_obb_mesh_plugin::gather_mesh_points(std::vector<Point_3>& points)
Scene_points_with_normal_item* point_set_item =
qobject_cast<Scene_points_with_normal_item*>(scene->item(index));
if(poly_item || selection_item)
if(item || selection_item)
{
typedef typename boost::property_map<FaceGraph, boost::vertex_point_t>::type PointPMap;
typedef typename boost::graph_traits<FaceGraph>::vertex_descriptor vertex_descriptor;
@ -115,30 +92,29 @@ void Create_obb_mesh_plugin::gather_mesh_points(std::vector<Point_3>& points)
std::vector<vertex_descriptor> selected_vertices;
if(poly_item != NULL)
if(item != NULL)
{
FaceGraph& pmesh = *poly_item->polyhedron();
FaceGraph& pmesh = *item->polyhedron();
selected_vertices.assign(vertices(pmesh).begin(), vertices(pmesh).end());
PointPMap pmap = get(CGAL::vertex_point, pmesh);
BOOST_FOREACH(vertex_descriptor v, selected_vertices)
for(vertex_descriptor v : selected_vertices)
points.push_back(get(pmap, v));
}
else if(selection_item != NULL) // using selection of faces
{
FaceGraph& pmesh = *selection_item->polyhedron();
BOOST_FOREACH(face_descriptor f, selection_item->selected_facets)
for(face_descriptor f : selection_item->selected_facets)
{
BOOST_FOREACH(vertex_descriptor v, vertices_around_face(halfedge(f, pmesh), pmesh))
{
for(vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh))
selected_vertices.push_back(v);
}
}
PointPMap pmap = get(CGAL::vertex_point, pmesh);
BOOST_FOREACH(vertex_descriptor v, selected_vertices)
for(vertex_descriptor v : selected_vertices)
points.push_back(get(pmap, v));
}
CGAL_assertion(points.size() >= 3);
}
@ -149,12 +125,9 @@ void Create_obb_mesh_plugin::gather_mesh_points(std::vector<Point_3>& points)
return;
std::cout << "points_set->size()= " << points_set->size() << std::endl;
BOOST_FOREACH(Point_3 p, points_set->points())
{
for(const Point_3& p : points_set->points())
points.push_back(p);
}
}
}
void Create_obb_mesh_plugin::obb()
@ -164,24 +137,14 @@ void Create_obb_mesh_plugin::obb()
gather_mesh_points(points);
// find obb
CGAL::Eigen_linear_algebra_traits la_traits;
std::vector<Point_3> obb_points(8);
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(points, obb_points, la_traits, true);
std::array<Point_3, 8> obb_points;
CGAL::optimal_bounding_box(points, obb_points);
Scene_item* item;
if(mw->property("is_polyhedorn_mode").toBool())
{
Polyhedron* p = new Polyhedron;
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
obb_points[4], obb_points[5], obb_points[6], obb_points[7], *p);
item = new Scene_polyhedron_item(p);
}
else {
SMesh* p = new SMesh;
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
obb_points[4], obb_points[5], obb_points[6], obb_points[7], *p);
item = new Scene_surface_mesh_item(p);
}
Scene_facegraph_item* item;
SMesh* p = new SMesh;
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
obb_points[4], obb_points[5], obb_points[6], obb_points[7], *p);
item = new Scene_facegraph_item(p);
item->setName("Optimal bbox mesh");
item->setRenderingMode(Wireframe);

322
Solver_interface/include/CGAL/Eigen_linear_algebra_traits.h
View File

@ -1,322 +0,0 @@
// Copyright (c) 2018 GeometryFactory (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$
// SPDX-License-Identifier: GPL-3.0+
//
//
// Author(s) : Konstantinos Katrioplas
#ifndef CGAL_EIGEN_LINEAR_ALGEBRA_TRAITS_H
#define CGAL_EIGEN_LINEAR_ALGEBRA_TRAITS_H
#include <CGAL/basic.h>
#include <Eigen/Dense>
#include <Eigen/QR>
namespace CGAL {
template<typename T, int D = Eigen::Dynamic>
class Eigen_dense_vector;
/*!
\ingroup PkgSolver
The class `Eigen_dense_matrix` is a wrapper around \ref thirdpartyEigen "Eigen" matrix type
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::DenseMatrix`</a>.
\tparam T Number type.
\tparam D1 Number of rows, or Dynamic.
\tparam D1 Number of cols, or Dynamic.
\sa `CGAL::Eigen_dense_vector<T, D>`
\sa `CGAL::Eigen_linear_algebra_traits`
*/
template<typename T, int D1 = Eigen::Dynamic, int D2 = Eigen::Dynamic>
class Eigen_dense_matrix
{
public:
/// The internal matrix type from \ref thirdpartyEigen "Eigen".
typedef Eigen::Matrix<T, D1, D2> EigenType;
/// Create a dense matrix
Eigen_dense_matrix(std::size_t nrows, std::size_t ncols)
: m_matrix(static_cast<int>(nrows), static_cast<int>(ncols))
{
CGAL_assertion(m_matrix.rows() > 0);
CGAL_assertion(m_matrix.cols() > 0);
}
/// Create a dense matrix
Eigen_dense_matrix(int nrows, int ncols)
: m_matrix(nrows, ncols)
{
CGAL_assertion(m_matrix.rows() > 0);
CGAL_assertion(m_matrix.cols() > 0);
}
/// Create a dense matrix out of a \ref thirdpartyEigen "Eigen" matrix type
Eigen_dense_matrix(const EigenType& eigen_mat)
: m_matrix(eigen_mat) {}
Eigen_dense_matrix() : m_matrix() {}
/// Read access to a matrix coefficient.
///
/// \pre 0 <= i < row_dimension().
/// \pre 0 <= j < column_dimension().
T& operator() (int i_, int j_)
{
return m_matrix(i_, j_);
}
/// Write access to a matrix coefficient: a_ij <- val
///
/// \pre 0 <= i < row_dimension().
/// \pre 0 <= j < column_dimension().
void set_coef(std::size_t i_, std::size_t j_, T val)
{
int i = static_cast<int>(i_);
int j = static_cast<int>(j_);
CGAL_precondition(i < m_matrix.rows());
CGAL_precondition(j < m_matrix.cols());
m_matrix.coeffRef(i,j) = val;
}
/// Return the matrix number of rows
std::size_t rows() const {return m_matrix.rows();}
/// Return the matrix number of cols
std::size_t cols() const {return m_matrix.cols();}
/// Resize to i rows and j cols
void resize(int i_, int j_) { m_matrix.resize(i_, j_);}
const T& coeff(int i_) const
{
return m_matrix.coeff(i_);
}
EigenType m_matrix;
};
/*!
\ingroup PkgSolver
The class `Eigen_vector` is a wrapper around \ref thirdpartyEigen "Eigen" dense 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`.
\tparam T Number type.
\tparam D Number of colums, or Dynamic.
\sa `CGAL::Eigen_dense_matrix<T>`
\sa `CGAL::Eigen_linear_algebra_traits<T>`
*/
template <typename T, int D>
class Eigen_dense_vector
{
private:
/// The internal vector type from \ref thirdpartyEigen "Eigen".
typedef Eigen::Matrix<T, D, 1> EigenType;
public:
/// Create a dense vector out of a \ref thirdpartyEigen "Eigen" vector type
Eigen_dense_vector(const EigenType& vec) : m_vector(vec) {}
/// Read and write and access to a vector coefficient: `a_i`
const T& coeff(std::size_t i)
{
CGAL_assertion(i >= 0);
CGAL_assertion(i < static_cast<std::size_t>(D));
return m_vector.coeff(i);
}
EigenType m_vector;
};
/*!
\ingroup PkgSolver
The class `Eigen_linear_algebra_traits` provides an interface to linear algebra functionalities of \ref thirdpartyEigen "Eigen".
\ref thirdpartyEigen "Eigen" version 3.1 (or later) must be available on the system.
\sa `CGAL::Eigen_dense_matrix<T, D1, D2>`
\sa `CGAL::Eigen_dense_vector<T, D>`
\sa http://eigen.tuxfamily.org
Example
--------------
\code{.cpp}
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
// dynamic matrix at run time to store a large amount of data
typedef Linear_algebra_traits::MatrixXd MatrixXd;
// preallocated 3x3 matrix at compile time
typedef Linear_algebra_traits::Matrix3d Matrix3d;
// preallocated 3-cols vector at compile time
typedef Linear_algebra_traits::Vector3d Vector3d;
\endcode
*/
class Eigen_linear_algebra_traits
{
public:
typedef double NT;
typedef int Index;
// dynamic size at run time
typedef CGAL::Eigen_dense_matrix<NT> MatrixXd;
// dynamic rows in run time, fixed cols in compile time
typedef CGAL::Eigen_dense_matrix<NT, Eigen::Dynamic, 3> MatrixX3d;
// fixed at compile time
typedef CGAL::Eigen_dense_matrix<NT, 3, 3> Matrix3d;
// fixed at compile time
typedef CGAL::Eigen_dense_vector<NT, 3> Vector3d;
/// Get the transpose of a `CGAL::Eigen_dense_matrix<T, D1, D2>` matrix
template <class Matrix>
static Matrix transpose(const Matrix& mat)
{
return Matrix(mat.m_matrix.transpose());
}
/// Get the determinant of a `CGAL::Eigen_dense_matrix<T, D1, D2>` matrix
template <class Matrix>
static NT determinant(const Matrix& mat)
{
return mat.m_matrix.determinant();
}
/// Performs QR decomposition of matrix A to a unitary matrix and an upper triagonal
/// and returns the unitary matrix.
template <class NT, int D1, int D2>
static CGAL::Eigen_dense_matrix<NT, D1, D2> qr_factorization(const CGAL::Eigen_dense_matrix<NT, D1, D2>& A)
{
Eigen::HouseholderQR<Eigen::Matrix<NT, D1, D2> > qr(A.m_matrix);
return CGAL::Eigen_dense_matrix<NT, D1, D2>(qr.householderQ());
}
template <class Matrix>
static void qr_factorization(std::vector<Matrix>& simplex)
{
for(std::size_t i = 0; i < simplex.size(); ++i)
{
Matrix mat = simplex[i].m_matrix;
simplex[i] = qr_factorization(mat);
}
}
// CGAL::Eigen_dense_vector<NT, D2> : the returned type with D2 may be -1 in compile time,
// and may not be equal to the expected type.
// Eigen manages to return a precompiled row out of a dynamic matrix but I don't know how.
/// Get the row vector out of a `CGAL::Eigen_dense_matrix<T, D1, D2>`. The result is stored in a
/// preallocated at compile time 3-column `CGAL::Eigen_dense_vector<T, 3>`
template <class NT, int D1, int D2>
static CGAL::Eigen_dense_vector<NT, 3> row3(const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
int i)
{
return CGAL::Eigen_dense_vector<NT, 3>(A.m_matrix.row(i));
}
};
/// Matrix multiplication. If the columns of A and the rows of B are equal at compile time,
/// the product is stored at a preallocated at compile time `CGAL::Eigen_dense_matrix`. Otherwise,
/// the product is stored in a dynamic at run time matrix.
template <class NT, int D1, int D2, int D3>
const CGAL::Eigen_dense_matrix<NT, D1, D3> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2 >& A,
const CGAL::Eigen_dense_matrix<NT, D2, D3 >& B)
{
return CGAL::Eigen_dense_matrix<NT, D1, D3>(A.m_matrix * B.m_matrix);
}
// D2 and D3 may not be equal at compile time, but equal at run time!
// This overload returns a dynamic matrix.
template <class NT, int D1, int D2, int D3, int D4>
const CGAL::Eigen_dense_matrix<NT> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2 >& A,
const CGAL::Eigen_dense_matrix<NT, D3, D4 >& B)
{
return CGAL::Eigen_dense_matrix<NT>(A.m_matrix * B.m_matrix);
}
// scalar - matrix multiplication
template <class NT, int D1, int D2>
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator* (const NT& scalar,
const CGAL::Eigen_dense_matrix<NT, D1, D2>& B)
{
return CGAL::Eigen_dense_matrix<NT, D1, D2>(scalar * B.m_matrix);
}
template <class NT, int D1, int D2>
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
const NT& scalar)
{
return CGAL::Eigen_dense_matrix<NT, D1, D2>(A.m_matrix * scalar);
}
template <class NT, int D1, int D2>
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator/ (const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
const double& scalar)
{
return CGAL::Eigen_dense_matrix<NT, D1, D2>(A.m_matrix / scalar);
}
template <class NT, int D1, int D2>
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator/ (const double& scalar,
const CGAL::Eigen_dense_matrix<NT, D1, D2> & A)
{
return CGAL::Eigen_dense_matrix<NT, D1, D2> (scalar / A.m_matrix);
}
// addition
template <class NT, int D1, int D2>
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator+ (const CGAL::Eigen_dense_matrix<NT, D1, D2> & A,
const CGAL::Eigen_dense_matrix<NT, D1, D2> & B)
{
return CGAL::Eigen_dense_matrix<NT, D1, D2> (A.m_matrix + B.m_matrix);
}
// vector - matrix multiplication
template <class NT, int D1, int D2>
const Eigen_dense_vector<NT, D1> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
const CGAL::Eigen_dense_vector<NT, D2>& V)
{
return Eigen_dense_vector<NT, D1>(A.m_matrix * V.m_vector);
}
} // end namespace
#endif // CGAL_EIGEN_LINEAR_ALGEBRA_TRAITS_H