[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的泛型机制真的很强大，拓展性很好。
代码语言：javascript复制// 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[])
{
//读取shp
const 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:%dn", nShapeType);
printf("Entities:%dn", 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:%dn", parts);
printf("nVertices:%dn", 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 polygon
mark_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, "plynformat %s 1.0n", "ascii");
fprintf(ply, "element vertex %dn",idx );
fprintf(ply, "property float xn");
fprintf(ply, "property float yn");
fprintf(ply, "property float zn");
fprintf(ply, "element face %dn", count);
fprintf(ply, "property list uint8 int32 vertex_indicesn");
fprintf(ply, "end_headern");
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 %fn", 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 %dn", 3, id0, id1, id2);
}
}
}
std::cout << "There are " << count << " facets in the domain." << std::endl;
system("pause");
return 0;
}效果图
