6

我从未使用过 CGAL，几乎没有 C/C++ 经验。但是按照谷歌，我设法在 Windows 7 64 位机器上使用 Visual Studio 2010 编译示例“Alpha_shapes_3”（\CGAL-4.1-beta1\examples\Alpha_shapes_3）。

在此处输入图像描述

现在，如果我们检查程序“ex_alpha_shapes_3”的源代码，我们会注意到一个名为“bunny_1000”的数据文件在 3d 点簇所在的位置是红色的。现在我的问题是如何更改源代码，以便在为给定点计算 alpha 形状后，将 alpha 形状的表面网格保存/写入外部文件中。它可以是简单的多边形列表和它们各自的 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
);

但是我仍然完全不知道如何在上面的代码中使用它！

4

2 回答 2

9

如此处所述，构面是一对 (Cell_handle c,int i)，定义为 c 中与索引 i 的顶点相对的构面。在此页面上，您可以了解单元格的顶点索引是如何的。

在下面的代码示例中，我添加了一个小输出，cout通过复制顶点来打印 OFF 文件。要做一些干净的事情，您可以使用 a为每个顶点关联一个唯一索引，或者像在这些示例std::map<Alpha_shape_3::Vertex_handle,int>中那样向顶点添加信息。

/// 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();
于 2013-04-10T05:42:00.633 回答
0

这是我的代码，它输出vtk文件以在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();
}
于 2017-12-08T21:10:07.507 回答