【CGAL】测试surface_mesh的使用

近几日研究了CGAL的surface_mesh生成和分割，

有些注意事项：

1.网格的三角网要一个接一个邻接着插入，不断验证定向（orient），不然最后可能就有的三角形缺失了。

2.采用四个点的面元插入，实际并不是三角形，无法进行Polygon_mesh_processing中的分割（split）

  

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/clip.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polygon_mesh_processing/transform.h>
#include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
#include <CGAL/boost/graph/Face_filtered_graph.h>
#include <CGAL/Polygon_mesh_processing/fair.h>
#include <boost/property_map/property_map.hpp>
#include <CGAL/convex_hull_2.h>

#include <iostream>
#include <fstream>
#include <sstream>

namespace PMP = CGAL::Polygon_mesh_processing;
namespace params = PMP::parameters;

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3> Surface_mesh;
typedef CGAL::Polyhedron_3<K> Polyhedron;

typedef Surface_mesh::Vertex_index vertex_descriptor;
typedef Surface_mesh::Face_index face_descriptor;

#define dot(u,v)   ((u).x() * (v).x() + (u).y() * (v).y() + (u).z() * (v).z())
#define norm(v)    sqrt(dot(v,v))  // norm = length of vector
#define d(u,v)     norm(u-v)       // distance = norm of difference

// pbase_Plane(): get base of perpendicular from point to a plane
//    Input:  P = a 3D point
//            PL = a plane with point V0 and normal n
//    Output: *B = base point on PL of perpendicular from P
//    Return: the distance from P to the plane PL
float pbase_Plane(K::Point_3 P, K::Plane_3 PL, K::Point_3* B)
{
    float    sb, sn, sd;
    K::Vector_3 n=PL.orthogonal_vector();
    K::Vector_3 t_p_v0 = (P - PL.point());
    sn = -dot(n, t_p_v0);
    sd = dot(n, n);
    sb = sn / sd;

    *B = P + sb * n;
    return d(P, *B);
}

int main()
{
    Surface_mesh tm1;
    std::ifstream input(CGAL::data_file_path("data/Box_stitched.off"));
    input >> tm1;

    if (!input)
    {
        std::cerr << "File not found. Aborting." << std::endl;
        assert(false);
        return 0;
    }
    input.close();
    //构造平面
    std::vector<K::Plane_3> m_planes;
    K::Plane_3 east(1, 0, 0, 0.5);
    K::Plane_3 north(0, 1, 0, -0.5);
    K::Plane_3 weast(1, 0, 0, -0.5);
    K::Plane_3 south(0, 1, 0, 0.5);
    K::Plane_3 up(0, 0, 1, -0.5);
    K::Plane_3 down(0, 0, 1, 0.5);

    K::Plane_3 cut_plane(0, 0, 1, 0);
    m_planes.push_back(east);
    m_planes.push_back(north);
    m_planes.push_back(down);
    m_planes.push_back(south);
    m_planes.push_back(up);
    m_planes.push_back(weast);
    m_planes.push_back(cut_plane);

    //PMP::stitch_boundary_cycle(mesh);
    PMP::stitch_borders(tm1, CGAL::parameters::vertex_point_map(get(CGAL::vertex_point, tm1)));
    PMP::split(tm1, cut_plane, params::use_compact_clipper(true));
    std::vector<Surface_mesh> meshes;
    PMP::split_connected_components(tm1, meshes, params::all_default());
    CGAL::clear(tm1);
    Surface_mesh mesh1 = meshes[0];//分割为2个多面体
    Surface_mesh mesh2 = meshes[1];
    Surface_mesh::Vertex_range r = mesh1.vertices();
    Surface_mesh m_NewMeshLeft;
    for (vertex_descriptor vd : mesh1.vertices())
    {
        K::Point_3 pt = mesh1.point(vd);
        m_NewMeshLeft.add_vertex(pt);
    }
    Surface_mesh::Vertex_range::iterator  vb, ve;
    vb = r.begin();
    ve = r.end();
    for (size_t id_p = 0; id_p < m_planes.size(); id_p++)
    {
        K::Plane_3 planeX= m_planes[id_p];
        std::list<K::Point_3> pts;
        std::size_t idx = 0;
        std::vector<K::Point_3> intersecting_points;
        std::vector<vertex_descriptor> intersecting_VertexIndexs;
        for (vertex_descriptor vd : m_NewMeshLeft.vertices())
        {
            K::Point_3 pt = m_NewMeshLeft.point(vd);
            K::Point_3 basePt;
            float dist = pbase_Plane(pt, planeX, &basePt);//靠近切平面的点
            if (dist < 0.001)
            {
                intersecting_points.push_back(pt);
                intersecting_VertexIndexs.push_back(vd);
                const K::Point_2& q = planeX.to_2d(pt);
                pts.push_back(K::Point_3(q.x(), q.y(), idx));
                idx++;
            }
        }
        if (pts.size()==0)
        {
            continue;
        }
        typedef CGAL::Projection_traits_xy_3<K>  Projection;
        std::list<K::Point_3> hull;
        //创建凸包围盒，平面与bbox边交点围成的凸包
        CGAL::convex_hull_2(pts.begin(), pts.end(), std::back_inserter(hull), Projection());
        std::vector<K::Point_3> polyline;
        std::vector<vertex_descriptor> ch_VertexIndexs;
        std::vector< std::set<const K::Plane_3*> > ch_source_planes;
        for (typename std::list<K::Point_3>::iterator it = hull.begin(); it != hull.end(); ++it)
        {
            std::size_t idx = std::size_t(it->z());
            polyline.push_back(intersecting_points[idx]);
            ch_VertexIndexs.push_back(intersecting_VertexIndexs[idx]);
            //ch_source_planes.push_back(intersecting_points_source_planes[idx]);
        }
        /*typename std::list<K::Point_3>::iterator it0 = hull.begin();
        std::size_t idx0 = std::size_t(it0->z());
        polyline.push_back(intersecting_points[idx0]);
        ch_VertexIndexs.push_back(intersecting_VertexIndexs[idx0]);*/
        //face_descriptor f = m_NewMeshLeft.add_face(ch_VertexIndexs[0], ch_VertexIndexs[1], ch_VertexIndexs[2], ch_VertexIndexs[3]);
        //if (f == Surface_mesh::null_face())
        //{
        //    //std::cerr << "The face could not be added because of an orientation error." << std::endl;
        //    f = m_NewMeshLeft.add_face(ch_VertexIndexs[3], ch_VertexIndexs[2], ch_VertexIndexs[1], ch_VertexIndexs[0]);
        //}
        // any type, having Type(int, int, int) constructor available, can be used to hold output triangles
        typedef CGAL::Triple<int, int, int> Triangle_int;
        std::vector<Triangle_int> patch;
        patch.reserve(polyline.size() - 2); // there will be exactly n-2 triangles in the patch
        CGAL::Polygon_mesh_processing::triangulate_hole_polyline(polyline, std::back_inserter(patch));
        for (std::size_t i = 0; i < patch.size(); ++i)
        {
            std::cout << "Triangle " << i << ": " << patch[i].first << " " << patch[i].second << " " << patch[i].third
                << std::endl;
            vertex_descriptor v1 = ch_VertexIndexs[patch[i].first];
            vertex_descriptor v2 = ch_VertexIndexs[patch[i].second];
            vertex_descriptor v3 = ch_VertexIndexs[patch[i].third];

            std::cout << "Triangle " << i << "OO: " << v1 << " " << v2 << " " << v3<< std::endl;
            face_descriptor f =m_NewMeshLeft.add_face( v1,v2, v3);
            if (f == Surface_mesh::null_face())
            {
                std::cerr << "The face could not be added because of an orientation error." << std::endl;
                f = m_NewMeshLeft.add_face(v1, v3, v2);
                //assert(f != Mesh::null_face());
            }
        }
    }
    
    //PMP::fair(meshCutPlane, ch_VertexIndexs);
    CGAL::IO::write_polygon_mesh("results/NewSurfaceMesh.off", m_NewMeshLeft, CGAL::parameters::stream_precision(17));
    meshes.clear();

    Surface_mesh mesh0000;
    const std::string filename000 = CGAL::data_file_path("results/NewSurfaceMesh.off");
    if (!PMP::IO::read_polygon_mesh(filename000, mesh0000))
    {
        std::cerr << "Invalid input." << std::endl;
        return 1;
    }
    PMP::stitch_borders(mesh0000);
    PMP::split(mesh0000, K::Plane_3(0, 1.2, 0, 0));

    std::vector<Surface_mesh> meshes2;
    PMP::split_connected_components(mesh0000, meshes2, params::all_default());
    Surface_mesh mesh122 = meshes2[0];
    CGAL::IO::write_polygon_mesh("results/HalfBox2.off", mesh122, CGAL::parameters::stream_precision(17));
    return EXIT_SUCCESS;
}