CGAL and the Boost Graph Library

265 阅读8分钟

CGAL and the Boost Graph Library

许多几何数据结构都可以解释为图,因为它们由顶点和边组成。对于halfedge数据结构、多面体曲面、arrangement以及二维三角剖分类来说,情况都是如此。利用对偶性,人们也可以将面解释为顶点,相邻面之间的边解释为对偶图的边。 CGAL的研究范围是几何算法而不是图算法。尽管如此,这个包提供了必要的类和函数,使Boost Graph Library(简称BGL)的算法能够使用CGAL数据结构。 此外,该软件包扩展了BGL,引入了 HalfedgeGraphFaceGraph 等概念,允许处理半边和面。这些概念反映了第6章中描述的半边数据结构的设计,即对偶半边和围绕顶点和面的圆形半边序列。

属性映射

BGL中广泛使用的另一个特性是Boost Property Map Library提供的属性映射。属性映射是一个通用的接口,用于将键对象映射到对应的值对象。 BGL使用属性映射将信息与顶点和边关联起来。该机制使用traits类( boost::property_traits )和自由函数来读取( get )和写入( put )信息,包括顶点、边以及半边和面。例如,BGL Dijksta的最短路径算法在这样的属性映射中写入每个顶点的前驱节点以及到源节点的距离。 一些默认的属性映射与图形类型相关联。它们被称为内部属性映射(internal property map),可以通过 get() 函数的重载来获取。

Surface_mesh类作为Boost Graph 概念的模型

Surface_mesh 类是BGLCGAL提供的大多数图概念的一个模型。完整的列表可以在boost::graph_traits的文档中找到。这些例子展示了如何使用 Surface_mesh 的一些BGL算法,以及如何使用CGAL提供的概念来实现一个简单的算法。

Surface_mesh的最小生成树MST

下面的示例程序计算曲面网格上的最小生成树。更多的例子可以在三角网格简化、三角网格分割和三角网格变形章节中找到。 surface mesh类使用整数索引来寻址顶点和边,它还带有一个内置的属性机制,可以很好地映射到BGL上。

最小生成树算法在属性映射中写入每个顶点的前驱节点。

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>

#include <iostream>
#include <fstream>

#include <CGAL/boost/graph/prim_minimum_spanning_tree.h>

typedef CGAL::Simple_cartesian<double>                       Kernel;
typedef Kernel::Point_3                                      Point;
typedef CGAL::Surface_mesh<Point>                            Mesh;

typedef boost::graph_traits<Mesh>::vertex_descriptor vertex_descriptor;

int main(int argc, char* argv[])
{
  const std::string filename = (argc>1) ? argv[1] : CGAL::data_file_path("meshes/prim.off");

  Mesh P;
  if(!CGAL::IO::read_polygon_mesh(filename, P))
  {
    std::cerr << "Invalid input." << std::endl;
    return 1;
  }
 //
  Mesh::Property_map<vertex_descriptor,vertex_descriptor> predecessor;
  predecessor = P.add_property_map<vertex_descriptor,vertex_descriptor>("v:predecessor").first;

  boost::prim_minimum_spanning_tree(P, predecessor, boost::root_vertex(*vertices(P).first));

  std::cout << "#VRML V2.0 utf8\n"
    "DirectionalLight {\n"
    "direction 0 -1 0\n"
    "}\n"
    "Shape {\n"
    "  appearance Appearance {\n"
    "    material Material { emissiveColor 1 0 0}}\n"
    "    geometry\n"
    "    IndexedLineSet {\n"
    "      coord Coordinate {\n"
    "        point [ \n";

  for(vertex_descriptor vd : vertices(P)){
    std::cout <<  "        " << P.point(vd) << "\n";
  }

  std::cout << "        ]\n"
    "     }\n"
    "      coordIndex [\n";
  for(vertex_descriptor vd : vertices(P)){
    if(predecessor[vd]!=vd){
      std::cout << "      " << std::size_t(vd) << ", " << std::size_t(predecessor[vd]) <<  ", -1\n";
    }
  }

  std::cout << "]\n"
    "  }#IndexedLineSet\n"
    "}# Shape\n";

  P.remove_property_map(predecessor);
  return 0;
}

meshes/prim.off网格曲面: prim 网格曲面的最小生成树,写入VRML格式的: 在这里插入图片描述 VRML结果文件使用three.js可视化: 在这里插入图片描述

Polyhedron_3类作为Boost Graph概念的模型

Polyhedron_3 类是BGLCGAL提供的大多数图概念的一个模型。完整的列表可以在boost::graph_traits的文档中找到。这些例子展示了如何使用 Polyhedron_3的一些BGL算法,以及如何使用CGAL提供的概念来实现一个简单的算法。

多面体表面Polyhedral Surface的最小生成树MST

下面的示例程序计算多面体表面上的最小生成树。更多的例子可以在“三角网格简化”一章中找到。

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <iostream>
#include <list>
#include <CGAL/boost/graph/kruskal_min_spanning_tree.h>
typedef CGAL::Simple_cartesian<double>                       Kernel;
//typedef Kernel::Vector_3                                     Vector;
typedef Kernel::Point_3                                      Point;
typedef CGAL::Polyhedron_3<Kernel>                           Polyhedron;

typedef boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Polyhedron>::vertex_iterator   vertex_iterator;
typedef boost::graph_traits<Polyhedron>::edge_descriptor   edge_descriptor;

// The BGL makes heavy use of indices associated to the vertices
// We use a std::map to store the index
typedef std::map<vertex_descriptor,int> Vertex_index_map;
Vertex_index_map vertex_index_map;
// A std::map is not a property map, because it is not lightweight
typedef boost::associative_property_map<Vertex_index_map> Vertex_index_pmap;
Vertex_index_pmap vertex_index_pmap(vertex_index_map);
void
kruskal(const Polyhedron& P)
{
  // associate indices to the vertices
  vertex_iterator vb, ve;
  int index = 0;
  // boost::tie assigns the first and second element of the std::pair
  // returned by boost::vertices to the variables vb and ve
  for(boost::tie(vb, ve)=vertices(P); vb!=ve; ++vb){
    vertex_index_pmap[*vb]= index++;
  }
  // We use the default edge weight which is the length of the edge
  // This property map is defined in graph_traits_Polyhedron_3.h
  // In the function call you can see a named parameter: vertex_index_map
  std::list<edge_descriptor> mst;
  boost::kruskal_minimum_spanning_tree(P,
                                       std::back_inserter(mst),
                                       boost::vertex_index_map(vertex_index_pmap));
  std::cout << "#VRML V2.0 utf8\n"
    "Shape {\n"
    "  appearance Appearance {\n"
    "    material Material { emissiveColor 1 0 0}}\n"
    "    geometry\n"
    "    IndexedLineSet {\n"
    "      coord Coordinate {\n"
    "        point [ \n";
  for(boost::tie(vb, ve) = vertices(P); vb!=ve; ++vb){
    std::cout <<  "        " << (*vb)->point() << "\n";
  }
  std::cout << "        ]\n"
               "     }\n"
    "      coordIndex [\n";
  for(std::list<edge_descriptor>::iterator it = mst.begin(); it != mst.end(); ++it)
  {
    edge_descriptor e = *it ;
    vertex_descriptor s = source(e,P);
    vertex_descriptor t = target(e,P);
    std::cout << "      " << vertex_index_pmap[s] << ", " << vertex_index_pmap[t] <<  ", -1\n";
  }
  std::cout << "]\n"
    "  }#IndexedLineSet\n"
    "}# Shape\n";
}
int main() {
  Polyhedron P;
  Point a(1,0,0);
  Point b(0,1,0);
  Point c(0,0,1);
  Point d(0,0,0);
  P.make_tetrahedron(a,b,c,d);
  kruskal(P);
  return 0;
}

四面体最小生成树: 在这里插入图片描述 在这里插入图片描述

使用带有ID的顶点和边

下面的示例程序调用了BGL Kruskal的最小生成树算法,访问了存储在多面体顶点中的 id() 字段。

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedron_items_with_id_3.h>
#include <iostream>
#include <list>
#include <CGAL/boost/graph/kruskal_min_spanning_tree.h>
typedef CGAL::Simple_cartesian<double>                               Kernel;
typedef Kernel::Point_3                                              Point;
typedef CGAL::Polyhedron_3<Kernel,CGAL::Polyhedron_items_with_id_3>  Polyhedron;
typedef boost::graph_traits<Polyhedron>::vertex_descriptor vertex_descriptor;
typedef boost::graph_traits<Polyhedron>::vertex_iterator   vertex_iterator;
typedef boost::graph_traits<Polyhedron>::edge_descriptor   edge_descriptor;
void
kruskal( const Polyhedron& P)
{
  // We use the default edge weight which is the length of the edge
  // This property map is defined in graph_traits_Polyhedron_3.h
  // This function call requires a vertex_index_map named parameter which
  // when  omitted defaults to "get(vertex_index,graph)".
  // That default works here because the vertex type has an "id()" method
  // field which is used by the vertex_index internal property.
  std::list<edge_descriptor> mst;
  boost::kruskal_minimum_spanning_tree(P,std::back_inserter(mst));
  std::cout << "#VRML V2.0 utf8\n"
    "Shape {\n"
    "appearance Appearance {\n"
    "material Material { emissiveColor 1 0 0}}\n"
    "geometry\n"
    "IndexedLineSet {\n"
    "coord Coordinate {\n"
    "point [ \n";
  vertex_iterator vb, ve;
  for(boost::tie(vb,ve) = vertices(P); vb!=ve; ++vb){
    std::cout << (*vb)->point() << "\n";
  }
  std::cout << "]\n"
    "}\n"
    "coordIndex [\n";
  for(std::list<edge_descriptor>::iterator it = mst.begin(); it != mst.end(); ++it){
    std::cout << source(*it,P)->id()
              << ", " << target(*it,P)->id() <<  ", -1\n";
  }
  std::cout << "]\n"
    "}#IndexedLineSet\n"
    "}# Shape\n";
}
int main() {
  Polyhedron P;
  Point a(1,0,0);
  Point b(0,1,0);
  Point c(0,0,1);
  Point d(0,0,0);
  P.make_tetrahedron(a,b,c,d);
  // associate indices to the vertices using the "id()" field of the vertex.
  vertex_iterator vb, ve;
  int index = 0;
  // boost::tie assigns the first and second element of the std::pair
  // returned by boost::vertices to the variables vit and ve
  for(boost::tie(vb,ve)=vertices(P); vb!=ve; ++vb ){
    vertex_descriptor  vd = *vb;
    vd->id() = index++;
  }
  kruskal(P);
  return 0;
}

在这里插入图片描述

三角剖分Triangulations 作为Boost Graph概念的模型

三角剖分有顶点和面,允许直接转换为图形。半边被定义为一对面柄和边缘的指数。完整的列表可以在boost::graph_traits的文档中找到。 结合计算几何和图论的一个经典例子是平面上点集的欧氏最小生成树。它可以通过在点集的Delaunay三角剖分上运行最小生成树算法来计算。

欧氏最小生成树

在下面的例子中,我们创建一个Delaunay三角网,并在其上运行Kruskal的最小生成树算法。因为三角网的顶点句柄不是数组中的索引,我们必须提供一个属性映射,将顶点句柄映射到 [0, t.number_of_vertices()) 范围内的整数。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h>

#include <CGAL/boost/graph/kruskal_min_spanning_tree.h>

#include <fstream>
#include <iostream>
#include <map>

typedef CGAL::Exact_predicates_inexact_constructions_kernel         K;
typedef K::Point_2                                                  Point;

typedef CGAL::Delaunay_triangulation_2<K>                           Triangulation;

typedef boost::graph_traits<Triangulation>::vertex_descriptor       vertex_descriptor;
typedef boost::graph_traits<Triangulation>::vertex_iterator         vertex_iterator;
typedef boost::graph_traits<Triangulation>::edge_descriptor         edge_descriptor;

// The BGL makes use of indices associated to the vertices
// We use a std::map to store the index
typedef std::map<vertex_descriptor,int>                             VertexIndexMap;

// A std::map is not a property map, because it is not lightweight
typedef boost::associative_property_map<VertexIndexMap>             VertexIdPropertyMap;

int main(int argc,char* argv[])
{
  const char* filename = (argc > 1) ? argv[1] : "data/points.xy";
  std::ifstream input(filename);
  Triangulation tr;

  Point p;
  while(input >> p)
    tr.insert(p);

  // Associate indices to the vertices
  VertexIndexMap vertex_id_map;
  VertexIdPropertyMap vertex_index_pmap(vertex_id_map);
  int index = 0;

  for(vertex_descriptor vd : vertices(tr))
    vertex_id_map[vd] = index++;

  // We use the default edge weight which is the squared length of the edge
  // This property map is defined in graph_traits_Triangulation_2.h

  // In the function call you can see a named parameter: vertex_index_map
  std::list<edge_descriptor> mst;
  boost::kruskal_minimum_spanning_tree(tr, std::back_inserter(mst),
                                       vertex_index_map(vertex_index_pmap));

  std::cout << "The edges of the Euclidean minimum spanning tree:" << std::endl;
  for(edge_descriptor ed : mst)
  {
    vertex_descriptor svd = source(ed, tr);
    vertex_descriptor tvd = target(ed, tr);
    Triangulation::Vertex_handle sv = svd;
    Triangulation::Vertex_handle tv = tvd;
    std::cout << "[ " << sv->point() << "  |  " << tv->point() << " ] " << std::endl;
  }

  return EXIT_SUCCESS;
}

在这里插入图片描述

在顶点中存储顶点ID

BGL算法广泛地使用了顶点索引。在前面的例子中,我们将索引存储在 std::map 中,并将该映射转换为属性映射。然后,这个属性映射被作为参数传递给最短路径函数。 如果用户没有显式传递一个属性映射,图算法会使用 get(boost::vertex_index,ft) 返回的属性映射。这个属性映射假定顶点有一个成员函数 id() ,该函数返回一个int类型的引用。因此,CGAL提供了一类 Triangulation_vertex_base_with_id_2 。正确设置索引是用户的责任。 实例进一步说明,图的特征也适用于Delaunay三角剖分。

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_vertex_base_with_id_2.h>

#include <CGAL/boost/graph/graph_traits_Delaunay_triangulation_2.h>
#include <CGAL/boost/graph/dijkstra_shortest_paths.h>

#include <fstream>

typedef CGAL::Exact_predicates_inexact_constructions_kernel             K;
typedef K::Point_2                                                      Point;

typedef CGAL::Triangulation_vertex_base_with_id_2<K>                    Tvb;
typedef CGAL::Triangulation_face_base_2<K>                              Tfb;
typedef CGAL::Triangulation_data_structure_2<Tvb, Tfb>                  Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds>                          Triangulation;

typedef boost::graph_traits<Triangulation>::vertex_descriptor           vertex_descriptor;
typedef boost::graph_traits<Triangulation>::vertex_iterator             vertex_iterator;

typedef boost::property_map<Triangulation, boost::vertex_index_t>::type VertexIdPropertyMap;

int main(int argc,char* argv[])
{
  const char* filename = (argc > 1) ? argv[1] : "data/points.xy";
  std::ifstream input(filename);
  Triangulation tr;

  Point p;
  while(input >> p)
    tr.insert(p);

  // associate indices to the vertices
  int index = 0;
  for(vertex_descriptor vd : vertices(tr))
    vd->id()= index++;

  VertexIdPropertyMap vertex_index_pmap = get(boost::vertex_index, tr);

  // Dijkstra's shortest path needs property maps for the predecessor and distance
  std::vector<vertex_descriptor> predecessor(num_vertices(tr));
  boost::iterator_property_map<std::vector<vertex_descriptor>::iterator, VertexIdPropertyMap>
    predecessor_pmap(predecessor.begin(), vertex_index_pmap);

  std::vector<double> distance(num_vertices(tr));
  boost::iterator_property_map<std::vector<double>::iterator, VertexIdPropertyMap>
    distance_pmap(distance.begin(), vertex_index_pmap);

  vertex_descriptor source = *vertices(tr).first;
  std::cout << "\nStart dijkstra_shortest_paths at " << source->point() << std::endl;

  boost::dijkstra_shortest_paths(tr, source, distance_map(distance_pmap)
                                            .predecessor_map(predecessor_pmap));

  for(vertex_descriptor vd : vertices(tr))
  {
    std::cout << vd->point() << " [" << vd->id() << "] ";
    std::cout << " has distance = "  << get(distance_pmap,vd) << " and predecessor ";

    vd = get(predecessor_pmap,vd);
    std::cout << vd->point() << " [" << vd->id() << "]\n";
  }

  return EXIT_SUCCESS;
}

在这里插入图片描述

参考

doc.cgal.org/latest/BGL/…