saving CGAL alpha shape surface mesh

Question

I have never used CGAL and have got almost no C/C++ experience. But following Google I have however managed to compile the example "Alpha_shapes_3" (\CGAL-4.1-beta1\examples\Alpha_shapes_3) on a Windows 7 64bit machine using visual studio 2010.

Now if we check the source code for the program "ex_alpha_shapes_3" we notice that a data file called "bunny_1000" is red where the 3d point cluster resides. Now my question is how can I change the source code so that after the alpha shape is computed for the given points, surface mesh of the alpha shape is saved/wrote in an external file. It can be simply the list of polygons and their respective 3D vertices. I guess these polygons will be defining the surface mesh of the alpha shape. If I can do that I can see the output of the alpha shape generation program in an external tool I am familiar with.

I know this is very straightforward but I could not figure this out with my limited knowledge of CGAL.

I know you gueys have the code but I am pasting it again for completion.

#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 searching a lot in the internet I found that probably we need to use something like

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

But I am still completely clueless how to use this in the above code!

Answer 1

Here is my code, which outputs vtk file for visualization in Paraview. Comparing with slorior's solutions, no duplicated points are saved in the file. But my code is just for the visualization, if you need to figure out the exterior or interior simplexes, you should modify 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();
}

Answer 2

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

In the following code sample, I added a small output that prints an OFF file on cout by duplicating the vertices. To do something clean, you can either use a std::map<Alpha_shape_3::Vertex_handle,int> to associate a unique index per vertex or add an info to the vertices like in those 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();