сохранение CGAL альфа-формы поверхности сетки

Вот мой код, который выводит ВТК файл для визуализации в ParaView. По сравнению с решениями 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();
}