Hello World 之 CGAL

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\n";

break;

case CGAL::LEFT_TURN:

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

break;

case CGAL::RIGHT_TURN:

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

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

}

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

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

}

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

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

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/.