CGAL通用多边形面积计算代码错误问题求助

问题背景 #

我尝试实现CGAL中通用多边形（含圆弧边界）的面积计算，参考相关解答写出的代码存在计算错误，以下是原代码、复现示例及测试结果。

原实现代码 #

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

最小复现示例 #

// 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 (){
  Polygon_set_2 S;
  const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 0.5*0.5));
  const auto circ2 = construct_polygon(Circle_2(Point_2(5, 0), 0.5*0.5));

  // Comment and uncomment these blocks to switch between examples

  // const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 2*2));
  // const auto circ2 = construct_polygon(Circle_2(Point_2(5, 0), 2*2));

  // const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 2*2));
  // const auto circ2 = construct_polygon(Circle_2(Point_2(2, 0), 2*2));

  // const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 1*1));
  // const auto circ2 = construct_polygon(Circle_2(Point_2(0.5, 0), 1*1));

  S.join(circ1);
  S.join(circ2);

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

注：修正了原main函数中Circle_2的构造参数（CGAL的Circle_2默认接受平方半径）

错误测试结果 #

• 测试1：两个圆心在(0,0)和(5,0)、半径0.5的不重叠圆，代码返回π，正确结果应为0.5π
• 测试2：两个圆心在(0,0)和(5,0)、半径2的不重叠圆，代码返回0
• 测试3：两个圆心在(0,0)和(2,0)、半径2的重叠圆，代码返回0
• 测试4：两个圆心在(0,0)和(0.5,0)、半径1的重叠圆，正确面积为4.13108，代码返回4.09047

问题分析 #

1. 多边形方向与符号处理错误：
   • CGAL的General_polygon_2边界方向（顺时针/逆时针）会影响面积计算的符号，原代码未区分方向导致面积相互抵消（如测试2中两个圆的面积一正一负相加为0）。
   • 孔洞的面积应该从外边界面积中减去，原代码错误地累加了孔洞面积。
2. 圆弧面积计算错误：
   • 原代码未考虑圆弧的方向（凸向/凹向），导致圆弧部分的面积符号错误。
   • 扇形面积公式有误，正确的扇形面积应基于圆心角计算，原公式仅适用于小角度近似，且未考虑方向。
3. 线性边面积公式错误：
   • 原线性边的面积累加公式(s.x()-t.x())*(s.y()+t.y())/2符号逻辑错误，正确的 shoelace 公式应为(s.x()*t.y() - t.x()*s.y())/2。

修正后的代码 #

#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 <cmath>

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;
typedef Kernel::FT                                        FT;

// 计算两点间的平方距离
FT squared_distance(const Point_2& P1, const Point_2& P2)
{
    return CGAL::squared_distance(P1, P2);
}

// 计算单个通用多边形的面积（考虑方向）
double area(const Polygon_2& P)
{
    double res = 0.0;
    const bool is_counterclockwise = P.is_counterclockwise_oriented();

    for (auto it = P.curves_begin(); it != P.curves_end(); ++it)
    {
        const auto& curve = *it;
        const Point_2 s = curve.source();
        const Point_2 t = curve.target();

        if (curve.is_linear())
        {
            // 线性边使用Shoelace公式
            double shoelace = CGAL::to_double(s.x() * t.y() - t.x() * s.y()) / 2.0;
            res += shoelace;
        }
        else if (curve.is_circular())
        {
            // 线性部分的Shoelace贡献
            double shoelace = CGAL::to_double(s.x() * t.y() - t.x() * s.y()) / 2.0;
            res += shoelace;

            // 计算圆弧部分的面积修正
            const Circle_2& circle = curve.supporting_circle();
            const Point_2 center = circle.center();
            const FT rs = circle.squared_radius();
            const FT ds = squared_distance(s, t);

            // 计算圆心角（用向量点积）
            Point_2 vec_s(s.x() - center.x(), s.y() - center.y());
            Point_2 vec_t(t.x() - center.x(), t.y() - center.y());
            FT dot = vec_s.x() * vec_t.x() + vec_s.y() * vec_t.y();
            FT cos_theta = dot / (std::sqrt(CGAL::to_double(rs)) * std::sqrt(CGAL::to_double(rs)));
            // 限制cos_theta在[-1,1]范围内，避免数值误差
            cos_theta = std::max(-1.0, std::min(1.0, CGAL::to_double(cos_theta)));
            double theta = std::acos(cos_theta);

            // 判断圆弧方向：顺时针则圆心角取负
            if ((is_counterclockwise && curve.is_clockwise()) || (!is_counterclockwise && !curve.is_clockwise()))
            {
                theta = -theta;
            }

            // 扇形面积：0.5 * r² * theta
            double sector_area = 0.5 * CGAL::to_double(rs) * theta;
            // 三角形面积：0.5 * |vec_s × vec_t|
            double triangle_area = 0.5 * std::abs(CGAL::to_double(vec_s.x() * vec_t.y() - vec_s.y() * vec_t.x()));

            // 圆弧部分的面积贡献 = 扇形面积 - 三角形面积
            res += (sector_area - triangle_area);
        }
    }

    // 返回绝对值，确保面积为正
    return std::abs(res);
}

// 计算带孔洞的多边形面积
double area(const Polygon_with_holes_2& P)
{
    double res = area(P.outer_boundary());
    // 孔洞面积需要减去（孔洞方向与外边界相反）
    for (auto it = P.holes_begin(); it != P.holes_end(); ++it)
    {
        res -= area(*it);
    }
    return std::abs(res);
}

// 从圆构造多边形
Polygon_2 construct_polygon(const Circle_2& circle)
{
    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);

    Polygon_2 pgn;
    X_monotone_curve_2 arc;
    for (const auto& obj : objects)
    {
        CGAL::assign(arc, obj);
        pgn.push_back(arc);
    }
    return pgn;
}

// 从矩形构造多边形
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;
}

int main()
{
    Polygon_set_2 S;
    // 使用平方半径构造Circle_2
    const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 0.5*0.5));
    const auto circ2 = construct_polygon(Circle_2(Point_2(5, 0), 0.5*0.5));

    // 切换测试用例
    // const auto circ1 = construct_polygon(Circle_2(Point_2(0, 0), 2*2));
    // const auto circ2 = construct_polygon(Circle_2(Point_2(5, 0), 2*2));

    // const auto circ1 = construct_polygon(Circle_2(