[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[])
{
    //读取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:%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 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, "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;
}