CGAL users discussion list

[Yaoyu Hu <>]

Hi Sebastien,

You are very helpful!

I tried to use CGAL::Polygon_mesh_processing::remove_isolated_vertices() and it worked. There are 4 isolated points removed by the remove_isolated_vertices() function. Then the downstream isotropic remeshing works without any problem. 

As suggested by you, I am sharing the point cloud and sample code to do the advancing front surface reconstruction with CGAL. You could find the point cloud at

https://drive.google.com/file/d/1t0cwW6CgPrqA30i_x4KP7T42d66-vvt0/view?usp=sharing

and the source code as the attachment.

I tested the sample code on my machine with CGAL 5.0.2, the reconstructed mesh will have 4 isolated points that get removed by remove_isolated_vertices(). Please note that, the advancing front surface reconstruction takes hundreds of seconds with the Debug build but only a couple of tens seconds if executed in the Release build. I am not sure if the efficiency difference indicating any other issues.

Hope this helps to identify any potential issues in the reconstruction.

Regards,

Yaoyu

On Tue, Sep 8, 2020 at 12:52 AM "Sebastien Loriot (GeometryFactory)" <> wrote:

I opened your file in the demo and did not see any issue.
The issue you mentioned that is detected by the code indicates that you
have isolated vertices (i.e. vertices with no edge/face attached).

If you have such an error for a mesh that is the output of the CGAL
advancing front algorithm, then this is a bug. If this is the case,
could you please try to come up with an example and a point cloud
showing the issue so that we can reproduce it and fix it?

In the meantime you can also try the function remove_isolated_vertices()

https://doc.cgal.org/latest/Polygon_mesh_processing/group__PMP__repairing__grp.html#ga91c02ef57e638faf2622eae93e7a25e2

Thanks,

Sebastien.


On 9/7/20 6:54 PM, Yaoyu Hu ( via cgal-discuss
Mailing List) wrote:
> Hi Sebastien,
>
> Thank you for your reply. I tried to use CGAL::is_valid_polygon_mesh(),
> it indeed returns false. However, I do not know how to visualize the
> content in the CGAL::Verbose_ostream, then I manually stepped into the
> source code. I found that the check fails
> inside is_valid_halfedge_graph() which is defined in
> cgal/BGL/include/CGAL/boost/graph/helpers.h . The actual place that
> triggers the failure is the following code section
>
> ========== Code section begins. ==========
>
> for(vertex_descriptor vbegin : vertices(g))
>    {
>      // Pointer integrity.
>      if(halfedge(vbegin, g) != boost::graph_traits<Graph>::null_halfedge())
>        valid = (target(halfedge(vbegin, g), g) == vbegin);
>      else
>        valid = false;
>
>      if(!valid)
>      {
>        verr << "vertex " << v << " halfedge incident to vertex is the
> null halfedge." << std::endl;
>        verr << "Halfedge Graph Structure is NOT VALID." << std::endl;
>        return false;
>      }
>
> .....other codes..........
> }
>
> ========== Code section ends. ==========
>
> The v value is 12718 for this particular mesh. (mesh link:
> https://drive.google.com/file/d/1Z9AMfRnbkJGmDuUOc7ZAdMAvbhhhx_Vn/view?usp=sharing
> )
>
> The input mesh, as I provided in the previous email and as the above
> link, is a resulting mesh from the advancing front surface
> reconstruction of CGAL (I wrote another separate piece of code to do
> that). The reconstructed mesh was saved to the file system
> with CGAL::write_ply() function. This mesh can be processed by the
> Polyhedron Demo program which performs isotropic remeshing without any
> warnings or errors.
>
> What should I do to make the mesh a valid polygon mesh? Is there a way
> to fix this situation? I tried to remove the self-intersection facets by
> first using CGAL::Polygon_mesh_processing::self_intersections() and then
> CGAL::Euler::remove_face(). Additionally, I
> tried CGAL::Polygon_mesh_processing::duplicate_non_manifold_vertices().
> No intersections and no non-manifold vertices are found. But the mesh is
> still not valid.
>
> I may try to look into the vertex 12718 ( if using 0-based indexing, it
> is 12717 I suppose ) and try to figure out what makes this vertex so
> special.
>
> Hoping to get more insights into the problem from your help.
>
> Regards,
>
> Yaoyu
>
> On Mon, Sep 7, 2020 at 12:50 AM "Sebastien Loriot (GeometryFactory)"
> < <mailto:>> wrote:
>
>     Could you check the value returned by
>     CGAL::is_valid_polygon_mesh(my_mesh) called on your mesh before calling
>     the CGAL function? It might be that the input is not a valid polygon
>     mesh (non-manifold vertices for example).
>
>     Thanks,
>
>     Sebastien.
>
>
>
>     On 9/6/20 6:47 AM, Yaoyu Hu (
>     <mailto:> via cgal-discuss
>     Mailing List) wrote:
>      > Hi,
>      >
>      > I was trying to
>     use CGAL::Polygon_mesh_processing::isotropic_remeshing()
>      > to remesh a CGAL::Surface_mesh. However, I got the following
>     error message:
>      >
>      > ========== Error message begins. ==========
>      >
>      > terminate called after throwing an instance of
>     'CGAL::Assertion_exception'
>      >    what():  CGAL ERROR: assertion violation!
>      > Expr: is_valid_polygon_mesh(mesh_)
>      > File:
>      >
>     /home/yaoyu/Libraries/cgal/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Isotropic_remeshing/remesh_impl.h
>      > Line: 1037
>      >
>      > ========== Error message ends. ==========
>      >
>      > I was trying the code listed at
>      >
>     https://doc.cgal.org/latest/Polygon_mesh_processing/Polygon_mesh_processing_2isotropic_remeshing_example_8cpp-example.html
>
>      > and I enclosed the source file. I am not sure about the reason
>     why CGAL
>      > thinks the Surface_mesh is not a valid polygon mesh.
>      >
>      > The input mesh is in PLY format and you could find the file at
>      >
>     https://drive.google.com/file/d/1Z9AMfRnbkJGmDuUOc7ZAdMAvbhhhx_Vn/view?usp=sharing
>      >
>      > One thing that might be interesting is that if I use the CGAL
>     Polyhedron
>      > Demo program to load and do the isotropic remeshing, no errors are
>      > encountered.
>      >
>      > Any comments are appreciated.
>      >
>      > Thank you!
>      >
>      > Yaoyu
>      >
>      > --
>      > You are currently subscribed to cgal-discuss.
>      > To unsubscribe or access the archives, go to
>      > https://sympa.inria.fr/sympa/info/cgal-discuss
>      >
>
>     --
>     You are currently subscribed to cgal-discuss.
>     To unsubscribe or access the archives, go to
>     https://sympa.inria.fr/sympa/info/cgal-discuss
>
>
>
> --
> You are currently subscribed to cgal-discuss.
> To unsubscribe or access the archives, go to
> https://sympa.inria.fr/sympa/info/cgal-discuss
>

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss

//
// Created by yaoyu on 9/8/20.
//
// Some of the codes are copied from
// https://doc.cgal.org/latest/Advancing_front_surface_reconstruction/Advancing_front_surface_reconstruction_2reconstruction_surface_mesh_8cpp-example.html
// Some of the codes are copied from
// cgal/Polyhedron/demo/Polyhedron/Plugins/Point_set/Surface_reconstruction_advancing_front_impl.cpp
//

#include <algorithm>
#include <iostream>
#include <string>

#include <boost/math/constants/constants.hpp>

#include <CGAL/Advancing_front_surface_reconstruction.h>
#include <CGAL/array.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/IO/read_ply_points.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <CGAL/Surface_mesh.h>

// Namespace.
namespace PMP = CGAL::Polygon_mesh_processing;

// Global constants.
const auto G_PI = boost::math::constants::pi<double>();

// Typedefs for CGAL.
typedef std::array< std::size_t, 3 > Facet_t;

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel_t;
typedef Kernel_t::Point_3                                   Point3_t;
typedef CGAL::Surface_mesh<Point3_t>                        SurfaceMesh_t;

struct MeshConstructorByPointIndices {
    typedef typename boost::property_map< SurfaceMesh_t, boost::vertex_point_t >::type VPMap_t;

    SurfaceMesh_t &mesh;
    VPMap_t       vpMap;
    std::vector< typename boost::graph_traits<SurfaceMesh_t>::vertex_descriptor > vertices;

    template < typename Point_T >
    MeshConstructorByPointIndices( SurfaceMesh_t &m, const std::vector<Point_T> &points )
            : mesh{m} {
        vpMap = get(boost::vertex_point, mesh);
        for ( const auto& p : points ) {
            boost::graph_traits<SurfaceMesh_t>::vertex_descriptor v;
            v = add_vertex( mesh );
            vertices.push_back(v);
            put( vpMap, v, p );
        }
    }

    MeshConstructorByPointIndices& operator = ( const Facet_t f ) {
        typedef boost::graph_traits<SurfaceMesh_t>::vertex_descriptor DescV_t;
        std::vector< DescV_t > facet(3);

        facet[0] = vertices[ f[0] ];
        facet[1] = vertices[ f[1] ];
        facet[2] = vertices[ f[2] ];

        CGAL::Euler::add_face( facet, mesh );

        return *this;
    }

    MeshConstructorByPointIndices& operator *  () { return *this; }
    MeshConstructorByPointIndices& operator ++ () { return *this; }
    MeshConstructorByPointIndices  operator ++ (int) { return *this; }
};

struct PriorityRadius {
    double bound;

    explicit PriorityRadius ( double b )
            : bound{b} {}

    template < typename AdvFront_T, typename CellHandle_T >
    double operator () ( const AdvFront_T &adv, CellHandle_T &c, const int &index ) const {
        if ( 0 == bound ) {
            return adv.smallest_radius_delaunay_sphere( c, index );
        }

        double d = 0;
        d = sqrt( CGAL::squared_distance( c->vertex( (index+1)%4 )->point(),
                                          c->vertex( (index+2)%4 )->point() ) );
        if ( d > bound ) return adv.infinity();

        d = sqrt( CGAL::squared_distance( c->vertex( (index+2)%4 )->point(),
                                          c->vertex( (index+3)%4 )->point() ) );
        if ( d > bound ) return adv.infinity();

        d = sqrt( CGAL::squared_distance( c->vertex( (index+1)%4 )->point(),
                                          c->vertex( (index+3)%4 )->point() ) );
        if ( d > bound ) return adv.infinity();

        return adv.smallest_radius_delaunay_sphere( c, index );
    }
};

static std::vector<Point3_t> read_points_from_ply( const std::string &fn ) {
    std::ifstream ifs { fn };
    if ( !ifs ) {
        std::stringstream ss;
        ss << "Cannot open " << fn << " for reading. ";
        throw std::runtime_error( ss.str() );
    }

    std::vector<Point3_t> cgalPoints;
    if ( !CGAL::read_ply_points( ifs, std::back_inserter( cgalPoints ) ) ) {
        std::stringstream ss;
        ss << "read_ply() from " << fn << " failed. ";
        throw std::runtime_error( ss.str() );
    }

    ifs.close();

    return cgalPoints;
}

template < typename PT >
void write_mesh_ply(const std::string &fn,
                    CGAL::Surface_mesh<PT> &sm,
                    bool flagBinary= true) {
    std::ofstream ofs;

    if ( flagBinary ) {
        ofs.open( fn, std::ios::binary );
    } else {
        ofs.open( fn );
    }

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

    CGAL::write_ply( ofs, sm );

    ofs.close();
}

static void reconstruct_surface_mesh(
        const std::vector<Point3_t> &points, SurfaceMesh_t &mesh,
        double longestEdge=0.1, double radiusRatioBound=5.0, double betaDeg=30.0 ) {
    std::cout << "Begin reconstruction by AFS. \n";
    PriorityRadius priority(longestEdge);
    MeshConstructorByPointIndices meshConstructor( mesh, points );
    std::cout << "AFS... \n";
    CGAL::advancing_front_surface_reconstruction(
            points.begin(), points.end(),
            meshConstructor, priority,
            radiusRatioBound, betaDeg/180.0*G_PI );
}

int main( int argc, char **argv ) {
    std::cout << "Hello, Mesh_Reconstruction! \n";

    if ( argc != 3 ) {
        std::stringstream ss;
        ss << "Expecting 3 input arguments. " << argc << " given. ";
        throw std::runtime_error( ss.str() );
    }

    std::string inPointsFn = argv[1];
    std::string outMeshFn  = argv[2];

    // Load the PLY point cloud by PCL and convert the points into CGAL point representation.
    auto points = read_points_from_ply( inPointsFn );

    SurfaceMesh_t mesh;
    reconstruct_surface_mesh( points, mesh );

    // Check the reconstructed mesh.
    auto res = CGAL::is_valid_polygon_mesh(mesh);
    std::cout << "rs = " << res << "\n";

    // Remove isolated points.
    auto nIsolated = PMP::remove_isolated_vertices( mesh );
    std::cout << "nIsolated = " << nIsolated << "\n";

    // Save the mesh.
    write_mesh_ply( outMeshFn, mesh );

    return 0;
}