在CGAL中查找三角形交叉點/兩個二維三角形的差異結果

Question

我試圖找到CGAL中兩個三角形之間的交集和差異的結果。

我目前的做法是將三角形轉換爲多邊形，計算相交/差異並最終對結果進行三角測量。

以下代碼給出了'在交集測試中初始化（或CGAL錯誤）之前使用的變量'。任何想法可能是錯的？

#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

您無法修改Triangle_2對象。以下是不正確的：

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); 

您應該使用3點的構造函數。