CGAL 入门基础

  本文是提供给了解c++和几何算法基础知识的CGAL新手的入门教程。第一部分展示了如何定义点和段类，以及如何在它们上应用几何谓词。本节进一步提醒大家，在使用浮点数作为坐标时，会出现严重的问题。第二部分展示了一个计算2D凸包的CGAL函数。第三部分展示了Traits类的含义，第四部分解释了概念和模型的概念。

一、点和线段基础操作

1、概述

  在第一个例子中，我们演示了如何构造一些点和一个线段，并对它们执行一些基本操作。所有CGAL头文件都在include/CGAL子目录中。所有CGAL类和函数都在命名空间CGAL中。类以大写字母开头，全局函数以小写字母开头，常量均为大写。对象的维度用后缀表示。与点类型一样，geometric primitives是在 kernel中定义的。第一个示例代码中，选择 kernel使用双精度浮点数作为该点的笛卡尔坐标。除了类型之外，我们还看到了谓词，如三点的方向测试，以及距离和中点计算等结构。谓词具有一组离散的可能结果，而构造则产生一个数字或另一个几何实体。

2、整型坐标

代码示例

#include 
#include  //笛卡尔坐标相关头文件

typedef CGAL::Simple_cartesian<double> Kernel; // 内核使用双精度浮点数作为该点的笛卡尔坐标
typedef Kernel::Point_2 Point_2;               // 二维点
typedef Kernel::Segment_2 Segment_2;           // 二维线段

int main()
{
    // 定义两个位于笛卡尔坐标系下的二维点坐标
    Point_2 p(1, 1), q(10, 10); 
    std::cout << "p = " << p << std::endl;
    std::cout << "q = " << q.x() << " " << q.y() << std::endl;
    // 计算两点之间的平方距离
    std::cout << "两点之间的平方距离：" << CGAL::squared_distance(p, q) << std::endl;
    // 计算m到线段pq的平方距离
    Segment_2 s(p, q); // p和q两点构成的线段
    Point_2 m(5, 9);   // 点坐标m

    std::cout << "m = " << m << std::endl;
    std::cout << "点m到线段pq的平方距离：" << CGAL::squared_distance(s, m) << std::endl;
    // 判断三点之间的位置关系
    std::cout << "p 到 q 再到 m 三点的关系为(与先后顺序有关)： ";
    switch (CGAL::orientation(p, q, m)) 
    {
    case CGAL::COLLINEAR:
        std::cout << "三点共线\n";
        break;
    case CGAL::LEFT_TURN:
        std::cout << "三点构成左转\n";
        break;
    case CGAL::RIGHT_TURN:
        std::cout << "三点构成右转\n";
        break;
    }

    std::cout << "p和q的中点为： " << CGAL::midpoint(p, q) << std::endl;
    return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

结果展示

p点的坐标为：1 1
q点的坐标为：x = 10, y = 10
p和q两点之间的平方距离：162
m点的坐标为：5 9
点m到线段pq的平方距离：8
p 到 q 再到 m 三点的关系为(与先后顺序有关)： 三点构成左转
p和q的中点为： 5.5 5.5
1
2
3
4
5
6
7

3、浮点型坐标

  使用浮点数进行几何运算，由于数值精度的问题，可能会出现意想不到的结果。
代码示例
surprising.cpp

#include 
#include 

typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::Point_2 Point_2;

int main()
{
    {
        Point_2 p(0, 0.3), q(1, 0.6), r(2, 0.9);
        std::cout << (CGAL::collinear(p, q, r) ? "collinear\n" : "not collinear\n");
    }
    {
        Point_2 p(0, 1.0 / 3.0), q(1, 2.0 / 3.0), r(2, 1);
        std::cout << (CGAL::collinear(p, q, r) ? "collinear\n" : "not collinear\n");
    }
    {
        Point_2 p(0, 0), q(1, 1), r(2, 2);
        std::cout << (CGAL::collinear(p, q, r) ? "collinear\n" : "not collinear\n");
    }
    return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

阅读代码时，我们可以假设它将打印三次“共线”。然而，实际输出如下:

not collinear
not collinear
collinear
1
2
3

  这是因为这些分数不能用双精度数来表示，而共线性检验将在内部计算一个$3\times 3$矩阵的行列式，它接近但不等于零，因此前两个检验是非共线性的。类似的情况也可能发生在执行左转的点上，但由于行列式计算期间的舍入误差，这些点似乎是共线的，或执行右转。如果一定要使用输入坐标的绝对浮点型精度，那么可以使用：CGAL::Exact_predicates_exact_constructions_kernel内核来执行精确的谓词和提取结构。
代码示例
exact.cpp

#include 
#include // 精确浮点型
#include 

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; // 精确浮点型
typedef Kernel::Point_2 Point_2;

int main()
{
    Point_2 p(0, 0.3), q, r(2, 0.9);
    {
        q = Point_2(1, 0.6);
        std::cout << (CGAL::collinear(p, q, r) ? "collinear\n" : "not collinear\n");
    }

    {
        std::istringstream input("0 0.3   1 0.6   2 0.9");
        input >> p >> q >> r;
        std::cout << (CGAL::collinear(p, q, r) ? "collinear\n" : "not collinear\n");
    }

    {
        q = CGAL::midpoint(p, r);
        std::cout << (CGAL::collinear(p, q, r) ? "collinear\n" : "not collinear\n");
    }

    return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

下面是输出结果，您可能仍然会感到惊讶

not collinear
collinear
collinear

1
2
3
4

  在第一个块中，这些点仍然不是共线的，原因很简单，您看到的文本坐标被转换为浮点数。当它们被转换为任意精确有理数时，它们完全代表浮点数，而不是文本!这与第二个块不同，它对应于从文件中读取数字。然后，从字符串直接构造任意精度有理数，以便它们准确地表示文本。在第三个块中，可以看到作为中点结构的构造是精确的，正如内核类型的名称所表示的那样。在许多情况下，您将拥有“精确”的浮点数，从某种意义上说，它们是由某些应用程序计算或从传感器获得的。它们不是字符串“0.1”或动态计算的“1.0/10.0”，而是一个全精度浮点数。如果它们被输入到一个没有构造的算法，那么可以使用一个提供精确谓词但不精确构造的内核。其中一个例子就是凸包算法，我们将在下一节中看到。输出是输入的子集，算法只比较坐标和执行方向。
  乍一看，执行精确谓词和构造的内核似乎是完美的选择，但考虑到性能要求或有限的内存资源，那么执行精确谓词和构造的内核也就未必是完美选择。此外，对于许多算法来说，精确构造是无关紧要的。例如一个曲面网格简化算法，它通过将边缘折叠到边缘的中点来迭代收缩边缘。
  大多数CGAL包都会说明它们应该使用或支持哪种内核。

二、点序列的凸包

  本节中的所有例子都计算一组点的二维凸包。我们展示了算法以表示点范围的开始/结束迭代器对的形式获取输入，并将结果(在示例中是凸包上的点)写入输出迭代器。

1、在数组Array中提取凸包点

  在第一个例子中，我们有一个5个点的数组作为输入。由于这些点的凸包是输入的子集，因此提供一个相同大小的数组来存储结果是安全的。
代码示例
array_convex_hull_2.cpp

#include 
#include 
#include  // 二维凸包头文件

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

int main()
{
    // 使用数组表示的二维点集
    Point_2 points[5] = { Point_2(0,0), Point_2(10,0), Point_2(10,10), Point_2(6,5), Point_2(4,1) };
    Point_2 result[5]; // 用于接收存储凸包点的数组
    // CGAL凸包算法
    Point_2* ptr = CGAL::convex_hull_2(points, points + 5, result);
    // 打印提取结果
    std::cout << ptr - result << " 个凸包上的点，分别为:" << std::endl;
    for (int i = 0; i < ptr - result; i++) {
        std::cout << result[i] << std::endl;
    }

    return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

输出结果

3 个凸包上的点，分别为:
0 0
10 0
10 10
1
2
3
4

  在上一节中，我们已经看到CGAL附带了几个内核。由于凸包算法只对输入点的坐标和方向进行比较，因此我们可以选择一个提供精确谓词的核，而不提供精确的几何结构。凸包函数convex_hull_2接受三个参数。前两个为输入参数，分别为输入数组points的起始指针和结束指针；第三个参数为输出参数，是结果数组result的起始指针。函数将指针返回到结果数组中，就在最后一个写入的凸包点后面，由ptr接收，因此指针差ptr-result告诉我们凸包上有多少个点。

2、在向量Vector中提取凸包点

  在第二个例子中，我们将内置数组替换为标准模板库的std::vector。

代码示例
vector_convex_hull_2.cpp

#include 
#include 

#include 

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef std::vector<Point_2> Points;

int main()
{
	Points points, result;
	points.push_back(Point_2(0, 0));
	points.push_back(Point_2(10, 0));
	points.push_back(Point_2(10, 10));
	points.push_back(Point_2(6, 5));
	points.push_back(Point_2(4, 1));


	CGAL::convex_hull_2(points.begin(), points.end(), std::back_inserter(result));
	std::cout << result.size() << " points on the convex hull" << std::endl;
	return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

  首先，调用std::vector类的push_back()方法在vector中放入一些点。然后，调用凸包函数，前两个参数points.begin()和points.end()是迭代器，是指针的泛化:它们可以解引用和递增，凸包函数是通用的，因为它接受任何可以解引用和递增的输入；
第三个参数是结果被写入的位置。在前面的例子中，我们提供了一个指向已分配内存的指针。这种指针的通用化是输出迭代器，它允许对解引用迭代器进行递增和赋值。在这个例子中，我们从一个空向量开始，它会根据需要增长。因此，我们不能简单地传递result.begin()，而是一个由辅助函数std::back_inserter(result)生成的输出迭代器。这个输出迭代器在加1时不做任何操作，而是对赋值调用result.push_back(..)。
  如果你了解STL(标准模板库)，上述内容就很有意义了，因为这是STL将算法与容器解耦的方式。如果您不了解STL，最好先熟悉一下它的基本思想。

三、关于kernel和Traits类

  在本节中，我们将展示如何表达必须满足的要求，这样才能将convex_hull_2()这样的函数用于任意点类型。如果查看convex_hull_2()函数和其他2D凸包算法的手册页，就会发现它们有两个版本。在我们目前看到的例子中，函数接受两个迭代器来表示输入点的范围，并接受一个输出迭代器来表示将结果写入其中。第二个版本有一个额外的模板参数Traits，以及一个这中类型的额外参数。

template<class InputIterator , class OutputIterator , class Traits >
OutputIterator
convex_hull_2(InputIterator first,
              InputIterator beyond,
              OutputIterator result,
              const Traits & ch_traits)

1
2
3
4
5
6
7

  典型的凸壳算法使用的几何基元（primitives）是什么？当然，这取决于算法，所以让我们考虑什么可能是最简单的高效算法，所谓的“Graham/Andrew Scan”。该算法首先对点从左到右进行排序，然后从排序后的列表中依次添加一个点，逐步构建凸包。要做到这一点，它必须至少知道一些点的类型，它应该知道如何排序这些点，它必须能够计算三组点的方向。
  这就是模板参数Traits的用武之地。ch_graham_andrew()必须提供以下嵌套类型：

• Traits::Point_2
• Traits::Less_xy_2
• Traits::Left_turn_2
• Traits::Equal_2

  如您所猜，Left_turn_2负责方向测试，Less_xy_2用于对点进行排序。这些类型必须满足的要求在ConvexHullTraits_2中有详细解释。
由于一个简单的原因，这些类型被重新组合。另一种选择是一个相当冗长的函数模板，以及一个更长的函数调用。

template <class InputIterator, class OutputIterator, class Point_2, class Less_xy_2, class Left_turn_2, class Equal_2>
OutputIterator
ch_graham_andrew( InputIterator  first,
                  InputIterator  beyond,
                  OutputIterator result);

1
2
3
4
5
6

  有两个明显的问题:什么可以用作模板参数的参数?为什么我们要有模板参数呢?
为了回答第一个问题，CGAL概念内核的任何模型都提供了概念ConvexHullTraits_2所需要的东西。
至于第二个问题，考虑一个应用我们要计算投影到yz平面的三维点的凸包。使用类Projection_traits_yz_3，这是对前面示例的一个小修改。
代码示例
convex_hull_yz.cpp

#include 
#include 
#include 
#include 
#include 

typedef CGAL::Exact_predicates_inexact_constructions_kernel K3;
typedef CGAL::Projection_traits_yz_3<K3> K;
typedef K::Point_2 Point_2;

int main()
{
	std::istream_iterator< Point_2 >  input_begin(std::cin);
	std::istream_iterator< Point_2 >  input_end;
	std::ostream_iterator< Point_2 >  output(std::cout, "\n");
	CGAL::convex_hull_2(input_begin, input_end, output, K());
	return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

  另一个示例是关于用户定义的点类型，或来自 CGAL 以外的第三方库的点类型。将点类型与此点类型所需的谓词（predicates）一起放在类的作用域中，就可以使用这些点运行convex_hull_2()。
  最后，让我们解释一下为什么一个特征对象(traits object)会传递给凸壳函数？它将允许使用更一般的投影特征对象来存储状态，例如，如果投影平面是由一个方向给出的，该方向在类Projection_traits_yz_3中是硬连线（hardwired）的。

四、概念与模型

  在上一节中，我们写道："CGAL概念内核的任何模型都提供了概念所需的内容ConvexHullTraits_2。概念(concept)是对类型(type)的一组要求，即它具有某些嵌套类型，某些成员函数，或者带有某些以类型为主的自由函数。概念的模型是满足概念要求的类(class)。
让我们看一下以下函数。

template <typename T>
T
duplicate(T t)
{
  return t;
}

1
2
3
4
5
6
7

如果你想用一个类C来实例化这个函数，类C至少必须提供一个复制构造函数，我们说这个类必须是可复制（CopyConstructible）的模型。单例类不满足此要求。
另一个例子是这个函数:

template <typename T>
T& std::min(const T& a, const T& b)
{
  return (a<b)?a:b;
}

1
2
3
4
5
6

  函数只有当operator<(…)在T中定义时才能编译，也就是说这个type必须是一个可比较（LessThanComparable）的模型。使用必需的自由函数概念的一个例子是CGAL包CGAL和Boost Graph Library中的HalfedgeListGraph。为了成为类G的HalfedgeListGraph模型，必须有一个全局函数half fedges(const G&)，等等。一个包含必需特征类的概念的例子是InputIterator。对于InputIterator的模型，std::iterator_traits类的特化必须存在(或者泛型模板必须适用)。