CGAL trouver des points intérieurs dans un maillage

Question

En essayant d'utiliser CGAL pour trouver les points qui sont à l'intérieur d'un maillage triangle. (c'est-à-dire l'ensemble des points qui ne sont pas sur sa frontière.) Il y a un similar example here pour un maillage 3D qui utilise la fonction CGAL :: Side_of_triangle_mesh <> Quelqu'un peut-il aider/conseiller sur la façon de moduler cela pour une triangulation 2D?

Mon code de test crée juste deux carrés, l'un dans l'autre, puis met une graine à l'origine (pour faire un trou) et effectue ensuite une triangulation 2D Delaunay Quand j'appelle la side_of_triangle_mesh <> classe.

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

Je reçois l'erreur:

/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 

Est-ce que cela signifie que side_of_triangle_mesh ne fonctionne que pour le maillage 3D Polyhedron_3? Si oui, quelqu'un peut-il suggérer un moyen de le faire pour un maillage 2D?

Mon code de test est ci-dessous: Merci

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

Answer 1

Merci à @sloriot je suis arrivé au bas de. La ligne du bas est que vous ne pouvez pas utiliser CGAL :: Side_of_triangle_mesh <> pour les maillages Delaunay 2D. Au lieu de cela, vous devez parcourir tous les sommets du maillage et tester chacun d'entre eux en regardant toutes les faces voisines autour du sommet. si l'une des faces est en dehors du domaine, le sommet ou le point est sur le bord du maillage.

Voici le code corrigé:

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

Lorsque vous exécutez vous devriez obtenir quelque chose comme: This

Answer 2

Vous devez utiliser la fonction locate. Il renverra un handle de visage contenant votre point de requête. Ensuite, vous devez vérifier si le visage est dans le domaine ou non en utilisant la fonction membre is_in_domain(). Notez que si vous avez plusieurs requêtes à faire, vous devez d'abord placer tous vos points de requête dans un conteneur, les trier en utilisant hilbert_sort et utiliser l'emplacement du point précédent comme deuxième paramètre de la fonction de localisation.

Voici un exemple de code de la façon d'obtenir les visages d'incidents dans les cas limites:

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. 

}