2D和3D线性几何内核

1 介绍

CGAL（ the Computational Geometry Algorithms Library），由C++编写。CGAL包含三个主要的部分：kernel（内核）、基本几何数据结构和算法的集合、非几何支持设施。

kernel（内核包含恒定大小的对象，例如点，向量，射线，四面体等等），每个类型都有其一系列的应用于类对象的函数（例如点的坐标、点相对于对象的位置的测试、返回边界框、对象的长度或面积的函数等等）。 CGAL 内核还包含基本操作，例如仿射变换、交叉点的检测和计算以及距离计算。

2 内核表示

笛卡尔核（Cartesian Kernels）

使用 Cartesian<FieldNumberType> 可以选择坐标的笛卡尔表示，必须同时声明坐标的类型（FieldNumberType，int不是FieldNumberType）。然鹅，对于某些使用笛卡尔表示的计算，不需要除法运算的话，RingNumberType 就足够了。

Cartesian<FieldNumberType> 在内部使用引用计数来节省复制成本。 CGAL 还提供了 Simple_cartesian<FieldNumberType>，这是一个使用笛卡尔表示但没有引用计数的内核。使用 Simple_cartesian<FieldNumberType> 进行调试更容易，因为坐标存储在类中，因此可以直接访问坐标。根据算法的不同，它也可能比 Cartesian<FieldNumberType> 的效率略高或略低。

齐次核（Homogeneous Kernels）

齐次坐标允许避免数值计算中的除法运算，因为附加坐标可以用作公分母。避免除法对于精确的几何计算很有用。使用 Homogeneous<RingNumberType> 可以为内核对象的坐标选择齐次表示。由于齐次表示不使用除法，因此与齐次表示类关联的数字类型必须是仅用于较弱概念RingNumberType的模型。

同样，CGAL 还提供了 Simple_homogeneous<RingNumberType>，这是一个使用齐次表示但没有引用计数的内核。特点同Simple_cartesian<FieldNumberType> 。

命名规则

几何对象的名称首字母大写，例如Point， Segment， Triangle。
下划线后跟对象的维度，例如 _2、_3 或 _d。
内核类作为参数，它本身使用数字类型进行参数化，例如 Cartesian<double> 或 Homogeneous<leda_integer>。

内核的选择和预定义

如果从积分笛卡尔坐标开始，许多几何计算将只涉及整数数值。尤其是对于仅评估谓词的几何计算而言，这是正确的，这相当于行列式计算。示例是点集的三角剖分和凸包计算。在这种情况下，笛卡尔表示可能是首选，即使是环类型。

如果要构建新点，例如两条线的交点，笛卡尔坐标的计算通常涉及除法。因此，需要使用具有笛卡尔表示的 FieldNumberType（此时不再使用RingNumberType），或者切换到齐次表示。如果计算的可靠性至关重要，那么正确的选择数字类型可能是保证精确计算的关键。

预定义内核

CGAL 为一般有用的内核提供了 5 个 typedef，它们都是笛卡尔表示的内核，且都支持double类型的笛卡尔坐标构造点，都提供了精确的几何谓词：

Exact_predicates_inexact_constructions_kernel：提供精确的几何谓词，但由于舍入误差，几何构造可能不精确。然而，它对于许多 CGAL 算法来说已经足够了，并且比下面精确构造的内核更快。
Exact_predicates_exact_constructions_kernel：除了精确的几何谓词之外，还提供精确的几何构造。
Exact_predicates_exact_constructions_kernel_with_sqrt：同Exact_predicates_exact_constructions_kernel，但数字类型是概念FieldWithSqrt的模型。
Exact_predicates_exact_constructions_kernel_with_kth_root：同 Exact_predicates_exact_constructions_kernel，但数字类型是概念 FieldWithKthRoot 的模型。
Exact_predicates_exact_constructions_kernel_with_root_of：同Exact_predicates_exact_constructions_kernel，但数字类型是概念FieldWithRootOf的模型。

3 内核几何

点和向量（Points and Vectors）

CGAL严格区分点、向量和方向。一个点是欧几里得空间中的一个点，一个向量是两个点差，表示向量空间中从一个点到另一个点的方向和距离，一个方向实质上是一个忽略长度的向量。

CGAL 定义了一个 Origin 类型的符号常量 ORIGIN，它表示原点处的点。该常数用于点和向量之间的转换。从点 p 处减去它会得到 p 的轨迹向量。

Cartesian<double>::Point_2 p(1.0, 1.0), q;
Cartesian<double>::Vector_2 v;
v = p - ORIGIN;
q = ORIGIN + v;
assert( p == q );

如果需要确定两点的中点，可以写成

q = p_1 + (p_2 - p_1) / 2.0;
// 或者直接调用函数midpoint(p_1,p_2)

内核对象

除了点 (Kernel::Point_2, Kernel::Point_3)、向量 (Kernel::Vector_2, Kernel::Vector_3) 和方向 (Kernel::Direction_2, Kernel::Direction_3)，CGAL 还提供线、射线、线段、平面、三角形、四面体、等长方形、等长方体、圆形和球体。

CGAL 中的线 (Kernel::Line_2, Kernel::Line_3) 是有方向的。在二维空间中，它们将平面划分为正侧和负侧。一条射线（Kernel::Ray_2, Kernel::Ray_3）是一条线上的半无限区间，这条线从这个区间的有限端点朝向这个区间的任何其他点。段 (Kernel::Segment_2, Kernel::Segment_3) 是有向线上的有界区间，端点是有序的，因此它们的方向与线的方向相同。

平面是环境空间 E3 中维数为 2 的仿射子空间，通过三个点，或者一个点和一条线、射线或线段。平面也是有方向的，并将空间划分为正面和负面。

关于多边形和多面体，内核提供三角形、等向矩形、等向长方体和四面体。更复杂的多边形和多面体或多面体表面可以从基本库（Polygon_2，Polyhedron_3）中获得。

方向和相对位置

CGAL 中的几何对象具有成员函数，用于测试点相对于对象的位置。

orient_side()，判断一个测试点是在正侧、负侧还是在有向边界上。这些函数返回 Oriented_side 类型的值。
bounded_side()，判断一个测试点相对于那些分成有界和无界部分的空间对象的位置，返回类型为 Bounded_side 。
has_on()，判断一个测试点相对于那些低维对象的位置，例如三维空间中的三角形或二维空间中的线段，返回类型为布尔值。

4 谓词和构造（Predicates and Constructions）

谓词

谓词是几何内核的核心。

CGAL 为点集的方向提供谓词（orientation()、left_turn()、right_turn()、collinear()、coplanar()），用于根据给定顺序比较点，特别是比较笛卡尔坐标（例如 lexicographically_xy_smaller() )、圆内和球内测试，以及比较距离的谓词。不仅返回布尔值的组件称为谓词，而且返回枚举类型（如比较结果或方向）的组件也称为谓词。

构造

生成既不是 bool 类型也不是 enum 类型的对象的函数和函数对象称为构造。构造涉及新数值的计算，并且由于舍入误差可能不精确，除非使用具有精确数字类型的内核。

示例

typedef Cartesian<double> K;
typedef K::Point_2 Point_2;
typedef K::Segment_2 Segment_2;
Segment_2 segment_1, segment_2;
std::cin >> segment_1 >> segment_2;
/* use auto */
auto v = intersection(segment_1, segment_2);
if (v) {
  /* not empty */
  if (const Point_2 *p = boost::get<Point_2>(&*v) ) {
    /* do something with *p */
  } else {
    const Segment_2 *s = boost::get<Segment_2>(&*v);
    /* do something with *s */
  }
} else {
  /* empty intersection */
}

建设性谓词（Constructive Predicates）

为了测试点 p 相对于由三个点 q、r 和 s 定义的平面的位置，可能会尝试构建平面 Kernel::Plane_3(q,r,s) 并使用方法 oriented_side(p)。然而，除非数字类型是精确的，否则构造的平面只是近似的，舍入误差可能导致oriented_side(p) 返回一个与p、q、r 和s 的真实方向不同的方向。在 CGAL 中提供了谓词可以使此类几何决策直接参考输入点 p、q、r、s 做出，而不需要引入像平面这样的中间对象。对于上述测试，获得结果的推荐方法是使用orientation(p,q,r,s)。对于精确数字类型，情况有所不同。如果要对同一平面进行多个测试，则构建平面并使用oriented_side(p) 是可行的。

5 可拓展内核（Extensible Kernel）

CGAL支持用户将自定义的几何类插入到现有的CGAL内核中，实现内核的拓展。

一个栗子

MyPointC2.h

#ifndef MY_POINTC2_H
#define MY_POINTC2_H

#include <CGAL/Origin.h>
#include <CGAL/Bbox_2.h>

class MyPointC2 {
private:
  double vec[2];
  int col;
public:
  MyPointC2()
    : col(0)
  {
    *vec = 0;
    *(vec+1) = 0;
  }
  MyPointC2(const double x, const double y, int c = 0)
    : col(c)
  {
    *vec = x;
    *(vec+1) = y;
  }
  const double& x() const  { return *vec; }
  const double& y() const { return *(vec+1); }
  double& x() { return *vec; }
  double& y() { return *(vec+1); }
  int color() const { return col; }
  int& color() { return col; }
  bool operator==(const MyPointC2 &p) const
  {
    return ( *vec == *(p.vec) )  && ( *(vec+1) == *(p.vec + 1) && ( col == p.col) );
  }
  bool operator!=(const MyPointC2 &p) const
  {
      return !(*this == p);
  }
};
#endif // MY_POINTC2_H

对于这个类，为了使它运行正常，必须提供Kernel::Construct_bbox_2这个仿函数。

MyConstruct_bbox_2.h

#ifndef MYCONSTRUCT_BBOX_2_H
#define MYCONSTRUCT_BBOX_2_H
#include "MyPointC2.h"
template <class ConstructBbox_2>
class MyConstruct_bbox_2 : public ConstructBbox_2 {
public:
  using ConstructBbox_2::operator();
  CGAL::Bbox_2 operator()(const MyPointC2& p) const {
    return CGAL::Bbox_2(p.x(), p.y(), p.x(), p.y());
  }
};
#endif //MYCONSTRUCT_BBOX_2_H

随机访问一个点的笛卡尔坐标也是类似的。由于坐标存储在双精度数组中，我们可以使用 double* 作为随机访问迭代器。

MyConstruct_coord_iterator.h

#ifndef MYCONSTRUCT_COORD_ITERATOR_H
#define MYCONSTRUCT_COORD_ITERATOR_H
#include "MyPointC2.h"
class MyConstruct_coord_iterator {
public:
  const double* operator()(const MyPointC2& p)
  {
    return &p.x();
  }
  const double* operator()(const MyPointC2& p, int)
  {
    const double* pyptr = &p.y();
    pyptr++;
    return pyptr;
  }
};
#endif //MYCONSTRUCT_COORD_ITERATOR_H

需要提供构造点的仿函数，不必将带有 Origin 作为参数的构造函数添加到类中，也不必强制添加具有齐次坐标的构造函数。

MyConstruct_point_2.h

#ifndef MYCONSTRUCT_POINT_2_H
#define MYCONSTRUCT_POINT_2_H
template <typename K, typename OldK>
class MyConstruct_point_2
{
  typedef typename K::RT         RT;
  typedef typename K::Point_2    Point_2;
  typedef typename K::Line_2     Line_2;
  typedef typename Point_2::Rep  Rep;
public:
  typedef Point_2                result_type;
  // Note : the CGAL::Return_base_tag is really internal CGAL stuff.
  // Unfortunately it is needed for optimizing away copy-constructions,
  // due to current lack of delegating constructors in the C++ standard.
  Rep // Point_2
  operator()(CGAL::Return_base_tag, CGAL::Origin o) const
  { return Rep(o); }
  Rep // Point_2
  operator()(CGAL::Return_base_tag, const RT& x, const RT& y) const
  { return Rep(x, y); }
  Rep // Point_2
  operator()(CGAL::Return_base_tag, const RT& x, const RT& y, const RT& w) const
  { return Rep(x, y, w); }
  Point_2
  operator()(const CGAL::Origin&) const
  { return MyPointC2(0, 0, 0); }
  Point_2
  operator()(const RT& x, const RT& y) const
  {
    return MyPointC2(x, y, 0);
  }
  const Point_2&
  operator()(const Point_2 & p) const
  {
    return p;
  }
  Point_2
  operator()(const Line_2& l) const
  {
    typename OldK::Construct_point_2 base_operator;
    Point_2 p = base_operator(l);
    return p;
  }
  Point_2
  operator()(const Line_2& l, int i) const
  {
    typename OldK::Construct_point_2 base_operator;
    return base_operator(l, i);
  }
  // We need this one, as such a functor is in the Filtered_kernel
  Point_2
  operator()(const RT& x, const RT& y, const RT& w) const
  {
    if(w != 1){
      return MyPointC2(x/w, y/w, 0);
    } else {
      return MyPointC2(x,y, 0);
    }
  }
};
#endif //MYCONSTRUCT_POINT_2_H

把程序合起来

MyKernel.h

#ifndef MYKERNEL_H
#define MYKERNEL_H
#include <CGAL/Cartesian.h>
#include "MyPointC2.h"
#include "MySegmentC2.h"
#include "MyConstruct_bbox_2.h"
#include "MyConstruct_coord_iterator.h"
#include "MyConstruct_point_2.h"
// K_ is the new kernel, and K_Base is the old kernel
template < typename K_, typename K_Base >
class MyCartesian_base
  : public K_Base::template Base<K_>::Type
{
  typedef typename K_Base::template Base<K_>::Type   OldK;
public:
  typedef K_                                Kernel;
  typedef MyPointC2                         Point_2;
  typedef MySegmentC2<Kernel>               Segment_2;
  typedef MyConstruct_point_2<Kernel, OldK>       Construct_point_2;
  typedef const double*                     Cartesian_const_iterator_2;
  typedef MyConstruct_coord_iterator        Construct_cartesian_const_iterator_2;
  typedef MyConstruct_bbox_2<typename OldK::Construct_bbox_2>
                                            Construct_bbox_2;
  Construct_point_2
  construct_point_2_object() const
  { return Construct_point_2(); }
  Construct_bbox_2
  construct_bbox_2_object() const
  { return Construct_bbox_2(); }
  Construct_cartesian_const_iterator_2
  construct_cartesian_const_iterator_2_object() const
  { return Construct_cartesian_const_iterator_2(); }
  template < typename Kernel2 >
  struct Base { typedef MyCartesian_base<Kernel2, K_Base>  Type; };
};
template < typename FT_ >
struct MyKernel
  : public CGAL::Type_equality_wrapper<
                MyCartesian_base<MyKernel<FT_>, CGAL::Cartesian<FT_> >,
                MyKernel<FT_> >
{};
#endif // MYKERNEL_H

最后，谓词和构造与新点一起使用，它们可用于构造线段和三角形，以及基本库中的数据结构，因为 Delaunay 三角剖分与它们一起使用。

MyKernel.cpp

#include <CGAL/basic.h>
#include <CGAL/Filtered_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/squared_distance_2.h>
#include <cassert>
#include "MyKernel.h"
#include "MyPointC2_iostream.h"
typedef MyKernel<double>                   MK;
typedef CGAL::Filtered_kernel_adaptor<MK>  K;
typedef CGAL::Delaunay_triangulation_2<K>  Delaunay_triangulation_2;
typedef K::Point_2         Point;
typedef K::Segment_2       Segment;
typedef K::Ray_2           Ray;
typedef K::Line_2          Line;
typedef K::Triangle_2      Triangle;
typedef K::Iso_rectangle_2 Iso_rectangle;
const int RED= 1;
const int BLACK=2;
int main()
{
  Point a(0,0), b(1,0), c(1,1), d(0,1);
  a.color()=RED;
  b.color()=BLACK;
  d.color()=RED;
  Delaunay_triangulation_2 dt;
  dt.insert(a);
  K::Orientation_2 orientation;
  orientation(a,b,c);
  Point p(1,2), q;
  p.color() = RED;
  q.color() = BLACK;
  std::cout << p << std::endl;
  K::Compute_squared_distance_2 squared_distance;
  std::cout << "squared_distance(a, b) == "
            << squared_distance(a, b) << std::endl;
  Segment s1(p,q), s2(a, c);
  K::Construct_midpoint_2 construct_midpoint_2;
  Point mp = construct_midpoint_2(p,q);
  std::cout << "midpoint(" << p << " , " << q << ") == " << mp << std::endl;
  assert(s1.source().color() == RED);
  K::Intersect_2 intersection;
  const auto intersect = intersection(s1, s2);
  K::Construct_cartesian_const_iterator_2 construct_it;
  K::Cartesian_const_iterator_2  cit = construct_it(a);
  assert(*cit == a.x());
  cit = construct_it(a,0);
  cit--;
  assert(*cit == a.y());
  Line l1(a,b), l2(p, q);
  intersection(l1, l2);
  intersection(s1, l1);
  Ray r1(d,b), r2(d,c);
  intersection(r1, r2);
  intersection(r1, l1);
  squared_distance(r1, r2);
  squared_distance(r1, l2);
  squared_distance(r1, s2);
  Triangle t1(a,b,c), t2(a,c,d);
  intersection(t1, t2);
  intersection(t1, l1);
  intersection(t1, s1);
  intersection(t1, r1);
  Iso_rectangle i1(a,c), i2(d,p);
  intersection(i1, i2);
  intersection(i1, s1);
  intersection(i1, r1);
  intersection(i1, l1);
  t1.orientation();
  std::cout << s1.source() << std::endl;
  std::cout << t1.bbox() << std::endl;
  std::cout << "done" << std::endl;
  return 0;
}

最后一个例子中官方文档里似乎缺少了几个头文件。