【CGAL】测试surface_mesh的使用

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

有些注意事项:

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

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

  

 

 

  1 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
  2 #include <CGAL/Polygon_mesh_processing/clip.h>
  3 #include <CGAL/Surface_mesh.h>
  4 #include <CGAL/Polyhedron_3.h>
  5 #include <CGAL/Polygon_mesh_processing/transform.h>
  6 #include <CGAL/Polygon_mesh_processing/IO/polygon_mesh_io.h>
  7 #include <CGAL/boost/graph/Face_filtered_graph.h>
  8 #include <CGAL/Polygon_mesh_processing/fair.h>
  9 #include <boost/property_map/property_map.hpp>
 10 #include <CGAL/convex_hull_2.h>
 11 
 12 #include <iostream>
 13 #include <fstream>
 14 #include <sstream>
 15 
 16 namespace PMP = CGAL::Polygon_mesh_processing;
 17 namespace params = PMP::parameters;
 18 
 19 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 20 typedef CGAL::Surface_mesh<K::Point_3> Surface_mesh;
 21 typedef CGAL::Polyhedron_3<K> Polyhedron;
 22 
 23 typedef Surface_mesh::Vertex_index vertex_descriptor;
 24 typedef Surface_mesh::Face_index face_descriptor;
 25 
 26 #define dot(u,v)   ((u).x() * (v).x() + (u).y() * (v).y() + (u).z() * (v).z())
 27 #define norm(v)    sqrt(dot(v,v))  // norm = length of vector
 28 #define d(u,v)     norm(u-v)       // distance = norm of difference
 29 
 30 // pbase_Plane(): get base of perpendicular from point to a plane
 31 //    Input:  P = a 3D point
 32 //            PL = a plane with point V0 and normal n
 33 //    Output: *B = base point on PL of perpendicular from P
 34 //    Return: the distance from P to the plane PL
 35 float pbase_Plane(K::Point_3 P, K::Plane_3 PL, K::Point_3* B)
 36 {
 37     float    sb, sn, sd;
 38     K::Vector_3 n=PL.orthogonal_vector();
 39     K::Vector_3 t_p_v0 = (P - PL.point());
 40     sn = -dot(n, t_p_v0);
 41     sd = dot(n, n);
 42     sb = sn / sd;
 43 
 44     *B = P + sb * n;
 45     return d(P, *B);
 46 }
 47 
 48 int main()
 49 {
 50     Surface_mesh tm1;
 51     std::ifstream input(CGAL::data_file_path("data/Box_stitched.off"));
 52     input >> tm1;
 53 
 54     if (!input)
 55     {
 56         std::cerr << "File not found. Aborting." << std::endl;
 57         assert(false);
 58         return 0;
 59     }
 60     input.close();
 61     //构造平面
 62     std::vector<K::Plane_3> m_planes;
 63     K::Plane_3 east(1, 0, 0, 0.5);
 64     K::Plane_3 north(0, 1, 0, -0.5);
 65     K::Plane_3 weast(1, 0, 0, -0.5);
 66     K::Plane_3 south(0, 1, 0, 0.5);
 67     K::Plane_3 up(0, 0, 1, -0.5);
 68     K::Plane_3 down(0, 0, 1, 0.5);
 69 
 70     K::Plane_3 cut_plane(0, 0, 1, 0);
 71     m_planes.push_back(east);
 72     m_planes.push_back(north);
 73     m_planes.push_back(down);
 74     m_planes.push_back(south);
 75     m_planes.push_back(up);
 76     m_planes.push_back(weast);
 77     m_planes.push_back(cut_plane);
 78 
 79     //PMP::stitch_boundary_cycle(mesh);
 80     PMP::stitch_borders(tm1, CGAL::parameters::vertex_point_map(get(CGAL::vertex_point, tm1)));
 81     PMP::split(tm1, cut_plane, params::use_compact_clipper(true));
 82     std::vector<Surface_mesh> meshes;
 83     PMP::split_connected_components(tm1, meshes, params::all_default());
 84     CGAL::clear(tm1);
 85     Surface_mesh mesh1 = meshes[0];//分割为2个多面体
 86     Surface_mesh mesh2 = meshes[1];
 87     Surface_mesh::Vertex_range r = mesh1.vertices();
 88     Surface_mesh m_NewMeshLeft;
 89     for (vertex_descriptor vd : mesh1.vertices())
 90     {
 91         K::Point_3 pt = mesh1.point(vd);
 92         m_NewMeshLeft.add_vertex(pt);
 93     }
 94     Surface_mesh::Vertex_range::iterator  vb, ve;
 95     vb = r.begin();
 96     ve = r.end();
 97     for (size_t id_p = 0; id_p < m_planes.size(); id_p++)
 98     {
 99         K::Plane_3 planeX= m_planes[id_p];
100         std::list<K::Point_3> pts;
101         std::size_t idx = 0;
102         std::vector<K::Point_3> intersecting_points;
103         std::vector<vertex_descriptor> intersecting_VertexIndexs;
104         for (vertex_descriptor vd : m_NewMeshLeft.vertices())
105         {
106             K::Point_3 pt = m_NewMeshLeft.point(vd);
107             K::Point_3 basePt;
108             float dist = pbase_Plane(pt, planeX, &basePt);//靠近切平面的点
109             if (dist < 0.001)
110             {
111                 intersecting_points.push_back(pt);
112                 intersecting_VertexIndexs.push_back(vd);
113                 const K::Point_2& q = planeX.to_2d(pt);
114                 pts.push_back(K::Point_3(q.x(), q.y(), idx));
115                 idx++;
116             }
117         }
118         if (pts.size()==0)
119         {
120             continue;
121         }
122         typedef CGAL::Projection_traits_xy_3<K>  Projection;
123         std::list<K::Point_3> hull;
124         //创建凸包围盒,平面与bbox边交点围成的凸包
125         CGAL::convex_hull_2(pts.begin(), pts.end(), std::back_inserter(hull), Projection());
126         std::vector<K::Point_3> polyline;
127         std::vector<vertex_descriptor> ch_VertexIndexs;
128         std::vector< std::set<const K::Plane_3*> > ch_source_planes;
129         for (typename std::list<K::Point_3>::iterator it = hull.begin(); it != hull.end(); ++it)
130         {
131             std::size_t idx = std::size_t(it->z());
132             polyline.push_back(intersecting_points[idx]);
133             ch_VertexIndexs.push_back(intersecting_VertexIndexs[idx]);
134             //ch_source_planes.push_back(intersecting_points_source_planes[idx]);
135         }
136         /*typename std::list<K::Point_3>::iterator it0 = hull.begin();
137         std::size_t idx0 = std::size_t(it0->z());
138         polyline.push_back(intersecting_points[idx0]);
139         ch_VertexIndexs.push_back(intersecting_VertexIndexs[idx0]);*/
140         //face_descriptor f = m_NewMeshLeft.add_face(ch_VertexIndexs[0], ch_VertexIndexs[1], ch_VertexIndexs[2], ch_VertexIndexs[3]);
141         //if (f == Surface_mesh::null_face())
142         //{
143         //    //std::cerr << "The face could not be added because of an orientation error." << std::endl;
144         //    f = m_NewMeshLeft.add_face(ch_VertexIndexs[3], ch_VertexIndexs[2], ch_VertexIndexs[1], ch_VertexIndexs[0]);
145         //}
146         // any type, having Type(int, int, int) constructor available, can be used to hold output triangles
147         typedef CGAL::Triple<int, int, int> Triangle_int;
148         std::vector<Triangle_int> patch;
149         patch.reserve(polyline.size() - 2); // there will be exactly n-2 triangles in the patch
150         CGAL::Polygon_mesh_processing::triangulate_hole_polyline(polyline, std::back_inserter(patch));
151         for (std::size_t i = 0; i < patch.size(); ++i)
152         {
153             std::cout << "Triangle " << i << ": " << patch[i].first << " " << patch[i].second << " " << patch[i].third
154                 << std::endl;
155             vertex_descriptor v1 = ch_VertexIndexs[patch[i].first];
156             vertex_descriptor v2 = ch_VertexIndexs[patch[i].second];
157             vertex_descriptor v3 = ch_VertexIndexs[patch[i].third];
158 
159             std::cout << "Triangle " << i << "OO: " << v1 << " " << v2 << " " << v3<< std::endl;
160             face_descriptor f =m_NewMeshLeft.add_face( v1,v2, v3);
161             if (f == Surface_mesh::null_face())
162             {
163                 std::cerr << "The face could not be added because of an orientation error." << std::endl;
164                 f = m_NewMeshLeft.add_face(v1, v3, v2);
165                 //assert(f != Mesh::null_face());
166             }
167         }
168     }
169     
170     //PMP::fair(meshCutPlane, ch_VertexIndexs);
171     CGAL::IO::write_polygon_mesh("results/NewSurfaceMesh.off", m_NewMeshLeft, CGAL::parameters::stream_precision(17));
172     meshes.clear();
173 
174     Surface_mesh mesh0000;
175     const std::string filename000 = CGAL::data_file_path("results/NewSurfaceMesh.off");
176     if (!PMP::IO::read_polygon_mesh(filename000, mesh0000))
177     {
178         std::cerr << "Invalid input." << std::endl;
179         return 1;
180     }
181     PMP::stitch_borders(mesh0000);
182     PMP::split(mesh0000, K::Plane_3(0, 1.2, 0, 0));
183 
184     std::vector<Surface_mesh> meshes2;
185     PMP::split_connected_components(mesh0000, meshes2, params::all_default());
186     Surface_mesh mesh122 = meshes2[0];
187     CGAL::IO::write_polygon_mesh("results/HalfBox2.off", mesh122, CGAL::parameters::stream_precision(17));
188     return EXIT_SUCCESS;
189 }

 

posted @ 2022-04-01 15:59  太一吾鱼水  阅读(1298)  评论(0编辑  收藏  举报