c++ - 在 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 点的构造函数。