如何获得 CGAL 中一般多边形集的面积？

Question

在以下取自 here 的 CGAL 示例中，我有一个圆形和矩形的并集。看起来是这样的：

如何使用 CGAL 获得合并结果的面积？

// Compile with: clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O2 -g -DNDEBUG -Wall -Wextra -pedantic -march=native -frounding-math bob.cpp -lgmpxx -lmpfr -lgmp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Gps_circle_segment_traits_2.h>
#include <CGAL/General_polygon_set_2.h>
#include <CGAL/Lazy_exact_nt.h>
#include <list>
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2                                   Point_2;
typedef Kernel::Circle_2                                  Circle_2;
typedef CGAL::Gps_circle_segment_traits_2<Kernel>         Traits_2;
typedef CGAL::General_polygon_set_2<Traits_2>             Polygon_set_2;
typedef Traits_2::General_polygon_2                       Polygon_2;
typedef Traits_2::General_polygon_with_holes_2            Polygon_with_holes_2;
typedef Traits_2::Curve_2                                 Curve_2;
typedef Traits_2::X_monotone_curve_2                      X_monotone_curve_2;
// Construct a polygon from a circle.
Polygon_2 construct_polygon (const Circle_2& circle)
{
  // Subdivide the circle into two x-monotone arcs.
  Traits_2 traits;
  Curve_2 curve (circle);
  std::list<CGAL::Object>  objects;
  traits.make_x_monotone_2_object() (curve, std::back_inserter(objects));
  CGAL_assertion (objects.size() == 2);
  // Construct the polygon.
  Polygon_2 pgn;
  X_monotone_curve_2 arc;
  std::list<CGAL::Object>::iterator iter;
  for (iter = objects.begin(); iter != objects.end(); ++iter) {
    CGAL::assign (arc, *iter);
    pgn.push_back (arc);
  }
  return pgn;
}
// Construct a polygon from a rectangle.
Polygon_2 construct_polygon (const Point_2& p1, const Point_2& p2,
                             const Point_2& p3, const Point_2& p4)
{
  Polygon_2 pgn;
  X_monotone_curve_2 s1(p1, p2);    pgn.push_back(s1);
  X_monotone_curve_2 s2(p2, p3);    pgn.push_back(s2);
  X_monotone_curve_2 s3(p3, p4);    pgn.push_back(s3);
  X_monotone_curve_2 s4(p4, p1);    pgn.push_back(s4);
  return pgn;
}
// The main program:
int main ()
{
  // Insert four non-intersecting circles.
  Polygon_set_2 S;
  Polygon_2 circ1, circ2, circ3, circ4;
  circ1 = construct_polygon(Circle_2(Point_2(1, 1), 1));  S.insert(circ1);
  circ2 = construct_polygon(Circle_2(Point_2(5, 1), 1));  S.insert(circ2);
  circ3 = construct_polygon(Circle_2(Point_2(5, 5), 1));  S.insert(circ3);
  circ4 = construct_polygon(Circle_2(Point_2(1, 5), 1));  S.insert(circ4);
  // Compute the union with four rectangles incrementally.
  Polygon_2 rect1, rect2, rect3, rect4;
  rect1 = construct_polygon(Point_2(1, 0), Point_2(5, 0),
                            Point_2(5, 2), Point_2(1, 2));
  S.join (rect1);
  rect2 = construct_polygon(Point_2(1, 4), Point_2(5, 4),
                            Point_2(5, 6), Point_2(1, 6));
  S.join (rect2);
  rect3 = construct_polygon(Point_2(0, 1), Point_2(2, 1),
                            Point_2(2, 5), Point_2(0, 5));
  S.join (rect3);
  rect4 = construct_polygon(Point_2(4, 1), Point_2(6, 1),
                            Point_2(6, 5), Point_2(4, 5));
  S.join (rect4);
  // Print the output.
  std::list<Polygon_with_holes_2> res;
  S.polygons_with_holes (std::back_inserter (res));
  std::copy (res.begin(), res.end(),
             std::ostream_iterator<Polygon_with_holes_2>(std::cout, "\n"));
  std::cout << std::endl;
  return 0;
}

Answer 1

可能具有圆形边缘的一般多边形没有 area 函数 - 因此我们需要深入了解它们的表示细节并自己编写此类函数（叹息！！）。这种 area 函数的示例如下：

auto squared_distance(const Traits_2::Point_2& P1, const Traits_2::Point_2& P2)
{
  const auto dx = P1.x() - P2.x();
  const auto dy = P1.y() - P2.y();
  return dx * dx + dy * dy;
}

auto area(const Polygon_2& P)
{
  auto res = 0.0;
  for (auto it = P.curves_begin(); it != P.curves_end(); ++it)
  {
    if (it->is_linear())
    {
      const auto s = it->source();
      const auto t = it->target();
      res += CGAL::to_double((s.x() - t.x()) * (s.y() + t.y()) / 2);
    }
    else if (it->is_circular())
    {
      const auto s = it->source();
      const auto t = it->target();
      res += CGAL::to_double((s.x() - t.x()) * (s.y() + t.y()) / 2);
      const auto ds = CGAL::to_double(squared_distance(s, t));
      const auto rs = CGAL::to_double(it->supporting_circle().squared_radius());
      const auto areaSector = rs * std::asin(std::sqrt(ds) / (std::sqrt(rs) * 2));
      const auto areaTriangle = std::sqrt(ds) * std::sqrt(rs * 4 - ds) / 4;
      res += (areaSector - areaTriangle);
    }
  }
  return res;
}

auto area(const Polygon_with_holes_2& P)
{
  auto res = area(P.outer_boundary());
  for (auto it = P.holes_begin(); it != P.holes_end(); ++it) res += area(*it);
  return res;
}

如您所见，我在这里使用 Polygon_with_holes_2 访问函数 outer_boundary、holes_begin 和 holes_end。对于 Polygon_2 我只使用了 curves_begin 和 curves_end 访问函数。处理此多边形的边缘需要更多的努力，这些边缘属于 X_monotone_curve_2 类型。对于您需要为一般多边形开发的任何其他功能，所有这些访问功能都是必需的。

包含您的示例的完整 MWE 如下：

// Compile with: clang++ -DBOOST_ALL_NO_LIB -DCGAL_USE_GMPXX=1 -O2 -g -DNDEBUG -Wall -Wextra -pedantic -march=native -frounding-math bob.cpp -lgmpxx -lmpfr -lgmp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Gps_circle_segment_traits_2.h>
#include <CGAL/General_polygon_set_2.h>
#include <CGAL/Lazy_exact_nt.h>

#include <list>

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2                                   Point_2;
typedef Kernel::Circle_2                                  Circle_2;
typedef CGAL::Gps_circle_segment_traits_2<Kernel>         Traits_2;
typedef CGAL::General_polygon_set_2<Traits_2>             Polygon_set_2;
typedef Traits_2::General_polygon_2                       Polygon_2;
typedef Traits_2::General_polygon_with_holes_2            Polygon_with_holes_2;
typedef Traits_2::Curve_2                                 Curve_2;
typedef Traits_2::X_monotone_curve_2                      X_monotone_curve_2;

auto squared_distance(const Traits_2::Point_2& P1, const Traits_2::Point_2& P2){
  const auto dx = P1.x() - P2.x();
  const auto dy = P1.y() - P2.y();
  return dx * dx + dy * dy;
}

auto area(const Polygon_2& P){
  double res = 0.0;
  for (auto it = P.curves_begin(); it != P.curves_end(); ++it){
    if (it->is_linear()){
      const auto s = it->source();
      const auto t = it->target();
      res += CGAL::to_double((s.x() - t.x()) * (s.y() + t.y()) / 2);
    } else if (it->is_circular()) {
      const auto s = it->source();
      const auto t = it->target();
      res += CGAL::to_double((s.x() - t.x()) * (s.y() + t.y()) / 2);
      const auto ds = CGAL::to_double(squared_distance(s, t));
      const auto rs = CGAL::to_double(it->supporting_circle().squared_radius());
      const auto areaSector = rs * std::asin(std::sqrt(ds) / (std::sqrt(rs) * 2));
      const auto areaTriangle = std::sqrt(ds) * std::sqrt(rs * 4 - ds) / 4;
      res += (areaSector - areaTriangle);
    }
  }
  return res;
}

auto area(const Polygon_with_holes_2& P){
  auto res = area(P.outer_boundary());
  for (auto it = P.holes_begin(); it != P.holes_end(); ++it) {
    res += area(*it);
  }
  return res;
}

// Construct a polygon from a circle.
Polygon_2 construct_polygon (const Circle_2& circle){
  // Subdivide the circle into two x-monotone arcs.
  Traits_2 traits;
  Curve_2 curve (circle);
  std::list<CGAL::Object>  objects;
  traits.make_x_monotone_2_object() (curve, std::back_inserter(objects));
  CGAL_assertion (objects.size() == 2);
  // Construct the polygon.
  Polygon_2 pgn;
  X_monotone_curve_2 arc;
  std::list<CGAL::Object>::iterator iter;
  for (iter = objects.begin(); iter != objects.end(); ++iter) {
    CGAL::assign (arc, *iter);
    pgn.push_back (arc);
  }
  return pgn;
}

// Construct a polygon from a rectangle.
Polygon_2 construct_polygon (
  const Point_2& p1,
  const Point_2& p2,
  const Point_2& p3,
  const Point_2& p4
){
  Polygon_2 pgn;
  X_monotone_curve_2 s1(p1, p2);    pgn.push_back(s1);
  X_monotone_curve_2 s2(p2, p3);    pgn.push_back(s2);
  X_monotone_curve_2 s3(p3, p4);    pgn.push_back(s3);
  X_monotone_curve_2 s4(p4, p1);    pgn.push_back(s4);
  return pgn;
}

// The main program:
int main (){
  // Insert four non-intersecting circles.
  Polygon_set_2 S;
  Polygon_2 circ1, circ2, circ3, circ4;
  circ1 = construct_polygon(Circle_2(Point_2(1, 1), 1));  S.insert(circ1);
  circ2 = construct_polygon(Circle_2(Point_2(5, 1), 1));  S.insert(circ2);
  circ3 = construct_polygon(Circle_2(Point_2(5, 5), 1));  S.insert(circ3);
  circ4 = construct_polygon(Circle_2(Point_2(1, 5), 1));  S.insert(circ4);
  // Compute the union with four rectangles incrementally.
  Polygon_2 rect1, rect2, rect3, rect4;
  rect1 = construct_polygon(Point_2(1, 0), Point_2(5, 0),
                            Point_2(5, 2), Point_2(1, 2));
  S.join (rect1);
  rect2 = construct_polygon(Point_2(1, 4), Point_2(5, 4),
                            Point_2(5, 6), Point_2(1, 6));
  S.join (rect2);
  rect3 = construct_polygon(Point_2(0, 1), Point_2(2, 1),
                            Point_2(2, 5), Point_2(0, 5));
  S.join (rect3);
  rect4 = construct_polygon(Point_2(4, 1), Point_2(6, 1),
                            Point_2(6, 5), Point_2(4, 5));
  S.join (rect4);

  // Print the output.
  std::list<Polygon_with_holes_2> res;
  S.polygons_with_holes (std::back_inserter (res));
  std::copy (res.begin(), res.end(),
             std::ostream_iterator<Polygon_with_holes_2>(std::cout, "\n"));
  std::cout << std::endl;

  for(const auto &x: res){
    std::cout << "Area = "<< area(x) <<std::endl;
  }

  return 0;
}

运行这给出了 31.1416 的输出。您有四个面积为 8 的矩形，每个矩形有四个面积为 1 的重叠区域和 4 个四分之一圆给出：4*8-4+π=31.141592653589793，因此结果与理论相符。

如何获得 CGAL 中一般多边形集的面积？

How does one obtain the area of a general polygon set in CGAL?

c++

cgal