Saving CGAL alpha surface shape mesh

Question

I never used CGAL and had little or no experience in C / C ++. But following Google, however, I managed to put together the example "Alpha_shapes_3" (\ CGAL-4.1-beta1 \ examples \ Alpha_shapes_3) on a 64-bit Windows 7 machine using visual studio 2010.

Now, if we check the source code of the program "ex_alpha_shapes_3", we notice that the data file called "bunny_1000" is red, where the 3D point is a cluster. Now my question is: how to change the source code so that after the alpha form is calculated for the given points, the surface mesh of the alpha form is saved / written in an external file. It can be just a list of polygons and their corresponding three-dimensional vertices. I think that these polygons will determine the alpha surface mesh. If I can do this, I can see the output of the alpha form generation program in an external tool that I am familiar with.

I know this is very simple, but I could not figure it out with my limited knowledge of CGAL. A.

I know you have gueys code, but I am pasting it again to complete.

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <CGAL/Delaunay_triangulation_3.h> #include <CGAL/Alpha_shape_3.h> #include <fstream> #include <list> #include <cassert> typedef CGAL::Exact_predicates_inexact_constructions_kernel Gt; typedef CGAL::Alpha_shape_vertex_base_3<Gt> Vb; typedef CGAL::Alpha_shape_cell_base_3<Gt> Fb; typedef CGAL::Triangulation_data_structure_3<Vb,Fb> Tds; typedef CGAL::Delaunay_triangulation_3<Gt,Tds> Triangulation_3; typedef CGAL::Alpha_shape_3<Triangulation_3> Alpha_shape_3; typedef Gt::Point_3 Point; typedef Alpha_shape_3::Alpha_iterator Alpha_iterator; int main() { std::list<Point> lp; //read input std::ifstream is("./data/bunny_1000"); int n; is >> n; std::cout << "Reading " << n << " points " << std::endl; Point p; for( ; n>0 ; n--) { is >> p; lp.push_back(p); } // compute alpha shape Alpha_shape_3 as(lp.begin(),lp.end()); std::cout << "Alpha shape computed in REGULARIZED mode by default" << std::endl; // find optimal alpha value Alpha_iterator opt = as.find_optimal_alpha(1); std::cout << "Optimal alpha value to get one connected component is " << *opt << std::endl; as.set_alpha(*opt); assert(as.number_of_solid_components() == 1); return 0; } 

After much searching on the Internet, I found that we probably need to use something like

 std::list<Facet> facets; alpha_shape.get_alpha_shape_facets ( std::back_inserter(facets),Alpha_shape::REGULAR ); 

But I still don’t know at all how to use this in the above code!

Answer 1

As described here , a cell is a pair (Cell_handle c, int i), defined as a cell in c, opposite the vertex of index i. On this page, you indicated how the vertex indices of the cell are indicated.

In the following code example, I added a little output that prints the OFF file to cout by duplicating the vertices. To do something clean, you can either use std::map<Alpha_shape_3::Vertex_handle,int> to associate a unique index to the vertex, or add information to the vertices, as in these examples .

 /// collect all regular facets std::vector<Alpha_shape_3::Facet> facets; as.get_alpha_shape_facets(std::back_inserter(facets), Alpha_shape_3::REGULAR); std::stringstream pts; std::stringstream ind; std::size_t nbf=facets.size(); for (std::size_t i=0;i<nbf;++i) { //To have a consistent orientation of the facet, always consider an exterior cell if ( as.classify( facets[i].first )!=Alpha_shape_3::EXTERIOR ) facets[i]=as.mirror_facet( facets[i] ); CGAL_assertion( as.classify( facets[i].first )==Alpha_shape_3::EXTERIOR ); int indices[3]={ (facets[i].second+1)%4, (facets[i].second+2)%4, (facets[i].second+3)%4, }; /// according to the encoding of vertex indices, this is needed to get /// a consistent orienation if ( facets[i].second%2==0 ) std::swap(indices[0], indices[1]); pts << facets[i].first->vertex(indices[0])->point() << "\n" << facets[i].first->vertex(indices[1])->point() << "\n" << facets[i].first->vertex(indices[2])->point() << "\n"; ind << "3 " << 3*i << " " << 3*i+1 << " " << 3*i+2 << "\n"; } std::cout << "OFF "<< 3*nbf << " " << nbf << " 0\n"; std::cout << pts.str(); std::cout << ind.str(); 

Answer 2

Here is my code that outputs a vtk file for rendering in Paraview . Compared to slorior solutions, duplicate points are not saved in the file. But my code is for visualization only, if you need to figure out external or internal simplexes, you must change the code to get these results.

 void writevtk(Alpha_shape_3 &as, const std::string &asfile) { // http://cgal-discuss.949826.n4.nabble.com/Help-with-filtration-and-filtration-with-alpha-values-td4659524.html#a4659549 std::cout << "Information of the Alpha_Complex:\n"; std::vector<Alpha_shape_3::Cell_handle> cells; std::vector<Alpha_shape_3::Facet> facets; std::vector<Alpha_shape_3::Edge> edges; // tetrahedron = cell, they should be the interior, it is inside the 3D space as.get_alpha_shape_cells(std::back_inserter(cells), Alpha_shape_3::INTERIOR); // triangles // for the visualiization, don't need regular because tetrahedron will show it //as.get_alpha_shape_facets(std::back_inserter(facets), Alpha_shape_3::REGULAR); as.get_alpha_shape_facets(std::back_inserter(facets), Alpha_shape_3::SINGULAR); // edges as.get_alpha_shape_edges(std::back_inserter(edges), Alpha_shape_3::SINGULAR); std::cout << "The alpha-complex has : " << std::endl; std::cout << cells.size() << " cells as tetrahedrons" << std::endl; std::cout << facets.size() << " triangles" << std::endl; std::cout << edges.size() << " edges" << std::endl; size_t tetra_num, tri_num, edge_num; tetra_num = cells.size(); tri_num = facets.size(); edge_num = edges.size(); // vertices: points <-> id std::map<Point, size_t> points; size_t index = 0; // finite_.. is from DT class for (auto v_it = as.finite_vertices_begin(); v_it != as.finite_vertices_end(); v_it++) { points[v_it->point()] = index; index++; } // write std::ofstream of(asfile); of << "# vtk DataFile Version 2.0\n\nASCII\nDATASET UNSTRUCTURED_GRID\n\n"; of << "POINTS " << index << " float\n"; for (auto v_it = as.finite_vertices_begin(); v_it != as.finite_vertices_end(); v_it++) { of << v_it->point() << std::endl; } of << std::endl; of << "CELLS " << tetra_num + tri_num + edge_num << " " << 5 * tetra_num + 4 * tri_num + 3 * edge_num << std::endl; for (auto cell:cells) { size_t v0 = points.find(cell->vertex(0)->point())->second; size_t v1 = points.find(cell->vertex(1)->point())->second; size_t v2 = points.find(cell->vertex(2)->point())->second; size_t v3 = points.find(cell->vertex(3)->point())->second; of << "4 " << v0 << " " << v1 << " " << v2 << " " << v3 << std::endl; } // https://doc.cgal.org/latest/TDS_3/classTriangulationDataStructure__3.html#ad6a20b45e66dfb690bfcdb8438e9fcae for (auto tri_it = facets.begin(); tri_it != facets.end(); ++tri_it) { of << "3 "; auto tmp_tetra = tri_it->first; for (int i = 0; i < 4; i++) { if (i != tri_it->second) { of << points.find(tmp_tetra->vertex(i)->point())->second << " "; } } of << std::endl; } // https://doc.cgal.org/latest/TDS_3/classTriangulationDataStructure__3.html#af31db7673a6d7d28c0bb90a3115ac695 for (auto e : edges) { of << "2 "; auto tmp_tetra = e.get<0>(); int p1, p2; p1 = e.get<1>(); p2 = e.get<2>(); of << points.find(tmp_tetra->vertex(p1)->point())->second << " " << points.find(tmp_tetra->vertex(p2)->point())->second << std::endl; } of << std::endl; of << "CELL_TYPES " << tetra_num + tri_num + edge_num << std::endl; for (int i = 0; i < tetra_num; i++) { of << "10 "; } for (int i = 0; i < tri_num; i++) { of << "5 "; } for (int i = 0; i < edge_num; i++) { of << "3 "; } of << std::endl; of.close(); }