main.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>

#include <CGAL/Octree.h>
#include <CGAL/Orthtree_traits_face_graph.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
#include <CGAL/Polygon_mesh_processing/border.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/autorefinement.h>
#include <CGAL/Random.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <boost/lexical_cast.hpp>
#include <CGAL/IO/polygon_soup_io.h>

#include <CGAL/boost/graph/properties.h>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>

#include <vector>

#include <CGAL/Polygon_mesh_processing/repair_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_mesh_to_polygon_soup.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/Polygon_mesh_processing/repair_self_intersections.h>

using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using Mesh = CGAL::Surface_mesh<K::Point_3>;

using OTraits = CGAL::Orthtree_traits_face_graph<Mesh, Mesh::Property_map<Mesh::Vertex_index, K::Point_3>>;

using namespace std;
namespace PMP = CGAL::Polygon_mesh_processing;
typedef boost::graph_traits<Mesh>::face_descriptor face_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor halfedge_descriptor;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;

Mesh makeHexahedron(K::Point_3 min_coord, K::Point_3 max_coord)
{
    // Declaring hexahedron to return
    Mesh hexahedron;

    // Declaring x, y, and z tmp coordinates
    auto x = min_coord.x(), y = min_coord.y(), z = min_coord.z();
    array<K::Point_3, 8> points;

    // Construct the hexhedron
    for (int i = 0; i < 8; i++)
    {
        // Change x coordinate
        if (i != 0 && (i + 1) % 2 == 0)
            x = x == min_coord.x() ? max_coord.x() : min_coord.x();

        // Change y coordinate
        if (i != 0 && i % 2 == 0)
            y = y == min_coord.y() ? max_coord.y() : min_coord.y();

        // Change z coordinate
        if (i != 0 && i % 4 == 0)
            z = max_coord.z();

        points[i] = K::Point_3(x, y, z);
    }

    // Create hexahedron and triangulate its faces
    CGAL::make_hexahedron(points[0], points[1], points[2], points[3], points[4], points[5], points[6], points[7], hexahedron);
    PMP::triangulate_faces(hexahedron);

    return hexahedron;
}

class Split_predicate_min_volume
{

    double m_min_vol;
    std::shared_ptr<CGAL::Side_of_triangle_mesh<Mesh, K>> m_side_of_ptr;

    double estimate_volume(const std::vector<K::Point_3> &samples, double cube_vol) const
    {
        std::size_t inside = 0;
        for (const K::Point_3 &p : samples)
            if ((*m_side_of_ptr)(p) == CGAL::ON_BOUNDED_SIDE)
                ++inside;
        return cube_vol / samples.size() * inside;
    }

    void fill_bbox(const OTraits::Bbox_d &bb, std::vector<K::Point_3> &samples, std::size_t nb_samples) const
    {
        samples.reserve(nb_samples);

        CGAL::Random &rg = CGAL::get_default_random();
        auto grid_dx = bb.xmax() - bb.xmin();
        auto grid_dy = bb.ymax() - bb.ymin();
        auto grid_dz = bb.zmax() - bb.zmin();

        for (std::size_t i = 0; i < nb_samples; ++i)
        {
            samples.push_back(
                K::Point_3(bb.xmin() + rg.get_double() * grid_dx,
                           bb.ymin() + rg.get_double() * grid_dy,
                           bb.zmin() + rg.get_double() * grid_dz));
        }
    }

public:
    Split_predicate_min_volume() {}

    Split_predicate_min_volume(const double &mv, const Mesh &tmesh)
        : m_min_vol(mv), m_side_of_ptr(new CGAL::Side_of_triangle_mesh<Mesh, K>(tmesh))
    {
    }

    /*!
      \brief returns `true` if `ni` should be split, `false` otherwise.
     */
    template <typename NodeIndex, typename Tree>
    bool operator()(NodeIndex ni, const Tree& tree) const
    {
        OTraits::Bbox_d bb = tree.bbox(ni);
        std::vector<K::Point_3> samples;
        fill_bbox(bb, samples, 100);
        double vol = estimate_volume(samples, CGAL::to_double(bb.volume()));
        const_cast<Tree&>(tree).template property<bool>("is_inside").value()[ni] = (vol!=0);

        if (tree.data(ni).empty())
            return false;



        return (vol > m_min_vol);
    }
};

using Octree = CGAL::Orthtree<OTraits>;

Mesh recreate_mesh(Octree &ot)
{
    int i = 0;
    Mesh result, to_increment;

    auto is_inside = ot.property<bool>("is_inside").value();

    for (Octree::Node_index node : ot.traverse(CGAL::Orthtrees::Leaves_traversal<Octree>(ot)))
    {

        if (!ot.is_leaf(node) || (/* ot.data(node).empty() && */ !is_inside[node]))
            continue;

        auto bb = ot.bbox(node);

        Mesh tmp;
        Mesh mesh = makeHexahedron(K::Point_3(bb.xmin(), bb.ymin(), bb.zmin()), K::Point_3(bb.xmax(), bb.ymax(), bb.zmax()));

        // PMP::corefine_and_compute_union(tmp, mesh, result, CGAL::parameters::throw_on_self_intersection(true));
        result.join(mesh);

        if (i % 100 == 0)
        {
            cout << "i = " << i << endl;
        }

        i++;
    }

    return result;
}

int main(int argc, char **argv)
{
    const std::string filename = "PipedPBase.stl";
    const double mv = argc > 2 ? atof(argv[2]) : 20;
    Mesh mesh;
    if (!CGAL::IO::read_polygon_mesh(filename, mesh))
    {
        std::cerr << "Error: cannot read file" << std::endl;
        return EXIT_FAILURE;
    }

    Octree tree(mesh, mesh.points());
    tree.add_property<bool>("is_inside",false);
    Split_predicate_min_volume sp(mv, mesh);

    cout << "REFINING THE TREE..." << endl;
    tree.refine(sp);
    cout << "TREE HAS BEEN REFINED SUCCESSFULLY !" << endl;

    cout << "RECREATING THE MESH..." << endl;
    Mesh mesh2 = recreate_mesh(tree);
    cout << "THE MESH HAS BEEN RECREATED SUCCESSFULLY !" << endl;

    std::vector<K::Point_3> points;
    std::vector<std::vector<std::size_t>> polygons;

    Mesh mesh3, stl_to_rewrite;

    // Repair the mesh to remove self intersections
    PMP::polygon_mesh_to_polygon_soup(mesh2, points, polygons);
    PMP::autorefine_triangle_soup(points, polygons);

    CGAL::IO::write_polygon_soup("output/result_before_repair.stl", points, polygons);

    PMP::repair_polygon_soup(points, polygons, CGAL::parameters::erase_all_duplicates(true));

    // Write the repaired soup into result.stl
    CGAL::IO::write_polygon_soup("output/result.stl", points, polygons);

    //~ // Read the repaired mesh (written just above)
    //~ CGAL::IO::read_polygon_mesh("output/result.stl", stl_to_rewrite);

    //~ // Write the mesh that we just read into a new STL file --> THIS DESTROYS THE MESH
    //~ CGAL::IO::write_polygon_mesh("output/rewrite_result.stl", stl_to_rewrite);

    return 0;
}