CGAL带洞剖分

先附上效果图

CGAL环境配置，自行百度。。。。。

完整代码如下：

#include"vtkPolyData.h"
#include"vtkPolyDataWriter.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 <CGAL/Polygon_2.h>
#include <iostream>
struct FaceInfo2
{
    FaceInfo2() {}
    int nesting_level;
    bool in_domain() {
        return nesting_level % 2 == 1;
    }
};
typedef CGAL::Exact_predicates_inexact_constructions_kernel       K;
typedef CGAL::Triangulation_vertex_base_2<K>                      Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K>    Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<K, Fbb>        Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb>               TDS;
typedef CGAL::Exact_predicates_tag                                Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag>  CDT;
typedef CDT::Point                                                Point;
typedef CGAL::Polygon_2<K>                                        Polygon_2;
typedef CDT::Face_handle                                          Face_handle;
void
mark_domains(CDT& ct,
    Face_handle start,
    int index,
    std::list<CDT::Edge>& border)
{
    if (start->info().nesting_level != -1) {
        return;
    }
    std::list<Face_handle> queue;
    queue.push_back(start);
    while (!queue.empty()) {
        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);
                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);
                }
            }
        }
    }
}
//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::Face_handle f : cdt.all_face_handles()) {
        f->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();
        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);
        }
    }
}

void ConstructPolyData(vtkPolyData*& p, const CDT& cdt)
{
    if (cdt.number_of_vertices() > 0) {

        vtkPoints* pp = vtkPoints::New();
        pp->SetDataTypeToDouble();
        p = vtkPolyData::New();
        p->SetPoints(pp);
        pp->Delete();


        std::map<CDT::Vertex_handle, int> V;
        int inum = 0;
        CDT::Finite_vertices_iterator v_iter;
        for (v_iter = cdt.finite_vertices_begin(); v_iter != cdt.finite_vertices_end(); ++v_iter) {
            Point p = v_iter->point();
            V[v_iter] = inum;
            inum++;
            pp->InsertNextPoint(p.x(), p.y(), 0.0);
        }



        if (cdt.number_of_faces() > 0)
        {
            vtkCellArray* tt = vtkCellArray::New();
            CDT::Finite_faces_iterator f_iter;

            for (f_iter = cdt.finite_faces_begin(); f_iter != cdt.finite_faces_end(); f_iter++) {
                if (f_iter->info().in_domain()) {
                    int v0 = V[f_iter->vertex(0)];
                    int v1 = V[f_iter->vertex(1)];
                    int v2 = V[f_iter->vertex(2)];

                    tt->InsertNextCell(3);
                    tt->InsertCellPoint(v0);
                    tt->InsertCellPoint(v1);
                    tt->InsertCellPoint(v2);
                }              
            }

            p->SetPolys(tt);
            tt->Delete();
        }
    }
}

int main()
{
    //construct two non-intersecting nested polygons
    Polygon_2 polygon1;
    polygon1.push_back(Point(0, 0));
    polygon1.push_back(Point(0, 2));
    polygon1.push_back(Point(2, 2));
    polygon1.push_back(Point(2, 0));

  
    Polygon_2 polygon2;
    polygon2.push_back(Point(0.5, 0.5));
    polygon2.push_back(Point(0.5, 1.5));
    polygon2.push_back(Point(1.5, 1.5));
    polygon2.push_back(Point(1.5, 0.5));
   
 
    //Insert the polygons into a constrained triangulation
    CDT cdt;
    cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);
    cdt.insert_constraint(polygon2.vertices_begin(), polygon2.vertices_end(), true);
    //Mark facets that are inside the domain bounded by the polygon
    mark_domains(cdt);
    int count = 0;
    for (Face_handle f : cdt.finite_face_handles())
    {
        if (f->info().in_domain()) ++count;
    }
    std::cout << "There are " << count << " facets in the domain." << std::endl;

    int face_num = cdt.number_of_faces();

    vtkPolyData* result = vtkPolyData::New();
   ConstructPolyData(result, cdt);

    vtkPolyDataWriter* w1 = vtkPolyDataWriter::New();
    w1->SetFileName("6666666666.vtk");
    w1->SetInputData(result);
    w1->Write();
    w1->Delete();
}