如何用CGAL通过大圆切割球体并计算球面多边形面积？

我需要计算球面上多个大圆分割出的球面多边形面积（实际只需最小面积，但预计不影响实现方式）。CGAL::Nef_polyhedron_S2看似完全适用，但遇到两个问题：

• 不确定如何用CGAL最优表述需求，以下是初步尝试的代码，显然不够理想；
• 使用精确核时速度过慢，无法满足场景需求；使用非精确核时，相同代码会出现段错误。

测试代码 #

#include <CGAL/Nef_polyhedron_S2.h>
#include <CGAL/Random.h>
#include <iostream>

#define EXACT 1

#if EXACT
// 速度过慢
#include <CGAL/Exact_integer.h>
#include <CGAL/Homogeneous.h>
typedef CGAL::Exact_integer RT;
typedef CGAL::Homogeneous<RT> Kernel;
#else
// 触发段错误
#include <CGAL/Cartesian.h>
typedef CGAL::Cartesian<double> Kernel;
#endif

typedef CGAL::Nef_polyhedron_S2<Kernel> Nef_polyhedron;
typedef Nef_polyhedron::Sphere_circle Sphere_circle;
typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Plane_3 Plane_3;


// 占位用的Marsaglia实现
Vector_3 create_random_unit_vector(CGAL::Random& rng) {
    double x1, x2, s;
    do {
        x1 = 2 * rng.get_double() - 1;
        x2 = 2 * rng.get_double() - 1;
        s = x1 * x1 + x2 * x2;
    } while (s >= 1); // 循环直到s小于1

    double z = 1 - 2 * s;
    double scale = 2 * std::sqrt(1 - z * z);
    double x = scale * x1;
    double y = scale * x2;

    const int N = 1000000;
    return Vector_3((int)(x * N), (int)(y * N), (int)(z * N));
}


int main()
{
  CGAL::Random rng(1);

  const int n = 10;

  std::cout << "开始构建" << std::endl;
  Nef_polyhedron N(Nef_polyhedron::EMPTY);
  for (int i=0; i<n; ++i) {
    Vector_3 v = create_random_unit_vector(rng);
    Plane_3 plane(CGAL::ORIGIN, v);
    Sphere_circle S(plane);
    N = N + Nef_polyhedron(S) * Nef_polyhedron(S.opposite());
  }

  std::cout << N.number_of_sfaces() << " 个面" << std::endl;
  std::cout << "理论值应为 " << n * n - n + 2 << std::endl;
  return 0;
}

问题分析与解决方案建议 #

1. 优化CGAL模型构建逻辑 #

当前代码的集合运算逻辑存在冗余：Nef_polyhedron(S)表示大圆分割出的一个半球，S.opposite()对应另一个半球，两者的交集是空集，这样的运算完全无法实现球面分割的预期效果。

正确的构建逻辑应从完整球面出发，逐个用大圆切割：

// 从完整球面初始化，而非空集
Nef_polyhedron N(Nef_polyhedron::COMPLETE);
for (int i=0; i<n; ++i) {
    Vector_3 v = create_random_unit_vector(rng);
    Plane_3 plane(CGAL::ORIGIN, v);
    Sphere_circle S(plane);
    // 直接用split操作分割当前多面体，效率远高于集合运算
    N.split(S);
}

2. 解决非精确核段错误问题 #

使用Cartesian<double>核出现段错误，核心原因是浮点数精度误差导致几何计算歧义（比如大圆相交时的共线、重合判定错误），可通过以下方式解决：

• 改用过滤核平衡精度与速度：使用带过滤的核，在保证精度的同时提升运算速度：

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Filtered_kernel.h>
typedef CGAL::Simple_cartesian<double> CK;
typedef CGAL::Filtered_kernel<CK> Kernel;

• 修正随机向量生成逻辑：当前代码将浮点数向量缩放转整数的操作完全多余，对于非精确核直接返回浮点数向量即可：

Vector_3 create_random_unit_vector(CGAL::Random& rng) {
    double x1, x2, s;
    do {
        x1 = 2 * rng.get_double() - 1;
        x2 = 2 * rng.get_double() - 1;
        s = x1 * x1 + x2 * x2;
    } while (s >= 1);

    double z = 1 - 2 * s;
    double scale = 2 * std::sqrt(1 - z * z);
    double x = scale * x1;
    double y = scale * x2;

    // 直接返回浮点数向量，无需缩放转整数
    return Vector_3(x, y, z);
}

3. 最小面积计算 #

遍历所有球面多边形并计算面积，最终筛选出最小值：

#include <CGAL/area.h>
#include <limits>

// 遍历所有面计算最小面积
double min_area = std::numeric_limits<double>::max();
for (auto sf = N.sfaces_begin(); sf != N.sfaces_end(); ++sf) {
    if (!sf->is_empty()) {
        double area = CGAL::area(*sf);
        if (area < min_area) {
            min_area = area;
        }
    }
}
std::cout << "最小面积: " << min_area << std::endl;

内容的提问来源于stack exchange，提问作者Dániel Varga