CGAL users discussion list

[Per Zetterlund <>]

I've created a minimal test case, see attached file.


MyMeshTrait::Construct_initial_points::operator() gets called
with n = 20, and if I feed it with 20 points I get

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: aspect_ratio >= 0 && aspect_ratio <= 1
File: /usr/include/CGAL/Surface_mesher/Standard_criteria.h
Line: 161
Aborted (core dumped)

If I instead feed it with 4 points it seems to work,
and of I input the center point of *all* triangles, I get

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: is_finite(d)
File: /usr/include/CGAL/GMP/Gmpq_type.h
Line: 132
Aborted (core dumped)

What is going on here? Please help me find what I'm doing wrong. I want
to input many initial points to cover all disjoint component in the mesh.

 /Per

Here is the input file:

https://dl.dropboxusercontent.com/u/28611985/in.off.zip

On Tue, Jul 22, 2014 at 11:10 AM, Per Zetterlund 
<>
 wrote:
> I've gotten a bit further with this.
>
> I'm using AABB Tree to speed up the intersects.
>
> However, my input mesh has several components and it
> is causing some trouble. The 3D Surface Mesh docs
> says that it handles disconnected components as long
> as there is an initial point on each component. This is
> not working for me; I get error messages such as
>
> terminate called after throwing an instance of 'CGAL::Assertion_exception'
>   what():  CGAL ERROR: assertion violation!
> Expr: is_finite(d)
> File: /usr/include/CGAL/GMP/Gmpq_type.h
> Line: 132
>
> But if I restrict the initial points to lay on one of the
> components I get a nice reconstruction of that component.
>
> Any ideas what I'm doing wrong?
>
>  /Per
>
>
> On Wed, Jul 16, 2014 at 2:57 PM, Per Zetterlund 
> <>
>  wrote:
>> Hi All,
>>
>> I have a triangular mesh that is ill formed (i.e. crosses
>> itself at vertices and edges) that I want to remesh to
>> a nice CGAL::Polyhedron.
>>
>> Right now I've written my own SurfaceMeshTraits class
>> that I use with make_surface_mesh (3D Surface Mesh Generation
>> package). It's brute force in that for each line/ray/segment I check
>> intersection with every triangle in the mesh. So it is very slow...
>> Is there a better way?
>>
>> I found the Surface_mesher/polyhedron_remesher.cpp example
>> but I can't find documentation for neither CGAL::Polyhedral_surface_3
>> nor CGAL::make_piecewise_smooth_surface_mesh so I'm guessing it's
>> outdated?
>>
>> What is the difference between the "3D Surface Mesh Generation"
>> package and the "3D Mesh Generation" package? Which one is
>> most suitable for my needs?
>>
>> Thanks,
>>
>>  Per
#include <fstream>
#include <set>
#include <stack>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_triangle_primitive.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

typedef Kernel::Point_3 Point;
typedef Kernel::Triangle_3 Triangle;

typedef std::vector<Triangle>::iterator Iterator;
typedef CGAL::AABB_triangle_primitive<Kernel,Iterator> Primitive;
typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;

struct MyMeshTrait {
    typedef Point Intersection_point;

    struct Surface_3 {
        const std::vector<Triangle>& mTriangles;
        Tree mTree;

        Surface_3(std::vector<Triangle> triangles) : mTriangles(triangles),
                                                     mTree(triangles.begin(),triangles.end()) {}
    };

    struct Construct_initial_points {
        template <class OutputIteratorPoints>
        OutputIteratorPoints operator() (const Surface_3& surf, OutputIteratorPoints pts, int n = 0) {

            assert(surf.mTriangles.size() >= 20);

            for(int i = 0; i < n; ++i) {
                *pts = CGAL::centroid(surf.mTriangles[i]);
                ++pts;
            }
            return pts;
        }
    };

    struct Intersect_3 {
        template <class T>
        CGAL::Object operator()(const Surface_3& surf, const T& s) const {
            auto res = surf.mTree.any_intersection(s);
            if(res) {
                return res->first;
            } else {
                return CGAL::Object();
            }
        }
    };

    Construct_initial_points construct_initial_points_object() const
    {
      return Construct_initial_points();
    }

    Intersect_3 intersect_3_object() const
    {
      return Intersect_3();
    }
};

std::vector<Triangle> importOffVectorOfTriangles(const std::string& filename) {
    std::ifstream in(filename);
    std::vector<Triangle> triangles;

    bool color = false;

    std::string header;
    in >> header;
    if(header == "COFF") {
        color = true;
    } else {
        assert(header == "OFF");
    }

    int num_verts;
    int num_triangles;
    int num_edges;

    in >> num_verts >> num_triangles >> num_edges;
    assert(num_edges == 0);


    std::vector<Point> verts(num_verts);
    for(int t = 0; t < num_verts; ++t) {
        float x,y,z;
        in >> x >> y >> z;
        verts[t] = Point(x,y,z);
        if(color) {
            int dummy;
            in >> dummy >> dummy >> dummy >> dummy;
        }
    }

    for(int t = 0; t < num_triangles; ++t) {
        std::vector<Point> tri(3);
        int num;
        in >> num;
        assert(num == 3);
        for(int i = 0; i < 3; ++i) {
            int idx;
            in >> idx;
            tri[i] = verts[idx];
        }
        triangles.emplace_back(tri[0], tri[1], tri[2]);
    }
    return triangles;
}

static void generateMesh(const std::vector<Triangle>& triangles) {
    CGAL::Surface_mesh_default_triangulation_3 tr;
    CGAL::Complex_2_in_triangulation_3<CGAL::Surface_mesh_default_triangulation_3> c2t3(tr);
    MyMeshTrait::Surface_3 surface(triangles);
    MyMeshTrait traits;
    CGAL::Surface_mesh_default_criteria_3<CGAL::Surface_mesh_default_triangulation_3> criteria(30, 0.06, 0.06);

    CGAL::make_surface_mesh(c2t3, surface, traits, criteria, CGAL::Manifold_with_boundary_tag(),3);

    CGAL::Polyhedron_3<Kernel> polyhedron;

    CGAL::output_surface_facets_to_polyhedron(c2t3, polyhedron);

    polyhedron.inside_out();

    std::ofstream out("out.off");
    out << polyhedron;
}

int main(int argc, char** argv) {
    generateMesh(importOffVectorOfTriangles("in.off"));
}