CGAL find interior points in a mesh

Question

Im trying to use CGAL to find the points that are in the interior of a triangle mesh. (i.e. the set of points that are not on its boundary. There is a similar example here for a 3D mesh which uses the function CGAL::Side_of_triangle_mesh<>. Can anyone help / advise on how to mod this for a 2D triangulation ?

My test code just creates two squares, one inside the other, then puts a seed at the origin (to make a hole) and then performs a 2D Delaunay triangulation. When I call the side_of_triangle_mesh<> class with:

Point_2 p = points2D[i];
CGAL::Bounded_side res = inside2D(p);

I get the error:

/usr/local/include/CGAL/Side_of_triangle_mesh.h:164:16: Candidate function not viable: no known conversion from 'Point_2' (aka 'Point_2<CGAL::Epick>') to 'const Point' (aka 'const Point_3<CGAL::Epick>') for 1st argument

Does this mean that side_of_triangle_mesh only works for 3D Polyhedron_3 mesh ? If so can anyone suggest a way to do this for a 2D mesh?

My test code is below: Thanks

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <vector>
#include <CGAL/Delaunay_mesher_2.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel     K;
typedef K::Point_2 Point_2;
typedef CGAL::Triangle_2<K> triangle;
typedef CGAL::Delaunay_mesh_vertex_base_2<K>                    Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K>                      Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb>            Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds>      CDT;
typedef CDT::Vertex_handle                                      Vertex_handle;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT>                Criteria;

int main(int argc, char* argv[])
{
    // Create a vector of the points
    //
    std::vector<Point_2> points2D ;
    points2D.push_back(Point_2(-5.0, -5.0));    // ----------
    points2D.push_back(Point_2( 5.0, -5.0));    // |   OUTER
    points2D.push_back(Point_2( 5.0,  5.0));    // |   SQUARE
    points2D.push_back(Point_2(-5.0,  5.0));    // ----------
    points2D.push_back(Point_2(-2.5, -2.5));    // ----------
    points2D.push_back(Point_2( 2.5, -2.5));    // |   INNER
    points2D.push_back(Point_2( 2.5,  2.5));    // |   SQUARE
    points2D.push_back(Point_2(-2.5,  2.5));    // ----------
    size_t numTestPoints = points2D.size();

    // Create a constrained delaunay triangulation and add the points
    //
    CDT cdt;
    std::vector<Vertex_handle> vhs;
    for (unsigned int i=0; i<numTestPoints; ++i){
        vhs.push_back(cdt.insert(points2D[i]));
    }

    // Creare constraints of the sides of the mesh
    //
    cdt.insert_constraint(vhs[0],vhs[1]);
    cdt.insert_constraint(vhs[1],vhs[2]);
    cdt.insert_constraint(vhs[2],vhs[3]);
    cdt.insert_constraint(vhs[3],vhs[0]);

    cdt.insert_constraint(vhs[4],vhs[5]);
    cdt.insert_constraint(vhs[5],vhs[6]);
    cdt.insert_constraint(vhs[6],vhs[7]);
    cdt.insert_constraint(vhs[7],vhs[4]);

    // Create a seed to make sure the inner square is a hole
    //
    std::list<Point_2> list_of_seeds;
    list_of_seeds.push_back(Point_2(0, 0));

    // Refine the mesh into triangles of max side length '1' and ensuring seeds are 'holes'
    //
    CGAL::refine_Delaunay_mesh_2(cdt, list_of_seeds.begin(), list_of_seeds.end(),
                                 Criteria(0.125, 1),false);


    // Call side_of_triangle_mesh
    //
    CGAL::Side_of_triangle_mesh<CDT, K> inside2D(cdt) ;

    int n_inside = 0;
    int n_boundary = 0;
    for(unsigned int i=0; i<numTestPoints; ++i)
    {
        Point_2 p = points2D[i];
        CGAL::Bounded_side res = inside2D(p);
        // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        // NO MATCHING FUNCTION Call

        if (res == CGAL::ON_BOUNDED_SIDE) { ++n_inside; }
        if (res == CGAL::ON_BOUNDARY) { ++n_boundary; }
    }
    std::cerr << "2D results for query size: " << cdt.number_of_vertices() << std::endl;
    std::cerr << "  " << n_inside << " points inside " << std::endl;
    std::cerr << "  " << n_boundary << " points on boundary " << std::endl;

    return 0 ;
}

CobaltGray · Accepted Answer · 2017-05-26T16:57:59.043

Thanks to @sloriot I got to the bottom of this. The bottom line is that you cannot use CGAL::Side_of_triangle_mesh<> for 2D Delaunay meshes. Instead what you have to do is loop through all the vertices in the mesh and test each one by looking at all the neighbouring faces around the vertex point. if any of the faces are outside of the domain then the vertex or point is on the edge of the mesh.

Here is the corrected code:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Delaunay_mesh_vertex_base_2.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <vector>
#include <CGAL/Delaunay_mesher_2.h>


typedef CGAL::Exact_predicates_inexact_constructions_kernel     K;
typedef K::Point_2                                              Point_2;
typedef CGAL::Delaunay_mesh_vertex_base_2<K>                    Vb;
typedef CGAL::Delaunay_mesh_face_base_2<K>                      Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb>            Tds;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, Tds>      CDT;
typedef CDT::Vertex_handle                                      Vertex_handle;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT>                Criteria;
typedef CDT::Vertex_iterator                                    Vertex_iterator;
typedef CDT::Vertex                                             Vertex;
typedef CDT::Face                                               Face ;
typedef Face::Face_handle                                       Face_handle ;
typedef CDT::Face_circulator                                    Face_circulator ;


int main(int argc, char* argv[])
{

    // Create a vector of the points
    //
    std::vector<Point_2> points2D ;
    points2D.push_back(Point_2(-5.0, -5.0));    // ----------
    points2D.push_back(Point_2( 5.0, -5.0));    // |   OUTER
    points2D.push_back(Point_2( 5.0,  5.0));    // |   SQUARE
    points2D.push_back(Point_2(-5.0,  5.0));    // ----------
    points2D.push_back(Point_2(-2.5, -2.5));    // ----------
    points2D.push_back(Point_2( 2.5, -2.5));    // |   INNER
    points2D.push_back(Point_2( 2.5,  2.5));    // |   SQUARE
    points2D.push_back(Point_2(-2.5,  2.5));    // ----------
    size_t numTestPoints = points2D.size();

    // Create a constrained delaunay triangulation and add the points
    //
    CDT cdt;
    std::vector<Vertex_handle> vhs;
    for (unsigned int i=0; i<numTestPoints; ++i){
        vhs.push_back(cdt.insert(points2D[i]));
    }

    // Creare constraints of the sides of the mesh
    //
    cdt.insert_constraint(vhs[0],vhs[1]);
    cdt.insert_constraint(vhs[1],vhs[2]);
    cdt.insert_constraint(vhs[2],vhs[3]);
    cdt.insert_constraint(vhs[3],vhs[0]);

    cdt.insert_constraint(vhs[4],vhs[5]);
    cdt.insert_constraint(vhs[5],vhs[6]);
    cdt.insert_constraint(vhs[6],vhs[7]);
    cdt.insert_constraint(vhs[7],vhs[4]);

    // Create a seed to make sure the inner square is a hole
    //
    std::list<Point_2> list_of_seeds;
    list_of_seeds.push_back(Point_2(0, 0));

    // Refine the mesh into triangles of max side length '1' and ensuring seeds are 'holes'
    //
    CGAL::refine_Delaunay_mesh_2(cdt, list_of_seeds.begin(), list_of_seeds.end(),
                                 Criteria(0.125, 1.5),false);

    // The mesh is now created. The next bit swings around each vertex point checking that
    // all faces around it are in the domain. If any are not then the vertex is on the
    // edge of the mesh.
    // thanks to @sloriot for this
    //

    Vertex v ;
    std::vector<Point_2> interior_points ;
    std::vector<Point_2> boundary_points ;
    CDT::Locate_type loc = CDT::Locate_type::VERTEX ;
    int li;

    Vertex_iterator vit = cdt.vertices_begin(), beyond = cdt.vertices_end() ;
    while (vit != beyond) {
        v = *vit ;
        Point_2 query = vit->point();
        Face_handle f =  cdt.locate(query, loc, li);

        bool current_face_in_domain ;
        Face_handle start = f ;
        do {
            current_face_in_domain = f->is_in_domain() ;
            Vertex_handle vh = f->vertex(li);
            f = f->neighbor(cdt.ccw(li));
            li = f->index(vh) ;
        }
        while(current_face_in_domain && f != start) ;

        if (f == start) {
            interior_points.push_back(query) ;
        }else{
            boundary_points.push_back(query) ;
        }
        ++vit ;
    }

    printf("Boundary points:\n");
    for(auto p = boundary_points.begin(); p != boundary_points.end(); ++p){
        printf("%f,%f\n",p->x(), p->y()) ;
    }

    printf("interior points:\n");
    for(auto p = interior_points.begin(); p != interior_points.end(); ++p){
        printf("%f,%f\n",p->x(), p->y()) ;
    }

    return 0 ;
}

When you run it you should get something like : This

sloriot · Answer 2 · 2017-05-26T14:04:42.577

0

You must use the locate function. It will return a face handle containing your query point. Then you need to check whether the face is in the domain or not using the is_in_domain() member function.

Note that if you have several queries to do, you should first put all your query points in a container, sort them using hilbert_sort and use the location of the previous point as second parameter to the locate function.

Here is a sample code of how to get the incident faces in boundary cases:

CDT_2 mesh;

CGAL::Locate_type loc;
int li;

Point_2 query;

Face_handle f =  mesh.locate(query, loc, li);

switch(loc)
{
  case FACE:
    f->is_in_domain();
  break;
  case EDGE:
  {
    bool face_1_in_domain = f->is_in_domain(); // first incident face status
    bool face_2_in_domain = f->neighbor(li)->is_in_domain(); // second incident face status
    break;
  }
  case VERTEX:
  {
    // turn around f->vertex(li)
    Face_handle start = f;
    do{
      bool current_face_in_domain = f->is_in_domain();
      Vertex_handle v = f->vertex(mesh.cw(li));
      f = f->neighbor(mesh.ccw(li));
      li = f->index(v);
    }
    while(f!=current);
    break;
  }
  default:
    // outside of the domain.

}

edited May 26 '17 at 14:04

answered May 26 '17 at 08:31

sloriot

6,070
18
27

Thanks @sloriot. What defines a 'Domain' and will a point on the edge of the hull be in or out of the domain? (I think 'in' from the manual which is not what I want. I want to isolate points on the boundary of the mesh from those in the interior.) You did give me a thought though. Could i just get the face handle associated with a vertex and from that get the other two vertices. Then find the neighbour face for each vertex and if any are 'infinite' then the test point must be on the boundary? There must be an easier way? Good tip on the hilbert_sort. – CobaltGray May 26 '17 at 08:57
Actually for the `locate` function I meant another overload that also include the location type (face, vertex or edge). I've updated my answer. The domain is defined by the connected components containing as least one seed point. – sloriot May 26 '17 at 12:18
Thanks for trying but that still doesn't answer my question - I need to discriminate between vertices that are on the boundary of the mesh and those that are in its interior. Ive implemented the is_in_domain() function as you suggested and it includes those vertices at the edge of the mesh. Unfortunately my own contrib didn't work either as a vertex point on the edge of a mesh can still have three non-infinite neighbours. – CobaltGray May 26 '17 at 12:47
If the point fall on an edge, look at the two incident faces and check the value of `is_in_domain()` for both faces. This will indicates if you are inside, outside, or on the boundary of the domain. It you fall on a vertex, it is equivalent but you need to consider the value of `is_in_domain()` for all incident faces. – sloriot May 26 '17 at 12:54
oooh K. I realise I'm starting to look dumb here but i'll carry on testing your patience.... I think I see what you are getting at but how do I find all the faces for a given vertex? I can loop through all vertices in the mesh and for each one, I can locate() its face but then I need to find the other faces also connected to this vertex. Once I have all the faces I can call is_in_domain() for each face. and if any face connecting to the vertex is outside the domain then its on the edge. Brill! SO how can I get the handle for each face associated with a vertex ? Thanks – CobaltGray May 26 '17 at 13:34
Ok thanks @sloriot. You had a bunch of bugs in your code which didn't help my understanding but you gave me just enough info to understand the API a bit better. I have edited your code and got mine working - i'll post my full answer for other. Thanks again - couldn't have done this without you. – CobaltGray May 26 '17 at 16:30

CGAL find interior points in a mesh

2 Answers2