mirror of https://github.com/CGAL/cgal
Merge branch 'CGAL-Optimal_bounding_box-GF-old' into CGAL-Optimal_bounding_box-GF
This commit is contained in:
Mael Rouxel-Labbé
commit
3e365756ca
|
|
@ -71,6 +71,7 @@ Nef_S2
|
|||
NewKernel_d
|
||||
Number_types
|
||||
OpenNL
|
||||
Optimal_bounding_box
|
||||
Optimal_transportation_reconstruction_2
|
||||
Optimisation_basic
|
||||
Partition_2
|
||||
|
|
|
|||
|
|
@ -131,6 +131,9 @@ CGAL_add_named_parameter(max_number_of_proxies_t, max_number_of_proxies, max_num
|
|||
CGAL_add_named_parameter(min_error_drop_t, min_error_drop, min_error_drop)
|
||||
CGAL_add_named_parameter(number_of_relaxations_t, number_of_relaxations, number_of_relaxations)
|
||||
|
||||
// List of named parameters used in Optimal_bounding_box package
|
||||
CGAL_add_named_parameter(use_convex_hull_t, use_convex_hull, use_convex_hull)
|
||||
|
||||
// meshing parameters
|
||||
CGAL_add_named_parameter(subdivision_ratio_t, subdivision_ratio, subdivision_ratio)
|
||||
CGAL_add_named_parameter(relative_to_chord_t, relative_to_chord, relative_to_chord)
|
||||
|
|
|
|||
|
|
@ -86,21 +86,24 @@ void test(const NamedParameters& np)
|
|||
assert(get_parameter(np, CGAL::internal_np::require_same_orientation).v == 49);
|
||||
assert(get_parameter(np, CGAL::internal_np::use_bool_op_to_clip_surface).v == 50);
|
||||
assert(get_parameter(np, CGAL::internal_np::face_size_map).v == 52);
|
||||
assert(get_parameter(np, CGAL::internal_np::snapping_tolerance).v == 57);
|
||||
assert(get_parameter(np, CGAL::internal_np::use_angle_smoothing).v == 53);
|
||||
assert(get_parameter(np, CGAL::internal_np::use_area_smoothing).v == 54);
|
||||
assert(get_parameter(np, CGAL::internal_np::use_Delaunay_flips).v == 55);
|
||||
assert(get_parameter(np, CGAL::internal_np::use_safety_constraints).v == 56);
|
||||
assert(get_parameter(np, CGAL::internal_np::area_threshold).v == 57);
|
||||
assert(get_parameter(np, CGAL::internal_np::volume_threshold).v == 58);
|
||||
assert(get_parameter(np, CGAL::internal_np::dry_run).v == 59);
|
||||
assert(get_parameter(np, CGAL::internal_np::do_lock_mesh).v == 60);
|
||||
assert(get_parameter(np, CGAL::internal_np::do_simplify_border).v == 61);
|
||||
assert(get_parameter(np, CGAL::internal_np::snapping_tolerance).v == 59);
|
||||
assert(get_parameter(np, CGAL::internal_np::dry_run).v == 60);
|
||||
assert(get_parameter(np, CGAL::internal_np::do_lock_mesh).v == 61);
|
||||
assert(get_parameter(np, CGAL::internal_np::do_simplify_border).v == 62);
|
||||
|
||||
// Named parameters that we use in the package 'Surface Mesh Simplification'
|
||||
assert(get_parameter(np, CGAL::internal_np::get_cost_policy).v == 34);
|
||||
assert(get_parameter(np, CGAL::internal_np::get_placement_policy).v == 35);
|
||||
|
||||
// Named parameters that we use in the package 'Optimal_bounding_box'
|
||||
assert(get_parameter(np, CGAL::internal_np::use_convex_hull).v == 63);
|
||||
|
||||
// To-be-documented named parameters
|
||||
assert(get_parameter(np, CGAL::internal_np::face_normal).v == 36);
|
||||
assert(get_parameter(np, CGAL::internal_np::random_seed).v == 37);
|
||||
|
|
@ -175,21 +178,24 @@ void test(const NamedParameters& np)
|
|||
check_same_type<49>(get_parameter(np, CGAL::internal_np::require_same_orientation));
|
||||
check_same_type<50>(get_parameter(np, CGAL::internal_np::use_bool_op_to_clip_surface));
|
||||
check_same_type<52>(get_parameter(np, CGAL::internal_np::face_size_map));
|
||||
check_same_type<57>(get_parameter(np, CGAL::internal_np::snapping_tolerance));
|
||||
check_same_type<53>(get_parameter(np, CGAL::internal_np::use_angle_smoothing));
|
||||
check_same_type<54>(get_parameter(np, CGAL::internal_np::use_area_smoothing));
|
||||
check_same_type<55>(get_parameter(np, CGAL::internal_np::use_Delaunay_flips));
|
||||
check_same_type<56>(get_parameter(np, CGAL::internal_np::use_safety_constraints));
|
||||
check_same_type<57>(get_parameter(np, CGAL::internal_np::area_threshold));
|
||||
check_same_type<58>(get_parameter(np, CGAL::internal_np::volume_threshold));
|
||||
check_same_type<59>(get_parameter(np, CGAL::internal_np::dry_run));
|
||||
check_same_type<60>(get_parameter(np, CGAL::internal_np::do_lock_mesh));
|
||||
check_same_type<61>(get_parameter(np, CGAL::internal_np::do_simplify_border));
|
||||
check_same_type<59>(get_parameter(np, CGAL::internal_np::snapping_tolerance));
|
||||
check_same_type<60>(get_parameter(np, CGAL::internal_np::dry_run));
|
||||
check_same_type<61>(get_parameter(np, CGAL::internal_np::do_lock_mesh));
|
||||
check_same_type<62>(get_parameter(np, CGAL::internal_np::do_simplify_border));
|
||||
|
||||
// Named parameters that we use in the package 'Surface Mesh Simplification'
|
||||
check_same_type<34>(get_parameter(np, CGAL::internal_np::get_cost_policy));
|
||||
check_same_type<35>(get_parameter(np, CGAL::internal_np::get_placement_policy));
|
||||
|
||||
// Named parameters that we use in the package 'Optimal_bounding_box'
|
||||
check_same_type<63>(get_parameter(np, CGAL::internal_np::use_convex_hull));
|
||||
|
||||
// To-be-documented named parameters
|
||||
check_same_type<36>(get_parameter(np, CGAL::internal_np::face_normal));
|
||||
check_same_type<37>(get_parameter(np, CGAL::internal_np::random_seed));
|
||||
|
|
@ -270,16 +276,17 @@ int main()
|
|||
.use_bool_op_to_clip_surface(A<50>(50))
|
||||
.use_binary_mode(A<51>(51))
|
||||
.face_size_map(A<52>(52))
|
||||
.snapping_tolerance(A<57>(57))
|
||||
.use_angle_smoothing(A<53>(53))
|
||||
.use_area_smoothing(A<54>(54))
|
||||
.use_Delaunay_flips(A<55>(55))
|
||||
.use_safety_constraints(A<56>(56))
|
||||
.use_safety_constraints(A<56>(56)
|
||||
.area_threshold(A<57>(57))
|
||||
.volume_threshold(A<58>(58))
|
||||
.dry_run(A<59>(59))
|
||||
.do_lock_mesh(A<60>(60))
|
||||
.do_simplify_border(A<61>(61))
|
||||
.snapping_tolerance(A<59>(59))
|
||||
.dry_run(A<60>(60))
|
||||
.do_lock_mesh(A<61>(61))
|
||||
.do_simplify_border(A<62>(62))
|
||||
.use_convex_hull(A<63>(63))
|
||||
);
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
|
|
|
|||
|
|
@ -130,6 +130,7 @@
|
|||
\package_listing{Box_intersection_d}
|
||||
\package_listing{AABB_tree}
|
||||
\package_listing{Spatial_sorting}
|
||||
\package_listing{Optimal_bounding_box}
|
||||
|
||||
\cgalPackageSection{PartGeometricOptimization,Geometric Optimization}
|
||||
|
||||
|
|
|
|||
|
|
@ -405,6 +405,17 @@ Boissonnat}
|
|||
,update = "97.08 kettner"
|
||||
}
|
||||
|
||||
@article{ cgal:cgm-fobbo-11,
|
||||
title={Fast oriented bounding box optimization on the rotation group SO (3, ℝ)},
|
||||
author={Chang, Chia-Tche and Gorissen, Bastien and Melchior, Samuel},
|
||||
journal={ACM Transactions on Graphics (TOG)},
|
||||
volume={30},
|
||||
number={5},
|
||||
pages={122},
|
||||
year={2011},
|
||||
publisher={ACM}
|
||||
}
|
||||
|
||||
@article{ cgal:cht-oacov-90
|
||||
,author = {M. Chang and N. Huang and C. Tang}
|
||||
,title = {An optimal algorithm for constructing oriented Voronoi diagrams
|
||||
|
|
|
|||
|
|
@ -0,0 +1,54 @@
|
|||
// Copyright (c) 2016 GeometryFactory SARL (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org)
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Andreas Fabri
|
||||
//
|
||||
// Warning: this file is generated, see include/CGAL/licence/README.md
|
||||
|
||||
#ifndef CGAL_LICENSE_OPTIMAL_BOUNDING_BOX_H
|
||||
#define CGAL_LICENSE_OPTIMAL_BOUNDING_BOX_H
|
||||
|
||||
#include <CGAL/config.h>
|
||||
#include <CGAL/license.h>
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE
|
||||
|
||||
# if CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE
|
||||
|
||||
# if defined(CGAL_LICENSE_WARNING)
|
||||
|
||||
CGAL_pragma_warning("Your commercial license for CGAL does not cover "
|
||||
"this release of the Optimal Bounding Box package.")
|
||||
# endif
|
||||
|
||||
# ifdef CGAL_LICENSE_ERROR
|
||||
# error "Your commercial license for CGAL does not cover this release \
|
||||
of the Optimal Bounding Box package. \
|
||||
You get this error, as you defined CGAL_LICENSE_ERROR."
|
||||
# endif // CGAL_LICENSE_ERROR
|
||||
|
||||
# endif // CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE < CGAL_RELEASE_DATE
|
||||
|
||||
#else // no CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE
|
||||
|
||||
# if defined(CGAL_LICENSE_WARNING)
|
||||
CGAL_pragma_warning("\nThe macro CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE is not defined."
|
||||
"\nYou use the CGAL Optimal Bounding Box package under "
|
||||
"the terms of the GPLv3+.")
|
||||
# endif // CGAL_LICENSE_WARNING
|
||||
|
||||
# ifdef CGAL_LICENSE_ERROR
|
||||
# error "The macro CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE is not defined.\
|
||||
You use the CGAL Optimal Bounding Box package under the terms of \
|
||||
the GPLv3+. You get this error, as you defined CGAL_LICENSE_ERROR."
|
||||
# endif // CGAL_LICENSE_ERROR
|
||||
|
||||
#endif // no CGAL_OPTIMAL_BOUNDING_BOX_COMMERCIAL_LICENSE
|
||||
|
||||
#endif // CGAL_LICENSE_OPTIMAL_BOUNDING_BOX_H
|
||||
|
|
@ -32,6 +32,7 @@ Minkowski_sum_3 3D Minkowski Sum of Polyhedra
|
|||
Nef_2 2D Boolean Operations on Nef Polygons
|
||||
Nef_3 3D Boolean Operations on Nef Polyhedra
|
||||
Nef_S2 2D Boolean Operations on Nef Polygons Embedded on the Sphere
|
||||
Optimal_bounding_box Optimal Bounding Box
|
||||
Optimal_transportation_reconstruction_2 Optimal Transportation Curve Reconstruction
|
||||
Partition_2 2D Polygon Partitioning
|
||||
Periodic_2_triangulation_2 2D Periodic Triangulations
|
||||
|
|
|
|||
|
|
@ -0,0 +1,29 @@
|
|||
# Created by the script cgal_create_CMakeLists
|
||||
# This is the CMake script for compiling a set of CGAL applications.
|
||||
|
||||
cmake_minimum_required(VERSION 3.1...3.15)
|
||||
project( Optimal_bounding_box_Benchmark )
|
||||
|
||||
# CGAL and its components
|
||||
find_package( CGAL QUIET )
|
||||
|
||||
if ( NOT CGAL_FOUND )
|
||||
message(STATUS "This project requires the CGAL library, and will not be compiled.")
|
||||
return()
|
||||
endif()
|
||||
|
||||
# include helper file
|
||||
include( ${CGAL_USE_FILE} )
|
||||
|
||||
find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater)
|
||||
if (NOT EIGEN3_FOUND)
|
||||
message(STATUS "This project requires the Eigen library, and will not be compiled.")
|
||||
return()
|
||||
else()
|
||||
include(${EIGEN3_USE_FILE})
|
||||
endif()
|
||||
|
||||
create_single_source_cgal_program("bench_obb.cpp")
|
||||
create_single_source_cgal_program("bench_perfomance.cpp")
|
||||
create_single_source_cgal_program("bench_custom.cpp")
|
||||
create_single_source_cgal_program("bench_fitness_function.cpp")
|
||||
|
|
@ -0,0 +1,64 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#include <CGAL/Optimal_bounding_box/optimal_bounding_box.h>
|
||||
#include <CGAL/subdivision_method_3.h>
|
||||
#include <CGAL/Timer.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
|
||||
//#define CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_BENCHMARK
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
void bench(const char* fname)
|
||||
{
|
||||
std::vector<K::Point_3> sm_points, obb_points;
|
||||
std::ifstream in(fname);
|
||||
|
||||
K::Point_3 p;
|
||||
int i = 0;
|
||||
while(in >> p)
|
||||
{
|
||||
if(i % 2 == 0) // avoid normals
|
||||
sm_points.push_back(p);
|
||||
|
||||
++i;
|
||||
}
|
||||
|
||||
std::cout << "input data (points + normals)= " << i << std::endl;
|
||||
std::cout << "number of points= " << sm_points.size() << std::endl;
|
||||
|
||||
CGAL::Eigen_linear_algebra_traits la_traits;
|
||||
|
||||
// use convex hull - > true
|
||||
// no convex hull - > false
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points, la_traits, false);
|
||||
|
||||
std::cout << "done" << '\n';
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_BENCHMARK
|
||||
std::cout.precision(17);
|
||||
for(int i =0; i < obb_points.size(); i++)
|
||||
std::cout << obb_points[i] << std::endl;
|
||||
|
||||
CGAL::Surface_mesh<K::Point_3> 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], mesh);
|
||||
|
||||
std::ofstream out("/tmp/result_obb.off");
|
||||
out << mesh;
|
||||
out.close();
|
||||
#endif
|
||||
}
|
||||
|
||||
int main(int argc, char* argv[])
|
||||
{
|
||||
bench(argv[1]);
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,64 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#include <CGAL/Optimal_bounding_box/fitness_function.h>
|
||||
#include <CGAL/Optimal_bounding_box/helper.h>
|
||||
#include <CGAL/Optimal_bounding_box/population.h>
|
||||
|
||||
#include <CGAL/subdivision_method_3.h>
|
||||
#include <CGAL/Timer.h>
|
||||
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
int main()
|
||||
{
|
||||
const char* fname = "data/elephant.off";
|
||||
|
||||
// 1) import a lot a mesh and subdivide it to create a big mesh
|
||||
std::ifstream input(fname);
|
||||
|
||||
CGAL::Surface_mesh<K::Point_3> mesh;
|
||||
if (!input || !(input >> mesh) || mesh.is_empty())
|
||||
{
|
||||
std::cerr << fname << " is not a valid off file.\n";
|
||||
std::exit(1);
|
||||
}
|
||||
|
||||
CGAL::Subdivision_method_3::CatmullClark_subdivision(mesh,
|
||||
CGAL::parameters::number_of_iterations(6));
|
||||
|
||||
int nb_points = static_cast<int>(vertices(mesh).size());
|
||||
std::cout << "number of points= " << nb_points << std::endl;
|
||||
|
||||
|
||||
// 2) fill a Matrix with them
|
||||
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
|
||||
typedef Linear_algebra_traits::MatrixX3d MatrixX3d;
|
||||
typedef Linear_algebra_traits::Matrix3d Matrix3d;
|
||||
|
||||
MatrixX3d points_mat(nb_points, 3);
|
||||
CGAL::Optimal_bounding_box::sm_to_matrix(mesh, points_mat);
|
||||
|
||||
// 3) create a population of simplices
|
||||
CGAL::Optimal_bounding_box::Population<Linear_algebra_traits> pop(50);
|
||||
|
||||
CGAL::Timer timer;
|
||||
timer.start();
|
||||
|
||||
// 4) compute fitness of population via the Fitness map
|
||||
CGAL::Optimal_bounding_box::Fitness_map<Linear_algebra_traits,
|
||||
Matrix3d, MatrixX3d> fit_map(pop, points_mat);
|
||||
double result = fit_map.get_best_fitness_value();
|
||||
|
||||
timer.stop();
|
||||
|
||||
std::cout << "took " << timer.time() << " to compute the fitness of all vertices.\n";
|
||||
std::cout << "value of fittest vertex= " << result << std::endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,106 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
#include <CGAL/Polyhedron_3.h>
|
||||
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#include <CGAL/Optimal_bounding_box/population.h>
|
||||
#include <CGAL/Optimal_bounding_box/optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/convex_hull_3.h>
|
||||
#include <CGAL/Timer.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
|
||||
//#define CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_BENCHMARK
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
bool assert_doubles(double d1, double d2, double epsilon)
|
||||
{
|
||||
return (d1 < d2 + epsilon && d1 > d2 - epsilon) ? true : false;
|
||||
}
|
||||
|
||||
template <typename SurfaceMesh, typename Point>
|
||||
void gather_mesh_points(SurfaceMesh& mesh, std::vector<Point>& points)
|
||||
{
|
||||
typedef typename boost::graph_traits<SurfaceMesh>::vertex_descriptor vertex_descriptor;
|
||||
typedef typename boost::property_map<SurfaceMesh, CGAL::vertex_point_t>::type PointPMap;
|
||||
PointPMap pmap = get(boost::vertex_point, mesh);
|
||||
for(vertex_descriptor v : vertices(mesh))
|
||||
points.push_back(get(pmap, v));
|
||||
}
|
||||
|
||||
void bench_finding_obb(std::string fname)
|
||||
{
|
||||
std::ifstream input(fname);
|
||||
|
||||
CGAL::Eigen_linear_algebra_traits la_traits;
|
||||
std::vector<K::Point_3> sm_points;
|
||||
|
||||
// import a mesh
|
||||
CGAL::Surface_mesh<K::Point_3> mesh;
|
||||
if (!input || !(input >> mesh) || mesh.is_empty())
|
||||
{
|
||||
std::cerr << fname << " is not a valid off file.\n";
|
||||
std::exit(1);
|
||||
}
|
||||
|
||||
// get mesh points
|
||||
gather_mesh_points(mesh, sm_points);
|
||||
|
||||
CGAL::Timer timer;
|
||||
|
||||
// 1) measure convex hull calculation
|
||||
timer.start();
|
||||
CGAL::Polyhedron_3<K> poly;
|
||||
convex_hull_3(sm_points.begin(), sm_points.end(), poly);
|
||||
std::vector<K::Point_3> ch_points(poly.points_begin(), poly.points_end());
|
||||
timer.stop();
|
||||
std::cout << "takes : " << timer.time() << " seconds to find the convex hull\n";
|
||||
|
||||
// 2) using convex hull
|
||||
timer.reset();
|
||||
timer.start();
|
||||
std::vector<K::Point_3> obb_points1;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points1, la_traits, true);
|
||||
timer.stop();
|
||||
std::cout << "found obb using convex hull: " << timer.time() << " seconds\n";
|
||||
|
||||
// 3) without convex hull
|
||||
timer.reset();
|
||||
timer.start();
|
||||
std::vector<K::Point_3> obb_points2;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points2, la_traits, false);
|
||||
timer.stop();
|
||||
std::cout << "found obb without convex hull: " << timer.time() << " seconds\n";
|
||||
timer.reset();
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_BENCHMARK
|
||||
CGAL::Surface_mesh<K::Point_3> result_mesh1;
|
||||
CGAL::make_hexahedron(obb_points1[0], obb_points1[1], obb_points1[2], obb_points1[3],
|
||||
obb_points1[4], obb_points1[5], obb_points1[6], obb_points1[7],
|
||||
result_mesh1);
|
||||
|
||||
CGAL::Surface_mesh<K::Point_3> result_mesh2;
|
||||
CGAL::make_hexahedron(obb_points2[0], obb_points2[1], obb_points2[2], obb_points2[3],
|
||||
obb_points2[4], obb_points2[5], obb_points2[6], obb_points2[7],
|
||||
result_mesh2);
|
||||
|
||||
std::ofstream out1("data/obb_result1.off");
|
||||
out1 << result_mesh1;
|
||||
out1.close();
|
||||
|
||||
std::ofstream out2("data/obb_result2.off");
|
||||
out2 << result_mesh2;
|
||||
out2.close();
|
||||
#endif
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
bench_finding_obb("data/elephant.off");
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,85 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#include <CGAL/Optimal_bounding_box/optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/subdivision_method_3.h>
|
||||
#include <CGAL/Timer.h>
|
||||
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
//#define CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_BENCHMARK
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
void bench_finding_obb(std::string fname)
|
||||
{
|
||||
std::ifstream input(fname);
|
||||
|
||||
// import a mesh
|
||||
CGAL::Surface_mesh<K::Point_3> mesh;
|
||||
if (!input || !(input >> mesh) || mesh.is_empty())
|
||||
{
|
||||
std::cerr << fname << " is not a valid off file.\n";
|
||||
std::exit(1);
|
||||
}
|
||||
|
||||
// export some times
|
||||
std::ofstream outt("data/times.txt");
|
||||
outt << "nb_vertices "<< "with_ch " << "without_ch" << std::endl;
|
||||
|
||||
CGAL::Timer timer;
|
||||
std::size_t measurements = 4;
|
||||
CGAL::Eigen_linear_algebra_traits la_traits;
|
||||
|
||||
for (std::size_t t = 0; t < measurements; ++t)
|
||||
{
|
||||
std::cout << "#vertices= " << vertices(mesh).size() << " |";
|
||||
|
||||
// 1) using convex hull
|
||||
timer.start();
|
||||
CGAL::Surface_mesh<K::Point_3> obb_mesh1;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(mesh, obb_mesh1, la_traits, true);
|
||||
timer.stop();
|
||||
double t_ch = timer.time();
|
||||
std::cout << " with ch: " << timer.time() << " s |";
|
||||
|
||||
// 2) without convex hull
|
||||
timer.reset();
|
||||
timer.start();
|
||||
CGAL::Surface_mesh<K::Point_3> obb_mesh2;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(mesh, obb_mesh2, la_traits, false);
|
||||
timer.stop();
|
||||
double t_no_ch = timer.time();
|
||||
std::cout << " without ch: " << timer.time() << " s\n";
|
||||
timer.reset();
|
||||
|
||||
outt << vertices(mesh).size() << " " << t_ch << " " << t_no_ch << std::endl;
|
||||
|
||||
// 3) subdivision
|
||||
CGAL::Subdivision_method_3::CatmullClark_subdivision(mesh,
|
||||
CGAL::parameters::number_of_iterations(1));
|
||||
}
|
||||
|
||||
outt.close();
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_BENCHMARK
|
||||
std::ofstream out1("data/obb_result1.off");
|
||||
out1 << obb_points1;
|
||||
out1.close();
|
||||
|
||||
std::ofstream out2("data/obb_result2.off");
|
||||
out2 << obb_points2;
|
||||
out2.close();
|
||||
#endif
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
bench_finding_obb("data/elephant.off");
|
||||
|
||||
return 0;
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
Binary file not shown.
|
After Width: | Height: | Size: 25 KiB |
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,5 @@
|
|||
nb_vertices with_ch without_ch
|
||||
2775 0.037331 0.327787
|
||||
16670 0.118251 1.96901
|
||||
66692 0.362664 7.75459
|
||||
266780 1.20984 32.166
|
||||
|
|
@ -0,0 +1,19 @@
|
|||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
#path-to-benchmarks
|
||||
benchmarkfile='data/times.txt'
|
||||
|
||||
data = np.loadtxt(benchmarkfile, skiprows = 1)
|
||||
|
||||
x = data[:, 0]
|
||||
y1 = data[:, 1]
|
||||
y2 = data[:, 2]
|
||||
|
||||
plt.plot(x, y1, 'go--', label='with convex hull')
|
||||
plt.plot(x, y2, 'ro--', label='without convex hull')
|
||||
legend = plt.legend(loc='best')
|
||||
plt.xlabel('#vertices')
|
||||
plt.ylabel('seconds')
|
||||
|
||||
plt.show()
|
||||
|
|
@ -0,0 +1 @@
|
|||
To draw a graph with the benchmark times set the path-to-the-measurments in draw_benchmark_times.py if necessary, and run the script with python 3.
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
/*!
|
||||
\ingroup PkgOptimalBoundingBoxConcepts
|
||||
\cgalConcept
|
||||
|
||||
The concept `OrientedBoundingBoxTraits_3` describes the requirements of the traits class
|
||||
used in the function `CGAL::oriented_bounding_box()`, and in particular the need for
|
||||
a 3x3 matrix type.
|
||||
|
||||
\cgalHasModel `CGAL::Oriented_bounding_box_traits_3`
|
||||
|
||||
*/
|
||||
class OrientedBoundingBoxTraits_3
|
||||
{
|
||||
public:
|
||||
/// The field number type; must be a model of the concept `FieldNumberType`
|
||||
typedef unspecified_type FT;
|
||||
|
||||
/// The 3D point type; must be model of `Kernel::Point_3`
|
||||
typedef unspecified_type Point_3;
|
||||
|
||||
/// The 3D affine transformation type; the template parameter `K` must be a model of `Kernel`
|
||||
/// and be compatible with the type `Point_3`
|
||||
typedef CGAL::Aff_transformation_3<K> Aff_transformation_3;
|
||||
|
||||
/// A construction object that must provide the function operator:
|
||||
/// `CGAL::Bbox_3 operator()(const Point_3&)`,
|
||||
/// which returns an axis-aligned bounding box that contains the point
|
||||
typedef unspecified_type Construct_bbox_3;
|
||||
|
||||
/// A 3x3 matrix type; model of the concept `SvdTraits::Matrix` and which supports
|
||||
/// matrix-matrix and scalar-matrix multiplication, as well as matrix-matrix addition
|
||||
typedef unspecified_type Matrix;
|
||||
|
||||
/// Returns the unitary matrix `Q` obtained in the QR-decomposition of the matrix `m`
|
||||
Matrix get_Q(const Matrix& m) const;
|
||||
|
||||
/// Returns the 3D box construction functor
|
||||
Construct_bbox_3 construct_bbox_3_object() const;
|
||||
};
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
|
||||
|
||||
PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Optimal Bounding Box"
|
||||
EXTRACT_ALL = false
|
||||
HIDE_UNDOC_CLASSES = true
|
||||
WARN_IF_UNDOCUMENTED = true
|
||||
|
||||
# macros to be used inside the code
|
||||
ALIASES += "cgalNamedParamsBegin=<dl class=\"params\"><dt>Named Parameters</dt><dd> <table class=\"params\">"
|
||||
ALIASES += "cgalNamedParamsEnd=</table> </dd> </dl>"
|
||||
ALIASES += "cgalParamBegin{1}=<tr><td class=\"paramname\">\ref OBB_\1 \"\1\"</td><td>"
|
||||
ALIASES += "cgalParamEnd=</td></tr>"
|
||||
|
||||
#macros for NamedParameters.txt
|
||||
ALIASES += "cgalNPTableBegin=<dl class=\"params\"><dt></dt><dd> <table class=\"params\">"
|
||||
ALIASES += "cgalNPTableEnd=</table> </dd> </dl>"
|
||||
ALIASES += "cgalNPBegin{1}=<tr><td class=\"paramname\">\1 </td><td>"
|
||||
ALIASES += "cgalNPEnd=</td></tr>"
|
||||
|
||||
MACRO_EXPANSION = YES
|
||||
EXPAND_ONLY_PREDEF = YES
|
||||
EXPAND_AS_DEFINED = CGAL_BGL_NP_TEMPLATE_PARAMETERS \
|
||||
CGAL_BGL_NP_CLASS
|
||||
EXCLUDE = ${CGAL_PACKAGE_INCLUDE_DIR}/CGAL/Optimal_bounding_box/internal
|
||||
EXCLUDE_SYMBOLS += experimental
|
||||
|
||||
|
|
@ -0,0 +1,65 @@
|
|||
/*!
|
||||
\defgroup obb_namedparameters Named Parameters for the package Optimal Bounding Box
|
||||
\ingroup PkgOptimalBoundingBoxRef
|
||||
|
||||
In this package, all functions optional parameters are implemented as BGL optional
|
||||
named parameters (see \ref BGLNamedParameters for more information on how to use them).
|
||||
Since the parameters of the various functions defined in this package are redundant,
|
||||
their long descriptions are centralized below.
|
||||
The sequence of named parameters should start with `CGAL::parameters::`.
|
||||
`CGAL::parameters::all_default()` can be used to indicate
|
||||
that default values of optional named parameters must be used.
|
||||
|
||||
In the following, we assume that the following types are provided as template parameters
|
||||
of functions and classes. Note that, for some of these functions,
|
||||
the type is more specific:
|
||||
<ul>
|
||||
<li>`PolygonMesh` is a model of the concept `FaceGraph`</li>.
|
||||
<li>`GeomTraits` a geometric traits class in which constructions are performed and
|
||||
predicates evaluated. Everywhere in this package, a \cgal `Kernel` fulfills the requirements.</li>
|
||||
</ul>
|
||||
|
||||
The following named parameters, offered by the package \ref PkgBGL
|
||||
(see \ref bgl_namedparameters), are used in this package:
|
||||
|
||||
\cgalNPTableBegin
|
||||
\cgalNPBegin{vertex_point_map} \anchor OBB_vertex_point_map
|
||||
is the property map with the points associated to the vertices of the polygon mesh `pmesh`.\n
|
||||
<b>Type:</b> a class model of `ReadablePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type and
|
||||
`GeomTraits::Point_3` as value type. \n
|
||||
<b>Default:</b> \code boost::get(CGAL::vertex_point, pmesh) \endcode
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{point_map} \anchor OBB_point_map
|
||||
is the property map containing the points associated to the elements of the point range `points`.\n
|
||||
<b>Type:</b> a class model of `ReadablePropertyMap` with `PointRange::iterator::value_type` as key type
|
||||
and `geom_traits::Point_3` as value type. \n
|
||||
<b>Default:</b> \code CGAL::Identity_property_map<geom_traits::Point_3>\endcode
|
||||
\cgalNPEnd
|
||||
\cgalNPTableEnd
|
||||
|
||||
In addition to these named parameters, this package offers the following named parameters:
|
||||
|
||||
\cgalNPTableBegin
|
||||
\cgalNPBegin{geom_traits} \anchor OBB_geom_traits
|
||||
is the geometric traits instance in which the mesh processing operation should be performed.\n
|
||||
<b>Type:</b> a Geometric traits class.\n
|
||||
<b>Default:</b>
|
||||
\code typename CGAL::Kernel_traits<
|
||||
typename boost::property_traits<
|
||||
typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type>::value_type>::Kernel \endcode
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{use_convex_hull} \anchor OBB_use_convex_hull
|
||||
Parameter used in the construction of oriented bounding box to indicate whether the algorithm should
|
||||
first extract the extreme points (points that are on the 3D convex hull) of the input data range
|
||||
to accelerate the computation of the bounding box.
|
||||
\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPTableEnd
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,16 @@
|
|||
namespace CGAL {
|
||||
|
||||
/*!
|
||||
\mainpage User Manual
|
||||
\anchor Chapter_Building_Optimal_Bounding_Box
|
||||
|
||||
\cgalAutoToc
|
||||
\authors Konstantinos Katrioplas, Mael Rouxel-Labbé
|
||||
|
||||
\section sectionSubIntro Introduction
|
||||
|
||||
AABB and stuff
|
||||
|
||||
\section secOBB Optimal Bounding Box
|
||||
*/
|
||||
} /* namespace CGAL */
|
||||
|
|
@ -0,0 +1,53 @@
|
|||
/// \defgroup PkgOptimalBoundingBoxRef Optimal Bounding Box Reference
|
||||
|
||||
/// \defgroup PkgOptimalBoundingBoxConcepts Concepts
|
||||
/// \ingroup PkgOptimalBoundingBoxRef
|
||||
|
||||
/// \defgroup PkgOptimalBoundingBoxClasses Optimal Bounding Box Classes
|
||||
/// \ingroup PkgOptimalBoundingBoxRef
|
||||
|
||||
/// \defgroup PkgOptimalBoundingBoxFunctions Optimal Bounding Box Functions
|
||||
/// \ingroup PkgOptimalBoundingBoxRef
|
||||
|
||||
/// \defgroup PkgOptimalBoundingBox_Oriented_bounding_box Oriented Bounding Box Functions
|
||||
/// \ingroup PkgOptimalBoundingBoxFunctions
|
||||
|
||||
/*!
|
||||
\addtogroup PkgOptimalBoundingBoxRef
|
||||
\cgalPkgDescriptionBegin{Optimal Bounding Box,PkgOptimalBoundingBox}
|
||||
\cgalPkgPicture{optimal_bounding_box.png}
|
||||
\cgalPkgSummaryBegin
|
||||
\cgalPkgAuthor{Konstantinos Katrioplas, Mael Rouxel-Labbé}
|
||||
\cgalPkgDesc{This package provides functions to compute tight oriented bounding boxes around a point set or a polygon mesh.}
|
||||
\cgalPkgManuals{Chapter_Building_Optimal_Bounding_Box,PkgOptimalBoundingBoxRef}
|
||||
\cgalPkgSummaryEnd
|
||||
\cgalPkgShortInfoBegin
|
||||
\cgalPkgSince{5.1}
|
||||
\cgalPkgDependsOn{\ref PkgConvexHull3}
|
||||
\cgalPkgBib{cgal:obb}
|
||||
\cgalPkgLicense{\ref licensesGPL "GPL"}
|
||||
\cgalPkgDemo{Polyhedron demo, polyhedron_3.zip}
|
||||
\cgalPkgShortInfoEnd
|
||||
\cgalPkgDescriptionEnd
|
||||
|
||||
\cgalClassifedRefPages
|
||||
|
||||
\cgalCRPSection{Parameters}
|
||||
|
||||
Optional parameters of the functions of this package are implemented as \ref BGLNamedParameters.
|
||||
The page \ref obb_namedparameters describes their usage and provides a list of the parameters
|
||||
that are used in this package.
|
||||
|
||||
\cgalCRPSection{Concepts}
|
||||
|
||||
- `OrientedBoundingBoxTraits_3`
|
||||
|
||||
\cgalCRPSection{Classes}
|
||||
|
||||
- `CGAL::Oriented_bounding_box_traits_3`
|
||||
|
||||
\cgalCRPSection{Methods}
|
||||
|
||||
- \link PkgOptimalBoundingBox_Oriented_bounding_box `CGAL::oriented_bounding_box()` \endlink
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,13 @@
|
|||
Manual
|
||||
Kernel_23
|
||||
STL_Extension
|
||||
Algebraic_foundations
|
||||
Circulator
|
||||
Stream_support
|
||||
Polyhedron
|
||||
BGL
|
||||
Solver_interface
|
||||
Surface_mesh
|
||||
AABB_tree
|
||||
Property_map
|
||||
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
/*!
|
||||
|
||||
\example Optimal_bounding_box/obb_example.cpp
|
||||
*/
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 11 KiB |
|
|
@ -0,0 +1,24 @@
|
|||
# Created by the script cgal_create_cmake_script
|
||||
# This is the CMake script for compiling a CGAL application.
|
||||
|
||||
cmake_minimum_required(VERSION 3.1...3.15)
|
||||
project( Optimal_bounding_box_Examples )
|
||||
|
||||
find_package(CGAL QUIET)
|
||||
|
||||
if (NOT CGAL_FOUND)
|
||||
message(STATUS "This project requires the CGAL library, and will not be compiled.")
|
||||
return()
|
||||
endif()
|
||||
|
||||
include( ${CGAL_USE_FILE} )
|
||||
|
||||
find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater)
|
||||
if (NOT EIGEN3_FOUND)
|
||||
message(STATUS "This project requires the Eigen library, and will not be compiled.")
|
||||
return()
|
||||
else()
|
||||
include(${EIGEN3_USE_FILE})
|
||||
endif()
|
||||
|
||||
create_single_source_cgal_program("obb_example.cpp")
|
||||
|
|
@ -0,0 +1,42 @@
|
|||
#define CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/optimal_bounding_box.h>
|
||||
#include <CGAL/Real_timer.h>
|
||||
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
typedef K::Point_3 Point;
|
||||
typedef CGAL::Surface_mesh<Point> Surface_mesh;
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
std::ifstream input(argv[1]);
|
||||
Surface_mesh sm;
|
||||
if (!input || !(input >> sm) || sm.is_empty())
|
||||
{
|
||||
std::cerr << argv[1] << " is not a valid off file.\n";
|
||||
return EXIT_FAILURE;
|
||||
}
|
||||
|
||||
CGAL::Real_timer timer;
|
||||
timer.start();
|
||||
|
||||
std::array<Point, 8> obb_points;
|
||||
CGAL::oriented_bounding_box(sm, obb_points);
|
||||
|
||||
std::cout << "Elapsed time: " << timer.time() << std::endl;
|
||||
|
||||
// Make a mesh out of the oriented bounding box
|
||||
Surface_mesh obb_sm;
|
||||
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], obb_sm);
|
||||
|
||||
std::ofstream("obb.off") << obb_sm;
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -0,0 +1,78 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
//
|
||||
// Author(s) : Mael Rouxel-Labbé
|
||||
//
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_BOUNDING_BOX_TRAITS_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_BOUNDING_BOX_TRAITS_H
|
||||
|
||||
#include <CGAL/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/Aff_transformation_3.h>
|
||||
|
||||
#ifdef CGAL_EIGEN3_ENABLED
|
||||
#include <CGAL/Eigen_matrix.h>
|
||||
#include <Eigen/QR>
|
||||
#endif
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
#if defined(CGAL_EIGEN3_ENABLED) || defined(DOXYGEN_RUNNING)
|
||||
|
||||
/// \ingroup PkgOptimalBoundingBoxClasses
|
||||
///
|
||||
/// The class `CGAL::Oriented_bounding_box_traits_3` is a traits type
|
||||
/// to be used with the overloads of the function `CGAL::oriented_bounding_box()`.
|
||||
///
|
||||
/// \attention This class requires the \ref thirdpartyEigen "Eigen" library.
|
||||
///
|
||||
/// \tparam K must be a model of `Kernel`
|
||||
///
|
||||
/// \cgalModels `OrientedBoundingBoxTraits_3`
|
||||
///
|
||||
template <typename K>
|
||||
class Oriented_bounding_box_traits_3
|
||||
{
|
||||
public:
|
||||
/// The field number type
|
||||
typedef typename K::FT FT;
|
||||
|
||||
/// The point type
|
||||
typedef typename K::Point_3 Point_3;
|
||||
|
||||
/// The affine transformation type
|
||||
typedef typename CGAL::Aff_transformation_3<K> Aff_transformation_3;
|
||||
|
||||
/// The axis-aligned bounding box construction object
|
||||
typedef typename K::Construct_bbox_3 Construct_bbox_3;
|
||||
|
||||
/// The matrix type
|
||||
typedef CGAL::Eigen_matrix<FT, 3, 3> Matrix;
|
||||
|
||||
private:
|
||||
typedef typename Matrix::EigenType EigenType;
|
||||
|
||||
public:
|
||||
/// Returns a default-constructed construction object
|
||||
Construct_bbox_3 construct_bbox_3_object() const { return Construct_bbox_3(); }
|
||||
|
||||
/// Performs the QR-decomposition of the matrix `m` to a unitary matrix and an upper triagonal
|
||||
/// and returns the unitary matrix
|
||||
Matrix get_Q(const Matrix& m) const
|
||||
{
|
||||
Eigen::HouseholderQR<EigenType> qr(m.eigen_object());
|
||||
return Matrix(EigenType(qr.householderQ()));
|
||||
}
|
||||
};
|
||||
#endif // defined(CGAL_EIGEN3_ENABLED) || defined(DOXYGEN_RUNNING)
|
||||
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_BOUNDING_BOX_TRAITS_H
|
||||
|
|
@ -0,0 +1,218 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// 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
|
||||
|
||||
#include <CGAL/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/internal/fitness_function.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/helper.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/nelder_mead_functions.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/population.h>
|
||||
|
||||
#include <CGAL/Random.h>
|
||||
#include <CGAL/number_utils.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Optimal_bounding_box {
|
||||
namespace internal {
|
||||
|
||||
template <typename PointRange, typename Traits>
|
||||
class Evolution
|
||||
{
|
||||
public:
|
||||
typedef typename Traits::FT FT;
|
||||
typedef typename Traits::Matrix Matrix;
|
||||
|
||||
typedef Population<Traits> Population;
|
||||
typedef typename Population::Simplex Simplex;
|
||||
|
||||
typedef Fitness_map<Population, PointRange, Traits> Fitness_map;
|
||||
|
||||
Evolution(const PointRange& points,
|
||||
CGAL::Random& rng,
|
||||
const Traits& traits)
|
||||
:
|
||||
m_population(traits),
|
||||
m_rng(rng),
|
||||
m_points(points),
|
||||
m_traits(traits)
|
||||
{ }
|
||||
|
||||
void genetic_algorithm(const std::size_t population_size = 50)
|
||||
{
|
||||
// random permutations
|
||||
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 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> group1(first_group_size), group2(first_group_size);
|
||||
std::vector<std::size_t> group3(second_group_size), group4(second_group_size);
|
||||
|
||||
int im = static_cast<int>(m);
|
||||
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
|
||||
FT bias = FT(0.1);
|
||||
|
||||
std::vector<Simplex> new_simplices(m);
|
||||
|
||||
for(std::size_t i=0; i<first_group_size; ++i)
|
||||
{
|
||||
std::array<Matrix, 4> offspring;
|
||||
for(int j=0; j<4; ++j)
|
||||
{
|
||||
const FT r = FT(m_rng.get_double());
|
||||
const FT fitnessA = compute_fitness<Traits>(m_population[group1[i]][j], m_points);
|
||||
const FT fitnessB = compute_fitness<Traits>(m_population[group2[i]][j], m_points);
|
||||
const FT threshold = (fitnessA < fitnessB) ? (0.5 + bias) : (0.5 - bias);
|
||||
|
||||
if(r < threshold)
|
||||
offspring[j] = m_population[group1[i]][j];
|
||||
else
|
||||
offspring[j] = m_population[group2[i]][j];
|
||||
}
|
||||
|
||||
new_simplices[i] = std::move(offspring);
|
||||
}
|
||||
|
||||
// crossover II
|
||||
bias = 0.1; // @fixme should the bias change? What should be the initial value?
|
||||
|
||||
for(std::size_t i=0; i<second_group_size; ++i)
|
||||
{
|
||||
std::array<Matrix, 4> offspring;
|
||||
for(int j=0; j<4; ++j)
|
||||
{
|
||||
const FT fitnessA = compute_fitness<Traits>(m_population[group3[i]][j], m_points);
|
||||
const FT fitnessB = compute_fitness<Traits>(m_population[group4[i]][j], m_points);
|
||||
const FT lambda = (fitnessA < fitnessB) ? (0.5 + bias) : (0.5 - bias);
|
||||
const FT rambda = 1 - lambda; // the 'l' in 'lambda' stands for left
|
||||
|
||||
// combine information from A and B
|
||||
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.get_Q(new_vertex);
|
||||
}
|
||||
|
||||
new_simplices[first_group_size + i] = std::move(offspring);
|
||||
}
|
||||
|
||||
m_population.simplices() = std::move(new_simplices);
|
||||
}
|
||||
|
||||
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 = 150;
|
||||
|
||||
// stopping criteria prameters
|
||||
FT prev_fit_value = 0.;
|
||||
FT new_fit_value = 0.;
|
||||
const FT tolerance = 1e-9;
|
||||
int stale = 0;
|
||||
|
||||
for(std::size_t t=0; t<generations; ++t)
|
||||
{
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
std::cout << "generation = " << t << "\n";
|
||||
#endif
|
||||
|
||||
genetic_algorithm(population_size);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
//std::cout << "pop after genetic" << std::endl;
|
||||
//pop.show_population();
|
||||
//std::cout << std::endl;
|
||||
#endif
|
||||
|
||||
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;
|
||||
#endif
|
||||
|
||||
new_fit_value = fitness_map.get_best_fitness_value();
|
||||
const FT difference = new_fit_value - prev_fit_value;
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
const Matrix& R_now = fitness_map.get_best();
|
||||
std::cout << "new best: " << R_now << std::endl;
|
||||
std::cout << "value difference with previous: " << difference << std::endl;
|
||||
#endif
|
||||
|
||||
if(CGAL::abs(difference) < tolerance * new_fit_value) // @todo should depend on input bbox diag
|
||||
++stale;
|
||||
|
||||
if(stale == 5)
|
||||
break;
|
||||
|
||||
prev_fit_value = new_fit_value;
|
||||
}
|
||||
}
|
||||
|
||||
const Matrix& get_best()
|
||||
{
|
||||
Fitness_map fitness_map(m_population, m_points);
|
||||
return fitness_map.get_best();
|
||||
}
|
||||
|
||||
private:
|
||||
Population m_population;
|
||||
CGAL::Random m_rng;
|
||||
const PointRange& m_points;
|
||||
const Traits& m_traits;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Optimal_bounding_box
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_EVOLUTION_H
|
||||
|
|
@ -0,0 +1,113 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// 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/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <limits>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Optimal_bounding_box {
|
||||
namespace internal {
|
||||
|
||||
template <typename Traits, typename PointRange>
|
||||
typename Traits::FT
|
||||
compute_fitness(const typename Traits::Matrix& R, // rotation matrix
|
||||
const PointRange& points)
|
||||
{
|
||||
typedef typename Traits::FT FT;
|
||||
typedef typename Traits::Point_3 Point;
|
||||
|
||||
CGAL_assertion(R.number_of_rows() == 3 && R.number_of_columns() == 3);
|
||||
CGAL_assertion(points.size() >= 3);
|
||||
|
||||
FT xmin, ymin, zmin, xmax, ymax, zmax;
|
||||
xmin = ymin = zmin = FT(std::numeric_limits<double>::max());
|
||||
xmax = ymax = zmax = FT(std::numeric_limits<double>::lowest());
|
||||
|
||||
for(const Point& pt : points)
|
||||
{
|
||||
const FT x = pt.x(), y = pt.y(), z = pt.z();
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
// volume
|
||||
return ((xmax - xmin) * (ymax - ymin) * (zmax - zmin));
|
||||
}
|
||||
|
||||
template <typename Population, typename PointRange, typename Traits>
|
||||
struct Fitness_map
|
||||
{
|
||||
typedef typename Traits::FT FT;
|
||||
typedef typename Population::Vertex Vertex;
|
||||
|
||||
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;
|
||||
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 = m_pop[i][j];
|
||||
const FT fitness = compute_fitness<Traits>(vertex, m_points);
|
||||
if(fitness < best_fitness)
|
||||
{
|
||||
simplex_id = i;
|
||||
vertex_id = j;
|
||||
best_fitness = fitness;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return m_pop[simplex_id][vertex_id];
|
||||
}
|
||||
|
||||
FT get_best_fitness_value() const
|
||||
{
|
||||
const Vertex& best_mat = get_best();
|
||||
return compute_fitness<Traits>(best_mat, m_points);
|
||||
}
|
||||
|
||||
private:
|
||||
const Population& m_pop;
|
||||
const PointRange& m_points;
|
||||
};
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Optimal_bounding_box
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_FITNESS_FUNCTION_H
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
// Mael Rouxel-Labbé
|
||||
//
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_INTERNAL_HELPER_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_INTERNAL_HELPER_H
|
||||
|
||||
#include <CGAL/license/Optimal_bounding_box.h>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Optimal_bounding_box {
|
||||
namespace internal {
|
||||
|
||||
template <typename Matrix>
|
||||
Matrix transpose(const Matrix& m)
|
||||
{
|
||||
Matrix tm;
|
||||
|
||||
tm.set(0, 0, m(0, 0)); tm.set(0, 1, m(1, 0)); tm.set(0, 2, m(2, 0));
|
||||
tm.set(1, 0, m(0, 1)); tm.set(1, 1, m(1, 1)); tm.set(1, 2, m(2, 1));
|
||||
tm.set(2, 0, m(0, 2)); tm.set(2, 1, m(1, 2)); tm.set(2, 2, m(2, 2));
|
||||
|
||||
return tm;
|
||||
}
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Optimal_bounding_box
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_INTERNAL_HELPER_H
|
||||
|
|
@ -0,0 +1,155 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// 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/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/internal/fitness_function.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/helper.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <array>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Optimal_bounding_box {
|
||||
namespace internal {
|
||||
|
||||
template <typename Matrix>
|
||||
Matrix reflection(const Matrix& S_centroid,
|
||||
const Matrix& S_worst)
|
||||
{
|
||||
return S_centroid * transpose(S_worst) * S_centroid;
|
||||
}
|
||||
|
||||
template <typename Matrix>
|
||||
Matrix expansion(const Matrix& S_centroid,
|
||||
const Matrix& S_worst,
|
||||
const Matrix& S_reflection)
|
||||
{
|
||||
return S_centroid * transpose(S_worst) * S_reflection;
|
||||
}
|
||||
|
||||
template <typename Matrix, typename Traits>
|
||||
Matrix mean(const Matrix& m1,
|
||||
const Matrix& m2,
|
||||
const Traits& traits)
|
||||
{
|
||||
typedef typename Traits::FT FT;
|
||||
|
||||
const Matrix reduction = FT(0.5) * (m1 + m2);
|
||||
|
||||
return traits.get_Q(reduction);
|
||||
}
|
||||
|
||||
template <typename Matrix, typename Traits>
|
||||
const Matrix nm_centroid(const Matrix& S1,
|
||||
const Matrix& S2,
|
||||
const Matrix& S3,
|
||||
const Traits& traits)
|
||||
{
|
||||
typedef typename Traits::FT FT;
|
||||
|
||||
const Matrix mean = (FT(1) / FT(3)) * (S1 + S2 + S3);
|
||||
|
||||
return traits.get_Q(mean);
|
||||
}
|
||||
|
||||
// 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)
|
||||
{
|
||||
typedef typename Traits::FT FT;
|
||||
typedef typename Traits::Matrix Matrix;
|
||||
|
||||
std::array<FT, 4> fitness;
|
||||
std::array<std::size_t, 4> indices = {{ 0, 1, 2, 3 }};
|
||||
|
||||
for(std::size_t t=0; t<nelder_mead_iterations; ++t)
|
||||
{
|
||||
for(std::size_t i=0; i<4; ++i)
|
||||
fitness[i] = compute_fitness<Traits>(simplex[i], points);
|
||||
|
||||
// get indices of sorted sequence
|
||||
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
|
||||
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 = std::move(s_simplex);
|
||||
fitness = std::move(s_fitness);
|
||||
|
||||
// centroid
|
||||
const Matrix v_centroid = nm_centroid(simplex[0], simplex[1], simplex[2], traits);
|
||||
|
||||
// find worst's vertex reflection
|
||||
const Matrix& v_worst = simplex[3];
|
||||
const Matrix v_refl = reflection(v_centroid, v_worst);
|
||||
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
|
||||
{
|
||||
// reflection
|
||||
simplex[3] = std::move(v_refl);
|
||||
}
|
||||
else
|
||||
{
|
||||
// expansion
|
||||
const Matrix v_expand = expansion(v_centroid, v_worst, v_refl);
|
||||
const FT f_expand = compute_fitness<Traits>(v_expand, points);
|
||||
if(f_expand < f_refl)
|
||||
simplex[3] = std::move(v_expand);
|
||||
else
|
||||
simplex[3] = std::move(v_refl);
|
||||
}
|
||||
}
|
||||
else // if reflected vertex is not better
|
||||
{
|
||||
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] = std::move(v_mean);
|
||||
}
|
||||
else
|
||||
{
|
||||
// reduction: move all vertices towards the best
|
||||
for(std::size_t i=1; i<4; ++i)
|
||||
simplex[i] = mean(simplex[i], simplex[0], traits);
|
||||
}
|
||||
}
|
||||
} // nelder mead iterations
|
||||
}
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Optimal_bounding_box
|
||||
} // namespace CGAL
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,109 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// 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/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
#include <CGAL/Random.h>
|
||||
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
namespace Optimal_bounding_box {
|
||||
namespace internal {
|
||||
|
||||
template<typename Traits>
|
||||
class Population
|
||||
{
|
||||
public:
|
||||
typedef typename Traits::FT FT;
|
||||
typedef typename Traits::Matrix Matrix;
|
||||
typedef Matrix Vertex;
|
||||
|
||||
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.get_Q(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
|
||||
|
||||
private:
|
||||
std::vector<Simplex> m_pop;
|
||||
const Traits& m_traits;
|
||||
};
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
template <typename Matrix>
|
||||
void Population<Matrix>::show_population()
|
||||
{
|
||||
std::size_t id = 0;
|
||||
for(const Simplex& simplex : m_pop)
|
||||
{
|
||||
std::cout << "Simplex: " << id++ << std::endl;
|
||||
for(const Matrix& R : simplex)
|
||||
std::cout << R << "\n\n";
|
||||
std:: cout << std:: endl;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
} // namespace internal
|
||||
} // namespace Optimal_bounding_box
|
||||
} // namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_POPULATION_H
|
||||
|
|
@ -0,0 +1,376 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
//
|
||||
// Author(s) : Mael Rouxel-Labbé
|
||||
//
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_ORIENTED_BOUNDING_BOX_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_ORIENTED_BOUNDING_BOX_H
|
||||
|
||||
#include <CGAL/license/Optimal_bounding_box.h>
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/internal/evolution.h>
|
||||
#include <CGAL/Optimal_bounding_box/internal/population.h>
|
||||
#include <CGAL/Optimal_bounding_box/Oriented_bounding_box_traits_3.h>
|
||||
|
||||
#include <CGAL/boost/graph/Named_function_parameters.h>
|
||||
#include <CGAL/boost/graph/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/Kernel_traits.h>
|
||||
#include <CGAL/Random.h>
|
||||
#include <CGAL/Simple_cartesian.h>
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
#include <CGAL/Real_timer.h>
|
||||
#endif
|
||||
|
||||
#include <array>
|
||||
#include <iostream>
|
||||
#include <iterator>
|
||||
#include <type_traits>
|
||||
#include <vector>
|
||||
|
||||
#ifdef DOXYGEN_RUNNING
|
||||
#define CGAL_BGL_NP_TEMPLATE_PARAMETERS NamedParameters
|
||||
#define CGAL_BGL_NP_CLASS NamedParameters
|
||||
#endif
|
||||
|
||||
namespace CGAL {
|
||||
namespace Optimal_bounding_box {
|
||||
namespace internal {
|
||||
|
||||
template <typename PointRange, typename Traits>
|
||||
void construct_oriented_bounding_box(std::array<typename Traits::Point_3, 8>& obb_points,
|
||||
const typename Traits::Aff_transformation_3& transformation,
|
||||
const typename Traits::Aff_transformation_3& inverse_transformation,
|
||||
const PointRange& points,
|
||||
const Traits& traits)
|
||||
{
|
||||
typedef typename Traits::Point_3 Point;
|
||||
|
||||
// Construct the bbox of the transformed point set
|
||||
CGAL::Bbox_3 bbox;
|
||||
for(const Point& pt : points)
|
||||
{
|
||||
const Point rotated_pt = transformation.transform(pt);
|
||||
bbox += traits.construct_bbox_3_object()(rotated_pt);
|
||||
}
|
||||
|
||||
obb_points[0] = Point(bbox.xmin(), bbox.ymin(), bbox.zmin());
|
||||
obb_points[1] = Point(bbox.xmax(), bbox.ymin(), bbox.zmin());
|
||||
obb_points[2] = Point(bbox.xmax(), bbox.ymax(), bbox.zmin());
|
||||
obb_points[3] = Point(bbox.xmin(), bbox.ymax(), bbox.zmin());
|
||||
|
||||
obb_points[4] = Point(bbox.xmin(), bbox.ymax(), bbox.zmax()); // see order in make_hexahedron()...
|
||||
obb_points[5] = Point(bbox.xmin(), bbox.ymin(), bbox.zmax());
|
||||
obb_points[6] = Point(bbox.xmax(), bbox.ymin(), bbox.zmax());
|
||||
obb_points[7] = Point(bbox.xmax(), bbox.ymax(), bbox.zmax());
|
||||
|
||||
// Apply the inverse rotation to the rotated axis aligned bounding box
|
||||
for(std::size_t i=0; i<8; ++i)
|
||||
obb_points[i] = inverse_transformation.transform(obb_points[i]);
|
||||
}
|
||||
|
||||
template <typename PointRange, typename Traits>
|
||||
void compute_best_transformation(typename Traits::Aff_transformation_3& transformation,
|
||||
typename Traits::Aff_transformation_3& inverse_transformation,
|
||||
const PointRange& points,
|
||||
CGAL::Random& rng,
|
||||
const Traits& traits)
|
||||
{
|
||||
typedef typename Traits::Matrix Matrix;
|
||||
typedef typename Traits::Aff_transformation_3 Aff_transformation_3;
|
||||
|
||||
const std::size_t max_generations = 50; // @todo hidden NP
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
CGAL::Real_timer timer;
|
||||
timer.start();
|
||||
#endif
|
||||
|
||||
Evolution<PointRange, Traits> search_solution(points, rng, traits);
|
||||
search_solution.evolve(max_generations);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
std::cout << "evolve: " << timer.time() << std::endl;
|
||||
timer.reset();
|
||||
#endif
|
||||
|
||||
const Matrix& rot = search_solution.get_best();
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
std::cout << "get best: " << timer.time() << std::endl;
|
||||
#endif
|
||||
|
||||
transformation = Aff_transformation_3(rot(0, 0), rot(0, 1), rot(0, 2),
|
||||
rot(1, 0), rot(1, 1), rot(1, 2),
|
||||
rot(2, 0), rot(2, 1), rot(2, 2));
|
||||
|
||||
// inverse transformation is simply the transposed since the matrix is unitary
|
||||
inverse_transformation = Aff_transformation_3(rot(0, 0), rot(1, 0), rot(2, 0),
|
||||
rot(0, 1), rot(1, 1), rot(2, 1),
|
||||
rot(0, 2), rot(1, 2), rot(2, 2));
|
||||
}
|
||||
|
||||
// Following two functions are overload to dispatch depending on return type
|
||||
template <typename PointRange, typename Traits>
|
||||
void construct_oriented_bounding_box(typename Traits::Aff_transformation_3& transformation,
|
||||
const PointRange& points,
|
||||
CGAL::Random& rng,
|
||||
const Traits& traits)
|
||||
{
|
||||
typename Traits::Aff_transformation_3 inverse_transformation;
|
||||
compute_best_transformation(transformation, inverse_transformation, points, rng, traits);
|
||||
}
|
||||
|
||||
template <typename PointRange, typename Traits>
|
||||
void construct_oriented_bounding_box(std::array<typename Traits::Point_3, 8>& obb_points,
|
||||
const PointRange& points,
|
||||
CGAL::Random& rng,
|
||||
const Traits& traits)
|
||||
{
|
||||
typename Traits::Aff_transformation_3 transformation, inverse_transformation;
|
||||
compute_best_transformation(transformation, inverse_transformation, points, rng, traits);
|
||||
|
||||
construct_oriented_bounding_box(obb_points, transformation, inverse_transformation, points, traits);
|
||||
}
|
||||
|
||||
// Entry point, decide whether to compute the CH_3 or not
|
||||
template <typename Output, typename PointRange, typename Traits>
|
||||
void construct_oriented_bounding_box(Output& output,
|
||||
const bool use_ch,
|
||||
const PointRange& points,
|
||||
CGAL::Random& rng,
|
||||
const Traits& traits)
|
||||
{
|
||||
typedef typename Traits::Matrix Matrix;
|
||||
typedef typename Traits::Point_3 Point;
|
||||
|
||||
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_oriented_bounding_box(output, ch_points, rng, traits);
|
||||
}
|
||||
else
|
||||
{
|
||||
return construct_oriented_bounding_box(output, points, rng, traits);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Optimal_bounding_box
|
||||
} // namespace internal
|
||||
|
||||
/// \addtogroup PkgOptimalBoundingBox_Oriented_bounding_box
|
||||
///
|
||||
/// The function `oriented_bounding_box` computes an approximation of the <i>optimal bounding box</i>,
|
||||
/// which is defined as the rectangular box with smallest volume of all the rectangular boxes containing
|
||||
/// the input points.
|
||||
///
|
||||
/// Internally, the algorithm uses an optimization process to compute a transformation (rotation)
|
||||
/// \f$ {\mathcal R}_b\f$ such that the axis-aligned box of the rotated input point set
|
||||
/// has a volume that is as small as possible given a fixed maximal number of optimization iterations.
|
||||
///
|
||||
/// \cgalHeading{Input}
|
||||
///
|
||||
/// The input can be either a range of points, or a polygon mesh.
|
||||
///
|
||||
/// \cgalHeading{Output}
|
||||
///
|
||||
/// The result of the algorithm can be retrieved as either:
|
||||
/// - the best affine transformation \f${\mathcal R}_b\f$ that the algorithm has found;
|
||||
/// - an array of eight points, representing the best oriented bounding box (\f${\mathcal B}_b\f$)
|
||||
/// that the algorithm has constructed, which is related to \f$ {\mathcal R}_b\f$ as it is
|
||||
/// the inverse transformation of the axis-aligned bounding box of the transformed point set.
|
||||
/// The order of the points in the array is the same as in the function
|
||||
/// \link PkgBGLHelperFct `CGAL::make_hexahedron()` \endlink,
|
||||
/// which can be used to construct a mesh from these points.
|
||||
///
|
||||
/// Note that when returning an array of points, these points are constructed from the axis-aligned
|
||||
/// bounding box and some precision loss should therefore be expected if a kernel not providing
|
||||
/// exact constructions is used.
|
||||
///
|
||||
/// The algorithm is based on a paper by Chang, Gorissen, and Melchior \cgalCite{cgal:cgm-fobbo-11}.
|
||||
|
||||
/// \ingroup PkgOptimalBoundingBox_Oriented_bounding_box
|
||||
///
|
||||
/// The function `oriented_bounding_box` computes an approximation of the <i>optimal bounding box</i>,
|
||||
/// which is defined as the rectangular box with smallest volume of all the rectangular boxes containing
|
||||
/// the input points.
|
||||
///
|
||||
/// See \ref PkgOptimalBoundingBox_Oriented_bounding_box for more information.
|
||||
///
|
||||
/// \tparam PointRange a model of `Range`. The value type may not be equal to the type `%Point_3` of the traits class
|
||||
/// if a point map is provided via named parameters (see below) to access points.
|
||||
/// \tparam Output either the type `Aff_transformation_3` of the traits class,
|
||||
/// or `std::array<Point, 8>` with `Point` being equivalent to the type `%Point_3` of the traits class
|
||||
/// \tparam NamedParameters a sequence of \ref obb_namedparameters "Named Parameters"
|
||||
///
|
||||
/// \param points the input range
|
||||
/// \param out the resulting array of points or affine transformation
|
||||
/// \param np an optional sequence of \ref obb_namedparameters "Named Parameters" among the ones listed below:
|
||||
///
|
||||
/// \cgalNamedParamsBegin
|
||||
/// \cgalParamBegin{point_map}
|
||||
/// a model of `ReadablePropertyMap` with value type the type `%Point_3` of the traits class.
|
||||
/// If this parameter is omitted, `CGAL::Identity_property_map<%Point_3>` is used.
|
||||
/// \cgalParamEnd
|
||||
/// \cgalParamBegin{geom_traits}
|
||||
/// a geometric traits class instance, model of the concept `OrientedBoundingBoxTraits_3`.
|
||||
/// %Default is a default construction object of type `CGAL::Oriented_bounding_box_traits_3<K>`
|
||||
/// where `K` is a kernel type deduced from the point type.
|
||||
/// \cgalParamEnd
|
||||
/// \cgalParamBegin{use_convex_hull}
|
||||
/// a Boolean value to indicate whether the algorithm should first extract the so-called extreme
|
||||
/// points of the data range (i.e. construct the convex hull) to reduce the input data range
|
||||
/// and accelerate the algorithm. %Default is `true`.
|
||||
/// \cgalParamEnd
|
||||
/// \cgalNamedParamsEnd
|
||||
///
|
||||
template <typename PointRange,
|
||||
typename Output,
|
||||
typename CGAL_BGL_NP_TEMPLATE_PARAMETERS>
|
||||
void oriented_bounding_box(const PointRange& points,
|
||||
Output& out,
|
||||
const CGAL_BGL_NP_CLASS& np
|
||||
#ifndef DOXYGEN_RUNNING
|
||||
, typename boost::enable_if<
|
||||
typename boost::has_range_iterator<PointRange>
|
||||
>::type* = 0
|
||||
#endif
|
||||
)
|
||||
{
|
||||
using CGAL::parameters::choose_parameter;
|
||||
using CGAL::parameters::get_parameter;
|
||||
|
||||
#if defined(CGAL_EIGEN3_ENABLED)
|
||||
typedef typename boost::range_value<PointRange>::type Point;
|
||||
typedef typename CGAL::Kernel_traits<Point>::type K;
|
||||
typedef Oriented_bounding_box_traits_3<K> Default_traits;
|
||||
#else
|
||||
typedef void Default_traits;
|
||||
#endif
|
||||
|
||||
typedef typename internal_np::Lookup_named_param_def<internal_np::geom_traits_t,
|
||||
CGAL_BGL_NP_CLASS,
|
||||
Default_traits>::type Geom_traits;
|
||||
|
||||
CGAL_static_assertion_msg(!(std::is_same<Geom_traits, void>::value),
|
||||
"You must provide a traits class or have Eigen enabled!");
|
||||
|
||||
Geom_traits traits = choose_parameter(get_parameter(np, internal_np::geom_traits), Geom_traits());
|
||||
|
||||
const bool use_ch = choose_parameter(get_parameter(np, internal_np::use_convex_hull), true);
|
||||
const unsigned int seed = choose_parameter(get_parameter(np, internal_np::random_seed), 0); // undocumented
|
||||
|
||||
CGAL::Random rng(seed);
|
||||
|
||||
// @todo handle those cases instead
|
||||
CGAL_assertion(points.size() >= 3);
|
||||
if(points.size() <= 3)
|
||||
{
|
||||
std::cerr << "The oriented bounding box cannot YET be computed for a mesh with fewer than 4 vertices!\n";
|
||||
return;
|
||||
}
|
||||
|
||||
return Optimal_bounding_box::internal::construct_oriented_bounding_box(out, use_ch, points, rng, traits);
|
||||
}
|
||||
|
||||
/// \ingroup PkgOptimalBoundingBox_Oriented_bounding_box
|
||||
///
|
||||
/// Extracts the vertices of the mesh as a point range and calls the overload using points as input.
|
||||
///
|
||||
/// \tparam PolygonMesh a model of `VertexListGraph`
|
||||
/// \tparam Output either the type `Aff_transformation_3` of the traits class,
|
||||
/// or `std::array<Point, 8>` with `Point` being equivalent to the type `%Point_3` of the traits class
|
||||
/// \tparam NamedParameters a sequence of \ref obb_namedparameters "Named Parameters"
|
||||
///
|
||||
/// \param pmesh the input mesh
|
||||
/// \param out the resulting array of points or affine transformation
|
||||
/// \param np an optional sequence of \ref obb_namedparameters "Named Parameters" among the ones listed below:
|
||||
///
|
||||
/// \cgalNamedParamsBegin
|
||||
/// \cgalParamBegin{vertex_point_map}
|
||||
/// the property map with the points associated to the vertices of `pmesh`.
|
||||
/// If this parameter is omitted, an internal property map for
|
||||
/// `CGAL::vertex_point_t` must be available in `PolygonMesh`
|
||||
/// \cgalParamEnd
|
||||
/// \cgalParamBegin{geom_traits}
|
||||
/// a geometric traits class instance, model of the concept `OrientedBoundingBoxTraits_3`.
|
||||
/// %Default is a default construction object of type `CGAL::Oriented_bounding_box_traits_3<K>`
|
||||
/// where `K` is a kernel type deduced from the point type.
|
||||
/// \cgalParamEnd
|
||||
/// \cgalParamBegin{use_convex_hull}
|
||||
/// a Boolean value to indicate whether the algorithm should first extract the so-called extreme
|
||||
/// points of the data range (i.e. construct the convex hull) to reduce the input data range
|
||||
/// and accelerate the algorithm. %Default is `true`.
|
||||
/// \cgalParamEnd
|
||||
/// \cgalNamedParamsEnd
|
||||
///
|
||||
template <typename PolygonMesh,
|
||||
typename Output,
|
||||
typename CGAL_BGL_NP_TEMPLATE_PARAMETERS>
|
||||
void oriented_bounding_box(const PolygonMesh& pmesh,
|
||||
Output& out,
|
||||
const CGAL_BGL_NP_CLASS& np
|
||||
#ifndef DOXYGEN_RUNNING
|
||||
, typename boost::disable_if<
|
||||
typename boost::has_range_iterator<PolygonMesh>
|
||||
>::type* = 0
|
||||
#endif
|
||||
)
|
||||
{
|
||||
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_BGL_NP_CLASS>::const_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;
|
||||
points.reserve(num_vertices(pmesh));
|
||||
|
||||
for(vertex_descriptor v : vertices(pmesh))
|
||||
points.push_back(get(vpm, v));
|
||||
|
||||
oriented_bounding_box(points, out, np);
|
||||
}
|
||||
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
/// Convenience overloads
|
||||
/////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
template <typename InputData /*range or mesh*/, typename OutputType /*array or transformation*/>
|
||||
void oriented_bounding_box(const InputData& data,
|
||||
OutputType& out)
|
||||
{
|
||||
return oriented_bounding_box(data, out, CGAL::parameters::all_default());
|
||||
}
|
||||
|
||||
/// \endcond
|
||||
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_ORIENTED_BOUNDING_BOX_H
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
// Copyright (c) 2018-2019 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// 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
|
||||
|
||||
#ifndef DOXYGEN_RUNNING
|
||||
#include <CGAL/license/Optimal_bounding_box.h>
|
||||
#endif
|
||||
|
||||
/**
|
||||
* \ingroup PkgOptimalBoundingBoxRef
|
||||
* \file CGAL/optimal_bounding_box.h
|
||||
* Convenience header file including the headers for all the classes and functions of this package.
|
||||
*/
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/Oriented_bounding_box_traits_3.h>
|
||||
#include <CGAL/Optimal_bounding_box/oriented_bounding_box.h>
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_OBB_H
|
||||
|
|
@ -0,0 +1 @@
|
|||
GeometryFactory (France)
|
||||
|
|
@ -0,0 +1,29 @@
|
|||
Algebraic_foundations
|
||||
Arithmetic_kernel
|
||||
BGL
|
||||
Cartesian_kernel
|
||||
Circulator
|
||||
Convex_hull_2
|
||||
Convex_hull_3
|
||||
Distance_2
|
||||
Distance_3
|
||||
Filtered_kernel
|
||||
Hash_map
|
||||
Homogeneous_kernel
|
||||
Installation
|
||||
Intersections_2
|
||||
Intersections_3
|
||||
Interval_support
|
||||
Kernel_23
|
||||
Kernel_d
|
||||
Modular_arithmetic
|
||||
Number_types
|
||||
Optimal_bounding_box
|
||||
Profiling_tools
|
||||
Property_map
|
||||
Random_numbers
|
||||
STL_Extension
|
||||
Solver_interface
|
||||
Stream_support
|
||||
TDS_2
|
||||
Triangulation_2
|
||||
|
|
@ -0,0 +1 @@
|
|||
GPL (v3 or later)
|
||||
|
|
@ -0,0 +1 @@
|
|||
GeometryFactory
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
# Created by the script cgal_create_CMakeLists
|
||||
# This is the CMake script for compiling a set of CGAL applications.
|
||||
|
||||
cmake_minimum_required(VERSION 3.1...3.15)
|
||||
project( Optimal_bounding_box_Tests )
|
||||
|
||||
find_package(CGAL QUIET)
|
||||
|
||||
if (NOT CGAL_FOUND)
|
||||
message(STATUS "This project requires the CGAL library, and will not be compiled.")
|
||||
return()
|
||||
endif()
|
||||
|
||||
include( ${CGAL_USE_FILE} )
|
||||
|
||||
find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater)
|
||||
if (NOT EIGEN3_FOUND)
|
||||
message(STATUS "This project requires the Eigen library, and will not be compiled.")
|
||||
return()
|
||||
else()
|
||||
include(${EIGEN3_USE_FILE})
|
||||
endif()
|
||||
|
||||
create_single_source_cgal_program("test_nelder_mead.cpp")
|
||||
create_single_source_cgal_program("test_OBB_traits.cpp")
|
||||
create_single_source_cgal_program("test_optimization_algorithms.cpp")
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,11 @@
|
|||
OFF
|
||||
4 4 0
|
||||
|
||||
-1 -0.1 0
|
||||
-1 0.1 0
|
||||
1 0 -0.1
|
||||
1 0 0.1
|
||||
3 0 1 2
|
||||
3 2 3 0
|
||||
3 1 3 2
|
||||
3 0 3 1
|
||||
|
|
@ -0,0 +1,11 @@
|
|||
OFF
|
||||
4 4 0
|
||||
|
||||
0.866802 0.740808 0.895304
|
||||
0.912651 0.761565 0.160330
|
||||
0.093661 0.892578 0.737412
|
||||
0.166461 0.149912 0.364944
|
||||
3 0 1 2
|
||||
3 2 3 0
|
||||
3 1 3 2
|
||||
3 0 3 1
|
||||
|
|
@ -0,0 +1,12 @@
|
|||
OFF
|
||||
4 4 0
|
||||
|
||||
0 1 0
|
||||
1 0 0
|
||||
0 0 0
|
||||
0 0 1
|
||||
3 0 1 2
|
||||
3 2 3 0
|
||||
3 1 3 2
|
||||
3 0 3 1
|
||||
|
||||
|
|
@ -0,0 +1,135 @@
|
|||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#include <CGAL/Optimal_bounding_box/nelder_mead_functions.h>
|
||||
#include <CGAL/Optimal_bounding_box/fitness_function.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
|
||||
bool assert_doubles(double d1, double d2, double epsilon)
|
||||
{
|
||||
return (d1 < d2 + epsilon && d1 > d2 - epsilon) ? true : false;
|
||||
}
|
||||
|
||||
void test_fitness_function()
|
||||
{
|
||||
typedef typename CGAL::Eigen_linear_algebra_traits::MatrixXd MatrixXd;
|
||||
MatrixXd data_points(4, 3);
|
||||
|
||||
data_points.set_coef(0, 0, 0.866802);
|
||||
data_points.set_coef(0, 1, 0.740808);
|
||||
data_points.set_coef(0, 2, 0.895304);
|
||||
|
||||
data_points.set_coef(1, 0, 0.912651);
|
||||
data_points.set_coef(1, 1, 0.761565);
|
||||
data_points.set_coef(1, 2, 0.160330);
|
||||
|
||||
data_points.set_coef(2, 0, 0.093661);
|
||||
data_points.set_coef(2, 1, 0.892578);
|
||||
data_points.set_coef(2, 2, 0.737412);
|
||||
|
||||
data_points.set_coef(3, 0, 0.166461);
|
||||
data_points.set_coef(3, 1, 0.149912);
|
||||
data_points.set_coef(3, 2, 0.364944);
|
||||
|
||||
typedef typename CGAL::Eigen_linear_algebra_traits::Matrix3d Matrix3d;
|
||||
Matrix3d rotation;
|
||||
rotation.set_coef(0, 0, -0.809204);
|
||||
rotation.set_coef(0, 1, 0.124296);
|
||||
rotation.set_coef(0, 2, 0.574230);
|
||||
rotation.set_coef(1, 0, -0.574694);
|
||||
rotation.set_coef(1, 1, 0.035719);
|
||||
rotation.set_coef(1, 2, -0.817589);
|
||||
rotation.set_coef(2, 0, -0.122134);
|
||||
rotation.set_coef(2, 1, -0.991602);
|
||||
rotation.set_coef(2, 2, 0.042528);
|
||||
|
||||
CGAL_assertion_code(typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits);
|
||||
CGAL_assertion_code(double fitness = CGAL::Optimal_bounding_box::
|
||||
compute_fitness<Linear_algebra_traits> (rotation, data_points));
|
||||
CGAL_assertion(assert_doubles(fitness, 0.58606, 1e-5));
|
||||
}
|
||||
|
||||
|
||||
|
||||
void test_eigen_matrix_interface()
|
||||
{
|
||||
CGAL_assertion_code(typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits);
|
||||
typedef CGAL::Eigen_dense_matrix<double, 3, 3> Matrix;
|
||||
|
||||
Matrix A(3, 3);
|
||||
A.set_coef(0, 0, 0.1);
|
||||
A.set_coef(0, 1, 0.2);
|
||||
A.set_coef(0, 2, 0.3);
|
||||
A.set_coef(1, 0, 0.4);
|
||||
A.set_coef(1, 1, 0.5);
|
||||
A.set_coef(1, 2, 0.6);
|
||||
A.set_coef(2, 0, 0.7);
|
||||
A.set_coef(2, 1, 0.8);
|
||||
A.set_coef(2, 2, 0.9);
|
||||
CGAL_assertion_code(Matrix B);
|
||||
CGAL_assertion_code(B = CGAL::Eigen_linear_algebra_traits::transpose(A));
|
||||
CGAL_assertion_code(Matrix S);
|
||||
CGAL_assertion_code(S = 0.5 * A);
|
||||
Matrix C(3,3);
|
||||
C.set_coef(0, 0, 0.3011944);
|
||||
C.set_coef(0, 1, 0.9932761);
|
||||
C.set_coef(0, 2, 0.5483701);
|
||||
C.set_coef(1, 0, 0.5149142);
|
||||
C.set_coef(1, 1, 0.5973263);
|
||||
C.set_coef(1, 2, 0.5162336);
|
||||
C.set_coef(2, 0, 0.0039213);
|
||||
C.set_coef(2, 1, 0.0202949);
|
||||
C.set_coef(2, 2, 0.9240308);
|
||||
|
||||
CGAL_assertion_code(Matrix Q = CGAL::Eigen_linear_algebra_traits::get_Q(C));
|
||||
CGAL_assertion_code(double epsilon = 1e-5);
|
||||
CGAL_assertion(assert_doubles(Q(0,0), -0.504895, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(0,1), 0.862834, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(0,2), -0.024447, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(1,0), -0.863156, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(1,1), -0.504894, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(1,2), 0.006687, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(2,0), -0.006573, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(2,1), 0.024478, epsilon));
|
||||
CGAL_assertion(assert_doubles(Q(2,2), 0.999679, epsilon));
|
||||
|
||||
Matrix D(3,3);
|
||||
D.set_coef(0, 0, -0.809204);
|
||||
D.set_coef(0, 1, 0.124296);
|
||||
D.set_coef(0, 2, 0.574230);
|
||||
D.set_coef(1, 0, -0.574694);
|
||||
D.set_coef(1, 1, 0.035719);
|
||||
D.set_coef(1, 2, -0.817589);
|
||||
D.set_coef(2, 0, -0.122134);
|
||||
D.set_coef(2, 1, -0.991602);
|
||||
D.set_coef(2, 2, 0.042528);
|
||||
|
||||
Matrix E(3,3);
|
||||
E.set_coef(0, 0, -0.45070);
|
||||
E.set_coef(0, 1, -0.32769);
|
||||
E.set_coef(0, 2, -0.83035);
|
||||
E.set_coef(1, 0, -0.13619);
|
||||
E.set_coef(1, 1, -0.89406);
|
||||
E.set_coef(1, 2, 0.42675);
|
||||
E.set_coef(2, 0, -0.88222);
|
||||
E.set_coef(2, 1, 0.30543);
|
||||
E.set_coef(2, 2, 0.35833);
|
||||
|
||||
CGAL_assertion_code(Matrix Sr = CGAL::Optimal_bounding_box::reflection<Linear_algebra_traits>(D, E));
|
||||
CGAL_assertion(assert_doubles(Sr(0,0), -0.13359, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(0,1), -0.95986, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(0,2), -0.24664, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(1,0), -0.60307, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(1,1), -0.11875, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(1,2), 0.78880, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(2,0), -0.78642, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(2,1), 0.25411, epsilon));
|
||||
CGAL_assertion(assert_doubles(Sr(2,2), -0.56300, epsilon));
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test_fitness_function();
|
||||
test_eigen_matrix_interface();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,196 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/optimal_bounding_box.h>
|
||||
|
||||
#include <iostream>
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
typedef K::Point_3 Point_3;
|
||||
|
||||
typedef CGAL::Oriented_bounding_box_traits_3<K> Traits;
|
||||
typedef Traits::Matrix Matrix;
|
||||
|
||||
bool are_equal(double d1, double d2, double epsilon)
|
||||
{
|
||||
return (d1 < d2 + epsilon && d1 > d2 - epsilon) ? true : false;
|
||||
}
|
||||
|
||||
void test_simplex_operations()
|
||||
{
|
||||
const double epsilon = 1e-5;
|
||||
|
||||
Matrix Sc(3, 3);
|
||||
Sc.set(0, 0, -0.809204); Sc.set(0, 1, 0.124296); Sc.set(0, 2, 0.574230);
|
||||
Sc.set(1, 0, -0.574694); Sc.set(1, 1, 0.035719); Sc.set(1, 2, -0.817589);
|
||||
Sc.set(2, 0, -0.122134); Sc.set(2, 1, -0.991602); Sc.set(2, 2, 0.042528);
|
||||
|
||||
Matrix S_worst(3, 3);
|
||||
S_worst.set(0, 0, -0.45070); S_worst.set(0, 1, -0.32769); S_worst.set(0, 2, -0.83035);
|
||||
S_worst.set(1, 0, -0.13619); S_worst.set(1, 1, -0.89406); S_worst.set(1, 2, 0.42675);
|
||||
S_worst.set(2, 0, -0.88222); S_worst.set(2, 1, 0.30543); S_worst.set(2, 2, 0.35833);
|
||||
|
||||
Matrix Sr = CGAL::Optimal_bounding_box::reflection<Linear_algebra_traits>(Sc, S_worst);
|
||||
assert(are_equal(Sr(0,0), -0.13359, epsilon));
|
||||
assert(are_equal(Sr(0,1), -0.95986, epsilon));
|
||||
assert(are_equal(Sr(0,2), -0.24664, epsilon));
|
||||
assert(are_equal(Sr(1,0), -0.60307, epsilon));
|
||||
assert(are_equal(Sr(1,1), -0.11875, epsilon));
|
||||
assert(are_equal(Sr(1,2), 0.78880, epsilon));
|
||||
assert(are_equal(Sr(2,0), -0.78642, epsilon));
|
||||
assert(are_equal(Sr(2,1), 0.25411, epsilon));
|
||||
assert(are_equal(Sr(2,2), -0.56300, epsilon));
|
||||
|
||||
Matrix Se = CGAL::Optimal_bounding_box::expansion<Linear_algebra_traits>(Sc, S_worst, Sr);
|
||||
assert(are_equal(Se(0,0), -0.87991, epsilon));
|
||||
assert(are_equal(Se(0,1), 0.36105, epsilon));
|
||||
assert(are_equal(Se(0,2), -0.30888, epsilon));
|
||||
assert(are_equal(Se(1,0), -0.11816, epsilon));
|
||||
assert(are_equal(Se(1,1), -0.79593, epsilon));
|
||||
assert(are_equal(Se(1,2), -0.59375, epsilon));
|
||||
assert(are_equal(Se(2,0), -0.460215, epsilon));
|
||||
assert(are_equal(Se(2,1), -0.48595, epsilon));
|
||||
assert(are_equal(Se(2,2), 0.74300, epsilon));
|
||||
|
||||
Matrix S_a(3, 3);
|
||||
S_a.set(0, 0, -0.277970); S_a.set(0, 1, 0.953559); S_a.set(0, 2, 0.116010);
|
||||
S_a.set(1, 0, -0.567497); S_a.set(1, 1, -0.065576); S_a.set(1, 2, -0.820760);
|
||||
S_a.set(2, 0, -0.775035); S_a.set(2, 1, -0.293982); S_a.set(2, 2, 0.559370);
|
||||
|
||||
Matrix S_b(3, 3);
|
||||
S_b.set(0, 0, -0.419979); S_b.set(0, 1, 0.301765); S_b.set(0, 2, -0.8558940);
|
||||
S_b.set(1, 0, -0.653011); S_b.set(1, 1, -0.755415); S_b.set(1, 2, 0.054087);
|
||||
S_b.set(2, 0, -0.630234); S_b.set(2, 1, 0.581624); S_b.set(2, 2, 0.514314);
|
||||
|
||||
Matrix S_c = CGAL::Optimal_bounding_box::mean<Linear_algebra_traits>(S_a, S_b);
|
||||
assert(are_equal(S_c(0,0), -0.35111, epsilon));
|
||||
assert(are_equal(S_c(0,1), 0.79308, epsilon));
|
||||
assert(are_equal(S_c(0,2), -0.49774, epsilon));
|
||||
assert(are_equal(S_c(1,0), -0.61398, epsilon));
|
||||
assert(are_equal(S_c(1,1), -0.59635, epsilon));
|
||||
assert(are_equal(S_c(1,2), -0.51710, epsilon));
|
||||
assert(are_equal(S_c(2,0), -0.70693, epsilon));
|
||||
assert(are_equal(S_c(2,1), 0.12405, epsilon));
|
||||
assert(are_equal(S_c(2,2), 0.69632, epsilon));
|
||||
}
|
||||
|
||||
void test_centroid()
|
||||
{
|
||||
const double epsilon = 1e-5;
|
||||
|
||||
Matrix S_a;
|
||||
S_a.set(0, 0, -0.588443); S_a.set(0, 1, 0.807140); S_a.set(0, 2, -0.047542);
|
||||
S_a.set(1, 0, -0.786228); S_a.set(1, 1, -0.584933); S_a.set(1, 2, -0.199246);
|
||||
S_a.set(2, 0, -0.188629); S_a.set(2, 1, -0.079867); S_a.set(2, 2, 0.978795);
|
||||
|
||||
Matrix S_b(3, 3);
|
||||
S_b.set(0, 0, -0.2192721); S_b.set(0, 1, 0.2792986); S_b.set(0, 2, -0.9348326);
|
||||
S_b.set(1, 0, -0.7772152); S_b.set(1, 1, -0.6292092); S_b.set(1, 2, -0.005686);
|
||||
S_b.set(2, 0, -0.5897934); S_b.set(2, 1, 0.7253193); S_b.set(2, 2, 0.3550431);
|
||||
|
||||
Matrix S_c(3, 3);
|
||||
S_c.set(0, 0, -0.32657); S_c.set(0, 1, -0.60013); S_c.set(0, 2, -0.730206);
|
||||
S_c.set(1, 0, -0.20022); S_c.set(1, 1, -0.71110); S_c.set(1, 2, 0.67398);
|
||||
S_c.set(2, 0, -0.92372); S_c.set(2, 1, 0.36630); S_c.set(2, 2, 0.11207);
|
||||
|
||||
Matrix S_centroid = CGAL::Optimal_bounding_box::nm_centroid<Linear_algebra_traits>(S_a, S_b, S_c);
|
||||
assert(are_equal(S_centroid(0,0), -0.419979, epsilon));
|
||||
assert(are_equal(S_centroid(0,1), 0.301765, epsilon));
|
||||
assert(are_equal(S_centroid(0,2), -0.855894, epsilon));
|
||||
assert(are_equal(S_centroid(1,0), -0.653011, epsilon));
|
||||
assert(are_equal(S_centroid(1,1), -0.755415, epsilon));
|
||||
assert(are_equal(S_centroid(1,2), 0.054087, epsilon));
|
||||
assert(are_equal(S_centroid(2,0), -0.630234, epsilon));
|
||||
assert(are_equal(S_centroid(2,1), 0.581624, epsilon));
|
||||
assert(are_equal(S_centroid(2,2), 0.514314, epsilon));
|
||||
}
|
||||
|
||||
void test_nelder_mead()
|
||||
{
|
||||
std::array<Point_3, 4> points;
|
||||
points[0] = Point_3(0.866802, 0.740808, 0.895304);
|
||||
points[1] = Point_3(0.912651, 0.761565, 0.160330);
|
||||
points[2] = Point_3(0.093661, 0.892578, 0.737412);
|
||||
points[3] = Point_3(0.166461, 0.149912, 0.364944);
|
||||
|
||||
// one simplex
|
||||
std::array<Matrix, 4> simplex;
|
||||
|
||||
Matrix v0, v1, v2, v3;
|
||||
v0(0,0) = -0.2192721; 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;
|
||||
v0(2,0) = -0.5897934; v0(2,1) = 0.7253193; v0(2,2) = 0.3550431;
|
||||
|
||||
v1(0,0) = -0.588443; v1(0,1) = 0.807140; v1(0,2) = -0.047542;
|
||||
v1(1,0) = -0.786228; v1(1,1) = -0.584933; v1(1,2) = -0.199246;
|
||||
v1(2,0) = -0.188629; v1(2,1) = -0.079867; v1(2,2) = 0.978795;
|
||||
|
||||
v2(0,0) = -0.277970; v2(0,1) = 0.953559; v2(0,2) = 0.116010;
|
||||
v2(1,0) = -0.567497; v2(1,1) = -0.065576; v2(1,2) = -0.820760;
|
||||
v2(2,0) = -0.775035; v2(2,1) = -0.293982; v2(2,2) = 0.559370;
|
||||
|
||||
v3(0,0) = -0.32657; v3(0,1) = -0.60013; v3(0,2) = -0.73020;
|
||||
v3(1,0) = -0.20022; v3(1,1) = -0.71110; v3(1,2) = 0.67398;
|
||||
v3(2,0) = -0.92372; v3(2,1) = 0.36630; v3(2,2) = 0.11207;
|
||||
|
||||
simplex[0] = v0;
|
||||
simplex[1] = v1;
|
||||
simplex[2] = v2;
|
||||
simplex[3] = v3;
|
||||
|
||||
std::size_t nm_iterations = 19;
|
||||
CGAL::Optimal_bounding_box::internal::nelder_mead<Linear_algebra_traits>(simplex, data_points, nm_iterations);
|
||||
|
||||
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));
|
||||
|
||||
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));
|
||||
|
||||
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));
|
||||
|
||||
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));
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test_simplex_operations();
|
||||
test_centroid();
|
||||
test_nelder_mead();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,305 @@
|
|||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#include <CGAL/Optimal_bounding_box/population.h>
|
||||
#include <CGAL/Optimal_bounding_box/optimal_bounding_box.h>
|
||||
#include <CGAL/Optimal_bounding_box/helper.h>
|
||||
#include <CGAL/Optimal_bounding_box/nelder_mead_functions.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
bool assert_doubles(double d1, double d2, double epsilon)
|
||||
{
|
||||
return (d1 < d2 + epsilon && d1 > d2 - epsilon) ? true : false;
|
||||
}
|
||||
|
||||
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(1,0) = 0.912651;
|
||||
data_points(1,1) = 0.761565;
|
||||
data_points(1,2) = 0.160330;
|
||||
data_points(2,0) = 0.093661;
|
||||
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,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();
|
||||
assert(pop.size() == 5);
|
||||
}
|
||||
|
||||
void test_random_unit_tetra()
|
||||
{
|
||||
// this is dynamic at run times
|
||||
CGAL::Eigen_dense_matrix<double, -1, -1> data_points(4, 3);
|
||||
|
||||
// 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(1,0) = 0.912651;
|
||||
data_points(1,1) = 0.761565;
|
||||
data_points(1,2) = 0.160330;
|
||||
data_points(2,0) = 0.093661;
|
||||
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,2) = 0.364944;
|
||||
|
||||
typedef CGAL::Simple_cartesian<double> K;
|
||||
typedef K::Point_3 Point;
|
||||
typedef CGAL::Surface_mesh<Point> Mesh;
|
||||
|
||||
// make a mesh and export it
|
||||
Mesh mesh;
|
||||
Point p1(0.866802, 0.740808, 0.895304);
|
||||
Point p2(0.912651, 0.761565, 0.160330);
|
||||
Point p3(0.093661, 0.892578, 0.737412);
|
||||
Point p4(0.166461, 0.149912, 0.364944);
|
||||
CGAL::make_tetrahedron(p1, p2, p3, p4, mesh);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
|
||||
std::ofstream out("data/random_unit_tetra.off");
|
||||
out << mesh;
|
||||
out.close();
|
||||
#endif
|
||||
|
||||
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);
|
||||
|
||||
Matrix R = evolution.get_best();
|
||||
const double epsilon = 1e-3;
|
||||
assert(assert_doubles(Linear_algebra_traits::compute_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);
|
||||
#endif
|
||||
}
|
||||
|
||||
void test_reference_tetrahedron(const char* fname)
|
||||
{
|
||||
std::ifstream input(fname);
|
||||
CGAL::Surface_mesh<K::Point_3> mesh;
|
||||
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::Matrix Matrix;
|
||||
|
||||
// points in a matrix
|
||||
Matrix points;
|
||||
CGAL::Optimal_bounding_box::sm_to_matrix(mesh, points);
|
||||
|
||||
std::size_t generations = 10;
|
||||
CGAL::Optimal_bounding_box::Population<Linear_algebra_traits> pop(50);
|
||||
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> experiment(pop, points);
|
||||
experiment.evolve(generations);
|
||||
|
||||
Matrix R = experiment.get_best();
|
||||
double epsilon = 1e-5;
|
||||
assert(assert_doubles(Linear_algebra_traits::compute_determinant(R), 1, epsilon));
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
|
||||
// postprocessing
|
||||
Matrix obb(8, 3);
|
||||
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
|
||||
#endif
|
||||
}
|
||||
|
||||
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())
|
||||
{
|
||||
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::Matrix Matrix;
|
||||
|
||||
// points in a matrix
|
||||
Matrix points;
|
||||
CGAL::Optimal_bounding_box::sm_to_matrix(mesh, points);
|
||||
|
||||
std::size_t max_generations = 10;
|
||||
CGAL::Optimal_bounding_box::Population<Linear_algebra_traits> pop(50);
|
||||
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> experiment(pop, points);
|
||||
experiment.evolve(max_generations);
|
||||
|
||||
Matrix R = experiment.get_best();
|
||||
double epsilon = 1e-3;
|
||||
assert(assert_doubles(Linear_algebra_traits::compute_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
|
||||
// postprocessing
|
||||
Matrix obb(8, 3);
|
||||
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
|
||||
#endif
|
||||
}
|
||||
|
||||
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())
|
||||
{
|
||||
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;
|
||||
|
||||
PointPMap pmap = get(boost::vertex_point, mesh);
|
||||
for(vertex_descriptor v : vertices(mesh))
|
||||
sm_points.push_back(get(pmap, v));
|
||||
|
||||
|
||||
CGAL::Eigen_linear_algebra_traits la_traits;
|
||||
std::vector<K::Point_3> obb_points;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points, la_traits, true);
|
||||
|
||||
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
|
||||
/*
|
||||
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);
|
||||
|
||||
std::ofstream out("data/obb_result.off");
|
||||
out << result_mesh;
|
||||
out.close();
|
||||
#endif
|
||||
}
|
||||
|
||||
void test_compute_obb_mesh(const std::string fname)
|
||||
{
|
||||
std::ifstream input(fname);
|
||||
CGAL::Surface_mesh<K::Point_3> mesh;
|
||||
if (!input || !(input >> mesh) || mesh.is_empty())
|
||||
{
|
||||
std::cerr << fname << " is not a valid off file.\n";
|
||||
std::exit(1);
|
||||
}
|
||||
|
||||
CGAL::Eigen_linear_algebra_traits la_traits;
|
||||
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
|
||||
std::ofstream out("/tmp/result_elephant.off");
|
||||
out << obbmesh;
|
||||
out.close();
|
||||
#endif
|
||||
}
|
||||
|
||||
void test_function_defaults_traits(const std::string fname1,
|
||||
const std::string fname2)
|
||||
{
|
||||
std::ifstream input1(fname1);
|
||||
CGAL::Surface_mesh<K::Point_3> mesh1;
|
||||
if (!input1 || !(input1 >> mesh1) || mesh1.is_empty())
|
||||
{
|
||||
std::cerr << fname1 << " is not a valid off file.\n";
|
||||
std::exit(1);
|
||||
}
|
||||
|
||||
std::ifstream input2(fname2);
|
||||
CGAL::Surface_mesh<K::Point_3> mesh2;
|
||||
if (!input2 || !(input2 >> mesh2) || mesh2.is_empty())
|
||||
{
|
||||
std::cerr << fname2 << " is not a valid off file.\n";
|
||||
std::exit(1);
|
||||
}
|
||||
|
||||
// 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);
|
||||
|
||||
for(vertex_descriptor v : vertices(mesh1))
|
||||
sm_points.push_back(get(pmap, v));
|
||||
|
||||
std::vector<K::Point_3> obb_points;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(sm_points, obb_points, true);
|
||||
|
||||
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::Optimal_bounding_box::compute_optimal_bounding_box(mesh2, obbmesh, true);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test_genetic_algorithm();
|
||||
|
||||
test_random_unit_tetra();
|
||||
test_reference_tetrahedron("data/reference_tetrahedron.off");
|
||||
test_long_tetrahedron("data/long_tetrahedron.off");
|
||||
test_compute_obb_evolution("data/random_unit_tetra.off");
|
||||
test_compute_obb_mesh("data/elephant.off");
|
||||
test_function_defaults_traits("data/random_unit_tetra.off", "data/elephant.off");
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -17,6 +17,9 @@ target_link_libraries(clipping_box_plugin PUBLIC scene_edit_box_item scene_basi
|
|||
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_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)
|
||||
target_link_libraries(basic_generator_plugin PUBLIC scene_surface_mesh_item scene_points_with_normal_item scene_polylines_item)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,154 @@
|
|||
#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 "Scene_polyhedron_selection_item.h"
|
||||
#include "Scene_points_with_normal_item.h"
|
||||
|
||||
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
|
||||
|
||||
#include <CGAL/optimal_bounding_box.h>
|
||||
|
||||
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
|
||||
{
|
||||
Q_OBJECT
|
||||
Q_INTERFACES(CGAL::Three::Polyhedron_demo_plugin_interface)
|
||||
Q_PLUGIN_METADATA(IID "com.geometryfactory.PolyhedronDemo.PluginInterface/1.0")
|
||||
|
||||
public:
|
||||
void init(QMainWindow* mainWindow, Scene_interface* scene_interface, Messages_interface*);
|
||||
QList<QAction*> actions() const;
|
||||
|
||||
bool applicable(QAction*) const
|
||||
{
|
||||
if(scene->mainSelectionIndex() != -1
|
||||
&& scene->item(scene->mainSelectionIndex())->isFinite())
|
||||
return true;
|
||||
return false;
|
||||
}
|
||||
|
||||
protected:
|
||||
void gather_mesh_points(std::vector<Point_3>& points);
|
||||
void obb();
|
||||
|
||||
public Q_SLOTS:
|
||||
void createObb()
|
||||
{
|
||||
QApplication::setOverrideCursor(Qt::WaitCursor);
|
||||
obb();
|
||||
QApplication::restoreOverrideCursor();
|
||||
}
|
||||
|
||||
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;
|
||||
mw = mainWindow;
|
||||
actionObb = new QAction(tr("Create &Optimal Bbox Mesh"), mainWindow);
|
||||
actionObb->setObjectName("createObbMeshAction");
|
||||
connect(actionObb, SIGNAL(triggered()), this, SLOT(createObb()));
|
||||
}
|
||||
|
||||
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* item = qobject_cast<Scene_facegraph_item*>(scene->item(index));
|
||||
|
||||
Scene_polyhedron_selection_item* selection_item =
|
||||
qobject_cast<Scene_polyhedron_selection_item*>(scene->item(index));
|
||||
|
||||
Scene_points_with_normal_item* point_set_item =
|
||||
qobject_cast<Scene_points_with_normal_item*>(scene->item(index));
|
||||
|
||||
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;
|
||||
typedef typename boost::graph_traits<FaceGraph>::face_descriptor face_descriptor;
|
||||
|
||||
std::vector<vertex_descriptor> selected_vertices;
|
||||
|
||||
if(item != NULL)
|
||||
{
|
||||
FaceGraph& pmesh = *item->polyhedron();
|
||||
selected_vertices.assign(vertices(pmesh).begin(), vertices(pmesh).end());
|
||||
PointPMap pmap = get(CGAL::vertex_point, pmesh);
|
||||
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();
|
||||
for(face_descriptor f : selection_item->selected_facets)
|
||||
{
|
||||
for(vertex_descriptor v : vertices_around_face(halfedge(f, pmesh), pmesh))
|
||||
selected_vertices.push_back(v);
|
||||
}
|
||||
|
||||
PointPMap pmap = get(CGAL::vertex_point, pmesh);
|
||||
for(vertex_descriptor v : selected_vertices)
|
||||
points.push_back(get(pmap, v));
|
||||
}
|
||||
|
||||
CGAL_assertion(points.size() >= 3);
|
||||
}
|
||||
|
||||
if(point_set_item)
|
||||
{
|
||||
Point_set* points_set = point_set_item->point_set();
|
||||
if(points_set == NULL)
|
||||
return;
|
||||
|
||||
std::cout << "points_set->size()= " << points_set->size() << std::endl;
|
||||
for(const Point_3& p : points_set->points())
|
||||
points.push_back(p);
|
||||
}
|
||||
}
|
||||
|
||||
void Create_obb_mesh_plugin::obb()
|
||||
{
|
||||
// gather point coordinates
|
||||
std::vector<Point_3> points;
|
||||
gather_mesh_points(points);
|
||||
|
||||
// find obb
|
||||
std::array<Point_3, 8> obb_points;
|
||||
CGAL::oriented_bounding_box(points, obb_points);
|
||||
|
||||
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);
|
||||
scene->addItem(item);
|
||||
}
|
||||
|
||||
#include "Create_obb_mesh_plugin.moc"
|
||||
|
|
@ -86,6 +86,9 @@ public:
|
|||
\cgalConcept
|
||||
Concept of matrix type used by the concept `SvdTraits`.
|
||||
|
||||
\cgalRefines `DefaultConstructible`
|
||||
\cgalRefines `Assignable`
|
||||
|
||||
\cgalHasModel `CGAL::Eigen_matrix<T>`
|
||||
*/
|
||||
class SvdTraits::Matrix
|
||||
|
|
|
|||
|
|
@ -8,377 +8,84 @@
|
|||
//
|
||||
// Author(s) : Gael Guennebaud
|
||||
|
||||
#ifndef CGAL_EIGEN_MATRIX_H
|
||||
#define CGAL_EIGEN_MATRIX_H
|
||||
#ifndef CGAL_SOLVER_INTERFACE_EIGEN_MATRIX_H
|
||||
#define CGAL_SOLVER_INTERFACE_EIGEN_MATRIX_H
|
||||
|
||||
#include <CGAL/basic.h> // include basic.h before testing #defines
|
||||
|
||||
#include <Eigen/Sparse>
|
||||
#include <CGAL/Eigen_sparse_matrix.h> // for backward compatibility
|
||||
|
||||
#include <Eigen/Dense>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolverInterfaceRef
|
||||
|
||||
The class `Eigen_sparse_matrix` is a wrapper around `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.
|
||||
|
||||
\cgalModels `SparseLinearAlgebraTraits_d::Matrix`
|
||||
|
||||
\tparam T Number type.
|
||||
|
||||
\sa `CGAL::Eigen_vector<T>`
|
||||
\sa `CGAL::Eigen_matrix<T>`
|
||||
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
|
||||
*/
|
||||
template<class T>
|
||||
struct Eigen_sparse_matrix
|
||||
{
|
||||
// Public types
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/// The internal matrix type from \ref thirdpartyEigen "Eigen".
|
||||
typedef Eigen::SparseMatrix<T> EigenType;
|
||||
|
||||
typedef T NT;
|
||||
/// @}
|
||||
|
||||
Eigen_sparse_matrix(const EigenType& et)
|
||||
: m_is_already_built(true), m_matrix(et), m_is_symmetric(false)
|
||||
{}
|
||||
|
||||
// Public operations
|
||||
public:
|
||||
Eigen_sparse_matrix() :
|
||||
m_is_already_built(false)
|
||||
{}
|
||||
|
||||
/// Create a square matrix initialized with zeros.
|
||||
Eigen_sparse_matrix(std::size_t dim, ///< Matrix dimension.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
:
|
||||
m_is_already_built(false),
|
||||
m_matrix(static_cast<int>(dim), static_cast<int>(dim))
|
||||
{
|
||||
CGAL_precondition(dim > 0);
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
m_triplets.reserve(dim); // reserve memory for a regular 3D grid
|
||||
}
|
||||
|
||||
/// Create a square matrix initialized with zeros.
|
||||
Eigen_sparse_matrix(int dim, ///< Matrix dimension.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_is_already_built(false),
|
||||
m_matrix(dim, dim)
|
||||
{
|
||||
CGAL_precondition(dim > 0);
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
// reserve memory for a regular 3D grid
|
||||
m_triplets.reserve(dim);
|
||||
}
|
||||
|
||||
/// Create a rectangular matrix initialized with zeros.
|
||||
///
|
||||
/// \pre rows == columns if `is_symmetric` is true.
|
||||
Eigen_sparse_matrix(std::size_t rows, ///< Number of rows.
|
||||
std::size_t columns, ///< Number of columns.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_is_already_built(false),
|
||||
m_matrix(static_cast<int>(rows), static_cast<int>(columns))
|
||||
{
|
||||
CGAL_precondition(rows > 0);
|
||||
CGAL_precondition(columns > 0);
|
||||
if(m_is_symmetric)
|
||||
{
|
||||
CGAL_precondition(rows == columns);
|
||||
}
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
// reserve memory for a regular 3D grid
|
||||
m_triplets.reserve(rows);
|
||||
}
|
||||
|
||||
void swap(Eigen_sparse_matrix& other)
|
||||
{
|
||||
std::swap(m_is_already_built, other.m_is_already_built);
|
||||
std::swap(m_is_symmetric, other.m_is_symmetric);
|
||||
m_matrix.swap(other.m_matrix);
|
||||
m_triplets.swap(other.m_triplets);
|
||||
}
|
||||
|
||||
|
||||
/// Delete this object and the wrapped matrix.
|
||||
~Eigen_sparse_matrix() { }
|
||||
|
||||
/// Create a rectangular matrix initialized with zeros.
|
||||
///
|
||||
/// \pre rows == columns if `is_symmetric` is true.
|
||||
Eigen_sparse_matrix(int rows, ///< Number of rows.
|
||||
int columns, ///< Number of columns.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_is_already_built(false),
|
||||
m_matrix(rows,columns)
|
||||
{
|
||||
CGAL_precondition(rows > 0);
|
||||
CGAL_precondition(columns > 0);
|
||||
if(is_symmetric)
|
||||
{
|
||||
CGAL_precondition(rows == columns);
|
||||
}
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
// reserve memory for a regular 3D grid
|
||||
m_triplets.reserve(rows);
|
||||
}
|
||||
|
||||
/// Return the matrix number of rows
|
||||
int row_dimension() const { return static_cast<int>(m_matrix.rows()); }
|
||||
/// Return the matrix number of columns
|
||||
int column_dimension() const { return static_cast<int>(m_matrix.cols()); }
|
||||
|
||||
/// Write access to a matrix coefficient: a_ij <- val.
|
||||
///
|
||||
/// Users can optimize calls to this function by setting 'new_coef' to `true`
|
||||
/// if the coefficient does not already exist in the matrix.
|
||||
///
|
||||
/// \warning For symmetric matrices, `Eigen_sparse_matrix` only stores the lower triangle
|
||||
/// and `set_coef()` does nothing if (i, j) belongs to the upper triangle.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
void set_coef(std::size_t i_, std::size_t j_, T val, bool new_coef = false)
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
CGAL_precondition(i < row_dimension());
|
||||
CGAL_precondition(j < column_dimension());
|
||||
|
||||
if (m_is_symmetric && (j > i))
|
||||
return;
|
||||
|
||||
if(m_is_already_built)
|
||||
{
|
||||
m_matrix.coeffRef(i,j) = val;
|
||||
}
|
||||
else
|
||||
{
|
||||
if(new_coef == false)
|
||||
{
|
||||
assemble_matrix();
|
||||
m_matrix.coeffRef(i,j) = val;
|
||||
}
|
||||
else
|
||||
{
|
||||
m_triplets.push_back(Triplet(i,j,val));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Write access to a matrix coefficient: a_ij <- a_ij + val.
|
||||
///
|
||||
/// \warning For symmetric matrices, Eigen_sparse_matrix only stores the lower triangle
|
||||
/// `add_coef()` does nothing if (i, j) belongs to the upper triangle.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
void add_coef(std::size_t i_, std::size_t j_, T val)
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
if(m_is_symmetric && (j > i))
|
||||
return;
|
||||
|
||||
if(m_is_already_built){
|
||||
CGAL_precondition(i < row_dimension());
|
||||
CGAL_precondition(j < column_dimension());
|
||||
m_matrix.coeffRef(i,j) += val;
|
||||
}else{
|
||||
m_triplets.push_back(Triplet(i,j,val));
|
||||
}
|
||||
}
|
||||
|
||||
/// Read access to a matrix coefficient.
|
||||
///
|
||||
/// \warning Complexity:
|
||||
/// - O(log(n)) if the matrix is already built.
|
||||
/// - O(n) if the matrix is not built.
|
||||
/// `n` being the number of entries in the matrix.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
NT get_coef (std::size_t i_, std::size_t j_) const
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
CGAL_precondition(i < row_dimension());
|
||||
CGAL_precondition(j < column_dimension());
|
||||
|
||||
if(m_is_symmetric && j > i)
|
||||
std::swap(i, j);
|
||||
|
||||
if (m_is_already_built)
|
||||
return m_matrix.coeffRef(i,j);
|
||||
else
|
||||
{
|
||||
NT val = 0;
|
||||
for(std::size_t t=0; t<m_triplets.size(); ++t)
|
||||
{
|
||||
if(m_triplets[t].col() == j &&
|
||||
m_triplets[t].row() == i)
|
||||
val += m_triplets[t].value();
|
||||
}
|
||||
return val;
|
||||
}
|
||||
}
|
||||
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
void assemble_matrix() const
|
||||
{
|
||||
m_matrix.setFromTriplets(m_triplets.begin(), m_triplets.end());
|
||||
m_is_already_built = true;
|
||||
m_triplets.clear(); // the matrix is built and will not be rebuilt
|
||||
}
|
||||
/// \endcond
|
||||
|
||||
/// Return the internal matrix, with type `EigenType`.
|
||||
const EigenType& eigen_object() const
|
||||
{
|
||||
if(!m_is_already_built)
|
||||
assemble_matrix();
|
||||
|
||||
// turns the matrix into compressed mode:
|
||||
// -> release some memory
|
||||
// -> required for some external solvers
|
||||
m_matrix.makeCompressed();
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
/// Return the internal matrix, with type `EigenType`.
|
||||
EigenType& eigen_object()
|
||||
{
|
||||
if(!m_is_already_built)
|
||||
assemble_matrix();
|
||||
|
||||
// turns the matrix into compressed mode:
|
||||
// -> release some memory
|
||||
// -> required for some external solvers
|
||||
m_matrix.makeCompressed();
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
|
||||
public:
|
||||
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
friend Eigen_sparse_matrix
|
||||
operator*(const T& c, const Eigen_sparse_matrix& M)
|
||||
{
|
||||
return Eigen_sparse_matrix(c* M.eigen_object());
|
||||
}
|
||||
|
||||
|
||||
friend Eigen_sparse_matrix
|
||||
operator+(const Eigen_sparse_matrix& M0, const Eigen_sparse_matrix& M1)
|
||||
{
|
||||
return Eigen_sparse_matrix(M0.eigen_object()+ M1.eigen_object());
|
||||
}
|
||||
/// \endcond
|
||||
|
||||
|
||||
// Fields
|
||||
private:
|
||||
mutable bool m_is_already_built;
|
||||
|
||||
typedef Eigen::Triplet<T,int> Triplet;
|
||||
mutable std::vector<Triplet> m_triplets;
|
||||
|
||||
mutable EigenType m_matrix;
|
||||
|
||||
// Symmetric/hermitian?
|
||||
bool m_is_symmetric;
|
||||
}; // Eigen_sparse_matrix
|
||||
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolverInterfaceRef
|
||||
|
||||
The class `Eigen_sparse_symmetric_matrix` is a wrapper around `Eigen` matrix type
|
||||
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::SparseMatrix` </a>
|
||||
|
||||
Since the matrix is symmetric, only the lower triangle part is stored.
|
||||
|
||||
\cgalModels `SparseLinearAlgebraTraits_d::Matrix`
|
||||
|
||||
\tparam T Number type.
|
||||
|
||||
\sa `CGAL::Eigen_vector<T>`
|
||||
\sa `CGAL::Eigen_sparse_matrix<T>`
|
||||
*/
|
||||
|
||||
template<class T>
|
||||
struct Eigen_sparse_symmetric_matrix
|
||||
: public Eigen_sparse_matrix<T>
|
||||
{
|
||||
/// Create a square *symmetric* matrix initialized with zeros.
|
||||
Eigen_sparse_symmetric_matrix(int dim) ///< Matrix dimension.
|
||||
: Eigen_sparse_matrix<T>(dim, true /* symmetric */)
|
||||
{
|
||||
}
|
||||
|
||||
/// Create a square *symmetric* matrix initialized with zeros.
|
||||
///
|
||||
/// \pre rows == columns.
|
||||
Eigen_sparse_symmetric_matrix(int rows, ///< Number of rows.
|
||||
int columns) ///< Number of columns.
|
||||
: Eigen_sparse_matrix<T>(rows, columns, true /* symmetric */)
|
||||
{
|
||||
}
|
||||
};
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolverInterfaceRef
|
||||
|
||||
The class `Eigen_matrix` is a wrapper around `Eigen` matrix type
|
||||
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html">`Eigen::Matrix`</a>.
|
||||
|
||||
\cgalModels `SvdTraits::Matrix`
|
||||
|
||||
\tparam T Number type.
|
||||
\tparam D1 Number of rows, or Dynamic
|
||||
\tparam D2 Number of columns, or Dynamic
|
||||
|
||||
\sa `CGAL::Eigen_vector<T>`
|
||||
\sa `CGAL::Eigen_sparse_matrix<T>`
|
||||
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
|
||||
*/
|
||||
template <class FT>
|
||||
template <class FT, int D1 = ::Eigen::Dynamic, int D2 = ::Eigen::Dynamic>
|
||||
struct Eigen_matrix
|
||||
: public ::Eigen::Matrix<FT, ::Eigen::Dynamic, ::Eigen::Dynamic>
|
||||
: public ::Eigen::Matrix<FT, D1, D2>
|
||||
{
|
||||
/// The internal matrix type from \ref thirdpartyEigen "Eigen".
|
||||
typedef ::Eigen::Matrix<FT, ::Eigen::Dynamic, ::Eigen::Dynamic> EigenType;
|
||||
typedef ::Eigen::Matrix<FT, D1, D2> EigenType;
|
||||
|
||||
/// Construct a matrix with `nr` rows and `nc` columns.
|
||||
/// Constructs a null matrix.
|
||||
Eigen_matrix() { }
|
||||
|
||||
/// Constructs an uninitialized matrix with `nr` rows and `nc` columns.
|
||||
/// This is useful for dynamic-size matrices.
|
||||
/// For fixed-size matrices, it is redundant to pass these parameters.
|
||||
Eigen_matrix(std::size_t nr, std::size_t nc) : EigenType(nr, nc) { }
|
||||
|
||||
/// Return the matrix number of rows.
|
||||
/// Constructs a matrix from an Eigen matrix.
|
||||
Eigen_matrix(const EigenType& b) : EigenType(b) { }
|
||||
|
||||
/// Returns the matrix number of rows.
|
||||
std::size_t number_of_rows() const { return this->rows(); }
|
||||
/// Return the matrix number of columns.
|
||||
/// Returns the matrix number of columns.
|
||||
std::size_t number_of_columns() const { return this->cols(); }
|
||||
|
||||
/// Return the value of the matrix at position (i,j).
|
||||
/// Returns the value of the matrix at position (i,j).
|
||||
FT operator()( std::size_t i , std::size_t j ) const { return EigenType::operator()(i,j); }
|
||||
|
||||
/// Write access to a matrix coefficient: `a_ij` <- `val`.
|
||||
/// Writes access to a matrix coefficient: `a_ij` <- `val`.
|
||||
void set(std::size_t i, std::size_t j, FT value) { this->coeffRef(i,j) = value; }
|
||||
|
||||
/// Return the internal matrix, with type `EigenType`.
|
||||
/// Returns the internal matrix, with type `EigenType`.
|
||||
const EigenType& eigen_object() const { return static_cast<const EigenType&>(*this); }
|
||||
|
||||
public:
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
friend Eigen_matrix operator*(const FT c, const Eigen_matrix& M)
|
||||
{
|
||||
return Eigen_matrix(c * M.eigen_object());
|
||||
}
|
||||
|
||||
friend Eigen_matrix operator*(const Eigen_matrix& M0, const Eigen_matrix& M1)
|
||||
{
|
||||
return Eigen_matrix(M0.eigen_object() * M1.eigen_object());
|
||||
}
|
||||
|
||||
friend Eigen_matrix operator+(const Eigen_matrix& M0, const Eigen_matrix& M1)
|
||||
{
|
||||
return Eigen_matrix(M0.eigen_object() + M1.eigen_object());
|
||||
}
|
||||
/// \endcond
|
||||
};
|
||||
|
||||
} //namespace CGAL
|
||||
|
||||
#endif // CGAL_EIGEN_MATRIX_H
|
||||
#endif // CGAL_SOLVER_INTERFACE_EIGEN_MATRIX_H
|
||||
|
|
|
|||
|
|
@ -0,0 +1,335 @@
|
|||
// Copyright (c) 2012 INRIA Bordeaux Sud-Ouest (France), All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org)
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
|
||||
//
|
||||
// Author(s) : Gael Guennebaud
|
||||
|
||||
#ifndef CGAL_SOLVER_INTERFACE_EIGEN_SPARSE_MATRIX_H
|
||||
#define CGAL_SOLVER_INTERFACE_EIGEN_SPARSE_MATRIX_H
|
||||
|
||||
#include <CGAL/basic.h> // include basic.h before testing #defines
|
||||
|
||||
#include <Eigen/Sparse>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolverInterfaceRef
|
||||
|
||||
The class `Eigen_sparse_matrix` is a wrapper around `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.
|
||||
|
||||
\cgalModels `SparseLinearAlgebraTraits_d::Matrix`
|
||||
|
||||
\tparam T Number type.
|
||||
|
||||
\sa `CGAL::Eigen_vector<T>`
|
||||
\sa `CGAL::Eigen_matrix<T>`
|
||||
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
|
||||
*/
|
||||
template<class T>
|
||||
struct Eigen_sparse_matrix
|
||||
{
|
||||
// Public types
|
||||
public:
|
||||
/// \name Types
|
||||
/// @{
|
||||
|
||||
/// The internal matrix type from \ref thirdpartyEigen "Eigen".
|
||||
typedef Eigen::SparseMatrix<T> EigenType;
|
||||
|
||||
typedef T NT;
|
||||
/// @}
|
||||
|
||||
Eigen_sparse_matrix(const EigenType& et)
|
||||
: m_is_already_built(true), m_matrix(et), m_is_symmetric(false)
|
||||
{}
|
||||
|
||||
// Public operations
|
||||
public:
|
||||
Eigen_sparse_matrix() : m_is_already_built(false) { }
|
||||
|
||||
/// Create a square matrix initialized with zeros.
|
||||
Eigen_sparse_matrix(std::size_t dim, ///< Matrix dimension.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
:
|
||||
m_is_already_built(false),
|
||||
m_matrix(static_cast<int>(dim), static_cast<int>(dim))
|
||||
{
|
||||
CGAL_precondition(dim > 0);
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
m_triplets.reserve(dim); // reserve memory for a regular 3D grid
|
||||
}
|
||||
|
||||
/// Create a square matrix initialized with zeros.
|
||||
Eigen_sparse_matrix(int dim, ///< Matrix dimension.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_is_already_built(false),
|
||||
m_matrix(dim, dim)
|
||||
{
|
||||
CGAL_precondition(dim > 0);
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
// reserve memory for a regular 3D grid
|
||||
m_triplets.reserve(dim);
|
||||
}
|
||||
|
||||
/// Create a rectangular matrix initialized with zeros.
|
||||
///
|
||||
/// \pre rows == columns if `is_symmetric` is true.
|
||||
Eigen_sparse_matrix(std::size_t rows, ///< Number of rows.
|
||||
std::size_t columns, ///< Number of columns.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_is_already_built(false),
|
||||
m_matrix(static_cast<int>(rows), static_cast<int>(columns))
|
||||
{
|
||||
CGAL_precondition(rows > 0);
|
||||
CGAL_precondition(columns > 0);
|
||||
if(m_is_symmetric)
|
||||
{
|
||||
CGAL_precondition(rows == columns);
|
||||
}
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
// reserve memory for a regular 3D grid
|
||||
m_triplets.reserve(rows);
|
||||
}
|
||||
|
||||
void swap(Eigen_sparse_matrix& other)
|
||||
{
|
||||
std::swap(m_is_already_built, other.m_is_already_built);
|
||||
std::swap(m_is_symmetric, other.m_is_symmetric);
|
||||
m_matrix.swap(other.m_matrix);
|
||||
m_triplets.swap(other.m_triplets);
|
||||
}
|
||||
|
||||
/// Delete this object and the wrapped matrix.
|
||||
~Eigen_sparse_matrix() { }
|
||||
|
||||
/// Create a rectangular matrix initialized with zeros.
|
||||
///
|
||||
/// \pre rows == columns if `is_symmetric` is true.
|
||||
Eigen_sparse_matrix(int rows, ///< Number of rows.
|
||||
int columns, ///< Number of columns.
|
||||
bool is_symmetric = false) ///< Symmetric/hermitian?
|
||||
: m_is_already_built(false),
|
||||
m_matrix(rows,columns)
|
||||
{
|
||||
CGAL_precondition(rows > 0);
|
||||
CGAL_precondition(columns > 0);
|
||||
if(is_symmetric)
|
||||
{
|
||||
CGAL_precondition(rows == columns);
|
||||
}
|
||||
|
||||
m_is_symmetric = is_symmetric;
|
||||
// reserve memory for a regular 3D grid
|
||||
m_triplets.reserve(rows);
|
||||
}
|
||||
|
||||
/// Return the matrix number of rows
|
||||
int row_dimension() const { return static_cast<int>(m_matrix.rows()); }
|
||||
/// Return the matrix number of columns
|
||||
int column_dimension() const { return static_cast<int>(m_matrix.cols()); }
|
||||
|
||||
/// Write access to a matrix coefficient: a_ij <- val.
|
||||
///
|
||||
/// Users can optimize calls to this function by setting 'new_coef' to `true`
|
||||
/// if the coefficient does not already exist in the matrix.
|
||||
///
|
||||
/// \warning For symmetric matrices, `Eigen_sparse_matrix` only stores the lower triangle
|
||||
/// and `set_coef()` does nothing if (i, j) belongs to the upper triangle.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
void set_coef(std::size_t i_, std::size_t j_, T val, bool new_coef = false)
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
CGAL_precondition(i < row_dimension());
|
||||
CGAL_precondition(j < column_dimension());
|
||||
|
||||
if (m_is_symmetric && (j > i))
|
||||
return;
|
||||
|
||||
if(m_is_already_built)
|
||||
{
|
||||
m_matrix.coeffRef(i,j) = val;
|
||||
}
|
||||
else
|
||||
{
|
||||
if(new_coef == false)
|
||||
{
|
||||
assemble_matrix();
|
||||
m_matrix.coeffRef(i,j) = val;
|
||||
}
|
||||
else
|
||||
{
|
||||
m_triplets.push_back(Triplet(i,j,val));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Write access to a matrix coefficient: a_ij <- a_ij + val.
|
||||
///
|
||||
/// \warning For symmetric matrices, Eigen_sparse_matrix only stores the lower triangle
|
||||
/// `add_coef()` does nothing if (i, j) belongs to the upper triangle.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
void add_coef(std::size_t i_, std::size_t j_, T val)
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
if(m_is_symmetric && (j > i))
|
||||
return;
|
||||
|
||||
if(m_is_already_built){
|
||||
CGAL_precondition(i < row_dimension());
|
||||
CGAL_precondition(j < column_dimension());
|
||||
m_matrix.coeffRef(i,j) += val;
|
||||
}else{
|
||||
m_triplets.push_back(Triplet(i,j,val));
|
||||
}
|
||||
}
|
||||
|
||||
/// Read access to a matrix coefficient.
|
||||
///
|
||||
/// \warning Complexity:
|
||||
/// - O(log(n)) if the matrix is already built.
|
||||
/// - O(n) if the matrix is not built.
|
||||
/// `n` being the number of entries in the matrix.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
NT get_coef (std::size_t i_, std::size_t j_) const
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
CGAL_precondition(i < row_dimension());
|
||||
CGAL_precondition(j < column_dimension());
|
||||
|
||||
if(m_is_symmetric && j > i)
|
||||
std::swap(i, j);
|
||||
|
||||
if (m_is_already_built)
|
||||
return m_matrix.coeffRef(i,j);
|
||||
else
|
||||
{
|
||||
NT val = 0;
|
||||
for(std::size_t t=0; t<m_triplets.size(); ++t)
|
||||
{
|
||||
if(m_triplets[t].col() == j &&
|
||||
m_triplets[t].row() == i)
|
||||
val += m_triplets[t].value();
|
||||
}
|
||||
return val;
|
||||
}
|
||||
}
|
||||
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
void assemble_matrix() const
|
||||
{
|
||||
m_matrix.setFromTriplets(m_triplets.begin(), m_triplets.end());
|
||||
m_is_already_built = true;
|
||||
m_triplets.clear(); // the matrix is built and will not be rebuilt
|
||||
}
|
||||
/// \endcond
|
||||
|
||||
/// Return the internal matrix, with type `EigenType`.
|
||||
const EigenType& eigen_object() const
|
||||
{
|
||||
if(!m_is_already_built)
|
||||
assemble_matrix();
|
||||
|
||||
// turns the matrix into compressed mode:
|
||||
// -> release some memory
|
||||
// -> required for some external solvers
|
||||
m_matrix.makeCompressed();
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
/// Return the internal matrix, with type `EigenType`.
|
||||
EigenType& eigen_object()
|
||||
{
|
||||
if(!m_is_already_built)
|
||||
assemble_matrix();
|
||||
|
||||
// turns the matrix into compressed mode:
|
||||
// -> release some memory
|
||||
// -> required for some external solvers
|
||||
m_matrix.makeCompressed();
|
||||
return m_matrix;
|
||||
}
|
||||
|
||||
public:
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
friend Eigen_sparse_matrix
|
||||
operator*(const T& c, const Eigen_sparse_matrix& M)
|
||||
{
|
||||
return Eigen_sparse_matrix(c* M.eigen_object());
|
||||
}
|
||||
|
||||
friend Eigen_sparse_matrix
|
||||
operator+(const Eigen_sparse_matrix& M0, const Eigen_sparse_matrix& M1)
|
||||
{
|
||||
return Eigen_sparse_matrix(M0.eigen_object()+ M1.eigen_object());
|
||||
}
|
||||
/// \endcond
|
||||
|
||||
// Fields
|
||||
private:
|
||||
mutable bool m_is_already_built;
|
||||
|
||||
typedef Eigen::Triplet<T,int> Triplet;
|
||||
mutable std::vector<Triplet> m_triplets;
|
||||
|
||||
mutable EigenType m_matrix;
|
||||
|
||||
// Symmetric/hermitian?
|
||||
bool m_is_symmetric;
|
||||
}; // Eigen_sparse_matrix
|
||||
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolverInterfaceRef
|
||||
|
||||
The class `Eigen_sparse_symmetric_matrix` is a wrapper around `Eigen` matrix type
|
||||
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::SparseMatrix` </a>
|
||||
|
||||
Since the matrix is symmetric, only the lower triangle part is stored.
|
||||
|
||||
\cgalModels `SparseLinearAlgebraTraits_d::Matrix`
|
||||
|
||||
\tparam T Number type.
|
||||
|
||||
\sa `CGAL::Eigen_vector<T>`
|
||||
\sa `CGAL::Eigen_sparse_matrix<T>`
|
||||
*/
|
||||
template<class T>
|
||||
struct Eigen_sparse_symmetric_matrix
|
||||
: public Eigen_sparse_matrix<T>
|
||||
{
|
||||
/// Create a square *symmetric* matrix initialized with zeros.
|
||||
Eigen_sparse_symmetric_matrix(int dim) ///< Matrix dimension.
|
||||
: Eigen_sparse_matrix<T>(dim, true /* symmetric */)
|
||||
{ }
|
||||
|
||||
/// Create a square *symmetric* matrix initialized with zeros.
|
||||
///
|
||||
/// \pre rows == columns.
|
||||
Eigen_sparse_symmetric_matrix(int rows, ///< Number of rows.
|
||||
int columns) ///< Number of columns.
|
||||
: Eigen_sparse_matrix<T>(rows, columns, true /* symmetric */)
|
||||
{ }
|
||||
};
|
||||
|
||||
} //namespace CGAL
|
||||
|
||||
#endif // CGAL_SOLVER_INTERFACE_EIGEN_SPARSE_MATRIX_H
|
||||
Loading…
Reference in New Issue