如何使用CGAL来使用大圆切割一个球体?

4
一堆在球面上的大圆将球面切割成球面多边形。我想知道每个多边形的面积。(实际上,我只需要最小的面积,但我不认为这会改变任何事情。)
CGAL::Nef_polyhedron_S2似乎非常适合这个任务,但我有两个问题。首先,我不确定如何最好地向CGAL提出问题。下面是一个初步尝试,但我确定它不是最佳选择。对于我的使用情况,使用精确核心太慢了。这也导致了第二个问题:相同的代码在不精确的核心下导致段错误。
#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;
}
1个回答

4
CGAL的2D Arrangement package(Arrangement_on_surface_2)从版本5.4开始支持由球面上的测地弧引起的2D排列。在名称中包含“spherical”一词的示例中寻找。阅读该软件包的用户手册,特别是第6章“曲面上的排列”(https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#aos_sec-curved_surfaces)和第7.2.4节嵌入在球体中的大圆弧(https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html#arr_ssectr_spherical)。
您可以构建、维护和操作这样的排列。您还可以独立使用特征类模板Arr_geodesic_arc_on_sphere_traits_2(就像软件包中的其他特征一样)。请注意,只使用有理算术。
这是一个示例,首先将3个(完整的)大圆聚合插入到一个排列中,然后使用增量插入方法插入另一个大圆。最终的排列由6个顶点,12条边和8个面组成。构造函数中的方向是包含大圆的平面的法线。有关更多信息,请参阅手册。它详细而全面。
#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;
}

谢谢 @Efi,这非常有帮助!这个软件包非常有前途,但也让人有些畏惧。如果您有时间的话,我真的很希望能得到更多的指导。比如说,如果我有一个表示n个大圆的排列,要如何添加第n+1个大圆到这个排列中呢?在我失败的Nef尝试中,这是一个单一的并集操作。 - undefined
太棒了,@Efi,简直太好用了!而且速度也非常快! - undefined

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接