CGAL user manual - 2D Polygons及常规算法介绍

CGAL user manual - 2D Polygons及常规算法介绍

CGAL 5.4 - 2D Polygons: User Manual

简介

此处的多边形是由闭合的线段连接起来的。该package由一些针对多边形的常用算法。对于其中一些算法，多边形必须是简单的，也就是说，多边形的边之间除了相邻的边之外不存在相交的情况。

目前支持的算法有：

• 找到最右边，最左边，最上边，最下边的顶点；
• 计算带符号面积；
• 判断多边形是不是简单的simple；
• 判断多边形是不是凸多边形；
• 判断多边形的朝向（顶点是顺时针的还是逆时针的）；只会考虑简单多边形
• 判断一个点是不是位于多边形的内部；

此处的所有的算法，都是输入前向迭代器进行计算的。Polygon_2用来表示多边形。

算法详解（补充项）

关于Cartesian坐标系的kernel内部的计算逻辑大多位于：CGAL-5.4\include\CGAL\Cartesian\function_objects.h

Polygon_2位于头文件：CGAL-5.4\include\CGAL\Polygon_2.h。

最右边的顶点

Polygon_2.right_vertex会调用到下面的函数。

template <typename Traits>
class Compare_vertices {
    typedef typename Traits::Less_xy_2 Less_xy_2;
    typedef typename Traits::Point_2 Point_2;
    Less_xy_2 less;
public:
    Compare_vertices(Less_xy_2 less) : less(less) {}

    // `Point_like` derives from `Point_2`
    template <typename Point_like>
    bool operator()(const Point_like& p1, const Point_like& p2) {
        return less(Point_2(p1), Point_2(p2));
    }
}; // end Compare_vertice

template <class ForwardIterator, class PolygonTraits>
ForwardIterator right_vertex_2(ForwardIterator first,
                                    ForwardIterator last,
                                    const PolygonTraits &traits)
{
    CGAL_polygon_precondition(first != last);
    internal::Polygon_2::Compare_vertices<PolygonTraits>
        less(traits.less_xy_2_object());
    return std::max_element(first, last, less);
}

即，该函数的重要实现是internal::Polygon_2::Compare_vertices。less_xy_2_object返回的是LessXY_2函数对象，这个函数先进行x比较，如果第一个参数的x小于第二个参数，返回true，第二个大于第一个，返回false，如果相等的时候比较y，第一个小于第二个的时候返回true。

其他最左边，最上边，最下边的顶点处理方式类似。

计算带符号面积

假设顶点是v0, v1, v2, ..., vn。那么对应的伪代码如下：

polygon_area(polygon):
	if polygon.n_vertices < 3:
		return 0
	first = v0
	second = v1
	third = v1
	while ++v1 != last:
		result += compute_area_2(first, second, third)
		second = third
	return result

Compute area的具体实现如下：

  template <typename K>
  class Compute_area_2
  {
    typedef typename K::FT                FT;
    typedef typename K::Iso_rectangle_2   Iso_rectangle_2;
    typedef typename K::Triangle_2        Triangle_2;
    typedef typename K::Point_2           Point_2;
  public:
    typedef FT               result_type;

    result_type
    operator()( const Point_2& p, const Point_2& q, const Point_2& r ) const
    {
      FT v1x = q.x() - p.x();
      FT v1y = q.y() - p.y();
      FT v2x = r.x() - p.x();
      FT v2y = r.y() - p.y();
      // 向量叉乘为平行四边形面积，三角形面积需要除以2
      return determinant(v1x, v1y, v2x, v2y)/2;
    }

    result_type
    operator()( const Iso_rectangle_2& r ) const
    { return (r.xmax()-r.xmin()) * (r.ymax()-r.ymin()); }

    result_type
    operator()( const Triangle_2& t ) const
    { return t.area(); }
  };

判断多边形是不是简单的simple

代码片段位于CGAL-5.4\include\CGAL\Polygon_2\Polygon_2_simplicity.h : is_simple_polygon

具体代码实现如下：

template <class Iterator, class PolygonTraits>
bool is_simple_polygon(Iterator points_begin, Iterator points_end,
                       const PolygonTraits& polygon_traits)
{
    typedef Iterator ForwardIterator;
    typedef i_polygon::Vertex_data<ForwardIterator, PolygonTraits> Vertex_data;
    typedef std::set<i_polygon::Vertex_index,
                     i_polygon::Less_segments<Vertex_data> >       Tree;

    // A temporary fix as the sweep in some cases doesn't discover vertices with degree > 2
    // Todo: fix the sweep code
    std::vector<typename PolygonTraits::Point_2> points(points_begin,points_end);
    std::sort(points.begin(), points.end(), polygon_traits.less_xy_2_object());

    typename std::vector<typename PolygonTraits::Point_2>::iterator
                                  succ(points.begin()) , it(succ++);
    for(;succ != points.end(); ++it,++succ){
      if(*it == *succ){
        return false;
      }
    }
    // end of fix
    Vertex_data   vertex_data(points_begin, points_end, polygon_traits);
    Tree tree(&vertex_data);
    vertex_data.init(&tree);
    vertex_data.sweep(&tree);
    return vertex_data.is_simple_result;
}

这个函数里面分两部分逻辑，第一部分是为了补充sweep的不足，处理顶点的度大于2的情况（即，有3条及以上的边共享顶点的情况）。对应的逻辑的伪代码如下：

sort points by Less_xy_2
for pt1 in points[0...n-2]
	if pt1 equals pt1.nextpt()
		return false

第二部分的主要逻辑位于sweep函数中。其中用到了如下kernel中的函数：

• orientation_2(p,q,r): 如果r位于直线pq的左边返回CGAL::LEFT_TURN，在右边返回CGAL::RIGHT_TURN，否则共线，返回CGAL::COLLINEAR

sweep函数如下：

template <class ForwardIterator, class PolygonTraits>
void Vertex_data<ForwardIterator, PolygonTraits>::
sweep(Tree *tree)
{
    if (this->m_size < 3)
            return;
    bool succes = true;
    // 从多边形最左边的顶点开始依次往右进行遍历
    for (Index_t i=0; i< this->m_size; ++i) {
        Vertex_index cur = index_at_rank(Vertex_order(i));
            Vertex_index prev_vt = prev(cur), next_vt = next(cur);
            if (ordered_left_to_right(cur, next_vt)) {
                if (ordered_left_to_right(cur, prev_vt))
	                // 最左边的一个顶点必然会走insertion逻辑
	                // 以逆时针的顺序插入prev_vt和cur到tree中
                    succes = insertion_event(tree, prev_vt, cur, next_vt);
                else
                    succes = replacement_event(tree, prev_vt, cur);
            } else {
                if (ordered_left_to_right(cur, prev_vt))
                    succes = replacement_event(tree, cur, prev_vt);
                else
                    succes = deletion_event(tree, prev_vt, cur);
            }
            if (!succes)
                break;
    }
    if (!succes)
            this->is_simple_result = false;
}

该函数的实现采用的是sweepline方法实现的。TODO：具体实现细节。

判断多边形是不是凸多边形

实现函数位于：CGAL-.4\include\CGAL\Polygon_2\Polygon_2_algorithms_impl.h : is_convex_2

函数的实现原理比较简单。关键代码如下：

// 此处不考虑共线点和重复点的存在

  bool HasClockwiseTriples = false; // 记录是否存在连续的三个点是顺时针的
  bool HasCounterClockwiseTriples = false; // 记录是否存在连续的三个点是逆时针的
  bool Order = less_xy_2(*previous, *current); // 前一个点和当前点的顺序是从左往右，还是从右往左
  int NumOrderChanges = 0;

  do {
  switch_orient:
    // 判断连续三个点的顺序
    switch (orientation(*previous, *current, *next)) {
      case CLOCKWISE:
        HasClockwiseTriples = true;
        break;
      case COUNTERCLOCKWISE:
        HasCounterClockwiseTriples = true;
        break;
      case ZERO:
        if(equal(*current, *next)) {
          if(next == first) {
            first = current;
          }
          ++next;
          if (next == last)
            next = first;
          goto switch_orient;
        }
        break;
    }

    bool NewOrder = less_xy_2(*current, *next);
    if (Order != NewOrder) NumOrderChanges++;

	// 从左往右变化为从右往左，或从右往左变化为从左往右，最多出现两次
    if (NumOrderChanges > 2) {
#ifdef CGAL_POLYGON_DEBUG
std::cout << "too many order changes: not convex!" << std::endl;
#endif
      return false;
    }

	// 如果同时出现逆时针和顺时针场景，则返回false
    if (HasClockwiseTriples && HasCounterClockwiseTriples) {
#ifdef CGAL_POLYGON_DEBUG
std::cout << "polygon not locally convex!" << std::endl;
#endif
      return false;
    }

    previous = current;
    current = next;
    ++next;
    if (next == last) next = first;
    Order = NewOrder;
  }
  while (previous != first);

  return true;

判断多边形的朝向

实现见：CGAL-5.4\include\CGAL\Polygon_2\Polygon_2_algorithms_impl.h : orientation_2

找到最左边的顶点，然后判断该点前一个点，和后一个点的朝向，如果是逆时针，那么多边形朝向就是逆时针，如果是顺时针，多边形朝向就是顺时针。

判断一个点是不是位于多边形的内部

实现见：CGAL-5.4\include\CGAL\Polygon_2\Polygon_2_algorithms_impl.h : bounded_side_2

具体思路是遍历多边形的线段（current,next)，并判断射线{(t, point.y()) | t >= point.x()}是不是与线段相交。如果是奇数次相交，那么位于多边形内部。如果是偶数次相交位于多边形外部。具体实现参见源代码。