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