Trouver le résultat triangulé de l'intersection/différence de deux triangles 2D dans CGAL

Question

J'essaie de trouver le résultat d'une intersection et la différence entre deux triangles dans CGAL.

Mon approche actuelle est de convertir les triangles en polygones, calculer l'intersection/différence et enfin trianguler le résultat.

Le code suivant me donne 'Variable utilisée avant d'être initialisée (ou bogue CGAL)' dans le test d'intersection. Des idées sur ce qui pourrait être faux?

#include <iostream> 
#include <CGAL/basic.h> 
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> 
#include <CGAL/Polygon_2.h> 
#include <CGAL/triangulate_polyhedron.h> 
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> 
#include <CGAL/Constrained_Delaunay_triangulation_2.h> 
#include <CGAL/Triangulation_face_base_with_info_2.h> 

#include <vector> 

#include <CGAL/Boolean_set_operations_2.h> 
#include <list> 


struct FaceInfo2 
{ 
    FaceInfo2(){} 
    int nesting_level; 
    bool in_domain(){ 
     return nesting_level%2 == 1; 
    } 
}; 


typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; 
typedef Kernel::Point_2         Point_2; 
typedef CGAL::Polygon_2<Kernel>       Polygon_2; 
typedef CGAL::Polygon_with_holes_2<Kernel>    Polygon_with_holes_2; 
typedef std::list<Polygon_with_holes_2>     Pwh_list_2; 
typedef CGAL::Triangulation_vertex_base_2<Kernel>      Vb; 
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2,Kernel> Fbb; 
typedef CGAL::Constrained_triangulation_face_base_2<Kernel,Fbb>  Fb; 
typedef CGAL::Triangulation_data_structure_2<Vb,Fb>    TDS; 
typedef CGAL::Exact_predicates_tag        Itag; 
typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel, TDS, Itag> CDT; 

/*typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; 
typedef Kernel::Point_2 Point_2; 
typedef Kernel::Segment_2 Segment_2; 

typedef CGAL::Polygon_2<Kernel> Polygon_2; 
*/ 
typedef Kernel::Triangle_2  Triangle; 

// from http://doc.cgal.org/latest/Triangulation_2/index.html#Subsection_37.8.2 
void mark_domains(CDT& ct, 
        CDT::Face_handle start, 
        int index, 
        std::list<CDT::Edge>& border) 
{ 
    if(start->info().nesting_level != -1){ 
     return; 
    } 
    std::list<CDT::Face_handle> queue; 
    queue.push_back(start); 
    while(! queue.empty()){ 
     CDT::Face_handle fh = queue.front(); 
     queue.pop_front(); 
     if(fh->info().nesting_level == -1){ 
      fh->info().nesting_level = index; 
      for(int i = 0; i < 3; i++){ 
       CDT::Edge e(fh,i); 
       CDT::Face_handle n = fh->neighbor(i); 
       if(n->info().nesting_level == -1){ 
        if(ct.is_constrained(e)) border.push_back(e); 
        else queue.push_back(n); 
       } 
      } 
     } 
    } 
} 

// from http://doc.cgal.org/latest/Triangulation_2/index.html#Subsection_37.8.2 
//explore set of facets connected with non constrained edges, 
//and attribute to each such set a nesting level. 
//We start from facets incident to the infinite vertex, with a nesting 
//level of 0. Then we recursively consider the non-explored facets incident 
//to constrained edges bounding the former set and increase the nesting level by 1. 
//Facets in the domain are those with an odd nesting level. 
void 
mark_domains(CDT& cdt) 
{ 
    for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){ 
     it->info().nesting_level = -1; 
    } 
    std::list<CDT::Edge> border; 
    mark_domains(cdt, cdt.infinite_face(), 0, border); 
    while(! border.empty()){ 
     CDT::Edge e = border.front(); 
     border.pop_front(); 
     CDT::Face_handle n = e.first->neighbor(e.second); 
     if(n->info().nesting_level == -1){ 
      mark_domains(cdt, n, e.first->info().nesting_level+1, border); 
     } 
    } 
} 

std::vector<Triangle> triangulate(Pwh_list_2 & polyList){ 
    std::vector<Triangle> res; 
    for (Polygon_with_holes_2 polygon1 : polyList) { 
     CDT cdt; 
     auto outer = polygon1.outer_boundary(); 

     cdt.insert_constraint(outer.vertices_begin(), outer.vertices_end(), true); 

     for (auto hole = polygon1.holes_begin(); hole != polygon1.holes_end(); hole++){ 
      cdt.insert_constraint(hole->vertices_begin(), hole->vertices_end(), true); 
     } 

     for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin(); 
      fit!=cdt.finite_faces_end();++fit) 
     { 
      std::cout <<"New triangle "<<std::endl; 
      if (fit->info().in_domain()) { 
       Triangle tri; 

       for (int i=0;i<3;i++){ 
        auto point = fit->vertex(i)->point(); 
        tri.vertex(i) =fit->vertex(i)->point(); 
        std::cout <<point.x()<<" "; 
       } 
       res.push_back(tri); 
       std::cout <<std::endl; 

      } 
     } 
    } 
    return res; 
} 

void ProjectOtherTriangle(Triangle thisTri, Triangle tri, 
          std::vector<Triangle>& differenceTriangle, std::vector<Triangle>& intersectionTriangles){ 

    Polygon_2 thisPoly; 
    for (int i=0;i<3;i++){ 
     thisPoly.push_back(thisTri[i]); 
    } 
    Polygon_2 poly; 
    for (int i=0;i<3;i++){ 
     poly.push_back(tri[i]); 
    } 

    // Compute the intersection of P and Q. 
    Pwh_list_2     intR; 
    CGAL::intersection (thisPoly, poly, std::back_inserter(intR)); 
    // add into newTriangles 
    auto tris = triangulate(intR); 
    intersectionTriangles.insert(intersectionTriangles.end(), tris.begin(), tris.end()); 

    // Find difference 
    CGAL::difference (thisPoly, poly, std::back_inserter(intR)); 
    // add into currentTriangle 
    tris = triangulate(intR); 
    differenceTriangle.insert(differenceTriangle.end(), tris.begin(), tris.end()); 
} 

void printTriangle(Triangle t){ 
    using namespace std; 
    for (int i=0;i<3;i++){ 
     cout << t[i] <<" "; 
    } 
    cout << endl; 
} 


int main() 
{ 
    // Construct the two input polygons. 
    Triangle P; 
    P[0] = Point_2 (0, 0); 
    P[1] = Point_2 (2, 0); 
    P[2] = Point_2 (1, 1.5); 
    printTriangle(P); 
    Triangle Q; 
    Q[0] = Point_2 (0, 2); 
    Q[1] = Point_2 (1, 0.5); 
    Q[2] = Point_2 (2, 2); 
    printTriangle(Q); 

    std::vector<Triangle> currentTriangle; 
    std::vector<Triangle> newTriangle; 

    ProjectOtherTriangle(P,Q, currentTriangle, newTriangle); 

    using namespace std; 

    cout << "currentTriangle"<<endl; 
    for (auto t : currentTriangle){ 
     printTriangle(t); 
    } 

    cout << "newTriangle"<<endl; 
    for (auto t : newTriangle){ 
     printTriangle(t); 
    } 


    return 0; 
} 

Answer 1

Vous ne pouvez pas modifier un objet Triangle_2. Ce qui suit n'est pas correct:

Triangle P; 
P[0] = Point_2 (0, 0); 
P[1] = Point_2 (2, 0); 
P[2] = Point_2 (1, 1.5); 
printTriangle(P); 
Triangle Q; 
Q[0] = Point_2 (0, 2); 
Q[1] = Point_2 (1, 0.5); 
Q[2] = Point_2 (2, 2); 
printTriangle(Q); 

Vous devez utiliser le constructeur à partir de 3 points.