Polygonal Surface Reconstruction

[petrasvestartas]

Issue Details

I am running polygonal surface reconstruction, the very first example, and I get just bunch of faces of a box, that is not the same as the input shape:

https://doc.cgal.org/latest/Polygonal_surface_reconstruction/index.html

There are these little planes:

Source Code

two_intersecting_boxes.zip file:

I check if the reconstruction failed or not, and it says there is no failure:

    // Reconstruct surface using mixed-integer programming solver
    if (!algo.reconstruct<MIP_Solver>(model))
    {
        std::cout << " Failed:...";
        // std::cerr << " Failed: " << algo.error_message() << std::endl;
    }else{
        std::cout << " Done. " << std::endl;
    }

I check if the reconstruction failed or not, and it says there is no failure:

    // Reconstruct surface using mixed-integer programming solver
    if (!algo.reconstruct<MIP_Solver>(model))
    {
        std::cout << " Failed:...";
        // std::cerr << " Failed: " << algo.error_message() << std::endl;
    }else{
        std::cout << " Done. " << std::endl;
    }

#include "polygonal_surface_reconstruction.h"

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/property_map.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Shape_detection/Efficient_RANSAC.h>
#include <CGAL/Polygonal_surface_reconstruction.h>

#define CGAL_USE_SCIP
#ifdef CGAL_USE_SCIP
#include <CGAL/SCIP_mixed_integer_program_traits.h>
typedef CGAL::SCIP_mixed_integer_program_traits<double> MIP_Solver;
#elif defined(CGAL_USE_GLPK)
#include <CGAL/GLPK_mixed_integer_program_traits.h>
typedef CGAL::GLPK_mixed_integer_program_traits<double> MIP_Solver;
#endif

#if defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)
#include <CGAL/Timer.h>
#include <fstream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::Point_3 Point;
typedef Kernel::Vector_3 Vector;

typedef boost::tuple<Point, Vector, int> PNI;
typedef std::vector<PNI> Point_vector;
typedef CGAL::Nth_of_tuple_property_map<0, PNI> Point_map;
typedef CGAL::Nth_of_tuple_property_map<1, PNI> Normal_map;
typedef CGAL::Nth_of_tuple_property_map<2, PNI> Plane_index_map;

typedef CGAL::Shape_detection::Efficient_RANSAC_traits<Kernel, Point_vector, Point_map, Normal_map> Traits;
typedef CGAL::Shape_detection::Efficient_RANSAC<Traits> Efficient_ransac;
typedef CGAL::Shape_detection::Plane<Traits> Plane;
typedef CGAL::Shape_detection::Point_to_shape_index_map<Traits> Point_to_shape_index_map;
typedef CGAL::Polygonal_surface_reconstruction<Kernel> Polygonal_surface_reconstruction;
typedef CGAL::Surface_mesh<Point> Surface_mesh;

std::tuple<compas::RowMatrixXd, compas::RowMatrixXi>
polygonal_surface_reconstruction_ransac(
    Eigen::Ref<const compas::RowMatrixXd> &P,
    Eigen::Ref<const compas::RowMatrixXd> &N)
{
    // Convert Eigen matrix to a vector of CGAL points.
    int p = P.rows();
    Point_vector points;

    for (int i = 0; i < p; i++)
    {
        points.push_back(boost::tuple<Point, Vector, int>(Point(P.row(i)[0], P.row(i)[1], P.row(i)[2]),
                                                          Vector(N.row(i)[0], N.row(i)[1], N.row(i)[2]), -1));
    }

    std::cout << " Done. " << points.size() << std::endl;

    Efficient_ransac ransac;
    ransac.set_input(points);
    ransac.add_shape_factory<Plane>();
    std::cout << "Extracting planes...";
    ransac.detect();
    Efficient_ransac::Plane_range planes = ransac.planes();
    std::size_t num_planes = planes.size();
    std::cout << " Done. " << num_planes << " planes extracted. " << std::endl;

    Point_to_shape_index_map shape_index_map(points, planes);
    for (std::size_t i = 0; i < points.size(); ++i)
    {
        int plane_index = get(shape_index_map, i);
        std::cout << "Point " << i << " is on plane " << plane_index << std::endl;
        points[i].get<2>() = plane_index;
    }

    std::cout << "Generating candidate faces...";
    Polygonal_surface_reconstruction algo(
        points,
        Point_map(),
        Normal_map(),
        Plane_index_map());

    Surface_mesh model;
    std::cout << "Reconstructing...";

    if (!algo.reconstruct<MIP_Solver>(model))
    {
        std::cerr << " Failed: " << algo.error_message() << std::endl;
    }

    Surface_mesh candidate_faces;
    algo.output_candidate_faces(candidate_faces);

    // Iterate faces of Surface_mesh and print vertices
    for (auto f : candidate_faces.faces())
    {
        std::cout << "face " << f << ": ";
        for (auto v : vertices_around_face(candidate_faces.halfedge(f), candidate_faces))
        {
            std::cout << v << " ";
        }
        std::cout << std::endl;
    }

    std::size_t v = candidate_faces.number_of_vertices();
    std::size_t f = candidate_faces.number_of_faces();

    compas::RowMatrixXd V(v, 3);
    compas::RowMatrixXi F(f, 4);

    compas::Mesh::Property_map<compas::Mesh::Vertex_index, compas::Kernel::Point_3> location = candidate_faces.points();

    for (compas::Mesh::Vertex_index vd : candidate_faces.vertices())
    {
        V(vd, 0) = (double)location[vd][0];
        V(vd, 1) = (double)location[vd][1];
        V(vd, 2) = (double)location[vd][2];
    }

    for (compas::Mesh::Face_index fd : candidate_faces.faces())
    {
        int i = 0;
        for (compas::Mesh::Vertex_index vd : vertices_around_face(candidate_faces.halfedge(fd), candidate_faces))
        {
            F(fd, i) = (int)vd;
            i++;
        }
    }

    std::tuple<compas::RowMatrixXd, compas::RowMatrixXi> result = std::make_tuple(V, F);

    return result;
}

void init_polygonal_surface_reconstruction(pybind11::module &m)
{
    pybind11::module submodule = m.def_submodule("polygonal_surface_reconstruction");

    submodule.def(
        "polygonal_surface_reconstruction_ransac",
        &polygonal_surface_reconstruction_ransac,
        pybind11::arg("P").noconvert(),
        pybind11::arg("N").noconvert());
};
#endif  // defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)

[petrasvestartas]

it was model: