How to use CGAL to slice up a sphere with great circles?

A bunch of great circles on the sphere slice up the sphere into spherical polygons. I'd like to know the areas of each of these polygons. (In fact, I only need the smallest of the areas, but I don't expect that to change anything.)

CGAL::Nef_polyhedron_S2 seems to be perfectly suitable for the task, but I have two issues. First, I'm not sure how best to formulate the question to CGAL. Below is an initial attempt, but I am sure it is suboptimal. With an exact kernel it is too slow for my use case. Which leads me to the second issue: the same code segfaults with an inexact kernel.

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

#define EXACT 1

#if EXACT
// TOO SLOW
#include <CGAL/Exact_integer.h>
#include <CGAL/Homogeneous.h>
typedef CGAL::Exact_integer RT;
typedef CGAL::Homogeneous<RT> Kernel;
#else
// SEGFAULTS
#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;


// ChatGPT-4's Marsaglia implementation as a placeholder
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); // Continue looping until s is less than 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 << "starting construction" << 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() << " faces" << std::endl;
  std::cout << "supposed to be " << n * n - n + 2 << std::endl;
  return 0;
}

Solution

The 2D Arrangement package (Arrangement_on_surface_2) of CGAL supports 2D arrangements induced by geodesic arcs on the sphere as of version 5.4. Look for examples that has the word "spherical" in their names. Read the user manual of the package, and in particular Chapter 6 "Arrangements on Curved Surfaces" (https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#aos_sec-curved_surfaces) and Section 7.2.4 Arcs of Great Circles Embedded in the Sphere (https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#arr_ssectr_spherical).

You can construct, maintain, and operate on such arrangements. You can also use the traits class template Arr_geodesic_arc_on_sphere_traits_2 independently (like any other traits in the package). Observe that only rational arithmetic is used.

Here is an example that first aggregately inserts 3 (full) great circles into an arrangement, and then uses the incremental insertion to insert another one. The final arrangement consists of 6 vertices, 12 edges, and 8 faces. The direction in the constructor of great circles is the normal to the plane that contains the great circle. For more info consult the manual. It is detailed and exhaustive.

#include <list>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Arrangement_on_surface_2.h>
#include <CGAL/Arr_geodesic_arc_on_sphere_traits_2.h>
#include <CGAL/Arr_spherical_topology_traits_2.h>
#include <CGAL/draw_arrangement_2.h>

using Kernel = CGAL::Exact_predicates_exact_constructions_kernel;
using Direction = Kernel::Direction_3;
using Geom_traits = CGAL::Arr_geodesic_arc_on_sphere_traits_2<Kernel>;
using Point = Geom_traits::Point_2;
using Curve = Geom_traits::Curve_2;
using Topol_traits = CGAL::Arr_spherical_topology_traits_2<Geom_traits>;
using Arrangement = CGAL::Arrangement_on_surface_2<Geom_traits, Topol_traits>;

int main() {
  Geom_traits traits;
  Arrangement arr(&traits);
  auto ctr_cv = traits.construct_curve_2_object();
  std::list<Curve> arcs;
  arcs.push_back(ctr_cv(Direction(0, 0, 1)));
  arcs.push_back(ctr_cv(Direction(0, 1, 0)));
  arcs.push_back(ctr_cv(Direction(1, 0, 0)));
  CGAL::insert(arr, arcs.begin(), arcs.end());
  std::cout << arr.number_of_vertices() << ", "
            << arr.number_of_edges() << ", "
            << arr.number_of_faces() << std::endl;
  CGAL::insert(arr, ctr_cv(Direction(0, 1, 1)));
  draw(arr, "aos", true);
  return 0;
}