Assertion failure doing collect_garbage() after removing isolated vertices from a surface mesh. CGAL.

test_remove.cpp

#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>

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

#include <CGAL/Surface_mesh/IO/PLY.h>

#include <CGAL/Polygon_mesh_processing/connected_components.h>
#include <CGAL/Polygon_mesh_processing/repair.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>

namespace PMP = CGAL::Polygon_mesh_processing;

typedef CGAL::Simple_cartesian<float> Kernel_t;
typedef Kernel_t::Point_3  Point_t;

typedef CGAL::Surface_mesh<Point_t> Mesh_t;

typedef Kernel_t::Compare_dihedral_angle_3 Compare_dihedral_angle_3;

template < typename MeshT >
void read_mesh(const std::string& fn, MeshT& mesh) {
    std::ifstream ifs(fn, std::ios::binary);

    if ( !ifs ) {
        std::stringstream ss;
        ss << "Cannot open " << fn;
        throw std::runtime_error(ss.str());
    }

    std::string comments;
    CGAL::IO::read_PLY(ifs, mesh, comments);
}

template < typename MeshT >
void duplicate_non_manifold_vertices(MeshT& mesh) {
    PMP::duplicate_non_manifold_vertices(mesh);
}

template < typename MeshT >
int remove_self_intersections( MeshT& mesh ) {
    typedef typename boost::graph_traits<MeshT>::face_descriptor Face_t;

    std::vector< std::pair< Face_t, Face_t > > intersections;
    PMP::self_intersections<CGAL::Parallel_if_available_tag>(
        faces(mesh), 
        mesh,
        std::back_inserter(intersections)
    );

    const int n_intersections = intersections.size();
    if ( n_intersections > 0 ) {
        for ( const auto& inter_pair : intersections ) {
            mesh.remove_face( inter_pair.first );
            mesh.remove_face( inter_pair.second );
        }

        mesh.collect_garbage();
    }

    return n_intersections;
}

template < typename MeshT >
int remove_isolated_vertices( MeshT& mesh ) {
    const int N = PMP::remove_isolated_vertices(mesh);
    return N;
}

// For computing the connected components.
template < typename MeshT >
struct ConnectionConstraint : public boost::put_get_helper< bool, ConnectionConstraint<MeshT> > {
    typedef typename boost::graph_traits<MeshT>::edge_descriptor EdgeT;

    typedef boost::readable_property_map_tag category;
    typedef bool                             value_type;
    typedef bool                             reference;
    typedef EdgeT                            key_type;

    ConnectionConstraint()
    : m_mesh(nullptr) {}

    ConnectionConstraint( MeshT& mesh, double bound )
    : m_mesh(&mesh), m_bound(bound) {}

    bool operator[] ( EdgeT edge ) const {
        const MeshT& mesh = *m_mesh;
        return m_compare(
            mesh.point( source( edge, mesh ) ),
            mesh.point( target( edge, mesh ) ),
            mesh.point( target( next( halfedge( edge, mesh ), mesh ), mesh ) ),
            mesh.point( target( next( opposite( halfedge( edge, mesh ), mesh), mesh ) ,mesh ) ),
            m_bound ) == CGAL::SMALLER;
    }

    const MeshT* m_mesh;
    Compare_dihedral_angle_3 m_compare;
    double m_bound;
};

template < typename MeshT >
void remove_small_connected_components( 
    MeshT& mesh, 
    int limit_num=10, 
    double limit_di_angle_rad=2.36 /* 3/4 pi */ ) {
    
    const double bound = std::cos( limit_di_angle_rad );

    PMP::keep_large_connected_components(
        mesh, 
        limit_num,
        PMP::parameters::edge_is_constrained_map( 
            ConnectionConstraint<MeshT>(mesh, bound) )
    );
}

int main(int argc, char** argv) {
    const std::string in_mesh_fn = 
        "Please download from https://drive.google.com/drive/folders/1AmQONiss-q26cgs4wNBp3fMSlVmo09Ej?usp=sharing";
    
    // Read the mesh.
    Mesh_t mesh;
    read_mesh( in_mesh_fn, mesh );

    // Duplicate non-manifold_vertices.
    duplicate_non_manifold_vertices(mesh);

    // Remove self-intersection faces.
    remove_self_intersections(mesh);

    // Remove isolated vertices.
    remove_isolated_vertices(mesh);

    // Remove small connected components.
    remove_small_connected_components(mesh);
    
    return 0;
}