【CGAL_几何内核】2D和3D线性几何内核

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

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

参考:CGAL 5.5 – 2D and 3D Linear Geometry Kernel: User Manual

© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

取消
昵称表情代码图片

    暂无评论内容