This commit is contained in:
Konstantinos Katrioplas 2018-05-03 12:12:51 +02:00
parent 319c8079ad
commit 37446e4672
7 changed files with 129 additions and 152 deletions

View File

@ -79,10 +79,10 @@ struct Fitness_map // -> a free function
double best_fitness = std::numeric_limits<int>::max();
for(std::size_t i = 0; i < pop.size(); ++i)
{
const std::vector<Matrix> simplex = pop[i];
//const std::vector<Matrix> simplex = pop[i];
for(std::size_t j =0; j < 4; ++j)
{
const Matrix vertex = simplex[j];
const Matrix vertex = pop[i][j];
//std::cout << "i= "<< i << " j=" << j<<"\n vertex= " << vertex << std::endl;
++count_vertices;
//std::cout << "vertex = " << vertex << std::endl;
@ -95,14 +95,13 @@ struct Fitness_map // -> a free function
vertex_id = j;
best_fitness = fitness;
}
}
}
std::vector<Matrix> best_simplex = pop[simplex_id];
Matrix temp = best_simplex[vertex_id];
//Matrix temp = best_simplex[vertex_id];
return temp;
return best_simplex[vertex_id];
}

View File

@ -74,7 +74,7 @@ const Matrix expansion(const Matrix& S_centroid, const Matrix& S_worst, const Ma
}
template <typename Matrix>
Matrix mean(const Matrix& m1, const Matrix& m2) // mean
Matrix mean(const Matrix& m1, const Matrix& m2)
{
// same API for reduction
CGAL_assertion(m1.rows() == 3);
@ -100,13 +100,6 @@ const Matrix centroid(Matrix& S1, Matrix& S2, Matrix& S3)
}
}} // end namespaces

View File

@ -33,8 +33,8 @@
#include <Eigen/Dense>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Surface_mesh.h> // used draw mesh
#include <CGAL/Polyhedron_3.h> // used to get the ch
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
@ -42,9 +42,6 @@ namespace CGAL {
namespace Optimal_bounding_box {
template <typename Matrix>
void evolution(Matrix& R, Matrix& points, std::size_t max_generations) // todo: points is const
{
@ -54,10 +51,8 @@ void evolution(Matrix& R, Matrix& points, std::size_t max_generations) // todo:
CGAL_assertion(R.rows() == 3);
CGAL_assertion(R.cols() == 3);
std::size_t nelder_mead_iterations = 20;
Population<Matrix> pop(50);
std::size_t nelder_mead_iterations = 20;
double prev_fit_value = 0;
double new_fit_value = 0;
@ -66,27 +61,32 @@ void evolution(Matrix& R, Matrix& points, std::size_t max_generations) // todo:
for(std::size_t t = 0; t < max_generations; ++t)
{
#ifdef OBB_DEBUG
std::cout << "generation= " << t << "\n";
#endif
genetic_algorithm(pop, points);
#ifdef OBB_DEBUG
//std::cout << "pop after genetic" << std::endl;
// pop.show_population();
// std::cout << std::endl;
//std::cin.get();
//pop.show_population();
//std::cout << std::endl;
#endif
for(std::size_t s = 0; s < pop.size(); ++s)
nelder_mead(pop[s], points, nelder_mead_iterations);
#ifdef OBB_DEBUG
//std::cout << "pop after nelder mead: " << std::endl;
//pop.show_population();
//std::cout << std::endl;
//std::cin.get();
// debugging
/*
Fitness_map<Matrix> fitness_map(pop, points);
Matrix R_now = fitness_map.get_best();
std::cout << "det= " << R_now.determinant() << std::endl;
*/
#endif
// stopping criteria
Fitness_map<Matrix> fitness_map(pop, points);
@ -102,7 +102,6 @@ void evolution(Matrix& R, Matrix& points, std::size_t max_generations) // todo:
prev_fit_value = new_fit_value;
}
// compute fitness of entire population
Fitness_map<Matrix> fitness_map(pop, points);
R = fitness_map.get_best();
}
@ -164,7 +163,7 @@ void fill_matrix(std::vector<Point>& v_points, Matrix& points_mat)
/// @param obb_points the 8 points of the obb.
/// @param use convex hull or not.
///
/// todo named parameters: max iterations
/// todo named parameters: max iterations, population size, tolerance.
template <typename Point>
void find_obb(std::vector<Point>& points, std::vector<Point>& obb_points, bool use_ch)
{
@ -175,8 +174,8 @@ void find_obb(std::vector<Point>& points, std::vector<Point>& obb_points, bool u
CGAL_assertion(obb_points.size() == 8);
typedef Eigen::MatrixXf MatrixXf; // using eigen internally
MatrixXf points_mat;
typedef Eigen::MatrixXd MatrixXd; // using eigen internally
MatrixXd points_mat;
// get the ch3
if(use_ch)
@ -196,11 +195,11 @@ void find_obb(std::vector<Point>& points, std::vector<Point>& obb_points, bool u
fill_matrix(points, points_mat);
}
MatrixXf R(3, 3);
MatrixXd R(3, 3);
std::size_t max_generations = 100;
CGAL::Optimal_bounding_box::evolution(R, points_mat, max_generations);
MatrixXf obb(8, 3);
MatrixXd obb(8, 3);
CGAL::Optimal_bounding_box::post_processing(points_mat, R, obb);
// matrix -> vector
@ -209,14 +208,12 @@ void find_obb(std::vector<Point>& points, std::vector<Point>& obb_points, bool u
Point p(obb(i, 0), obb(i, 1), obb(i, 2));
obb_points[i] = p;
}
}
// it is called after post processing
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;
@ -230,7 +227,6 @@ void matrix_to_mesh_and_draw(Matrix& data_points, std::string filename)
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);
@ -238,7 +234,6 @@ void matrix_to_mesh_and_draw(Matrix& data_points, std::string filename)
std::ofstream out(filename);
out << mesh;
out.close();
}

View File

@ -154,9 +154,7 @@ void nelder_mead(std::vector<Matrix>& simplex, Matrix& points, std::size_t nb_it
CGAL_assertion(simplex.size() == 4); // tetrahedron
} // iterations
}
struct Random_int_generator
@ -224,15 +222,12 @@ void genetic_algorithm(Population<Simplex>& pop, Simplex& points)
for(std::size_t i = 0; i < ids4.size(); ++i)
group4[i] = pop[ids4[i]];
/*
#ifdef OBB_DEBUG
check_det(group1);
check_det(group2);
check_det(group3);
check_det(group4);
*/
#endif
// crossover I
Population<Simplex> offspringsA(size_first_group);
@ -259,14 +254,11 @@ void genetic_algorithm(Population<Simplex>& pop, Simplex& points)
offspringsA[i] = offspring;
}
/*
#ifdef OBB_DEBUG
std::cout << "offspringsA: \n" ;
check_det(offspringsA);
std::cin.get();
*/
#endif
// crossover II
Population<Simplex> offspringsB(size_second_group);
@ -294,21 +286,15 @@ void genetic_algorithm(Population<Simplex>& pop, Simplex& points)
offspringsB[i] = offspring;
}
/*
#ifdef OBB_DEBUG
std::cout << "offspringsB: \n" ;
check_det(offspringsB);
std::cin.get();
*/
#endif
CGAL_assertion(offspringsA.size() == size_first_group);
CGAL_assertion(offspringsB.size() == size_second_group);
CGAL_assertion(offspringsA.size() + offspringsB.size() == pop.size());
// next generatrion
for(std::size_t i = 0; i < size_first_group; ++i)
{
@ -320,10 +306,9 @@ void genetic_algorithm(Population<Simplex>& pop, Simplex& points)
pop[size_first_group + i] = offspringsB[i];
}
}
#ifdef OBB_DEBUG
template <typename Simplex>
void check_det(Population<Simplex>& pop)
{
@ -331,20 +316,12 @@ void check_det(Population<Simplex>& pop)
{
for(int j = 0; j < 4; ++j)
{
auto A = pop[i][j];
auto A = pop[i][j]; // Simplex
std::cout << A.determinant() << std::endl;
}
}
}
#endif
} } // end namespaces

View File

@ -129,11 +129,6 @@ private:
}
}
Individual somebody; //temp
CGAL::Random random_generator;
std::size_t n;
std::vector<Individual> pop;

View File

@ -11,7 +11,7 @@
#include <Eigen/Dense>
typedef Eigen::MatrixXf MatrixXf;
typedef Eigen::MatrixXd MatrixXd;
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
@ -51,43 +51,49 @@ void gather_mesh_points(SurfaceMesh& mesh, std::vector<Point>& points)
points.push_back(get(pmap, v));
}
template <typename Point>
double calculate_volume(std::vector<Point> points)
{
CGAL::Bbox_3 bbox;
bbox = bbox_3(points.begin(), points.end());
K::Iso_cuboid_3 ic(bbox);
return ic.volume();
}
void test_population()
{
CGAL::Optimal_bounding_box::Population<MatrixXf> pop(5);
CGAL::Optimal_bounding_box::Population<MatrixXd> pop(5);
//pop.show_population();
CGAL_assertion(pop.size() == 5);
}
void test_nelder_mead()
{
MatrixXf data_points(4, 3);
MatrixXd data_points(4, 3);
data_points << 0.866802, 0.740808, 0.895304,
0.912651, 0.761565, 0.160330,
0.093661, 0.892578, 0.737412,
0.166461, 0.149912, 0.364944;
// one simplex
std::vector<MatrixXf> simplex(4);
std::vector<MatrixXd> simplex(4);
MatrixXf v0(3, 3);
MatrixXd v0(3, 3);
v0 << -0.2192721, 0.2792986, -0.9348326,
-0.7772152, -0.6292092, -0.0056861,
-0.5897934, 0.7253193, 0.3550431;
MatrixXf v1(3, 3);
MatrixXd v1(3, 3);
v1 << -0.588443, 0.807140, -0.047542,
-0.786228, -0.584933, -0.199246,
-0.188629, -0.079867, 0.978795;
MatrixXf v2(3, 3);
MatrixXd v2(3, 3);
v2 << -0.277970, 0.953559, 0.116010,
-0.567497, -0.065576, -0.820760,
-0.775035, -0.293982, 0.559370;
MatrixXf v3(3, 3);
MatrixXd v3(3, 3);
v3 << -0.32657, -0.60013, -0.73020,
-0.20022, -0.71110, 0.67398,
-0.92372, 0.36630, 0.11207;
@ -102,7 +108,7 @@ void test_nelder_mead()
double epsilon = 1e-5;
MatrixXf v0_new = simplex[0];
MatrixXd 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));
@ -113,7 +119,7 @@ void test_nelder_mead()
CGAL_assertion(assert_doubles(v0_new(2,1), 0.5111260, epsilon));
CGAL_assertion(assert_doubles(v0_new(2,2), 0.84094, epsilon));
MatrixXf v1_new = simplex[1];
MatrixXd 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));
@ -124,7 +130,7 @@ void test_nelder_mead()
CGAL_assertion(assert_doubles(v1_new(2,1), 0.338040, epsilon));
CGAL_assertion(assert_doubles(v1_new(2,2), 0.937987, epsilon));
MatrixXf v2_new = simplex[2];
MatrixXd 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));
@ -135,7 +141,7 @@ void test_nelder_mead()
CGAL_assertion(assert_doubles(v2_new(2,1), 0.334057, epsilon));
CGAL_assertion(assert_doubles(v2_new(2,2), 0.941423, epsilon));
MatrixXf v3_new = simplex[3];
MatrixXd 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));
@ -145,25 +151,24 @@ void test_nelder_mead()
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()
{
MatrixXf data_points(4, 3);
MatrixXd data_points(4, 3);
data_points << 0.866802, 0.740808, 0.895304,
0.912651, 0.761565, 0.160330,
0.093661, 0.892578, 0.737412,
0.166461, 0.149912, 0.364944;
CGAL::Optimal_bounding_box::Population<MatrixXf> pop(5);
CGAL::Optimal_bounding_box::Population<MatrixXd> pop(5);
CGAL::Optimal_bounding_box::genetic_algorithm(pop, data_points);
CGAL_assertion(pop.size() == 5);
}
void test_random_unit_tetra()
{
MatrixXf data_points(4, 3); // points on their convex hull
MatrixXd data_points(4, 3); // points on their convex hull
data_points << 0.866802, 0.740808, 0.895304,
0.912651, 0.761565, 0.160330,
0.093661, 0.892578, 0.737412,
@ -180,17 +185,18 @@ void test_random_unit_tetra()
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 OBB_DEBUG_TEST
std::ofstream out("data/random_unit_tetra.off");
out << mesh;
out.close();
#endif
MatrixXf R(3, 3);
MatrixXd R(3, 3);
std::size_t generations = 10;
CGAL::Optimal_bounding_box::evolution(R, data_points, generations);
std::cout << R << std::endl;
double epsilon = 1e-5;
double epsilon = 1e-3;
CGAL_assertion(assert_doubles(R.determinant(), 1, epsilon));
CGAL_assertion(assert_doubles(R(0,0), -0.25791, epsilon));
CGAL_assertion(assert_doubles(R(0,1), 0.796512, epsilon));
@ -202,9 +208,12 @@ void test_random_unit_tetra()
CGAL_assertion(assert_doubles(R(2,1), 0.512847, epsilon));
CGAL_assertion(assert_doubles(R(2,2), 0.836992, epsilon));
MatrixXf obb(8, 3);
#ifdef OBB_DEBUG_TEST
// postprocessing
MatrixXd 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)
@ -217,21 +226,21 @@ void test_reference_tetrahedron(const char* fname)
}
// points in a matrix
MatrixXf points;
MatrixXd points;
sm_to_matrix(mesh, points);
MatrixXf R(3, 3);
MatrixXd R(3, 3);
std::size_t generations = 10;
CGAL::Optimal_bounding_box::evolution(R, points, generations);
double epsilon = 1e-5;
std::cout << R << std::endl;
CGAL_assertion(assert_doubles(R.determinant(), 1, epsilon));
#ifdef OBB_DEBUG_TEST
// postprocessing
MatrixXf obb(8, 3);
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)
@ -244,14 +253,13 @@ void test_long_tetrahedron(std::string fname)
}
// points in a matrix
MatrixXf points;
MatrixXd points;
sm_to_matrix(mesh, points);
MatrixXf R(3, 3);
MatrixXd R(3, 3);
std::size_t generations = 10;
CGAL::Optimal_bounding_box::evolution(R, points, generations);
double epsilon = 1e-5;
std::cout << R << std::endl;
double epsilon = 1e-3;
CGAL_assertion(assert_doubles(R.determinant(), 1, epsilon));
CGAL_assertion(assert_doubles(R(0,0), -1, epsilon));
CGAL_assertion(assert_doubles(R(0,1), 0, epsilon));
@ -265,44 +273,12 @@ void test_long_tetrahedron(std::string fname)
assert_doubles(R(1,2), -0.707106, epsilon));
CGAL_assertion(assert_doubles(R(2,2), 0.707107, epsilon));
#ifdef OBB_DEBUG_TEST
// postprocessing
MatrixXf obb(8, 3);
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");
}
void rotate_tetrahedron(const char* fname, const char* Rname)
{
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";
exit(1);
}
MatrixXf R(3, 3);
std::ifstream input_R(Rname);
double x, y, z;
std::size_t i = 0;
while (input_R >> x >> y >> z) // not safe, but it's going away
{
R(i, 0) = x;
R(i, 1) = y;
R(i, 2) = z;
++i;
}
std:: cout << "det(benchmark)= " << R.determinant() << std::endl;
// points in a matrix
MatrixXf points;
sm_to_matrix(mesh, points);
// postprocessing
MatrixXf obb(8, 3);
CGAL::Optimal_bounding_box::post_processing(points, R, obb);
CGAL::Optimal_bounding_box::matrix_to_mesh_and_draw(obb, "data/inverse_rotated.off");
#endif
}
void test_find_obb(std::string fname)
@ -321,7 +297,8 @@ void test_find_obb(std::string fname)
std::vector<K::Point_3> obb_points;
CGAL::Optimal_bounding_box::find_obb(sm_points, obb_points, true);
double epsilon = 1e-4;
double epsilon = 1e-3;
/*
CGAL_assertion(assert_doubles(obb_points[0].x(), 1.01752, epsilon));
CGAL_assertion(assert_doubles(obb_points[0].y(), 0.437675, epsilon));
CGAL_assertion(assert_doubles(obb_points[0].z(), 0.382697, epsilon));
@ -346,8 +323,12 @@ void test_find_obb(std::string fname)
CGAL_assertion(assert_doubles(obb_points[7].x(), -0.0732647, epsilon));
CGAL_assertion(assert_doubles(obb_points[7].y(), 0.836134, epsilon));
CGAL_assertion(assert_doubles(obb_points[7].z(), 0.733927, epsilon));
*/
/* debug
double vol = calculate_volume(obb_points);
CGAL_assertion(assert_doubles(vol, 0.883371, epsilon));
#ifdef OBB_DEBUG_TEST
for(int i = 0; i < 8; ++i)
{
std::cout << obb_points[i].x() << " " << obb_points[i].y() << " " << obb_points[i].z() << "\n" ;
@ -359,27 +340,65 @@ void test_find_obb(std::string fname)
std::ofstream out("data/obb_result.off");
out << result_mesh;
out.close();
*/
#endif
}
void
bench(const char* fname)
{
std::vector<K::Point_3> sm_points, obb_points;
std::ifstream in(fname);
//double x,y,z;
K::Point_3 p;
int i = 0;
while(in >> p){
if(i % 2 == 0)
{
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;
std::cout.precision(17);
CGAL::Optimal_bounding_box::find_obb(sm_points, obb_points, true);
/*
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/octopus_result.off");
out << mesh;
out.close();
}
int main()
int main(int argc, char* argv[])
{
/*
test_population();
test_nelder_mead();
test_genetic_algorithm();
test_random_unit_tetra();
test_reference_tetrahedron("data/reference_tetrahedron.off");
test_long_tetrahedron("data/long_tetrahedron.off");
rotate_tetrahedron("data/random_unit_tetra.off", "data/rotation.dat");
test_find_obb("data/random_unit_tetra.off");
*/
bench(argv[1]);
return 0;
}

View File

@ -139,7 +139,6 @@ void Create_obb_mesh_plugin::gather_mesh_points(std::vector<Point_3>& points)
}
}
void Create_obb_mesh_plugin::obb()
{
// gather point coordinates