CGAL users discussion list

[Yaoyu Hu <>]

Hi,

I am trying to fix and clear some issues from a surface mesh. I do the following three steps in sequence: 

(1) Duplicate non-manifold vertices.

(2) Remove self-intersections. And,

(3) Remove isolated vertices.

(4) Remove small connected components.

I always get an assertion error in (4). In my code, (1), (3), and (4) are directly calling CGAL functions to do the job. For (2) I have a simple procedure where I first identify all the self-intersection pairs and then remove all the faces that have self-intersection. I have set up a sample code and a sample mesh. The code is also available as a GitHub Gist 

https://gist.github.com/huyaoyu/c10ff67a2fe67c141af77142e475a83d

The sample mesh is saved in my personal Google Drive.

https://drive.google.com/drive/folders/1AmQONiss-q26cgs4wNBp3fMSlVmo09Ej?usp=sharing

I am using CGAL 5.3 on Ubuntu 20.04.

=== Code begins. ===

#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;
}

=== Code ends. ===

Step (2) is the remove_self_intersections() function. The assertion error shows in my terminal is

===  Terminal output begins. ===

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: _idx < data_.size()
File: /usr/local/include/CGAL/Surface_mesh/Properties.h
Line: 193

=== Terminal output ends. ===

I have also saved the backtrace from GDB in the Google Drive folder (a txt and a screenshot). please download from https://drive.google.com/drive/folders/1AmQONiss-q26cgs4wNBp3fMSlVmo09Ej?usp=sharing

From the backtrace, it looks like when iterating through the faces of the mesh for computing the connected components, a face descriptor is out of the valid range.

Another thing that I found during the debugging is that, if you check the output of mesh.has_garbage() after calling PMP::remove_isolated_vertices() in my remove_isolated_vertices() function, mesh.has_garbage() is true. I suppose  PMP::remove_isolated_vertices() should not leave the mesh in a state where it has uncollected garbage.

Any help is appreciated!

Thank you!

Yaoyu