Unexpected result when merging several meshes into one using corefine_and_compute_union()

[AlbertDeTerre]

Issue

Hello,

I have numerous different meshes that I would like to merge to create a single unified mesh.
To do so, I'm using the corefine_and_compute_union function.

However, the issue is that once all the meshes are assembled, the result doesn't form a single cohesive mesh but rather a mesh with separations between each "child mesh" (see image below).

Is there a way to avoid that ?

Picture of the resulting mesh

My code

Mesh final_mesh, tmp;
for (Mesh mesh : meshes){
    PMP::corefine_and_compute_union(tmp, mesh, final_mesh);
}

Additional informations

• Using Exact_predicates_inexact_constructions_kernel
• Using Surface_mesh
• I tried to compute union like this: PMP::corefine_and_compute_union(tmp, to_append, final_mesh) but it ended with a self intersection error
• I tried using copy_face_graph and it ended with the same result.

[sloriot]

iteration variable is mesh and it is not used in the function call within the loop. Could it be the issue?

[AlbertDeTerre]

Sorry I just retyped the code instead of copy/paste it for some reasons so I made a mistake while retyping. Fixed that but it's not the issue.

[AlbertDeTerre]

It seems like all the child meshes are not fully merged but just fitted together like bricks. Is there a way to avoid that @sloriot ?

[sloriot]

My guess is that you have cubic cells that you are trying to union. But at some intermediate step you have two cubic cells sharing a vertex and it becomes non-manifold. While it could be only "temporary" as more cell will be added, I guess you could also have that in the output, right?

[AlbertDeTerre]

Yes it is what I think too.
Some of the small cubes share the same edges and vertices and are not deleted when doing the union so this end having too much vertices in the file (I could see that because the STL file of generated mesh is abnormally large).
What could I do to solve that ?

Here is the code if you want to try it yourself:

#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/boost/graph/copy_face_graph.h>
#include <CGAL/Random.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <CGAL/IO/polygon_soup_io.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;

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
    double 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();
        double grid_dx = bb.xmax() - bb.xmin();
        double grid_dy = bb.ymax() - bb.ymin();
        double 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
    {
        if (tree.data(ni).empty())
            return false;

        OTraits::Bbox_d bb = tree.bbox(ni);
        std::vector<K::Point_3> samples;
        fill_bbox(bb, samples, 100);

        return (estimate_volume(samples, bb.volume()) > m_min_vol);
    }
};

using Octree = CGAL::Orthtree<OTraits>;

Mesh dump_as_polylines(const Octree &ot)
{
    Mesh output, tmp;
    for (Octree::Node_index node : ot.traverse(CGAL::Orthtrees::Leaves_traversal<Octree>(ot)))
    {
        
        if (!ot.is_leaf(node) || ot.data(node).empty())
            continue;

        auto bb = ot.bbox(node);

        tmp = makeHexahedron(K::Point_3(bb.xmin(), bb.ymin(), bb.zmin()), K::Point_3(bb.xmax(), bb.ymax(), bb.zmax()));
        CGAL::copy_face_graph(tmp, output);
    }

    return output;
}

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());
    Split_predicate_min_volume sp(mv, mesh);
    tree.refine(sp);

    Mesh mesh2 = dump_as_polylines(tree);
    
    if (!CGAL::IO::write_STL("../../output/output.stl", mesh3))
    {
        cout << "FAIL" << endl;
    }
    cout << "AFTER" << endl;
}

[AlbertDeTerre]

Update !!

I found that doing this:

Mesh final_mesh, tmp;
for (Mesh mesh : meshes){
    PMP::corefine_and_compute_union(tmp, mesh, final_mesh);
}

creates a mesh that has self-intersections without throwing any error.

So, to solve this I tried to do it like this:

Surface_mesh result;
for (auto mesh : meshes)
{
    Surface_mesh output;
    PMP::corefine_and_compute_union(result, mesh, output);

    result = output;
}
CGAL::IO::write_STL("./result.stl", result)

This threw a core dumped error due to self intersections...

What I found

By trying new things I found that it doesn't throw any error if the vector that contains all the meshes is sorted.

For example I made an algorithm that creates a 3D stairs-like mesh:

vector<Surface_mesh> meshes;

double offset = 2;

for (int i = 0; i < 20; i++)
{
    for (int j = 0; j < i + 1; j++)
    {
        for (int k = 0; k < j + 1; k++)
        {
            meshes.push_back(makeHexahedron(Point(offset * i, offset * j, offset * k), Point(offset * (i + 1), offset * (j + 1), offset * (k + 1))));
        }
    }
}

If I shuffle the vector before merging all the mesh, it throws a core dumped error because of self intersections when trying to merge the meshes to create the final mesh.

BUT, if the vector isn't shuffled before merging all the mesh, the final mesh is well created and without any self-intersections.

This was used to merge all the meshes:

Surface_mesh result;
for (auto mesh : meshes)
{
    Surface_mesh output;
    PMP::corefine_and_compute_union(result, mesh, output);

    result = output;
}
CGAL::IO::write_STL("./result.stl", result)

Conclusion

So, to conclude, I think the issue happened in the code I gave because the meshes that I merge are not in a right order... So how could I solve that ?

Full code

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

#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>

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

#include <algorithm>
#include <random>

namespace PMP = CGAL::Polygon_mesh_processing;

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_3 Point;
typedef CGAL::Surface_mesh<Point> Surface_mesh;
typedef boost::graph_traits<Surface_mesh>::face_descriptor face_descriptor;

using namespace std;

Surface_mesh makeHexahedron(Point min_coord, Point max_coord)
{
    // Declaring hexahedron to return
    Surface_mesh hexahedron;

    // Declaring x, y, and z tmp coordinates
    auto x = min_coord.x(), y = min_coord.y(), z = min_coord.z();
    array<Point, 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] = Point(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;
}

int main()
{
    vector<Surface_mesh> meshes;

    double offset = 2;

    for (int i = 0; i < 20; i++)
    {
        for (int j = 0; j < i + 1; j++)
        {
            for (int k = 0; k < j + 1; k++)
            {
                meshes.push_back(makeHexahedron(Point(offset * i, offset * j, offset * k), Point(offset * (i + 1), offset * (j + 1), offset * (k + 1))));
            }
        }
    }

    // IF ENABLED --> Self intersections error will be thrown 
    // auto rng = std::default_random_engine{};
    // std::shuffle(std::begin(meshes), std::end(meshes), rng);

    Surface_mesh result, tmp;

    int i = 0;
    for (auto mesh : meshes)
    {

        Surface_mesh output;
        PMP::corefine_and_compute_union(tmp, mesh, output);

        tmp = output;

        cout << "Does self intersect [ " << i << "]" << PMP::does_self_intersect(output) << endl;
        i++;
    }

    cout << CGAL::IO::write_STL("./result1.stl", tmp);

    return 0;
}

[sloriot]

The order change prevent non-manifold features to appear, which is why you could compute the whole union. Unfortunately I don't see an easy solution for scheduling operations. What I'm not sure is if you really need to use that function. To me it seems that you are only doing union of iso-cuboids with intersection being of volume 0. There might be a most efficient way to handle that. From the top of my head, I'd use autorefine_triangle_soup() if isocuboids are not conformal, followed by a call to repair_polygon_soup() with the option to remove a face if it is present twice.

[AlbertDeTerre]

Thanks for your answer @sloriot !
I tried using autorefine_triangle_soup(), then repair_polygon_soup() and the result was quite ok but it seems that some triangles are missing that cause holes in the mesh (see image below).

My question

Is there a way to fill these holes ?

What I already tried

• Hole filling using this:

    unsigned int nb_holes = 0;
    vector<halfedge_descriptor> border_cycles;

    PMP::extract_boundary_cycles(mesh, back_inserter(border_cycles));

    for (halfedge_descriptor h : border_cycles)
    {

        vector<face_descriptor> patch_facets;
        vector<vertex_descriptor> patch_vertices;

        bool success = std::get<0>(PMP::triangulate_refine_and_fair_hole(mesh,
                                                                         h,
                                                                         CGAL::parameters::face_output_iterator(back_inserter(patch_facets))
                                                                             .vertex_output_iterator(back_inserter(patch_vertices))));
        ++nb_holes;
    }

It didn't do anything.

• Triangulate faces after autorefine and repair -> Didn't work

Image

[AlbertDeTerre]

Ok I just found that I have to write the mesh using write_polygon_soup instead of transforming soup into mesh then write the mesh.
Is there a reason to that @sloriot ?

[AlbertDeTerre]

I also found that if I write the mesh and then read the file again, it destroys my repaired mesh... why ?

(see images below)

Code

    Mesh mesh, stl_to_rewrite;

    // Repair the mesh
    PMP::polygon_mesh_to_polygon_soup(mesh, points, polygons);
    PMP::autorefine_triangle_soup(points, polygons);
    PMP::repair_polygon_soup(points, polygons, CGAL::parameters::erase_all_duplicates(true));

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

    // Read the STL file containing the repaired mesh and write the new mesh into rewrite_result.stl
    CGAL::IO::read_polygon_mesh("../../output/result.stl", stl_to_rewrite);
    CGAL::IO::write_polygon_mesh("../../output/rewrite_result.stl", stl_to_rewrite);

Images

result.stl

rewrite_result.stl

[sloriot]

If you provide a (small) self-contained example I can compile and run it and have a look at the output. No need to give the code for the first step, you can dump meshes using maximum stream precision (using std::setprecision(17) for example).

[AlbertDeTerre]

Sorry but I'm not sure to understand what I have to give you (I'm french speaker so english is not my first language)... Which part of the code should I give ?

[sloriot]

The one that you think is not working correctly.

[AlbertDeTerre]

The part of the code that I think is not working correctly is this:

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

    // Read the STL file containing the repaired mesh and write the new mesh into rewrite_result.stl
    CGAL::IO::read_polygon_mesh("../../output/result.stl", stl_to_rewrite);
    CGAL::IO::write_polygon_mesh("../../output/rewrite_result.stl", stl_to_rewrite);

It's when I try to read the repaired mesh from the file result.stl and then rewrite it into another file (rewrite_result.stl)

I'm sorry if this is not what you need I'm doing my best ^^'

[AlbertDeTerre]

Full code with the PipedPBase.stl file

File: PipedPBase.zip

#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
    {
        if (tree.data(ni).empty())
            return false;

        OTraits::Bbox_d bb = tree.bbox(ni);
        std::vector<K::Point_3> samples;
        fill_bbox(bb, samples, 100);

        return (estimate_volume(samples, CGAL::to_double(bb.volume())) > m_min_vol);
    }
};

using Octree = CGAL::Orthtree<OTraits>;

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

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

        if (!ot.is_leaf(node) || ot.data(node).empty())
            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));

        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());
    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);
    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;
}

[sloriot]

As I said you don't need corefinement_and_compute_union().
I modified a bit your example and got what I think you expected.

https://gist.github.com/sloriot/50f6c9b1ae9f1d2e0b0a45900703e45a

However, the output contains non-manifold edges so the output, if oriented to be put in a triangle mesh will end up being self-intersecting.

[AlbertDeTerre]

WOW ! Thank's for everything you're a HERO !
This solved my issue, and using join() instead of corefine_and_compute_union() made my program runs a lot faster !
Sorry for bothering you those past days and thanks a lot for your patience :)

[sloriot]

If you want to further improve the thing, the is_inside property for empty leaves could to filled using side_of and the centroid of the bbox of the node. Then the split predicate could be restored as in the previous version.

[sloriot]

And you could also probably extract the boundary faces directly from the octree using navigation with neighbors (a bit more trickier to write I guess)

[AlbertDeTerre]

Thanks, I'll probably try that later if I need better performances but it works already quite well.