c++ - guardar la malla de superficie de forma alfa de CGAL

download cgal

Question

Nunca he usado CGAL y casi no tengo experiencia en C / C ++. Pero después de Google, sin embargo, pude compilar el ejemplo "Alpha_shapes_3" (/ CGAL-4.1-beta1 / examples / Alpha_shapes_3) en una máquina con Windows 7 de 64 bits utilizando Visual Studio 2010.

Ahora, si comprobamos el código fuente del programa "ex_alpha_shapes_3", notamos que un archivo de datos llamado "bunny_1000" es rojo donde reside el clúster de puntos 3d. Ahora mi pregunta es cómo puedo cambiar el código fuente para que después de que la forma alfa se compute para los puntos dados, la malla de la superficie de la forma alfa se guarde / escriba en un archivo externo. Puede ser simplemente la lista de polígonos y sus respectivos vértices 3D. Supongo que estos polígonos definirán la malla de superficie de la forma alfa. Si puedo hacer eso, puedo ver el resultado del programa de generación de forma alfa en una herramienta externa con la que estoy familiarizado.

Sé que esto es muy sencillo, pero no pude resolver esto con mi conocimiento limitado de CGAL.

Sé que los gueys tienen el código, pero lo vuelvo a pegar para completarlo.

#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;
} 

Después de buscar mucho en internet, descubrí que probablemente necesitábamos algo como

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

¡Pero todavía no tengo ni idea de cómo usar esto en el código anterior!

Answer 1

Como se documenta aquí , una faceta es un par (Cell_handle c, int i) definido como la faceta en c opuesta al vértice del índice i. En esta página , tiene la descripción de cómo son los índices de vértice de una celda.

En el siguiente ejemplo de código, agregué un pequeño resultado que imprime un archivo OFF en cout duplicando los vértices. Para hacer algo limpio, puede usar std::map<Alpha_shape_3::Vertex_handle,int> para asociar un índice único por vértice o agregar una información a los vértices como en esos ejemplos .

/// 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

Aquí está mi código, que muestra el archivo vtk para visualización en Paraview . Comparando con las soluciones de slorior, no se guardan puntos duplicados en el archivo. Pero mi código es solo para la visualización, si necesita averiguar los símplex exteriores o interiores, debe modificar el código para obtener estos resultados.

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();
}