保存CGAL alpha形状表面网格

Question

我从来没有使用过CGAL，几乎没有C/C++的经验。但是，随着谷歌 我已经然而设法编译示例 “Alpha_shapes_3” （\ CGAL-4.1-β1\例子\ Alpha_shapes_3）使用 Visual Studio 2010中

现在Windows 7的64位计算机上如果我们检查程序“ex_alpha_shapes_3”的源代码，我们 注意到名为“bunny_1000”的数据文件在红点位于集群所在的3d点 处。 现在我的问题是如何更改源代码，以便在为给定点计算阿尔法 形状后，阿尔法形状的表面网格将在外部文件中保存/写入 。它可以是简单的多边形列表和它们各自的3D顶点。我想这些多边形将定义alpha形状的表面网格。如果我可以这样做，我可以在我熟悉的外部工具中看到 alpha形状生成程序的输出。

我知道这很直接，但我不知道我的 有限的CGAL知识。

我知道你有代码，但我粘贴它再次完成。

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

在网上搜索了很多后，我发现，也许我们需要使用像

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

但是，我还是完全无能如何在上面的代码中使用此！

Answer 1

如文献here所述，小平面是一对（Cell_handle c，int i），其定义为c中与索引i的顶点相对的小平面。 在this page上，您已经描述了单元格的顶点索引。

在以下代码示例中，我添加了一个小型输出，通过复制顶点在cout上打印OFF文件。要做一些干净的事情，你可以使用std::map<Alpha_shape_3::Vertex_handle,int>来关联每个顶点的唯一索引或者向顶点添加一个信息，如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(); 

Answer 2

这是我的代码，其输出在Paraviewvtk文件用于可视化。与slorior的解决方案相比，文件中没有重复的点。但是我的代码只是为了可视化，如果你需要弄清楚外部或内部单形，你应该修改代码来获得这些结果。

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