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
edefc616f3
|
|
@ -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,115 @@
|
|||
#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);
|
||||
BOOST_FOREACH(vertex_descriptor v, vertices(mesh))
|
||||
points.push_back(get(pmap, v));
|
||||
}
|
||||
|
||||
template <typename Point>
|
||||
double calculate_volume(std::vector<Point> points)
|
||||
{
|
||||
CGAL::Bbox_3 bbox;
|
||||
bbox = bbox_3(points.begin(), points.end());
|
||||
K::Iso_cuboid_3 ic(bbox);
|
||||
return ic.volume();
|
||||
}
|
||||
|
||||
void bench_finding_obb(std::string fname)
|
||||
{
|
||||
std::ifstream input(fname);
|
||||
|
||||
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,6 @@
|
|||
@INCLUDE = ${CGAL_DOC_PACKAGE_DEFAULTS}
|
||||
|
||||
PROJECT_NAME = "CGAL ${CGAL_DOC_VERSION} - Optimal Bounding Box"
|
||||
EXTRACT_ALL = false
|
||||
HIDE_UNDOC_CLASSES = true
|
||||
WARN_IF_UNDOCUMENTED = false
|
||||
|
|
@ -0,0 +1,315 @@
|
|||
/*!
|
||||
\defgroup pmp_namedparameters Named Parameters for Polygon Mesh Processing
|
||||
\ingroup PkgPolygonMeshProcessing
|
||||
|
||||
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 polygon mesh processing 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 polygon mesh processing 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 PkgBGLSummary
|
||||
(see \ref bgl_namedparameters), are used in this package:
|
||||
|
||||
\cgalNPTableBegin
|
||||
\cgalNPBegin{vertex_point_map} \anchor PMP_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{vertex_index_map} \anchor PMP_vertex_index_map
|
||||
is the property map containing the index of each vertex of the input polygon mesh.\n
|
||||
<b>Type:</b> a class model of `ReadablePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type and the value type
|
||||
\code typename boost::property_traits<typename boost::property_map<PolygonMesh, CGAL::vertex_index_t>::type>::value_type \endcode
|
||||
<b>Default:</b> \code boost::get(CGAL::vertex_index, pmesh)\endcode
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{face_index_map} \anchor PMP_face_index_map
|
||||
is the property map containing the index of each face of the input polygon mesh.\n
|
||||
<b>Type:</b> a class model of `ReadablePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%face_descriptor` as key type and the value type:
|
||||
\code typename boost::property_traits<typename boost::property_map<PolygonMesh, CGAL::face_index_t>::type>::value_type \endcode
|
||||
<b>Default:</b> \code boost::get(CGAL::face_index, pmesh)\endcode
|
||||
If this internal property map exists, its values should be initialized.
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{edge_is_constrained_map} \anchor PMP_edge_is_constrained_map
|
||||
is the property map containing information about edges of the input polygon mesh
|
||||
being marked or not. In `isotropic_remeshing()` and `connected_components()`,
|
||||
the marked edges are constrained.\n
|
||||
<b>Type:</b> a class model of `ReadWritePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%edge_descriptor` as key type and
|
||||
`bool` as value type. It should be default constructible.\n
|
||||
<b>Default:</b> a default property map where no edge is constrained
|
||||
\cgalNPEnd
|
||||
\cgalNPTableEnd
|
||||
|
||||
In addition to these named parameters, this package offers the following named parameters:
|
||||
|
||||
\cgalNPTableBegin
|
||||
\cgalNPBegin{geom_traits} \anchor PMP_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{vertex_incident_patches_map} \anchor PMP_vertex_incident_patches_map
|
||||
is the property map containing the surface patches incident to each vertex of the input polygon mesh.\n
|
||||
<b>Type:</b> a class model of `LvaluePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type. Its value type
|
||||
must be a container of `boost::property_traits<PatchIdMap>::%value_type` and have a function `insert()`.
|
||||
A `std::set` or a `boost::unordered_set` are recommended, as a patch index may be
|
||||
inserted several times.\n
|
||||
<b>Default:</b> \code boost::get(CGAL::vertex_incident_patches_t, pmesh)\endcode
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{vertex_feature_degree_map} \anchor PMP_vertex_feature_degree_map
|
||||
is the property map containing the number of feature edges being incident to the vertices of the polygon mesh `pmesh`.\n
|
||||
<b>Type:</b> a class model of `ReadWritePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type and
|
||||
`int` as value type. It should be default constructible.\n
|
||||
<b>Default:</b> \code boost::get(CGAL::vertex_feature_degree_t(), pmesh) \endcode
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{vertex_is_constrained_map} \anchor PMP_vertex_is_constrained_map
|
||||
is the property map containing information about vertices of the input polygon mesh being constrained or not.
|
||||
Constrained vertices may be replaced by new vertices, but the number and location
|
||||
of vertices remain unchanged.\n
|
||||
<b>Type:</b> a class model of `ReadWritePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%vertex_descriptor` as key type and
|
||||
`bool` as value type. It should be default constructible.\n
|
||||
<b>Default:</b> a default property map where no vertex is constrained is provided.
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{face_patch_map} \anchor PMP_face_patch_map
|
||||
is a property map containing information about faces.
|
||||
It is particularly well-suited for preserving surface patch IDs,
|
||||
or face colors.
|
||||
The edges at the interface between surface patches are treated similarly
|
||||
to the ones of `edge_is_constrained_map`.\n
|
||||
<b>Type:</b> a class model of `ReadWritePropertyMap` with
|
||||
`boost::graph_traits<PolygonMesh>::%face_descriptor` as key type and
|
||||
the desired property, model of `CopyConstructible` as value type.\n
|
||||
<b>Default:</b> a default property map where each face is associated with the ID of
|
||||
the connected component it belongs to. Connected components are
|
||||
computed with respect to the constrained edges listed in the property map
|
||||
`edge_is_constrained_map`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{first_index} \anchor PMP_first_index
|
||||
is the index of the first surface patch.\n
|
||||
<b>Type:</b> `std::size_t`\n
|
||||
<b>Default:</b> 1
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{density_control_factor} \anchor PMP_density_control_factor
|
||||
controls the density of the mesh generated by refinement, with larger values causing denser refinements.
|
||||
The density of vertices in the refined region is this factor times higher than before refinement.\n
|
||||
<b>Type:</b> floating scalar value\n
|
||||
<b>Default:</b> `CGAL::sqrt(2)`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{fairing_continuity} \anchor PMP_fairing_continuity
|
||||
controls the tangential continuity of the output surface in `fair()`.
|
||||
The possible values are 0, 1 and 2, refering to the C<sup>0</sup>, C<sup>1</sup>
|
||||
and C<sup>2</sup> continuity.\n
|
||||
<b>Type:</b> \c unsigned \c int between 0 and 2\n
|
||||
<b>Default:</b> `1`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{sparse_linear_solver} \anchor PMP_sparse_linear_solver
|
||||
is the solver used in `fair()`.\n
|
||||
<b>Type:</b> a class model of `SparseLinearAlgebraWithFactorTraits_d`.\n
|
||||
<b>Default:</b> if \ref thirdpartyEigen "Eigen" 3.2 (or greater) is available and
|
||||
`CGAL_EIGEN3_ENABLED` is defined, then the following overload of `Eigen_solver_traits`
|
||||
is provided as default value:\n
|
||||
\code CGAL::Eigen_solver_traits<Eigen::SparseLU<CGAL::Eigen_sparse_matrix<double>::EigenType, Eigen::COLAMDOrdering<int> > > \endcode
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_iterations} \anchor PMP_number_of_iterations
|
||||
is the number of iterations of the sequence of iterations performed in `isotropic_remeshing()`.\n
|
||||
<b>Type:</b> \c unsigned \c int \n
|
||||
<b>Default:</b> `1`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{protect_constraints} \anchor PMP_protect_constraints
|
||||
enables the protection of constraints listed by \ref PMP_edge_is_constrained_map
|
||||
"edge_is_constrained_map" and boundary edges
|
||||
in `isotropic_remeshing()`. If `true`, constraint edges cannot be modified at all
|
||||
during the remeshing process.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `false`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{relax_constraints} \anchor PMP_relax_constraints
|
||||
enables the tangential relaxation step in `isotropic_remeshing()`
|
||||
to be performed on vertices that are endpoints of constraints listed
|
||||
by \ref PMP_edge_is_constrained_map "edge_is_constrained_map", and boundary edges.
|
||||
The vertices move along the constrained polylines they belong to.
|
||||
Corners (i.e. vertices incident to more than 2 constraints, and vertices listed in
|
||||
\ref PMP_vertex_is_constrained_map "vertex_is_constrained_map") are not allowed
|
||||
to move at all.
|
||||
If \ref PMP_protect_constraints "protect_constraints" is
|
||||
set to `true`, this parameter is ignored.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_relaxation_steps} \anchor PMP_number_of_relaxation_steps
|
||||
is the number of iterations of tangential relaxation that are performed at each iteration
|
||||
of `isotropic_remeshing()`. A larger number of relaxation steps lead to
|
||||
a more isotropic mesh.\n
|
||||
<b>Type:</b> \c unsigned \c int \n
|
||||
<b>Default:</b> `1`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{use_delaunay_triangulation} \anchor PMP_use_delaunay_triangulation
|
||||
enables the use of the Delaunay triangulation facet search space for hole filling functions.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{use_random_uniform_sampling} \anchor PMP_use_random_uniform_sampling
|
||||
is a parameter used in `sample_triangle_mesh()` to indicate if points should be picked
|
||||
in a random uniform way.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{use_grid_sampling} \anchor PMP_use_grid_sampling
|
||||
is a parameter used in `sample_triangle_mesh()` to indicate if points should be picked
|
||||
in on a grid in each face.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `false`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{use_monte_carlo_sampling} \anchor PMP_use_monte_carlo_sampling
|
||||
is a parameter used in `sample_triangle_mesh()` to indicate if points should be picked
|
||||
using a Monte-Carlo approach.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `false`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{sample_edges} \anchor PMP_sample_edges
|
||||
is a parameter used in `sample_triangle_mesh()` to indicate if a dedicated sampling
|
||||
of edges should be done.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{sample_vertices} \anchor PMP_sample_vertices
|
||||
is a parameter used in `sample_triangle_mesh()` to indicate if triangle vertices should
|
||||
be copied in the output iterator.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{sample_faces} \anchor PMP_sample_faces
|
||||
is a parameter used in `sample_triangle_mesh()` to indicate if the interior of faces
|
||||
should be considered for the sampling.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_points_on_faces} \anchor PMP_number_of_points_on_faces
|
||||
is a parameter used in `sample_triangle_mesh()` to set the number of points picked
|
||||
using the random uniform method on faces.\n
|
||||
<b>Type:</b> `std::size_t` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_points_on_edges} \anchor PMP_number_of_points_on_edges
|
||||
is a parameter used in `sample_triangle_mesh()` to set the number of points picked
|
||||
using the random uniform method on edges.\n
|
||||
<b>Type:</b> `std::size_t` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_points_per_face} \anchor PMP_number_of_points_per_face
|
||||
is a parameter used in `sample_triangle_mesh()` to set the number of points picked
|
||||
per face using the Monte-Carlo method.\n
|
||||
<b>Type:</b> `std::size_t` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_points_per_edge} \anchor PMP_number_of_points_per_edge
|
||||
is a parameter used in `sample_triangle_mesh()` to set the number of points picked
|
||||
per edge using the Monte-Carlo method.\n
|
||||
<b>Type:</b> `std::size_t` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{grid_spacing} \anchor PMP_grid_spacing
|
||||
is a parameter used in `sample_triangle_mesh()` to set the grid spacing when using
|
||||
the grid sampling method.\n
|
||||
<b>Type:</b> `double` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_points_per_area_unit} \anchor PMP_number_of_points_per_area_unit
|
||||
is a parameter used in `sample_triangle_mesh()` to set the number of points per
|
||||
area unit to be picked up in faces for the random uniform sampling and
|
||||
Monte-Carlo methods.\n
|
||||
<b>Type:</b> `double` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{number_of_points_per_distance_unit} \anchor PMP_number_of_points_per_distance_unit
|
||||
is a parameter used in `sample_triangle_mesh()` to set the number of points per
|
||||
distance unit to be picked up on edges for the random uniform sampling and
|
||||
Monte-Carlo methods.\n
|
||||
<b>Type:</b> `double` \n
|
||||
<b>Default:</b> `0`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{do_project} \anchor PMP_do_project
|
||||
is a parameter used in `random_perturbation()` to set whether vertices should be re-projected
|
||||
to the input surface after their geometric perturbation.\n
|
||||
<b>Type:</b> `bool` \n
|
||||
<b>Default:</b> `true`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{random_seed} \anchor PMP_random_seed
|
||||
is a parameter used in `random_perturbation()` to choose a seed to initialize
|
||||
the random number generator `CGAL::Random()`.
|
||||
If this parameter is not provided, the perturbation is not deterministic
|
||||
(i.e. not reproducible from one run to the other).\n
|
||||
<b>Type:</b> `unsigned int` \n
|
||||
<b>Default:</b> the random number generator is initialized with `CGAL::Random()`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPBegin{outward_orientation} \anchor PMP_outward_orientation
|
||||
Parameter used in orientation functions to choose between an outward or inward orientation.
|
||||
\n
|
||||
\b Type : `bool` \n
|
||||
\b Default value is `true`
|
||||
|
||||
\cgalNPBegin{do_overlap_test_of_bounded_sides} \anchor PMP_do_overlap_test_of_bounded_sides
|
||||
Parameter used in intersection test functions to indicate whether overlapping tests of bounded sides
|
||||
of close meshes should be done in addition to surface intersection tests.
|
||||
\n
|
||||
\b Type : `bool` \n
|
||||
\b Default value is `false`
|
||||
\cgalNPEnd
|
||||
|
||||
\cgalNPTableEnd
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,9 @@
|
|||
namespace CGAL {
|
||||
/*!
|
||||
\mainpage User Manual
|
||||
\anchor Chapter_OptimalBoundingBox
|
||||
|
||||
\cgalAutoToc
|
||||
\authors
|
||||
|
||||
} /* namespace CGAL */
|
||||
|
|
@ -0,0 +1,30 @@
|
|||
/// \defgroup PkgOptimalBoundingBox Optimal Bounding Box Reference
|
||||
/// \defgroup PkgOptimalBoundingBoxConcepts Concepts
|
||||
/// \ingroup PkgOptimalBoundingBox
|
||||
|
||||
\cgalPkgDescriptionBegin{Optimal Bounding Box, PkgOptimalBoundingBoxSummary}
|
||||
\cgalPkgPicture{}
|
||||
|
||||
\cgalPkgSummaryBegin
|
||||
\cgalPkgAuthor{}
|
||||
\cgalPkgDesc{This package provides stuff.}
|
||||
\cgalPkgManuals{Chapter_OptimalBoundingBox,PkgOptimalBoundingBox}
|
||||
\cgalPkgSummaryEnd
|
||||
|
||||
\cgalPkgShortInfoBegin
|
||||
\cgalPkgSince{5.2}
|
||||
\cgalPkgDependsOn{documented for each function;}
|
||||
\cgalPkgBib{cgal:lty-pmp}
|
||||
\cgalPkgLicense{\ref licensesGPL "GPL"}
|
||||
\cgalPkgDemo{Polyhedron demo,polyhedron_3.zip}
|
||||
\cgalPkgShortInfoEnd
|
||||
|
||||
\cgalPkgDescriptionEnd
|
||||
|
||||
\cgalClassifedRefPages
|
||||
|
||||
## Parameters ##
|
||||
|
||||
## Functions ##
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,14 @@
|
|||
Manual
|
||||
Kernel_23
|
||||
STL_Extension
|
||||
Algebraic_foundations
|
||||
Circulator
|
||||
Stream_support
|
||||
Polyhedron
|
||||
BGL
|
||||
Solver_interface
|
||||
Surface_mesh
|
||||
Surface_mesh_deformation
|
||||
AABB_tree
|
||||
Triangulation_2
|
||||
Spatial_sorting
|
||||
|
|
@ -0,0 +1,4 @@
|
|||
/*!
|
||||
|
||||
\example Optimal_bounding_box/example_todo.cpp
|
||||
*/
|
||||
|
|
@ -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("example.cpp")
|
||||
|
|
@ -0,0 +1,252 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_EVOLUTION_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_EVOLUTION_H
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/population.h>
|
||||
#include <CGAL/Optimal_bounding_box/nelder_mead_functions.h>
|
||||
#include <CGAL/Optimal_bounding_box/fitness_function.h>
|
||||
|
||||
#include <CGAL/Random.h>
|
||||
#include <CGAL/number_utils.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Optimal_bounding_box {
|
||||
|
||||
template <typename Linear_algebra_traits>
|
||||
class Evolution
|
||||
{
|
||||
typedef typename Linear_algebra_traits::MatrixXd MatrixXd;
|
||||
typedef typename Linear_algebra_traits::Matrix3d Matrix3d;
|
||||
typedef typename Linear_algebra_traits::Vector3d Vector3d;
|
||||
|
||||
public:
|
||||
Evolution(Population<Linear_algebra_traits>& pop, MatrixXd& points)
|
||||
: population(pop), point_data(points)
|
||||
{}
|
||||
|
||||
void genetic_algorithm()
|
||||
{
|
||||
// random permutations
|
||||
std::size_t m = population.size();
|
||||
|
||||
//groups 1,2 : size m/2 groups 3,4 : size (m - m/2). m/2 is floored
|
||||
std::size_t size_first_group = m/2;
|
||||
std::size_t size_second_group = m - m/2;
|
||||
|
||||
std::vector<std::size_t> ids1(m/2), ids2(m/2);
|
||||
std::vector<std::size_t> ids3(m - m/2), ids4(m - m/2);
|
||||
|
||||
CGAL::Random rng;
|
||||
int im = static_cast<int>(m);
|
||||
std::generate(ids1.begin(), ids1.end(),
|
||||
[&rng, &im] () { return rng.get_int(0, im); });
|
||||
std::generate(ids2.begin(), ids2.end(),
|
||||
[&rng, &im] () { return rng.get_int(0, im); });
|
||||
std::generate(ids3.begin(), ids3.end(),
|
||||
[&rng, &im] () { return rng.get_int(0, im); });
|
||||
std::generate(ids4.begin(), ids4.end(),
|
||||
[&rng, &im] () { return rng.get_int(0, im); });
|
||||
|
||||
Population<Linear_algebra_traits> group1(m/2), group2(m/2);
|
||||
Population<Linear_algebra_traits> group3(m - m/2), group4(m - m/2);
|
||||
|
||||
for(std::size_t i = 0; i < ids1.size(); ++i)
|
||||
group1[i] = population[ids1[i]];
|
||||
|
||||
for(std::size_t i = 0; i < ids2.size(); ++i)
|
||||
group2[i] = population[ids2[i]];
|
||||
|
||||
for(std::size_t i = 0; i < ids3.size(); ++i)
|
||||
group3[i] = population[ids3[i]];
|
||||
|
||||
for(std::size_t i = 0; i < ids4.size(); ++i)
|
||||
group4[i] = population[ids4[i]];
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
check_det(group1);
|
||||
check_det(group2);
|
||||
check_det(group3);
|
||||
check_det(group4);
|
||||
#endif
|
||||
|
||||
// crossover I
|
||||
Population<Linear_algebra_traits> offspringsA(size_first_group);
|
||||
double bias = 0.1;
|
||||
|
||||
for(std::size_t i = 0; i < size_first_group; ++i)
|
||||
{
|
||||
std::vector<Matrix3d> offspring(4);
|
||||
for(int j = 0; j < 4; ++j)
|
||||
{
|
||||
double r = rng.get_double();
|
||||
double fitnessA = compute_fitness<Linear_algebra_traits>(group1[i][j], point_data);
|
||||
double fitnessB = compute_fitness<Linear_algebra_traits>(group2[i][j], point_data);
|
||||
double threshold;
|
||||
|
||||
if(fitnessA < fitnessB)
|
||||
threshold = 0.5 + bias;
|
||||
else
|
||||
threshold = 0.5 - bias;
|
||||
|
||||
if(r < threshold) // choose A
|
||||
offspring[j] = group1[i][j];
|
||||
else // choose B
|
||||
offspring[j] = group2[i][j];
|
||||
}
|
||||
offspringsA[i] = offspring;
|
||||
}
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
std::cout << "offspringsA: \n" ;
|
||||
check_det(offspringsA);
|
||||
#endif
|
||||
|
||||
// crossover II
|
||||
Population<Linear_algebra_traits> offspringsB(size_second_group);
|
||||
bias = 0.1;
|
||||
|
||||
for(std::size_t i = 0; i < size_second_group; ++i)
|
||||
{
|
||||
std::vector<Matrix3d> offspring(4);
|
||||
for(int j = 0; j < 4; ++j)
|
||||
{
|
||||
double fitnessA = compute_fitness<Linear_algebra_traits>(group3[i][j], point_data);
|
||||
double fitnessB = compute_fitness<Linear_algebra_traits>(group4[i][j], point_data);
|
||||
double lambda;
|
||||
if(fitnessA < fitnessB)
|
||||
lambda = 0.5 + bias;
|
||||
else
|
||||
lambda = 0.5 - bias;
|
||||
// combine information from A and B
|
||||
offspring[j] = lambda * group3[i][j] + lambda * group4[i][j];
|
||||
}
|
||||
|
||||
// qr factorization of the offspring
|
||||
Linear_algebra_traits::qr_factorization(offspring);
|
||||
offspringsB[i] = offspring;
|
||||
}
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
std::cout << "offspringsB: \n" ;
|
||||
check_det(offspringsB);
|
||||
#endif
|
||||
|
||||
CGAL_assertion(offspringsA.size() == size_first_group);
|
||||
CGAL_assertion(offspringsB.size() == size_second_group);
|
||||
CGAL_assertion(offspringsA.size() + offspringsB.size() == population.size());
|
||||
|
||||
// next generatrion
|
||||
for(std::size_t i = 0; i < size_first_group; ++i)
|
||||
population[i] = offspringsA[i];
|
||||
|
||||
for(std::size_t i = 0; i < size_second_group; ++i)
|
||||
population[size_first_group + i] = offspringsB[i];
|
||||
}
|
||||
|
||||
void evolve(std::size_t generations)
|
||||
{
|
||||
// hardcoded nelder_mead_iterations
|
||||
std::size_t nelder_mead_iterations = 20;
|
||||
|
||||
// stopping criteria prameters
|
||||
double prev_fit_value = 0;
|
||||
double new_fit_value = 0;
|
||||
double tolerance = 1e-2;
|
||||
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();
|
||||
|
||||
#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<Linear_algebra_traits>(population[s], point_data, nelder_mead_iterations);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
//std::cout << "pop after nelder mead: " << std::endl;
|
||||
//pop.show_population();
|
||||
//std::cout << std::endl;
|
||||
Fitness_map<Linear_algebra_traits, Matrix3d, MatrixXd> fitness_map_debug(population, point_data);
|
||||
Matrix3d R_now = fitness_map_debug.get_best();
|
||||
std::cout << "det= " << Linear_algebra_traits::determinant(R_now) << std::endl;
|
||||
#endif
|
||||
|
||||
// stopping criteria
|
||||
Fitness_map<Linear_algebra_traits, Matrix3d, MatrixXd> fitness_map(population, point_data);
|
||||
new_fit_value = fitness_map.get_best_fitness_value();
|
||||
double difference = new_fit_value - prev_fit_value;
|
||||
|
||||
if(CGAL::abs(difference) < tolerance * new_fit_value)
|
||||
stale++;
|
||||
|
||||
if(stale == 5)
|
||||
break;
|
||||
|
||||
prev_fit_value = new_fit_value;
|
||||
}
|
||||
}
|
||||
|
||||
const Matrix3d get_best()
|
||||
{
|
||||
Fitness_map<Linear_algebra_traits, Matrix3d, MatrixXd> fitness_map(population, point_data);
|
||||
return fitness_map.get_best();
|
||||
}
|
||||
|
||||
private:
|
||||
// data
|
||||
Population<Linear_algebra_traits> population;
|
||||
MatrixXd point_data;
|
||||
};
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
template <typename Simplex>
|
||||
void check_det(Population<Simplex>& pop)
|
||||
{
|
||||
for(std::size_t i = 0; i < pop.size(); ++i)
|
||||
{
|
||||
for(std::size_t j = 0; j < 4; ++j)
|
||||
{
|
||||
auto A = pop[i][j]; // Simplex
|
||||
std::cout << Linear_algebra_traits::determinant(A) << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
} // end namespace Optimal_bounding_box
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_EVOLUTION_H
|
||||
|
|
@ -0,0 +1,117 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_FITNESS_FUNCTION_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_FITNESS_FUNCTION_H
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/population.h>
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
|
||||
#include <limits>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Optimal_bounding_box {
|
||||
|
||||
template <typename Linear_algebra_traits, typename Vertex, typename Matrix>
|
||||
double compute_fitness(const Vertex& R, const Matrix& data)
|
||||
{
|
||||
// R: rotation matrix
|
||||
CGAL_assertion(R.cols() == 3);
|
||||
CGAL_assertion(R.rows() == 3);
|
||||
// data: points
|
||||
CGAL_assertion(data.cols() == 3);
|
||||
CGAL_assertion(data.rows() >= 3);
|
||||
|
||||
typedef typename Linear_algebra_traits::Vector3d Vector3d;
|
||||
typedef typename Linear_algebra_traits::Index Index;
|
||||
|
||||
double xmin, xmax, ymin, ymax, zmin, zmax;
|
||||
for(Index i = 0; i < static_cast<Index>(data.rows()); ++i){
|
||||
|
||||
Vector3d vec = Linear_algebra_traits::row3(data, i);
|
||||
vec = R * vec;
|
||||
|
||||
if(i == 0){
|
||||
xmin = xmax = vec.coeff(0);
|
||||
ymin = ymax = vec.coeff(1);
|
||||
zmin = zmax = vec.coeff(2);
|
||||
}else {
|
||||
if(vec.coeff(0) < xmin) xmin = vec.coeff(0);
|
||||
if(vec.coeff(1) < ymin) ymin = vec.coeff(1);
|
||||
if(vec.coeff(2) < zmin) zmin = vec.coeff(2);
|
||||
if(vec.coeff(0) > xmax) xmax = vec.coeff(0);
|
||||
if(vec.coeff(1) > ymax) ymax = vec.coeff(1);
|
||||
if(vec.coeff(2) > zmax) zmax = vec.coeff(2);
|
||||
}
|
||||
}
|
||||
|
||||
CGAL_assertion(xmax > xmin);
|
||||
CGAL_assertion(ymax > ymin);
|
||||
CGAL_assertion(zmax > zmin);
|
||||
|
||||
// volume
|
||||
return ((xmax - xmin) * (ymax - ymin) * (zmax - zmin));
|
||||
}
|
||||
|
||||
template <typename Linear_algebra_traits, typename Vertex, typename Matrix>
|
||||
struct Fitness_map
|
||||
{
|
||||
Fitness_map(Population<Linear_algebra_traits>& p, Matrix& points)
|
||||
: pop(p), points(points)
|
||||
{}
|
||||
|
||||
const Vertex get_best()
|
||||
{
|
||||
std::size_t simplex_id, vertex_id;
|
||||
double best_fitness = std::numeric_limits<int>::max();
|
||||
for(std::size_t i = 0; i < pop.size(); ++i)
|
||||
{
|
||||
for(std::size_t j =0; j < 4; ++j)
|
||||
{
|
||||
const Vertex vertex = pop[i][j];
|
||||
const double fitness = compute_fitness<Linear_algebra_traits>(vertex, points);
|
||||
if (fitness < best_fitness)
|
||||
{
|
||||
simplex_id = i;
|
||||
vertex_id = j;
|
||||
best_fitness = fitness;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return pop[simplex_id][vertex_id];
|
||||
}
|
||||
|
||||
double get_best_fitness_value()
|
||||
{
|
||||
const Vertex best_mat = get_best();
|
||||
return compute_fitness<Linear_algebra_traits>(best_mat, points);
|
||||
}
|
||||
|
||||
Population<Linear_algebra_traits>& pop;
|
||||
const Matrix& points;
|
||||
};
|
||||
|
||||
} // end namespace Optimal_bounding_box
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_FITNESS_FUNCTION_H
|
||||
|
|
@ -0,0 +1,113 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_HELPER_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_HELPER_H
|
||||
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
|
||||
#include <CGAL/Bbox_3.h>
|
||||
#include <CGAL/boost/graph/helpers.h>
|
||||
#include <CGAL/number_utils.h>
|
||||
#include <CGAL/Simple_cartesian.h>
|
||||
#include <CGAL/Surface_mesh.h>
|
||||
|
||||
#include <fstream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Optimal_bounding_box {
|
||||
|
||||
template <typename Matrix, typename Point>
|
||||
void fill_matrix(const std::vector<Point>& v_points, Matrix& points_mat)
|
||||
{
|
||||
points_mat.resize(v_points.size(), 3);
|
||||
for(std::size_t i = 0; i < v_points.size(); ++i)
|
||||
{
|
||||
Point p = v_points[i];
|
||||
points_mat.set_coef(i, 0, CGAL::to_double(p.x()));
|
||||
points_mat.set_coef(i, 1, CGAL::to_double(p.y()));
|
||||
points_mat.set_coef(i, 2, CGAL::to_double(p.z()));
|
||||
}
|
||||
}
|
||||
|
||||
template <typename SurfaceMesh, typename Matrix>
|
||||
void sm_to_matrix(SurfaceMesh& sm, Matrix& mat)
|
||||
{
|
||||
typedef typename boost::property_map<SurfaceMesh, boost::vertex_point_t>::const_type Vpm;
|
||||
typedef typename boost::property_traits<Vpm>::reference Point_ref;
|
||||
typedef typename boost::graph_traits<SurfaceMesh>::vertex_descriptor vertex_descriptor;
|
||||
Vpm vpm = get(boost::vertex_point, sm);
|
||||
|
||||
mat.resize(vertices(sm).size(), 3);
|
||||
std::size_t i = 0;
|
||||
for(vertex_descriptor v : vertices(sm))
|
||||
{
|
||||
Point_ref p = get(vpm, v);
|
||||
mat.set_coef(i, 0, CGAL::to_double(p.x()));
|
||||
mat.set_coef(i, 1, CGAL::to_double(p.y()));
|
||||
mat.set_coef(i, 2, CGAL::to_double(p.z()));
|
||||
++i;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename Point>
|
||||
double calculate_volume(const std::vector<Point>& points)
|
||||
{
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
|
||||
CGAL::Bbox_3 bbox = bbox_3(points.begin(), points.end());
|
||||
K::Iso_cuboid_3 ic(bbox);
|
||||
return ic.volume();
|
||||
}
|
||||
|
||||
// it is called after post processing in debug only
|
||||
template <typename Matrix>
|
||||
void matrix_to_mesh_and_draw(Matrix& data_points, std::string filename)
|
||||
{
|
||||
typedef CGAL::Simple_cartesian<double> K;
|
||||
typedef K::Point_3 Point;
|
||||
typedef CGAL::Surface_mesh<Point> Mesh;
|
||||
|
||||
// Simplex -> std::vector
|
||||
std::vector<Point> points;
|
||||
|
||||
for(int i = 0; i < data_points.rows(); ++i)
|
||||
{
|
||||
Point p(data_points(i, 0), data_points(i, 1), data_points(i, 2));
|
||||
points.push_back(p);
|
||||
}
|
||||
|
||||
Mesh mesh;
|
||||
CGAL::make_hexahedron(points[0], points[1], points[2], points[3],
|
||||
points[4], points[5], points[6], points[7], mesh);
|
||||
|
||||
std::ofstream out(filename);
|
||||
out << mesh;
|
||||
out.close();
|
||||
}
|
||||
|
||||
} // end namespace Optimal_bounding_box
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_HELPER_H
|
||||
|
|
@ -0,0 +1,184 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_NEALDER_MEAD_FUNCTIONS_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_NEALDER_MEAD_FUNCTIONS_H
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
#include <boost/iterator/counting_iterator.hpp>
|
||||
#include <CGAL/Optimal_bounding_box/fitness_function.h>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Optimal_bounding_box {
|
||||
|
||||
template <typename Linear_algebra_traits, typename Matrix>
|
||||
const Matrix reflection(const Matrix& S_centroid, const Matrix& S_worst)
|
||||
{
|
||||
CGAL_assertion(S_centroid.rows() == 3);
|
||||
CGAL_assertion(S_centroid.rows() == 3);
|
||||
CGAL_assertion(S_worst.cols() == 3);
|
||||
CGAL_assertion(S_worst.cols() == 3);
|
||||
|
||||
return S_centroid * Linear_algebra_traits::transpose(S_worst) * S_centroid;
|
||||
}
|
||||
|
||||
template <typename Linear_algebra_traits, typename Matrix>
|
||||
const Matrix expansion(const Matrix& S_centroid, const Matrix& S_worst, const Matrix& S_reflection)
|
||||
{
|
||||
CGAL_assertion(S_centroid.rows() == 3);
|
||||
CGAL_assertion(S_centroid.rows() == 3);
|
||||
CGAL_assertion(S_worst.cols() == 3);
|
||||
CGAL_assertion(S_worst.cols() == 3);
|
||||
CGAL_assertion(S_reflection.cols() == 3);
|
||||
CGAL_assertion(S_reflection.cols() == 3);
|
||||
|
||||
return S_centroid * Linear_algebra_traits::transpose(S_worst) * S_reflection;
|
||||
}
|
||||
|
||||
template <typename Linear_algebra_traits, typename Matrix>
|
||||
Matrix mean(const Matrix& m1, const Matrix& m2)
|
||||
{
|
||||
// same API for reduction
|
||||
CGAL_assertion(m1.rows() == 3);
|
||||
CGAL_assertion(m1.rows() == 3);
|
||||
CGAL_assertion(m2.cols() == 3);
|
||||
CGAL_assertion(m2.cols() == 3);
|
||||
|
||||
Matrix reduction = 0.5 * m1 + 0.5 * m2;
|
||||
Matrix Q = Linear_algebra_traits::qr_factorization(reduction);
|
||||
double det = Linear_algebra_traits::determinant(Q);
|
||||
return Q / det;
|
||||
}
|
||||
|
||||
template <typename Linear_algebra_traits, typename Matrix>
|
||||
const Matrix nm_centroid(const Matrix& S1, const Matrix& S2, const Matrix& S3)
|
||||
{
|
||||
Matrix mean = (S1 + S2 + S3) / 3.0;
|
||||
Matrix Q = Linear_algebra_traits::qr_factorization(mean);
|
||||
double det = Linear_algebra_traits::determinant(Q);
|
||||
return Q / det;
|
||||
}
|
||||
|
||||
// needed in nelder mead algorithm
|
||||
struct Comparator
|
||||
{
|
||||
Comparator(const std::vector<double>& in) : fitness(in) {}
|
||||
|
||||
inline bool operator() (std::size_t& i, std::size_t& j) {
|
||||
return fitness[i] < fitness[j];
|
||||
}
|
||||
|
||||
const std::vector<double>& fitness;
|
||||
};
|
||||
|
||||
// simplex: 4 rotation matrices are its vertices
|
||||
template <typename Linear_algebra_traits>
|
||||
void nelder_mead(std::vector<typename Linear_algebra_traits::Matrix3d>& simplex,
|
||||
const typename Linear_algebra_traits::MatrixXd& point_data,
|
||||
std::size_t nelder_mead_iterations)
|
||||
{
|
||||
CGAL_assertion(simplex.size() == 4); // tetrahedron
|
||||
|
||||
|
||||
typedef typename Linear_algebra_traits::Matrix3d Matrix3d;
|
||||
|
||||
std::vector<double> fitness(4);
|
||||
std::vector<std::size_t> indices(boost::counting_iterator<std::size_t>(0),
|
||||
boost::counting_iterator<std::size_t>(simplex.size()));
|
||||
|
||||
for(std::size_t t = 0; t < nelder_mead_iterations; ++t)
|
||||
{
|
||||
for(std::size_t i = 0; i < 4; ++i)
|
||||
{
|
||||
fitness[i] = compute_fitness<Linear_algebra_traits>(simplex[i], point_data);
|
||||
}
|
||||
|
||||
CGAL_assertion(fitness.size() == 4);
|
||||
CGAL_assertion(indices.size() == 4);
|
||||
|
||||
// get indices of sorted sequence
|
||||
Comparator compare_indices(fitness);
|
||||
std::sort(indices.begin(), indices.end(), compare_indices);
|
||||
|
||||
// new sorted simplex & fitness
|
||||
std::vector<Matrix3d> s_simplex(4);
|
||||
std::vector<double> s_fitness(4);
|
||||
for(int i = 0; i < 4; ++i)
|
||||
{
|
||||
s_simplex[i] = simplex[indices[i]];
|
||||
s_fitness[i] = fitness[indices[i]];
|
||||
}
|
||||
|
||||
simplex = s_simplex;
|
||||
fitness = s_fitness;
|
||||
|
||||
// centroid
|
||||
const Matrix3d v_centroid = nm_centroid<Linear_algebra_traits>(simplex[0], simplex[1], simplex[2]);
|
||||
|
||||
// find worst's vertex reflection
|
||||
const Matrix3d v_worst = simplex[3];
|
||||
const Matrix3d v_refl = reflection<Linear_algebra_traits>(v_centroid, v_worst);
|
||||
const double f_refl = compute_fitness<Linear_algebra_traits>(v_refl, point_data);
|
||||
|
||||
if(f_refl < fitness[2])
|
||||
{
|
||||
if(f_refl >= fitness[0]) // if reflected point is not better than the best
|
||||
{
|
||||
// do reflection
|
||||
simplex[3] = v_refl;
|
||||
}
|
||||
else
|
||||
{
|
||||
// expansion
|
||||
const Matrix3d v_expand = expansion<Linear_algebra_traits>(v_centroid, v_worst, v_refl);
|
||||
const double f_expand = compute_fitness<Linear_algebra_traits>(v_expand, point_data);
|
||||
if(f_expand < f_refl)
|
||||
simplex[3] = v_expand;
|
||||
else
|
||||
simplex[3] = v_refl;
|
||||
}
|
||||
}
|
||||
else // if reflected vertex is not better
|
||||
{
|
||||
const Matrix3d v_mean = mean<Linear_algebra_traits>(v_centroid, v_worst);
|
||||
const double f_mean = compute_fitness<Linear_algebra_traits>(v_mean, point_data);
|
||||
if(f_mean <= fitness[3])
|
||||
// contraction of worst
|
||||
simplex[3] = v_mean;
|
||||
else
|
||||
{
|
||||
// reduction: move all vertices towards the best
|
||||
for(std::size_t i=1; i < 4; ++i)
|
||||
{
|
||||
simplex[i] = mean<Linear_algebra_traits>(simplex[i], simplex[0]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
CGAL_assertion(simplex.size() == 4); // tetrahedron
|
||||
} // iterations
|
||||
}
|
||||
|
||||
} // end namespace Optimal_bounding_box
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,273 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_OBB_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_OBB_H
|
||||
|
||||
#include <CGAL/Optimal_bounding_box/population.h>
|
||||
#include <CGAL/Optimal_bounding_box/evolution.h>
|
||||
#include <CGAL/Optimal_bounding_box/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/Simple_cartesian.h>
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
#include <CGAL/Timer.h>
|
||||
#endif
|
||||
|
||||
#if defined(CGAL_EIGEN3_ENABLED)
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
#endif
|
||||
|
||||
#include <iostream>
|
||||
#include <iterator>
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Optimal_bounding_box {
|
||||
|
||||
// works on matrices only
|
||||
/// \cond SKIP_IN_MANUAL
|
||||
template <typename Linear_algebra_traits, typename Vertex, typename Matrix>
|
||||
void post_processing(const Matrix& points, Vertex& R, Matrix& obb)
|
||||
{
|
||||
CGAL_assertion(points.cols() == 3);
|
||||
CGAL_assertion(R.rows() == 3);
|
||||
CGAL_assertion(R.cols() == 3);
|
||||
CGAL_assertion(obb.rows() == 8);
|
||||
CGAL_assertion(obb.cols() == 3);
|
||||
|
||||
// 1) rotate points with R
|
||||
Matrix rotated_points(points.rows(), points.cols());
|
||||
rotated_points = points * Linear_algebra_traits::transpose(R);
|
||||
|
||||
// 2) get AABB from rotated points
|
||||
typedef CGAL::Simple_cartesian<double> K;
|
||||
typedef K::Point_3 Point;
|
||||
typedef typename Linear_algebra_traits::Index index;
|
||||
|
||||
// Simplex -> std::vector
|
||||
std::vector<Point> v_points;
|
||||
for(index i = 0; i < static_cast<index>(rotated_points.rows()); ++i)
|
||||
{
|
||||
Point p(rotated_points(i, 0), rotated_points(i, 1), rotated_points(i, 2));
|
||||
v_points.push_back(p);
|
||||
}
|
||||
CGAL::Bbox_3 bbox;
|
||||
bbox = bbox_3(v_points.begin(), v_points.end());
|
||||
K::Iso_cuboid_3 ic(bbox);
|
||||
|
||||
Matrix aabb(8, 3);
|
||||
for(std::size_t i = 0; i < 8; ++i)
|
||||
{
|
||||
aabb.set_coef(i, 0, ic[i].x());
|
||||
aabb.set_coef(i, 1, ic[i].y());
|
||||
aabb.set_coef(i, 2, ic[i].z());
|
||||
}
|
||||
|
||||
// 3) apply inverse rotation to rotated AABB
|
||||
obb = aabb * R;
|
||||
}
|
||||
/// \endcond
|
||||
|
||||
/// \ingroup OBB_grp
|
||||
/// calculates the optimal bounding box.
|
||||
///
|
||||
/// @tparam Point the point type
|
||||
/// @tparam LinearAlgebraTraits a model of `LinearAlgebraTraits`. If no instance of `LinearAlgebraTraits`
|
||||
/// is provided, then `CGAL::Eigen_linear_algebra_traits` is used.
|
||||
///
|
||||
/// @param points the input points that are included in the optimal bounding box.
|
||||
/// @param obb_points the eight points of the optimal bounding box to be calculated.
|
||||
/// @param use_ch a bool flag to indicating whether to use the convex hull of the input points
|
||||
/// as an optimization step.
|
||||
template <typename Point, typename LinearAlgebraTraits>
|
||||
void compute_optimal_bounding_box(const std::vector<Point>& points,
|
||||
std::vector<Point>& obb_points,
|
||||
LinearAlgebraTraits&,
|
||||
bool use_ch)
|
||||
{
|
||||
CGAL_assertion(points.size() >= 3);
|
||||
|
||||
if(obb_points.size() != 8)
|
||||
obb_points.resize(8);
|
||||
CGAL_assertion(obb_points.size() == 8);
|
||||
|
||||
// eigen linear algebra traits
|
||||
typedef typename LinearAlgebraTraits::MatrixXd MatrixXd;
|
||||
typedef typename LinearAlgebraTraits::Matrix3d Matrix3d;
|
||||
MatrixXd points_mat;
|
||||
|
||||
if(use_ch) // get the ch3
|
||||
{
|
||||
std::vector<Point> ch_points;
|
||||
CGAL::extreme_points_3(points, std::back_inserter(ch_points));
|
||||
CGAL::Optimal_bounding_box::fill_matrix(ch_points, points_mat);
|
||||
}
|
||||
else
|
||||
{
|
||||
CGAL::Optimal_bounding_box::fill_matrix(points, points_mat);
|
||||
}
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
CGAL::Timer timer;
|
||||
#endif
|
||||
|
||||
std::size_t max_generations = 100;
|
||||
Population<LinearAlgebraTraits> pop(50);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
timer.start();
|
||||
#endif
|
||||
|
||||
CGAL::Optimal_bounding_box::Evolution<LinearAlgebraTraits> search_solution(pop, points_mat);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
timer.stop();
|
||||
std::cout << "constructor: " << timer.time() << std::endl;
|
||||
timer.reset();
|
||||
timer.start();
|
||||
#endif
|
||||
|
||||
search_solution.evolve(max_generations);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
timer.stop();
|
||||
std::cout << "evolve: " << timer.time() << std::endl;
|
||||
timer.reset();
|
||||
timer.start();
|
||||
#endif
|
||||
|
||||
Matrix3d rotation = search_solution.get_best();
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
timer.stop();
|
||||
std::cout << "get best: " << timer.time() << std::endl;
|
||||
#endif
|
||||
|
||||
MatrixXd obb; // may be preallocated at compile time
|
||||
obb.resize(8, 3);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
timer.reset();
|
||||
timer.start();
|
||||
#endif
|
||||
|
||||
post_processing<LinearAlgebraTraits>(points_mat, rotation, obb);
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_BENCHMARKS
|
||||
timer.stop();
|
||||
std::cout << "post porcessing: " << timer.time() << std::endl;
|
||||
#endif
|
||||
|
||||
// matrix -> vector
|
||||
for(std::size_t i = 0; i < 8; ++i)
|
||||
{
|
||||
Point p(obb(i, 0), obb(i, 1), obb(i, 2));
|
||||
obb_points[i] = p;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename Point>
|
||||
void compute_optimal_bounding_box(const std::vector<Point>& points,
|
||||
std::vector<Point>& obb_points,
|
||||
bool use_ch)
|
||||
{
|
||||
#if defined(CGAL_EIGEN3_ENABLED)
|
||||
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
|
||||
#else
|
||||
#pragma message("Error: You must either provide linear traits or link CGAL with the Eigen library")
|
||||
Linear_algebra_traits; // no parameter provided, and Eigen is not enabled --> don't compile!
|
||||
#endif
|
||||
|
||||
Linear_algebra_traits la_traits;
|
||||
compute_optimal_bounding_box(points, obb_points, la_traits, use_ch);
|
||||
}
|
||||
|
||||
/// \ingroup OBB_grp
|
||||
/// calculates the optimal bounding box.
|
||||
///
|
||||
/// @tparam PolygonMesh a model of `FaceListGraph`
|
||||
/// @tparam LinearAlgebraTraits a model of `LinearAlgebraTraits`. If no instance of `LinearAlgebraTraits`
|
||||
/// is provided, then `CGAL::Eigen_linear_algebra_traits` is used.
|
||||
///
|
||||
/// @param pmesh the input mesh.
|
||||
/// @param obbmesh the hexaedron mesh to be built out of the optimal bounding box.
|
||||
/// @param la_traits an instance of the linear algebra traits.
|
||||
/// @param use_ch a bool flag to indicating whether to use the convex hull of the input points
|
||||
/// as an optimization step.
|
||||
template <typename PolygonMesh, typename LinearAlgebraTraits>
|
||||
void compute_optimal_bounding_box(const PolygonMesh& pmesh,
|
||||
PolygonMesh& obbmesh,
|
||||
LinearAlgebraTraits& la_traits,
|
||||
bool use_ch)
|
||||
{
|
||||
CGAL_assertion(vertices(pmesh).size() >= 3);
|
||||
|
||||
if(vertices(pmesh).size() <= 3)
|
||||
{
|
||||
std::cerr << "Not enough points in the mesh!\n";
|
||||
return;
|
||||
}
|
||||
|
||||
typedef typename boost::graph_traits<PolygonMesh>::vertex_descriptor vertex_descriptor;
|
||||
typedef typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type Vpm;
|
||||
typedef typename boost::property_traits<Vpm>::value_type Point;
|
||||
|
||||
std::vector<Point> points;
|
||||
Vpm pmap = get(boost::vertex_point, pmesh);
|
||||
BOOST_FOREACH(vertex_descriptor v, vertices(pmesh))
|
||||
points.push_back(get(pmap, v));
|
||||
|
||||
|
||||
std::vector<Point> obb_points;
|
||||
compute_optimal_bounding_box(points, obb_points, la_traits, use_ch);
|
||||
|
||||
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
|
||||
obb_points[4], obb_points[5], obb_points[6], obb_points[7], obbmesh);
|
||||
}
|
||||
|
||||
template <typename PolygonMesh>
|
||||
void compute_optimal_bounding_box(const PolygonMesh& pmesh,
|
||||
PolygonMesh& obbmesh,
|
||||
bool use_ch)
|
||||
{
|
||||
#if defined(CGAL_EIGEN3_ENABLED)
|
||||
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
|
||||
#else
|
||||
#pragma message("Error: You must either provide linear traits or link CGAL with the Eigen library")
|
||||
Linear_algebra_traits; // no parameter provided, and Eigen is not enabled --> don't compile!
|
||||
#endif
|
||||
|
||||
Linear_algebra_traits la_traits;
|
||||
compute_optimal_bounding_box(pmesh, obbmesh, la_traits, use_ch);
|
||||
}
|
||||
|
||||
} // end namespace Optimal_bounding_box
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_OBB_H
|
||||
|
|
@ -0,0 +1,136 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_OPTIMAL_BOUNDING_BOX_POPULATION_H
|
||||
#define CGAL_OPTIMAL_BOUNDING_BOX_POPULATION_H
|
||||
|
||||
#include <CGAL/assertions.h>
|
||||
#include <CGAL/Random.h>
|
||||
|
||||
#include <vector>
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
namespace Optimal_bounding_box {
|
||||
|
||||
template<typename Linear_algebra_traits>
|
||||
class Population
|
||||
{
|
||||
typedef typename Linear_algebra_traits::Matrix3d Matrix;
|
||||
typedef std::vector<Matrix> Simplex;
|
||||
|
||||
public:
|
||||
Population(std::size_t size)
|
||||
: n(size), random_generator(CGAL::Random())
|
||||
{
|
||||
// reserve pop space
|
||||
pop.reserve(n);
|
||||
|
||||
// create simplices
|
||||
for(std::size_t i = 0 ; i < n; ++i)
|
||||
{
|
||||
Simplex simplex(4);
|
||||
create_simplex(simplex);
|
||||
CGAL_assertion(simplex.size() == 4);
|
||||
pop.push_back(simplex);
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
void show_population();
|
||||
#endif
|
||||
|
||||
std::size_t size(){return n;}
|
||||
|
||||
// access simplex
|
||||
Simplex& operator[](std::size_t i)
|
||||
{
|
||||
CGAL_assertion(i < n);
|
||||
return pop[i];
|
||||
}
|
||||
|
||||
const Simplex& operator[](std::size_t i) const
|
||||
{
|
||||
CGAL_assertion(i < n);
|
||||
return pop[i];
|
||||
}
|
||||
|
||||
private:
|
||||
// create random population
|
||||
void create_simplex(Simplex& simplex)
|
||||
{
|
||||
CGAL_assertion(simplex.size() == 4);
|
||||
for(std::size_t i = 0; i < 4; ++i)
|
||||
{
|
||||
Matrix R;
|
||||
if(R.cols() == 0 || R.rows() == 0)
|
||||
R.resize(3, 3);
|
||||
|
||||
create_vertex(R);
|
||||
Matrix Q = Linear_algebra_traits::qr_factorization(R);
|
||||
|
||||
simplex[i] = Q;
|
||||
}
|
||||
CGAL_assertion(simplex.size() == 4);
|
||||
}
|
||||
|
||||
void create_vertex(Matrix& R)
|
||||
{
|
||||
CGAL_assertion(R.rows() == 3);
|
||||
CGAL_assertion(R.cols() == 3);
|
||||
|
||||
for(std::size_t i = 0; i < 3; ++i)
|
||||
{
|
||||
for(std::size_t j = 0; j < 3; ++j)
|
||||
{
|
||||
R.set_coef(i, j, random_generator.get_double());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
std::size_t n;
|
||||
CGAL::Random random_generator;
|
||||
std::vector<Simplex> pop;
|
||||
};
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG
|
||||
template <typename Matrix>
|
||||
void Population<Matrix>::show_population()
|
||||
{
|
||||
std::size_t id = 0;
|
||||
for(const Simplex i : pop)
|
||||
{
|
||||
CGAL_assertion(i.size() == 4);
|
||||
std:: cout << "Simplex: "<< id++ << std::endl;
|
||||
for(const Matrix R : i)
|
||||
{
|
||||
std::cout << R; // eigen out
|
||||
std::cout << "\n\n";
|
||||
}
|
||||
std:: cout << std:: endl;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
} // end namespace Optimal_bounding_box
|
||||
} // end namespace CGAL
|
||||
|
||||
#endif // CGAL_OPTIMAL_BOUNDING_BOX_POPULATION_H
|
||||
|
|
@ -0,0 +1 @@
|
|||
GeometryFactory (France)
|
||||
|
|
@ -0,0 +1 @@
|
|||
GPL (v3 or later)
|
||||
|
|
@ -0,0 +1 @@
|
|||
GeometryFactory
|
||||
|
|
@ -0,0 +1,26 @@
|
|||
# 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_linear_algebra_functions.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,301 @@
|
|||
#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_qr_factorization()
|
||||
{
|
||||
typedef CGAL::Eigen_dense_matrix<double, 3, 3> Mat;
|
||||
Mat A(3, 3);
|
||||
A.set_coef(0, 0, 0.3011944);
|
||||
A.set_coef(0, 1, 0.9932761);
|
||||
A.set_coef(0, 2, 0.5483701);
|
||||
A.set_coef(1, 0, 0.5149142);
|
||||
A.set_coef(1, 1, 0.5973263);
|
||||
A.set_coef(1, 2, 0.5162336);
|
||||
A.set_coef(2, 0, 0.0039213);
|
||||
A.set_coef(2, 1, 0.0202949);
|
||||
A.set_coef(2, 2, 0.9240308);
|
||||
|
||||
CGAL_assertion_code(Mat Q = CGAL::Eigen_linear_algebra_traits::qr_factorization(A));
|
||||
CGAL_assertion_code(double epsilon = 1e-6);
|
||||
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));
|
||||
}
|
||||
|
||||
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_simplex_operations()
|
||||
{
|
||||
CGAL_assertion_code(typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits);
|
||||
typedef CGAL::Eigen_dense_matrix<double, 3, 3> Matrix;
|
||||
|
||||
Matrix Sc(3, 3);
|
||||
Sc.set_coef(0, 0, -0.809204);
|
||||
Sc.set_coef(0, 1, 0.124296);
|
||||
Sc.set_coef(0, 2, 0.574230);
|
||||
Sc.set_coef(1, 0, -0.574694);
|
||||
Sc.set_coef(1, 1, 0.035719);
|
||||
Sc.set_coef(1, 2, -0.817589);
|
||||
Sc.set_coef(2, 0, -0.122134);
|
||||
Sc.set_coef(2, 1, -0.991602);
|
||||
Sc.set_coef(2, 2, 0.042528);
|
||||
|
||||
Matrix S_worst(3, 3);
|
||||
S_worst.set_coef(0, 0, -0.45070);
|
||||
S_worst.set_coef(0, 1, -0.32769);
|
||||
S_worst.set_coef(0, 2, -0.83035);
|
||||
S_worst.set_coef(1, 0, -0.13619);
|
||||
S_worst.set_coef(1, 1, -0.89406);
|
||||
S_worst.set_coef(1, 2, 0.42675);
|
||||
S_worst.set_coef(2, 0, -0.88222);
|
||||
S_worst.set_coef(2, 1, 0.30543);
|
||||
S_worst.set_coef(2, 2, 0.35833);
|
||||
|
||||
CGAL_assertion_code(Matrix Sr = CGAL::Optimal_bounding_box::reflection<Linear_algebra_traits>(Sc, S_worst));
|
||||
CGAL_assertion_code(double epsilon = 1e-5);
|
||||
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));
|
||||
|
||||
CGAL_assertion_code(Matrix Se =
|
||||
CGAL::Optimal_bounding_box::expansion<Linear_algebra_traits>(Sc, S_worst, Sr));
|
||||
CGAL_assertion(assert_doubles(Se(0,0), -0.87991, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(0,1), 0.36105, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(0,2), -0.30888, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(1,0), -0.11816, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(1,1), -0.79593, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(1,2), -0.59375, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(2,0), -0.460215, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(2,1), -0.48595, epsilon));
|
||||
CGAL_assertion(assert_doubles(Se(2,2), 0.74300, epsilon));
|
||||
|
||||
Matrix S_a(3, 3);
|
||||
S_a.set_coef(0, 0, -0.277970);
|
||||
S_a.set_coef(0, 1, 0.953559);
|
||||
S_a.set_coef(0, 2, 0.116010);
|
||||
S_a.set_coef(1, 0, -0.567497);
|
||||
S_a.set_coef(1, 1, -0.065576);
|
||||
S_a.set_coef(1, 2, -0.820760);
|
||||
S_a.set_coef(2, 0, -0.775035);
|
||||
S_a.set_coef(2, 1, -0.293982);
|
||||
S_a.set_coef(2, 2, 0.559370);
|
||||
|
||||
Matrix S_b(3, 3);
|
||||
S_b.set_coef(0, 0, -0.419979);
|
||||
S_b.set_coef(0, 1, 0.301765);
|
||||
S_b.set_coef(0, 2, -0.8558940);
|
||||
S_b.set_coef(1, 0, -0.653011);
|
||||
S_b.set_coef(1, 1, -0.755415);
|
||||
S_b.set_coef(1, 2, 0.054087);
|
||||
S_b.set_coef(2, 0, -0.630234);
|
||||
S_b.set_coef(2, 1, 0.581624);
|
||||
S_b.set_coef(2, 2, 0.514314);
|
||||
|
||||
CGAL_assertion_code(Matrix S_c =
|
||||
CGAL::Optimal_bounding_box::mean<Linear_algebra_traits>(S_a, S_b));
|
||||
CGAL_assertion(assert_doubles(S_c(0,0), -0.35111, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(0,1), 0.79308, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(0,2), -0.49774, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(1,0), -0.61398, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(1,1), -0.59635, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(1,2), -0.51710, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(2,0), -0.70693, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(2,1), 0.12405, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_c(2,2), 0.69632, epsilon));
|
||||
}
|
||||
|
||||
void test_centroid()
|
||||
{
|
||||
CGAL_assertion_code(typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits);
|
||||
typedef CGAL::Eigen_dense_matrix<double, 3, 3> Matrix;
|
||||
|
||||
Matrix S_a;
|
||||
S_a.set_coef(0, 0, -0.588443);
|
||||
S_a.set_coef(0, 1, 0.807140);
|
||||
S_a.set_coef(0, 2, -0.047542);
|
||||
S_a.set_coef(1, 0, -0.786228);
|
||||
S_a.set_coef(1, 1, -0.584933);
|
||||
S_a.set_coef(1, 2, -0.199246);
|
||||
S_a.set_coef(2, 0, -0.188629);
|
||||
S_a.set_coef(2, 1, -0.079867);
|
||||
S_a.set_coef(2, 2, 0.978795);
|
||||
|
||||
Matrix S_b(3, 3);
|
||||
S_b.set_coef(0, 0, -0.2192721);
|
||||
S_b.set_coef(0, 1, 0.2792986);
|
||||
S_b.set_coef(0, 2, -0.9348326);
|
||||
S_b.set_coef(1, 0, -0.7772152);
|
||||
S_b.set_coef(1, 1, -0.6292092);
|
||||
S_b.set_coef(1, 2, -0.005686);
|
||||
S_b.set_coef(2, 0, -0.5897934);
|
||||
S_b.set_coef(2, 1, 0.7253193);
|
||||
S_b.set_coef(2, 2, 0.3550431);
|
||||
|
||||
Matrix S_c(3, 3);
|
||||
S_c.set_coef(0, 0, -0.32657);
|
||||
S_c.set_coef(0, 1, -0.60013);
|
||||
S_c.set_coef(0, 2, -0.730206);
|
||||
S_c.set_coef(1, 0, -0.20022);
|
||||
S_c.set_coef(1, 1, -0.71110);
|
||||
S_c.set_coef(1, 2, 0.67398);
|
||||
S_c.set_coef(2, 0, -0.92372);
|
||||
S_c.set_coef(2, 1, 0.36630);
|
||||
S_c.set_coef(2, 2, 0.11207);
|
||||
|
||||
CGAL_assertion_code(Matrix S_centroid =
|
||||
CGAL::Optimal_bounding_box::nm_centroid<Linear_algebra_traits>(S_a, S_b, S_c));
|
||||
CGAL_assertion_code(double epsilon = 1e-5);
|
||||
CGAL_assertion(assert_doubles(S_centroid(0,0), -0.419979, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(0,1), 0.301765, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(0,2), -0.855894, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(1,0), -0.653011, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(1,1), -0.755415, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(1,2), 0.054087, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(2,0), -0.630234, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(2,1), 0.581624, epsilon));
|
||||
CGAL_assertion(assert_doubles(S_centroid(2,2), 0.514314, epsilon));
|
||||
}
|
||||
|
||||
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::qr_factorization(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_qr_factorization();
|
||||
test_fitness_function();
|
||||
test_simplex_operations();
|
||||
test_centroid();
|
||||
test_eigen_matrix_interface();
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
|
@ -0,0 +1,423 @@
|
|||
#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_nelder_mead()
|
||||
{
|
||||
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
|
||||
typedef Linear_algebra_traits::Matrix3d Matrix3d;
|
||||
typedef Linear_algebra_traits::MatrixXd MatrixXd;
|
||||
|
||||
MatrixXd data_points(4,3);
|
||||
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;
|
||||
|
||||
// one simplex
|
||||
std::vector<Matrix3d> simplex(4);
|
||||
Matrix3d v0(3,3);
|
||||
Matrix3d v1(3,3);
|
||||
Matrix3d v2(3,3);
|
||||
Matrix3d v3(3,3);
|
||||
|
||||
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::nelder_mead<Linear_algebra_traits>(simplex, data_points, nm_iterations);
|
||||
|
||||
CGAL_assertion_code(double epsilon = 1e-5);
|
||||
CGAL_assertion_code(Matrix3d v0_new = simplex[0]);
|
||||
CGAL_assertion(assert_doubles(v0_new(0,0), -0.288975, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(0,1), 0.7897657, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(0,2), -0.541076, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(1,0), -0.9407046, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(1,1), -0.3391466, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(1,2), 0.0073817, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(2,0), -0.1776743, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(2,1), 0.5111260, epsilon));
|
||||
CGAL_assertion(assert_doubles(v0_new(2,2), 0.84094, epsilon));
|
||||
|
||||
CGAL_assertion_code(Matrix3d v1_new = simplex[1]);
|
||||
CGAL_assertion(assert_doubles(v1_new(0,0), -0.458749, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(0,1), 0.823283, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(0,2), -0.334296, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(1,0), -0.885235, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(1,1), -0.455997, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(1,2), 0.091794, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(2,0), -0.076866, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(2,1), 0.338040, epsilon));
|
||||
CGAL_assertion(assert_doubles(v1_new(2,2), 0.937987, epsilon));
|
||||
|
||||
CGAL_assertion_code(Matrix3d v2_new = simplex[2]);
|
||||
CGAL_assertion(assert_doubles(v2_new(0,0), -0.346582, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(0,1), 0.878534, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(0,2), -0.328724, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(1,0), -0.936885, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(1,1), -0.341445, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(1,2), 0.075251, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(2,0), -0.046131, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(2,1), 0.334057, epsilon));
|
||||
CGAL_assertion(assert_doubles(v2_new(2,2), 0.941423, epsilon));
|
||||
|
||||
CGAL_assertion_code(Matrix3d v3_new = simplex[3]);
|
||||
CGAL_assertion(assert_doubles(v3_new(0,0), -0.394713, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(0,1), 0.791782, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(0,2), -0.466136, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(1,0), -0.912112, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(1,1), -0.398788, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(1,2), 0.094972, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(2,0), -0.110692, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(2,1), 0.462655, epsilon));
|
||||
CGAL_assertion(assert_doubles(v3_new(2,2), 0.879601, epsilon));
|
||||
}
|
||||
|
||||
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();
|
||||
CGAL_assertion(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;
|
||||
CGAL_assertion_code(typedef Linear_algebra_traits::Matrix3d Matrix3d);
|
||||
std::size_t generations = 10;
|
||||
CGAL::Optimal_bounding_box::Population<Linear_algebra_traits> pop(50);
|
||||
CGAL::Optimal_bounding_box::Evolution<Linear_algebra_traits> evolution(pop, data_points);
|
||||
evolution.evolve(generations);
|
||||
|
||||
CGAL_assertion_code(Matrix3d R = evolution.get_best());
|
||||
CGAL_assertion_code(double epsilon = 1e-3);
|
||||
CGAL_assertion(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(0,0), -0.25791, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(0,1), 0.796512, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(0,2), -0.546855, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(1,0), -0.947128, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(1,1), -0.320242, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(1,2), -0.0197553, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(2,0), -0.190861, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(2,1), 0.512847, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(2,2), 0.836992, epsilon));
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
|
||||
// postprocessing
|
||||
CGAL::Eigen_dense_matrix<double> obb(8, 3);
|
||||
CGAL::Optimal_bounding_box::post_processing(data_points, R, obb);
|
||||
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, "data/random_unit_tetra_result.off");
|
||||
#endif
|
||||
}
|
||||
|
||||
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::MatrixXd MatrixXd;
|
||||
CGAL_assertion_code(typedef Linear_algebra_traits::Matrix3d Matrix3d);
|
||||
|
||||
// points in a matrix
|
||||
MatrixXd 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);
|
||||
|
||||
CGAL_assertion_code(Matrix3d R = experiment.get_best());
|
||||
CGAL_assertion_code(double epsilon = 1e-5);
|
||||
CGAL_assertion(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
|
||||
// postprocessing
|
||||
MatrixXd obb(8, 3);
|
||||
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
|
||||
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, "data/OBB.off");
|
||||
#endif
|
||||
}
|
||||
|
||||
void test_long_tetrahedron(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::MatrixXd MatrixXd;
|
||||
CGAL_assertion_code(typedef Linear_algebra_traits::Matrix3d Matrix3d);
|
||||
|
||||
// points in a matrix
|
||||
MatrixXd 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);
|
||||
|
||||
CGAL_assertion_code(Matrix3d R = experiment.get_best());
|
||||
CGAL_assertion_code(double epsilon = 1e-3);
|
||||
CGAL_assertion(assert_doubles(Linear_algebra_traits::determinant(R), 1, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(0,0), -1, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(0,1), 0, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(0,2), 0, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(1,0), 0, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(1,1), -0.707107, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(1,2), 0.707106, epsilon) ||
|
||||
assert_doubles(R(1,2), -0.707106, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(2,0), 0, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(2,1), 0.707106, epsilon) ||
|
||||
assert_doubles(R(1,2), -0.707106, epsilon));
|
||||
CGAL_assertion(assert_doubles(R(2,2), 0.707107, epsilon));
|
||||
|
||||
#ifdef CGAL_OPTIMAL_BOUNDING_BOX_DEBUG_TEST
|
||||
// postprocessing
|
||||
MatrixXd obb(8, 3);
|
||||
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
|
||||
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, fname + "result.off");
|
||||
#endif
|
||||
}
|
||||
|
||||
void test_compute_obb_evolution(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);
|
||||
BOOST_FOREACH(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);
|
||||
|
||||
CGAL_assertion_code(double epsilon = 1e-3);
|
||||
CGAL_assertion_code(double vol = CGAL::Optimal_bounding_box::calculate_volume(obb_points));
|
||||
CGAL_assertion(assert_doubles(vol, 0.883371, epsilon));
|
||||
|
||||
#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(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(std::string fname1, 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);
|
||||
|
||||
BOOST_FOREACH(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);
|
||||
|
||||
CGAL_assertion_code(double epsilon = 1e-3);
|
||||
CGAL_assertion_code(double vol = CGAL::Optimal_bounding_box::calculate_volume(obb_points));
|
||||
CGAL_assertion(assert_doubles(vol, 0.883371, epsilon));
|
||||
|
||||
// test two
|
||||
CGAL::Surface_mesh< K::Point_3> obbmesh;
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(mesh2, obbmesh, true);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test_nelder_mead();
|
||||
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_polyhedron_item scene_polyhedron_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,191 @@
|
|||
#include <QtCore/qglobal.h>
|
||||
#include <CGAL/Three/Scene_item.h>
|
||||
#include <CGAL/Three/Scene_interface.h>
|
||||
|
||||
#include <QAction>
|
||||
#include <QMainWindow>
|
||||
#include <QApplication>
|
||||
|
||||
#include "Scene_surface_mesh_item.h"
|
||||
#include "Polyhedron_type.h"
|
||||
#include "Scene_polyhedron_item.h"
|
||||
#include <CGAL/Three/Polyhedron_demo_plugin_interface.h>
|
||||
#include "Scene_polyhedron_selection_item.h"
|
||||
#include "Scene_points_with_normal_item.h"
|
||||
|
||||
#include <boost/graph/graph_traits.hpp>
|
||||
#include <CGAL/Optimal_bounding_box/optimal_bounding_box.h>
|
||||
#include <CGAL/Eigen_linear_algebra_traits.h>
|
||||
|
||||
//typedef Scene_surface_mesh_item Scene_facegraph_item;
|
||||
//typedef Scene_facegraph_item Scene_facegraph_item;
|
||||
|
||||
typedef Scene_facegraph_item::Face_graph FaceGraph;
|
||||
typedef Polyhedron::Point_3 Point_3;
|
||||
using namespace CGAL::Three;
|
||||
|
||||
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->selectionIndices().size() == 1)
|
||||
{
|
||||
return qobject_cast<Scene_facegraph_item*>(scene->item(scene->mainSelectionIndex()))
|
||||
|| qobject_cast<Scene_polyhedron_selection_item*>(scene->item(scene->mainSelectionIndex()));
|
||||
}
|
||||
|
||||
Q_FOREACH(int index, scene->selectionIndices())
|
||||
{
|
||||
if (qobject_cast<Scene_facegraph_item*>(scene->item(index)))
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
*/
|
||||
|
||||
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* poly_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(poly_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(poly_item != NULL)
|
||||
{
|
||||
FaceGraph& pmesh = *poly_item->polyhedron();
|
||||
selected_vertices.assign(vertices(pmesh).begin(), vertices(pmesh).end());
|
||||
PointPMap pmap = get(CGAL::vertex_point, pmesh);
|
||||
BOOST_FOREACH(vertex_descriptor v, selected_vertices)
|
||||
points.push_back(get(pmap, v));
|
||||
|
||||
}
|
||||
else if(selection_item != NULL) // using selection of faces
|
||||
{
|
||||
FaceGraph& pmesh = *selection_item->polyhedron();
|
||||
BOOST_FOREACH(face_descriptor f, selection_item->selected_facets)
|
||||
{
|
||||
BOOST_FOREACH(vertex_descriptor v, vertices_around_face(halfedge(f, pmesh), pmesh))
|
||||
{
|
||||
selected_vertices.push_back(v);
|
||||
}
|
||||
}
|
||||
|
||||
PointPMap pmap = get(CGAL::vertex_point, pmesh);
|
||||
BOOST_FOREACH(vertex_descriptor v, selected_vertices)
|
||||
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;
|
||||
BOOST_FOREACH(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
|
||||
CGAL::Eigen_linear_algebra_traits la_traits;
|
||||
std::vector<Point_3> obb_points(8);
|
||||
CGAL::Optimal_bounding_box::compute_optimal_bounding_box(points, obb_points, la_traits, true);
|
||||
|
||||
Scene_item* item;
|
||||
if(mw->property("is_polyhedorn_mode").toBool())
|
||||
{
|
||||
Polyhedron* p = new Polyhedron;
|
||||
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
|
||||
obb_points[4], obb_points[5], obb_points[6], obb_points[7], *p);
|
||||
item = new Scene_polyhedron_item(p);
|
||||
}
|
||||
else {
|
||||
SMesh* p = new SMesh;
|
||||
CGAL::make_hexahedron(obb_points[0], obb_points[1], obb_points[2], obb_points[3],
|
||||
obb_points[4], obb_points[5], obb_points[6], obb_points[7], *p);
|
||||
item = new Scene_surface_mesh_item(p);
|
||||
}
|
||||
|
||||
item->setName("Optimal bbox mesh");
|
||||
item->setRenderingMode(Wireframe);
|
||||
scene->addItem(item);
|
||||
}
|
||||
|
||||
#include "Create_obb_mesh_plugin.moc"
|
||||
|
|
@ -0,0 +1,326 @@
|
|||
// Copyright (c) 2018 GeometryFactory (France).
|
||||
// All rights reserved.
|
||||
//
|
||||
// This file is part of CGAL (www.cgal.org).
|
||||
// You can redistribute it and/or modify it under the terms of the GNU
|
||||
// General Public License as published by the Free Software Foundation,
|
||||
// either version 3 of the License, or (at your option) any later version.
|
||||
//
|
||||
// Licensees holding a valid commercial license may use this file in
|
||||
// accordance with the commercial license agreement provided with the software.
|
||||
//
|
||||
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
|
||||
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// $URL$
|
||||
// $Id$
|
||||
// SPDX-License-Identifier: GPL-3.0+
|
||||
//
|
||||
//
|
||||
// Author(s) : Konstantinos Katrioplas
|
||||
|
||||
#ifndef CGAL_EIGEN_LINEAR_ALGEBRA_TRAITS_H
|
||||
#define CGAL_EIGEN_LINEAR_ALGEBRA_TRAITS_H
|
||||
|
||||
#include <CGAL/basic.h>
|
||||
#include <Eigen/Dense>
|
||||
#include <Eigen/QR>
|
||||
|
||||
|
||||
namespace CGAL {
|
||||
|
||||
|
||||
template<typename T, int D = Eigen::Dynamic>
|
||||
class Eigen_dense_vector;
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolver
|
||||
|
||||
The class `Eigen_dense_matrix` is a wrapper around \ref thirdpartyEigen "Eigen" matrix type
|
||||
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::DenseMatrix`</a>.
|
||||
|
||||
\tparam T Number type.
|
||||
\tparam D1 Number of rows, or Dynamic.
|
||||
\tparam D1 Number of cols, or Dynamic.
|
||||
|
||||
\sa `CGAL::Eigen_dense_vector<T, D>`
|
||||
\sa `CGAL::Eigen_linear_algebra_traits`
|
||||
*/
|
||||
template<typename T, int D1 = Eigen::Dynamic, int D2 = Eigen::Dynamic>
|
||||
class Eigen_dense_matrix
|
||||
{
|
||||
public:
|
||||
|
||||
/// The internal matrix type from \ref thirdpartyEigen "Eigen".
|
||||
typedef Eigen::Matrix<T, D1, D2> EigenType;
|
||||
|
||||
/// Create a dense matrix
|
||||
Eigen_dense_matrix(std::size_t nrows, std::size_t ncols)
|
||||
: m_matrix(static_cast<int>(nrows), static_cast<int>(ncols))
|
||||
{
|
||||
CGAL_assertion(m_matrix.rows() > 0);
|
||||
CGAL_assertion(m_matrix.cols() > 0);
|
||||
}
|
||||
|
||||
/// Create a dense matrix
|
||||
Eigen_dense_matrix(int nrows, int ncols)
|
||||
: m_matrix(nrows, ncols)
|
||||
{
|
||||
CGAL_assertion(m_matrix.rows() > 0);
|
||||
CGAL_assertion(m_matrix.cols() > 0);
|
||||
}
|
||||
|
||||
/// Create a dense matrix out of a \ref thirdpartyEigen "Eigen" matrix type
|
||||
Eigen_dense_matrix(const EigenType& eigen_mat)
|
||||
: m_matrix(eigen_mat) {}
|
||||
|
||||
Eigen_dense_matrix() : m_matrix() {}
|
||||
|
||||
/// Read access to a matrix coefficient.
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
T& operator() (int i_, int j_)
|
||||
{
|
||||
return m_matrix(i_, j_);
|
||||
}
|
||||
|
||||
/// Write access to a matrix coefficient: a_ij <- val
|
||||
///
|
||||
/// \pre 0 <= i < row_dimension().
|
||||
/// \pre 0 <= j < column_dimension().
|
||||
void set_coef(std::size_t i_, std::size_t j_, T val)
|
||||
{
|
||||
int i = static_cast<int>(i_);
|
||||
int j = static_cast<int>(j_);
|
||||
CGAL_precondition(i < m_matrix.rows());
|
||||
CGAL_precondition(j < m_matrix.cols());
|
||||
|
||||
m_matrix.coeffRef(i,j) = val;
|
||||
}
|
||||
|
||||
/// Return the matrix number of rows
|
||||
std::size_t rows() const {return m_matrix.rows();}
|
||||
/// Return the matrix number of cols
|
||||
std::size_t cols() const {return m_matrix.cols();}
|
||||
|
||||
/// Resize to i rows and j cols
|
||||
void resize(int i_, int j_) { m_matrix.resize(i_, j_);}
|
||||
|
||||
const T& coeff(int i_) const
|
||||
{
|
||||
return m_matrix.coeff(i_);
|
||||
}
|
||||
|
||||
EigenType m_matrix;
|
||||
};
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolver
|
||||
|
||||
The class `Eigen_vector` is a wrapper around \ref thirdpartyEigen "Eigen" dense vector
|
||||
type <a href="http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html"> </a>,
|
||||
which is a simple array of numbers.
|
||||
|
||||
\cgalModels `SvdTraits::Vector`
|
||||
\cgalModels `SparseLinearAlgebraTraits_d::Vector`.
|
||||
|
||||
\tparam T Number type.
|
||||
\tparam D Number of colums, or Dynamic.
|
||||
|
||||
\sa `CGAL::Eigen_dense_matrix<T>`
|
||||
\sa `CGAL::Eigen_linear_algebra_traits<T>`
|
||||
*/
|
||||
template <typename T, int D>
|
||||
class Eigen_dense_vector
|
||||
{
|
||||
private:
|
||||
|
||||
/// The internal vector type from \ref thirdpartyEigen "Eigen".
|
||||
typedef Eigen::Matrix<T, D, 1> EigenType;
|
||||
|
||||
public:
|
||||
|
||||
/// Create a dense vector out of a \ref thirdpartyEigen "Eigen" vector type
|
||||
Eigen_dense_vector(const EigenType& vec) : m_vector(vec) {}
|
||||
|
||||
/// Read and write and access to a vector coefficient: `a_i`
|
||||
const T& coeff(std::size_t i)
|
||||
{
|
||||
CGAL_assertion(i >= 0);
|
||||
CGAL_assertion(i < static_cast<std::size_t>(D));
|
||||
return m_vector.coeff(i);
|
||||
}
|
||||
|
||||
EigenType m_vector;
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
\ingroup PkgSolver
|
||||
|
||||
The class `Eigen_linear_algebra_traits` provides an interface to linear algebra functionalities of \ref thirdpartyEigen "Eigen".
|
||||
\ref thirdpartyEigen "Eigen" version 3.1 (or later) must be available on the system.
|
||||
|
||||
\sa `CGAL::Eigen_dense_matrix<T, D1, D2>`
|
||||
\sa `CGAL::Eigen_dense_vector<T, D>`
|
||||
\sa http://eigen.tuxfamily.org
|
||||
|
||||
Example
|
||||
--------------
|
||||
|
||||
\code{.cpp}
|
||||
|
||||
typedef CGAL::Eigen_linear_algebra_traits Linear_algebra_traits;
|
||||
|
||||
// dynamic matrix at run time to store a large amount of data
|
||||
typedef Linear_algebra_traits::MatrixXd MatrixXd;
|
||||
|
||||
// preallocated 3x3 matrix at compile time
|
||||
typedef Linear_algebra_traits::Matrix3d Matrix3d;
|
||||
|
||||
// preallocated 3-cols vector at compile time
|
||||
typedef Linear_algebra_traits::Vector3d Vector3d;
|
||||
|
||||
\endcode
|
||||
*/
|
||||
|
||||
class Eigen_linear_algebra_traits
|
||||
{
|
||||
public:
|
||||
typedef double NT;
|
||||
typedef int Index;
|
||||
|
||||
// dynamic size at run time
|
||||
typedef CGAL::Eigen_dense_matrix<NT> MatrixXd;
|
||||
|
||||
// dynamic rows in run time, fixed cols in compile time
|
||||
typedef CGAL::Eigen_dense_matrix<NT, Eigen::Dynamic, 3> MatrixX3d;
|
||||
|
||||
// fixed at compile time
|
||||
typedef CGAL::Eigen_dense_matrix<NT, 3, 3> Matrix3d;
|
||||
|
||||
// fixed at compile time
|
||||
typedef CGAL::Eigen_dense_vector<NT, 3> Vector3d;
|
||||
|
||||
/// Get the transpose of a `CGAL::Eigen_dense_matrix<T, D1, D2>` matrix
|
||||
template <class Matrix>
|
||||
static Matrix transpose(const Matrix& mat)
|
||||
{
|
||||
return Matrix(mat.m_matrix.transpose());
|
||||
}
|
||||
|
||||
/// Get the determinant of a `CGAL::Eigen_dense_matrix<T, D1, D2>` matrix
|
||||
template <class Matrix>
|
||||
static NT determinant(const Matrix& mat)
|
||||
{
|
||||
return mat.m_matrix.determinant();
|
||||
}
|
||||
|
||||
/// Performs QR decomposition of matrix A to a unitary matrix and an upper triagonal
|
||||
/// and returns the unitary matrix.
|
||||
template <class NT, int D1, int D2>
|
||||
static CGAL::Eigen_dense_matrix<NT, D1, D2> qr_factorization(const CGAL::Eigen_dense_matrix<NT, D1, D2>& A)
|
||||
{
|
||||
Eigen::HouseholderQR<Eigen::Matrix<NT, D1, D2> > qr(A.m_matrix);
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D2>(qr.householderQ());
|
||||
}
|
||||
|
||||
template <class Matrix>
|
||||
static void qr_factorization(std::vector<Matrix>& simplex)
|
||||
{
|
||||
for(std::size_t i = 0; i < simplex.size(); ++i)
|
||||
{
|
||||
Matrix mat = simplex[i].m_matrix;
|
||||
simplex[i] = qr_factorization(mat);
|
||||
}
|
||||
}
|
||||
|
||||
// CGAL::Eigen_dense_vector<NT, D2> : the returned type with D2 may be -1 in compile time,
|
||||
// and may not be equal to the expected type.
|
||||
// Eigen manages to return a precompiled row out of a dynamic matrix but I don't know how.
|
||||
|
||||
/// Get the row vector out of a `CGAL::Eigen_dense_matrix<T, D1, D2>`. The result is stored in a
|
||||
/// preallocated at compile time 3-column `CGAL::Eigen_dense_vector<T, 3>`
|
||||
template <class NT, int D1, int D2>
|
||||
static CGAL::Eigen_dense_vector<NT, 3> row3(const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
|
||||
int i)
|
||||
{
|
||||
return CGAL::Eigen_dense_vector<NT, 3>(A.m_matrix.row(i));
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
/// Matrix multiplication. If the columns of A and the rows of B are equal at compile time,
|
||||
/// the product is stored at a preallocated at compile time `CGAL::Eigen_dense_matrix`. Otherwise,
|
||||
/// the product is stored in a dynamic at run time matrix.
|
||||
template <class NT, int D1, int D2, int D3>
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D3> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2 >& A,
|
||||
const CGAL::Eigen_dense_matrix<NT, D2, D3 >& B)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D3>(A.m_matrix * B.m_matrix);
|
||||
}
|
||||
|
||||
// D2 and D3 may not be equal at compile time, but equal at run time!
|
||||
// This overload returns a dynamic matrix.
|
||||
template <class NT, int D1, int D2, int D3, int D4>
|
||||
const CGAL::Eigen_dense_matrix<NT> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2 >& A,
|
||||
const CGAL::Eigen_dense_matrix<NT, D3, D4 >& B)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT>(A.m_matrix * B.m_matrix);
|
||||
}
|
||||
|
||||
// scalar - matrix multiplication
|
||||
template <class NT, int D1, int D2>
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator* (const NT& scalar,
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2>& B)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D2>(scalar * B.m_matrix);
|
||||
}
|
||||
|
||||
template <class NT, int D1, int D2>
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
|
||||
const NT& scalar)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D2>(A.m_matrix * scalar);
|
||||
}
|
||||
|
||||
template <class NT, int D1, int D2>
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator/ (const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
|
||||
const double& scalar)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D2>(A.m_matrix / scalar);
|
||||
}
|
||||
|
||||
template <class NT, int D1, int D2>
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator/ (const double& scalar,
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> & A)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D2> (scalar / A.m_matrix);
|
||||
}
|
||||
|
||||
// addition
|
||||
template <class NT, int D1, int D2>
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> operator+ (const CGAL::Eigen_dense_matrix<NT, D1, D2> & A,
|
||||
const CGAL::Eigen_dense_matrix<NT, D1, D2> & B)
|
||||
{
|
||||
return CGAL::Eigen_dense_matrix<NT, D1, D2> (A.m_matrix + B.m_matrix);
|
||||
}
|
||||
|
||||
// vector - matrix multiplication
|
||||
template <class NT, int D1, int D2>
|
||||
const Eigen_dense_vector<NT, D1> operator* (const CGAL::Eigen_dense_matrix<NT, D1, D2>& A,
|
||||
const CGAL::Eigen_dense_vector<NT, D2>& V)
|
||||
{
|
||||
return Eigen_dense_vector<NT, D1>(A.m_matrix * V.m_vector);
|
||||
}
|
||||
|
||||
|
||||
|
||||
}
|
||||
// end namespace
|
||||
|
||||
|
||||
|
||||
#endif // CGAL_EIGEN_LINEAR_ALGEBRA_TRAITS_H
|
||||
Loading…
Reference in New Issue