Can I use isotropic_remeshing on a non-manifold mesh

[citystrawman]

Please use the following template to help us solving your issue.

Issue Details

I am imitating Polygon_mesh_processing/corefinement_difference_remeshed.cpp to remesh a non-manifold mesh (its background is here)

In general, what I have done is:

1. clip the slope with slide;
2. merge the slope and slide and remove the duplicated (intersecting) vertices; (after merging, the merged geometry becomes non-manifold)
3. get the vertices and select faces that have these vertices, then expand the face by 2 iteration;
4. use isotropic remesh

However after isotropic remesh, the off file has some problem. When I open the off file, it pops out warnings: after clicking OK, the mesh looks weird(I dont know if it is original mesh or it is corrected by meshlab):

and the console output is as follows:

I am not sure what is the reason for this problem: is it because I do not code it correctly, or is it because isotropic_remeshing only supports manifold geometry?

Source Code

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Surface_mesh/Surface_mesh.h>
#include <CGAL/draw_surface_mesh.h>
#include <CGAL/boost/graph/selection.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>

#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <CGAL/Polygon_mesh_processing/clip.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>

#include <iostream>
#include <string>
#include <CGAL/boost/graph/graph_traits_Surface_mesh.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_predicates_exact_constructions_kernel EK;
typedef CGAL::Surface_mesh<K::Point_3> Mesh;
typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;
typedef Mesh::Property_map<vertex_descriptor, EK::Point_3> Exact_point_map;
typedef boost::graph_traits<Mesh>::face_descriptor            face_descriptor;
typedef boost::graph_traits<Mesh>::edge_descriptor            edge_descriptor;
typedef boost::graph_traits<Mesh>::halfedge_descriptor        halfedge_descriptor;

namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = CGAL::parameters;

struct Vector_pmap_wrapper
{
    std::vector<bool>& vect;
    Vector_pmap_wrapper(std::vector<bool>& v) : vect(v) {}
    friend bool get(const Vector_pmap_wrapper& m, face_descriptor f)
    {
        return m.vect[f];
    }
    friend void put(const Vector_pmap_wrapper& m, face_descriptor f, bool b)
    {
        m.vect[f] = b;
    }
};

class Insert_iterator
{
    typedef std::unordered_map<face_descriptor, face_descriptor> Container;
    Container& container;
public:

    Insert_iterator(Container& c)
        : container(c) {}

    Insert_iterator&
        operator=(const std::pair<face_descriptor, face_descriptor>& p)
    {
        container[p.second] = p.first;
        return *this;
    }

    Insert_iterator&
        operator*() { return *this; }

    Insert_iterator
        operator++(int) { return *this; }

};
struct Visitor : public CGAL::Polygon_mesh_processing::Corefinement::Default_visitor<Mesh>
{
    typedef std::unordered_map<face_descriptor, face_descriptor> Container;

    Container& container;
    face_descriptor qfd;
    //std::vector<face_descriptor>& createdNewFace;
    std::vector<vertex_descriptor>& tmNewVertex;
    std::vector<vertex_descriptor>& splitterNewVertex;
    std::vector<halfedge_descriptor>& tmNewHalfEdge;
    std::vector<halfedge_descriptor>& splitterNewHalfEdge;

    Mesh& tm;
    Mesh& splitter;

    Visitor(Container& container, Mesh& tm, Mesh& splitter, std::vector<vertex_descriptor>& tmNewVertex, std::vector<vertex_descriptor>& splitterNewVertex, std::vector<halfedge_descriptor>& tmNewHalfEdge, std::vector<halfedge_descriptor>& splitterNewHalfEdge)
        : container(container), tm(tm), splitter(splitter), tmNewVertex(tmNewVertex), splitterNewVertex(splitterNewVertex),
        tmNewHalfEdge(tmNewHalfEdge), splitterNewHalfEdge(splitterNewHalfEdge)
    {}

    //==========================FACE==========================
    void before_subface_creations(face_descriptor f_old, Mesh&)
    {

    }

    void after_subface_created(face_descriptor f_new, Mesh&)
    {
    }

    void after_subface_creations(Mesh&)
    {
    }

    //==========================VERTEX==========================

    void intersection_point_detected(std::size_t  node_id,
        int  sdim,
        halfedge_descriptor  principal_edge,
        halfedge_descriptor  additional_edge,
        const Mesh& tm1,
        const Mesh& tm2,
        bool  is_target_coplanar,
        bool  is_source_coplanar)
    {
    }

    void new_vertex_added(std::size_t  node_id,
        vertex_descriptor  vh,
        const Mesh& tm)
    {
        std::cout << "node_id: " << node_id << "; vertex_descriptor: " << vh << "; Mesh: " << tm.number_of_vertices() << std::endl;
        if (&this->tm == std::addressof(tm)) {
            tmNewVertex.push_back(vh);
        }
        else if (&this->splitter == std::addressof(tm)) {
            splitterNewVertex.push_back(vh);
        }
    }

    void before_vertex_copy(vertex_descriptor /*v_src*/, const Mesh& /*tm_src*/, Mesh& /*tm_tgt*/) {}

    void after_vertex_copy(vertex_descriptor v_src, const Mesh& tm_src,
        vertex_descriptor  v_tgt, Mesh& tm_tgt)
    {
    }

    //==========================halfedge==========================
    void after_edge_copy(halfedge_descriptor h_old, const Mesh& m1, halfedge_descriptor  h_new, Mesh& m2)
    {
        std::cout << "after_edge_copy called: h_old: " << h_old << ", h_new: " << h_new << std::endl;
    }

    void after_edge_duplicated(halfedge_descriptor h_old, halfedge_descriptor  h_new, Mesh& m)
    {
        std::cout << "after_edge_duplicated called: h_old: " << h_old << ", h_new: " << h_new << std::endl;
    }

    void intersection_edge_copy(halfedge_descriptor  h_old1, const Mesh& tm1,
        halfedge_descriptor  h_old2, const Mesh& tm2,
        halfedge_descriptor  h_new, Mesh& tm_new)
    {
        std::cout << "intersection_edge_copy called: h_old1: " << h_old1 << "h_old2: " << h_old2 << ", h_new : " << h_new << std::endl;
    }

    void add_retriangulation_edge(halfedge_descriptor  h, Mesh& tm)
    {
        std::cout << "add_retriangulation_edge called: h: " << h << std::endl;
    }

    void edge_split(halfedge_descriptor  hnew, Mesh& tm) //新产生的vertex两两相连产生的halfedge
    {
        std::cout << "edge_split called: hnew: " << hnew << ", hnew's incident vertex: " << tm.target(hnew) << std::endl;
        if (&this->tm == std::addressof(tm)) {
            tmNewHalfEdge.push_back(hnew);
        }
        else if (&this->splitter == std::addressof(tm)) {
            splitterNewHalfEdge.push_back(hnew);
        }
    }
};

int main(int argc, char* argv[])
{
    const std::string filename1 = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/slide.off"); //slideRes-flip.off
    const std::string filename2 = (argc > 2) ? argv[2] : CGAL::data_file_path("meshes/slope.off");  //slopeRes.off
    double target_edge_length = (argc > 3) ? std::stod(std::string(argv[3])) : 15;
    unsigned int nb_iter = (argc > 4) ? std::stoi(std::string(argv[4])) : 10;

    Mesh mesh1, mesh2;
    if (!PMP::IO::read_polygon_mesh(filename1, mesh1) || !PMP::IO::read_polygon_mesh(filename2, mesh2))
    {
        std::cerr << "Invalid input." << std::endl;
        return 1;
    }

    if (PMP::does_self_intersect(mesh1) || PMP::does_self_intersect(mesh2))
    {
        std::cerr << "self intersect input." << std::endl;
        return 1;
    }

    Mesh mesh1_clip, mesh2_clip;

    CGAL::copy_face_graph(mesh1, mesh1_clip);
    CGAL::copy_face_graph(mesh2, mesh2_clip);

    std::unordered_map<face_descriptor, face_descriptor> t2q;
    std::vector<vertex_descriptor> tmNewVertex;
    std::vector<vertex_descriptor> splitterNewVertex;
    std::vector<halfedge_descriptor> tmNewHalfEdge;
    std::vector<halfedge_descriptor> splitterNewHalfEdge;

    Visitor v_m1(t2q, mesh1_clip, mesh2_clip, tmNewVertex, splitterNewVertex, tmNewHalfEdge, splitterNewHalfEdge);

    PMP::clip(mesh1_clip, mesh2_clip, params::visitor(v_m1), params::do_not_modify(false));

    CGAL::IO::write_polygon_mesh("slide_clip.off", mesh1_clip, CGAL::parameters::stream_precision(17));

    CGAL::IO::write_polygon_mesh("slope_clip.off", mesh2_clip, CGAL::parameters::stream_precision(17));

    std::unordered_map<vertex_descriptor, vertex_descriptor> v2v;
    std::unordered_map<halfedge_descriptor, halfedge_descriptor> h2h;
    std::unordered_map<face_descriptor, face_descriptor> f2f;

    CGAL::copy_face_graph(mesh1_clip, mesh2_clip, CGAL::parameters::vertex_to_vertex_map(boost::make_assoc_property_map(v2v))
        .halfedge_to_halfedge_output_iterator(std::inserter(h2h, h2h.end()))
        .face_to_face_map(boost::make_assoc_property_map(f2f)));

    mesh1_clip.collect_garbage();
    mesh2_clip.collect_garbage();

    int num_remove = 0;
    std::vector<vertex_descriptor> intersectVertices;
    for each (vertex_descriptor v in tmNewVertex) {
        vertex_descriptor v_merged = v2v[v];
        for each (vertex_descriptor find in mesh2_clip.vertices()) {
            if ((mesh2_clip.point(find).cartesian(1) == mesh2_clip.point(v_merged).cartesian(1)
                && mesh2_clip.point(find).cartesian(2) == mesh2_clip.point(v_merged).cartesian(2)
                && mesh2_clip.point(find).cartesian(0) == mesh2_clip.point(v_merged).cartesian(0)) && find != v_merged)
            {
                for (halfedge_descriptor h : CGAL::halfedges_around_target(find, mesh2_clip)) {
                    set_target(h, v_merged, mesh2_clip);
                }
                mesh2_clip.remove_vertex(find);
                intersectVertices.push_back(v_merged);
                num_remove++;
            }
        }
    }

    CGAL::draw(mesh2_clip);

    CGAL::IO::write_polygon_mesh("slope_slide_clip_merge.off", mesh2_clip, CGAL::parameters::stream_precision(17));

    //仿照corefinement_difference_remeshed尝试对模型进行remesh
    std::vector<face_descriptor> selected_faces;
    std::vector<bool> is_selected(num_faces(mesh2_clip), false);
    for each (vertex_descriptor v in intersectVertices) {
        for (halfedge_descriptor h : halfedges_around_target(v, mesh2_clip)) {
            if (!is_border(h, mesh2_clip)) {
                face_descriptor f = face(h, mesh2_clip);
                if (!is_selected[f])
                {
                    selected_faces.push_back(f);
                    is_selected[f] = true;
                }
            }
        }
    }

    // increase the face selection
    CGAL::expand_face_selection(selected_faces, mesh2_clip, 2,
        Vector_pmap_wrapper(is_selected), std::back_inserter(selected_faces));
    std::cout << selected_faces.size()
        << " faces were selected for the remeshing step\n";

    // remesh the region around the intersection polylines
    PMP::isotropic_remeshing(selected_faces, target_edge_length, mesh2_clip);
    //mesh2_clip.collect_garbage();
    CGAL::IO::write_polygon_mesh("slope_slide_clip_merge_remeshed.off", mesh2_clip, CGAL::parameters::stream_precision(17));

    return 0;
}

Environment

• Operating system (Windows/Mac/Linux, 32/64 bits): windows
• Compiler: VS2022
• Release or debug mode:
• Specific flags used (if any):
• CGAL version:
• Boost version:
• Other libraries versions if used (Eigen, TBB, etc.):

[janetournois]

Indeed, PMP::isotropic_remeshing() can only remesh surface patches with manifold geometry, which belong to a MutableFaceGraph (as mentioned in its doc). If your input mesh has non-manifold geometry, you can remesh your surface patch-per-patch.

However, from your log, it seems that your input mesh is not a valid polygon mesh at all, not only to be remeshed by isotropic_remeshing(). Non-manifold geometry is allowed, but not non-manifold connectivity. If it is the issue, you can use PMP::polygon_soup_to_polygon_mesh() which may duplicate non-manifold vertices or edges to have a manifold connectivity.

[citystrawman]

Indeed, PMP::isotropic_remeshing() can only remesh surface patches with manifold geometry, which belong to a MutableFaceGraph (as mentioned in its doc). If your input mesh has non-manifold geometry, you can remesh your surface patch-per-patch.

However, from your log, it seems that your input mesh is not a valid polygon mesh at all, not only to be remeshed by isotropic_remeshing(). Non-manifold geometry is allowed, but not non-manifold connectivity. If it is the issue, you can use PMP::polygon_soup_to_polygon_mesh() which may duplicate non-manifold vertices or edges to have a manifold connectivity.

May I know what is the difference between Non-manifold geometry and non-manifold connectivity?

[janetournois]

Sure Think of a vertex with two umbrellas sharing this vertex as a tip point (like a sand hourglass). Both the geometry and connectivity are non-manifold. Now if you duplicate this vertex to separate the 2 umbrellas, the geometry is still non-manifold (unless you change the coordinates of one of the duplicates), but now the connectivity is manifold, as 2 manifold surfaces.

PMP::isotropic_remeshing() can handle the double-with-2-tips-umbrella, not the double-with-1-tip-umbrella. However, the tip vertices can be constrained to make sure they remain exactly on both sides.

[citystrawman]

Sure Think of a vertex with two umbrellas sharing this vertex as a tip point (like a sand hourglass). Both the geometry and connectivity are non-manifold. Now if you duplicate this vertex to separate the 2 umbrellas, the geometry is still non-manifold (unless you change the coordinates of one of the duplicates), but now the connectivity is manifold, as 2 manifold surfaces.

PMP::isotropic_remeshing() can handle the double-with-2-tips-umbrella, not the double-with-1-tip-umbrella. However, the tip vertices can be constrained to make sure they remain exactly on both sides.

as you said, the tip vertices can be constrained to make sure they remain exactly on both sides, is this achieved by setting vertex_is_constrained_map in PMP::isotropic_remeshing()?

[janetournois]

exactly

[citystrawman]

exactly

Then here comes my latest issue #8388 : in this issue actually what I did was trying to use vertex_is_constrained_map on some of the tip vertices (I only need some of the vertices because some of the tip vertices are too close and if keep all these vertices the remeshed triangle should be ill shaped), and you can see that vertex_is_constrained_map caused an unwanted effect: the constrained vertices are not at its original position. Is my idea correct or wrong? if it is wrong, then do you have any idea? Thank you!