[CGAL]带岛多边形三角化
CGAL带岛多边形三角化，并输出(*.ply)格式的模型

模型输出的关键是节点和索引

#include <CGAL/Triangulation_vertex_base_with_id_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>

因此注意这两个泛型，对比不带信息的

#include <CGAL/Triangulation_vertex_base_2.h>
#include <CGAL/Triangulation_face_base_2.h>，这两个增加了部分信息作为拓展。

这样Vertex_handle就可以读取这部分拓展的信息。

心得：CGAL的泛型机制真的很强大，拓展性很好。

// AxModelDelaunay.cpp : 定义控制台应用程序的入口点。
//#include "stdafx.h"
#include "shapefil.h"#include "CGAL/exceptions.h"
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_id_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_with_id_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::Vertex_handle    Vertex_handle;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);}}}}
}
//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);}}
}int _tmain(int argc, _TCHAR* argv[])
{//读取shpconst char * pszShapeFile = "data\\walltest.shp";SHPHandle hShp = SHPOpen(pszShapeFile, "r");int nShapeType, nVertices;int nEntities = 0;double* minB = new double[4];double* maxB = new double[4];SHPGetInfo(hShp, &nEntities, &nShapeType, minB, maxB);printf("ShapeType:%d\n", nShapeType);printf("Entities:%d\n", nEntities);CDT cdt;for (int i = 0; i < nEntities; i++){int iShape = i;SHPObject *obj = SHPReadObject(hShp, iShape);printf("--------------Feature:%d------------\n", iShape);int parts = obj->nParts;int* partStart = obj->panPartStart;int verts = obj->nVertices;printf("nParts:%d\n", parts);printf("nVertices:%d\n", verts);for (int k = 0; k < parts; k++){Polygon_2 polygon1;int fromIdx = partStart[k];int toIdx = fromIdx;if (k<parts-1){toIdx = partStart[k + 1];}else{toIdx = verts;}for (size_t j = fromIdx; j < toIdx; j++){double x = obj->padfX[j];double y = obj->padfY[j];//Point_2 pt(x, y);printf("%f,%f;", x, y);polygon1.push_back(Point(x, y));}cdt.insert_constraint(polygon1.vertices_begin(), polygon1.vertices_end(), true);}printf("\n");}//construct two non-intersecting nested polygons  //Polygon_2 polygon1;//polygon1.push_back(Point(0, 0));//polygon1.push_back(Point(2, 0));//polygon1.push_back(Point(2, 2));//polygon1.push_back(Point(0, 2));//Polygon_2 polygon2;//polygon2.push_back(Point(0.5, 0.5));//polygon2.push_back(Point(1.5, 0.5));//polygon2.push_back(Point(1.5, 1.5));//polygon2.push_back(Point(0.5, 1.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 polygonmark_domains(cdt);FILE *ply = fopen("data\\floorpeint.ply", "w");int idx = 0;for (CDT::Vertex_iterator v = cdt.vertices_begin(); v != cdt.vertices_end(); ++v){Vertex_handle vv = v->handle();vv->id() = idx;idx++;}int count = 0;for (CDT::Finite_faces_iterator fit = cdt.finite_faces_begin();fit != cdt.finite_faces_end(); ++fit){if (fit->info().in_domain()){++count;for (int i = 0; i < 3; i++){Vertex_handle vert = fit->vertex(i);int x=vert->id();std::cout << "The Id is " << x << std::endl;CDT::Edge ed(fit, i);ed.second;}}}if (ply){fprintf(ply, "ply\nformat %s 1.0\n", "ascii");fprintf(ply, "element vertex %d\n",idx );fprintf(ply, "property float x\n");fprintf(ply, "property float y\n");fprintf(ply, "property float z\n");fprintf(ply, "element face %d\n", count);fprintf(ply, "property list uint8 int32 vertex_indices\n");fprintf(ply, "end_header\n");for (CDT::Vertex_iterator v = cdt.vertices_begin(); v != cdt.vertices_end(); ++v){Vertex_handle vv = v->handle();double x = vv->point().x();double y = vv->point().y();fprintf(ply, "%f %f %f\n", x, y, 0.0);}for (CDT::Finite_faces_iterator fit = cdt.finite_faces_begin();    fit != cdt.finite_faces_end(); ++fit){if (fit->info().in_domain()){Vertex_handle vertId0 = fit->vertex(0);Vertex_handle vertId1 = fit->vertex(1);Vertex_handle vertId2 = fit->vertex(2);int id0 = vertId0->id();int id1 = vertId1->id();int id2 = vertId2->id();fprintf(ply, "%d %d %d %d\n", 3, id0, id1, id2);}}}std::cout << "There are " << count << " facets in the domain." << std::endl;system("pause");return 0;
}

效果图