Hello World 之 CGAL

Posted 2021-01-19 myboat

CGAL有神秘的面纱，让我不断想看清其真面目。开始吧！

 

1 Three Points and One Segment

第一个例子是创建3个点和一条线段，并且在其上进行一些操作。

所有的CGAL头文件都在CGAL目录下。所有的CGAL类和函数都在CGAL的命名空间。类以大写字母开头，常量全大写，全局函数名小写。对象的空间维度由后缀给出。

几何元，如点，在一个kernel中定义。第一个例子中我们选择的kernel采用double精度的浮点数作为笛卡尔空间坐标。

另外，我们有predicate（断言），如位置测试断言，我们有construction（构建），如距离和中点的计算，都是construction。一个predicate的结果是一个离散集，一个construction产生一个值，也可能产生一个新的几何实体。

File Kernel_23/points_and_segment.cpp

#include <iostream>

#include <CGAL/Simple_cartesian.h>

 

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 << "sqdist(p,q) = "

<< CGAL::squared_distance(p,q) << std::endl;

 

Segment_2 s(p,q);

Point_2 m(5, 9);

 

std::cout << "m = " << m << std::endl;

std::cout << "sqdist(Segment_2(p,q), m) = "

<< CGAL::squared_distance(s,m) << std::endl;

 

std::cout << "p, q, and m ";

switch (CGAL::orientation(p,q,m)){

case CGAL::COLLINEAR:

std::cout << "are collinear ";

break;

case CGAL::LEFT_TURN:

std::cout << "make a left turn ";

break;

case CGAL::RIGHT_TURN:

std::cout << "make a right turn ";

break;

}

 

std::cout << " midpoint(p,q) = " << CGAL::midpoint(p,q) << std::endl;

return 0;

}

 

下面采用浮点数进行的几何运行让我们吃惊。

File Kernel_23/surprising.cpp

#include <iostream>

#include <CGAL/Simple_cartesian.h>

 

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 " : "not collinear ");

}

{

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 " : "not collinear ");

}

{

Point_2 p(0,0), q(1, 1), r(2, 2);

std::cout << (CGAL::collinear(p,q,r) ? "collinear " : "not collinear ");

}

return 0;

}

如果只从代码上看，我们会发现前两种情况也应当是共线的，但实际的结果是：

not collinear
not collinear
collinear

因为分数作为双精度数是不可被描述的，共线测试内部的计算是一个3X3行列式（determinant），它可以得到近似值，但不能得到误差为0的精确值。所以得出前两种情况为不花线的结论。

其他的predicate也会有同样的问题，如CGAL::orientation(p,q,m)运算可能会由于舍入误差，可能得出不同实际的结论。

如果你需要使数被全精度解析，你可以使用精确断言和精确构建的CGAL kernel。 

File Kernel_23/exact.cpp

#include <iostream>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>

#include <sstream>

 

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 " : "not collinear ");

}

 

{

std::istringstream input("0 0.3 1 0.6 2 0.9");

input >> p >> q >> r;

std::cout << (CGAL::collinear(p,q,r) ? "collinear " : "not collinear ");

}

 

{

q = CGAL::midpoint(p,r);

std::cout << (CGAL::collinear(p,q,r) ? "collinear " : "not collinear ");

}

 

return 0;

}

这里的结果仍然让我们吃惊：

not collinear
collinear
collinear

第一个结果仍然是错的，原因与前面相同，它们仍然是浮点运算。第二个结果不同，它由字符串生成（construct），则精确地代表了字符串所表示的数。第三个结果通过构建（construct)中点得到第三个点，构建操作是精确的，所以结果也是正确的。

在很多情况下，你操作“精确”浮点数据，认为它们是由应用计算得到或由传感器得到的。它们不是象“0.1”这样的字符串，也不是象"1.0/10.0"这样动态（ on the fly）生成的，它是一个全精度的浮点数。如果它们只是被传递入某个算法并且没有构建（construct）操作时，你可以使用支持精确断言（predicate）和非精确构建（construct）的kernel。这样的例子包括下一节我们看到的“凸包”算法。它的输出是输入的一个子集，这个算法只进行坐标比较和位置测试。

由于高精度的计算需要消耗比普通计算多的资源，内存、时间等，所以使用时需要考虑。

大部分CGAL包说明它们需要或支持哪种kernel。

2 The Convex Hull of a Sequence of Points

 本节所有例子都是凸包算法。输入一个点序列，输出所有凸包边界上的点序列。

下面的例子输入和输出的都是一个坐标数组。

File Convex_hull_2/array_convex_hull_2.cpp

#include <iostream>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/convex_hull_2.h>

 

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];

 

Point_2 *ptr = CGAL::convex_hull_2( points, points+5, result );

std::cout << ptr - result << " points on the convex hull:" << std::endl;

for(int i = 0; i < ptr - result; i++){

std::cout << result[i] << std::endl;

}

return 0;

}

如同上节所说，我们采用了精确断言和非精确构建的kernel来生成程序。

下面的例子则采用了标准库中的vector类来进行输入和输出。

File Convex_hull_2/vector_convex_hull_2.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/convex_hull_2.h>

 

#include <vector>

 

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;

}

3 About Kernels and Traits Classes

 本节简介traits的内涵和意义.

在上个例子中，如果我们阅读convex_hull_2()的手册时，会发现它及其他2D convex_hull_2()算法都有两个版本。其中一个版本包含了traits参数。

template<class InputIterator , class OutputIterator , class Traits >

OutputIterator

convex_hull_2(InputIterator first,

InputIterator beyond,

OutputIterator result,

const Traits & ch_traits)

什么原因要我们使用traits呢？泛型编程需要使用抽象元语对算法进行抽象，而在实现中将元语变为实际的类型和算法。那么convex hull算法需要哪些元语（primitive）呢？最简单的"Graham/Andrew Scan"算法过程是：（1)将所有输入的点进行从左到右排序；（2）从左向右顺序加入，逐步形成convex hull。这个过程(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）我们为什么要用这些模板参数？

对第一个问题：ConvexHullTraits_2所要求的任何模型，这些模型由CGAL概念Kernel提供。

对第二个问题：如果我们将来需要计算投影到yz平面上的的3D点集的convex hull时，我们设计一个新的traits——Projection_traits_yz_3，这样前面的例子就不需要进行大的修改。

File Convex_hull_2/convex_hull_yz.cpp

#include <iostream>

#include <iterator>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Projection_traits_yz_3.h>

#include <CGAL/convex_hull_2.h>

 

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, " " );

CGAL::convex_hull_2( input_begin, input_end, output, K() );

return 0;

}

另一个例子是关于使用已经定义的空间点类型，或者来自非CGAL库中的点类型，将这些点类型及其相应的断言（predicates）加入类范围，然后你就可以基于新的点类型运行convex_hull_2。

最后，为什么需要将一个traits对象作为参数传入该方法呢？主要原因在于我们可以用一个更加一般的投影特征对象（projection trait）来保存状态。例子：如果这个投影平面是由一个向量给出的方向，而且是通过硬编码的方式加入Projection_traits_yz_3。

 

4 Concepts and Models

一个概念（concept）是一个类型的一个需求集（requirment set），它包括一些内嵌的类型，成员函数或一些处理该类型自由函数。

一个概念的模型（model of a concept）是一个用于实现概念所有需求的一个类。

下面有一个函数：

template <typename T>

T

duplicate(T t)

{

return t;

}

如果你用一个类C来实例化该函数，则C必须提供一个复制构造函数（copy constructor），我们称类C必须是“可复制构造的”（CopyConstructible）。

另一个例子：

template <typename T>

T& std::min(const T& a, const T& b)

{

return (a<b)?a:b;

}

这个函数只有在类型T的operator<(..)有定义时才能编译。我们称类C必须是“小于关系可比较的”（LessThanComparable）

关于自由函数的一个例子：CGAL包和Boost Graph库中的HalfedgeListGraph概念。如果一个类G要成为HalfedgeListGraph的一个模型，必须有一个全局函数halfedges(const G&)处理该类。

关于带有traits需求的概念的一个例子是InputIterator。对于一个InputIterator的模型，一个特化的std::iterator_traits类必须存在（或其通用的模板必须可用）。

5 Further Reading

阅读：

 "The C++ Standard Library, A Tutorial and Reference" by Nicolai M. Josuttis from Addison-Wesley, or "Generic Programming and the STL" by Matthew H. Austern for the STL and its notion of concepts and models.

Other resources for CGAL are the rest of the tutorials and the user support page at https://www.cgal.org/.