Блог

Я новичок в CGAL, и в настоящее время мне нужно вычислить площадь сферических многогранников с помощью CGAL. Однако я застрял на первом шаге, то есть на построении многогранников сферы из точек сферы. Вот мой код

 typedef CGAL::Simple_cartesian<double> K;
typedef CGAL::Nef_polyhedron_S2<K> Nef_polyhedron;
typedef Nef_polyhedron::Sphere_point Sphere_point;
typedef Nef_polyhedron::Sphere_segment Sphere_segment;
typedef Nef_polyhedron::SVertex_const_iterator SVertex_const_iterator;

Sphere_point p1(1, 0, 0),  p2(0, 0, 1), p3(1, -1, 0);
Sphere_segment s1(p1, p2), s2(p2, p3), s3(p3, p1);
std::vector<Sphere_segment> tri({ s1, s2, s3 });
Nef_polyhedron S(tri.begin(), tri.end());

int i = 0;
for (SVertex_const_iterator it = S.svertices_begin(); it != S.svertices_end();   it)
{
    Sphere_point p = it->point();
    std::cout << "Point " << i << ":" << p.x() << " " << p.y() << " " << p.z() << std::endl;
    i  = 1;
}
 

Я ожидаю, что результат будет иметь размер три, который соответствует трем вершинам треугольника сферы. Тем не менее, я получаю 6 вершин:

 Point 0:0 -1 0
Point 1:-1 -0 -0
Point 2:0 0 1
Point 3:1 -1 0
Point 4:1 0 0
Point 5:0 1 0
 

Мой вопрос в том, как я могу построить многогранник сферы, окруженный входными сегментами сферы, с многогранником в левой части первого входного сегмента?

обновление Я попытался построить многогранник с помощью другого подхода, то есть путем вычисления пересечения нескольких полушарий. Это, по-видимому, гарантирует, что количество вершин в построенном сферическом многограннике соответствует входным данным, учитывая, что входные координаты являются целыми числами. Но код выдает ошибку, когда вводится double номер общего типа. Вот код

 typedef CGAL::Simple_cartesian<double> Ker;
typedef CGAL::Nef_polyhedron_S2<Ker> Nef_polyhedron;
typedef Nef_polyhedron::Sphere_circle Sphere_circle;
typedef Nef_polyhedron::Sphere_point Sphere_point;

//CPoint is my custom Cartesian point class that 
//support vector cross product, which I used to 
//compute the normal. Its Cartesian coordinates 
//can be accessed by [] operator.
Nef_polyhedron bd_vert_to_sph_poly(const std::vector<CPoint>amp; pts)
{
    size_t n = pts.size();
    //! 1. compute half-spheres using normal
    std::cout << "Input Cartesian Coordinates:" << std::endl;
    std::vector<Sphere_circle> circles(n);
    for (size_t i = 0; i < n; i  )
    {
        CPoint p1 = pts[i];
        CPoint p2 = pts[(i   1) % n];
        CPoint normal = p1 ^ p2;
        circles[i] = Sphere_circle(normal[0], normal[1], normal[2]);
        std::cout << p1[0] << " " << p1[1] << " " << p1[2] << std::endl;
    }
    //! 2. intersect S
    Nef_polyhedron S = Nef_polyhedron(Nef_polyhedron::COMPLETE);
    for (size_t i = 0; i < n; i  )
    {
        S *= Nef_polyhedron(circles[i]);
    }
    //debugging
    std::cout << "Consturct vertices Cartesian:" << std::endl;
    for (SVertex_const_iterator it = S.svertices_begin(); it != S.svertices_end();   it)
    {
        CPoint p(CGAL::to_double(it->point().x()),
            CGAL::to_double(it->point().y()),
            CGAL::to_double(it->point().z()));
        std::cout << p[0] << " " << p[1] << " " << p[2] << std::endl;
    }
    return S;
}
 

Входные pts данные состоят из 7 вершин с double декартовыми координатами типа. Этот код выдает мне ошибку при первом вызове в строке

 S *= Nef_polyhedron(circles[i]); // i==0 when error occurs
 

Здесь, во время выполнения, normal[0] = -0.0076270755275557253, normal[1] = 0.0062103059131449548, normal[2] = 0.016149138504648830 (цифры усекаются. Вы могли бы восстановить их, нормализуя вектор нормали, вызвав функцию sqrt). И вывод оболочки

 CGAL error: assertion violation!
Expression : c.has_on(p1)amp;amp;c.has_on(p2)
File       : C:devCGAL-5.1.1includeCGAL/Nef_S2/Sphere_segment.h
Line       : 53
Explanation:
Refer to the bug-reporting instructions at https://www.cgal.org/bug_report.html
 

При просмотре стека вызовов ошибка возникает из-за предварительного условия, которое проверяет, находится ли точка на сфере. Итак, мой вопрос здесь в том, как я могу заставить приведенный выше код работать?

Я попытался обойти эту проблему, изменив ядро на CGAL::Exact_predicates_exact_constructions_kernel , что позволяет выполнять точные вычисления. Теперь ошибка не выдается, но результат вычисления неверен, что

 Input Cartesian Coordinates:
0.733412 0.674178 0.0871214
0.722935 0.686566 0.0774091
0.721097 0.688361 0.0785988
0.710427 0.695728 0.106097
0.726191 0.675485 0.127931
0.739532 0.663825 0.111483
0.739734 0.66373 0.110705

Consturct vertices Cartesian:
1.15104e-274 1.07067e-274 2.02774e-275
0 0 0
0 0 0
1.40006e-136 1.37109e-136 2.09089e-137
6.67726e-235 6.13797e-235 7.93185e-236
5.9606e-63 5.69e-63 6.49699e-64
4.11434e-27 3.90736e-27 4.40548e-28
 

Координаты вершин выходной сферы не совпадают с входными. Я также пытался преобразовать тип normal[i] из double в CGAL::Exact_predicates_exact_constructions_kernel::FT , но результат остается тем же. Есть ли способ сделать эти координаты совпадающими или, по крайней мере, близкими?

После создания сферического многоугольника с тремя сегментами у вас будет шесть вершин, потому что вы исследуете не только свой многоугольник, но и весь DCEL, представляющий разбиение сферы на грани. Каждая грань ограничена внешним циклом полукругов (в простейшем случае), и она может содержать отверстия, также ограниченные циклами полукругов. Приведенный ниже код показывает, как исследовать грани такого раздела одну за другой:

 #include <iostream>
#include <sstream>
#include <string>
#include <vector>

#include <CGAL/Exact_rational.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Nef_polyhedron_S2.h>

using Kernel = CGAL::Simple_cartesian<CGAL::Exact_rational>;
using SPolygon = CGAL::Nef_polyhedron_S2<Kernel>;
using SPoint = SPolygon::Sphere_point;
using SSegment = SPolygon::Sphere_segment;

using std::cout;
using std::endl;

std::string pointAsString(const SPointamp; P)
{
  std::stringstream ss;
  ss << "(" << P.x() << "," << P.y() << "," << P.z() << ")";
  return ss.str();
}

std::string halfEdgeCycleAsString(const SPolygon::SHalfedge_const_iterator HEIT)
{
  std::stringstream ss;
  auto heit = HEIT; 
  while (true)
  {
    ss << pointAsString(heit->source()->point()) << " => " << pointAsString(heit->target()->point());
    heit = heit->snext();
    if (heit != HEIT)
    {
      ss << ", ";
    }
    else
    {
      break;
    }
  }
  return ss.str();
}

void printFaces(const SPolygonamp; P)
{
  unsigned i = 0;
  for (auto fit = P.sfaces_begin(); fit != P.sfaces_end();   fit)
  {
    cout << "face #" << i   << ":" << endl;
    for (auto fcit = fit->sface_cycles_begin(); fcit != fit->sface_cycles_end();   fcit)
    {
      cout << "  face boundary:" << endl;
      if (fcit.is_svertex())
      {
        const SPolygon::SVertex_const_handle vit = fcit;
        cout << "    vertex " << pointAsString(vit->point()) << endl;
      }
      else if (fcit.is_shalfedge())
      {
        const SPolygon::SHalfedge_const_handle heit = fcit;
        cout << "    halfedge cycle [" << halfEdgeCycleAsString(heit) << "]" << endl;
      }
      else if (fcit.is_shalfloop())
      {
        const SPolygon::SHalfloop_const_handle hlit = fcit;
        cout << "    halfloop (" << hlit->circle() << ")" << endl;
      }
    }
  }
}

int main()
{
  const SPoint p1(1, 0, 0),  p2(0, 0, 1), p3(1, -1, 0);
  const std::vector<SSegment> v{{p1, p2}, {p2, p3}, {p3, p1}};
  const SPolygon p(v.cbegin(), v.cend());
  p.print_statistics();
  printFaces(p);
}
 

Вывод этой программы:

 Sphere Map - Statistics
|V| = 6
|E| = 14
|L| = 0
|F| = 3
|Fcs| = 3

face #0:
  face boundary:
    halfedge cycle [(-1,0,0) => (0,-1,0), (0,-1,0) => (1,-1,0), (1,-1,0) => (0,0,1), (0,0,1) => (1,0,0), (1,0,0) => (0,1,0), (0,1,0) => (-1,0,0)]
face #1:
  face boundary:
    halfedge cycle [(1,0,0) => (0,0,1), (0,0,1) => (1,-1,0), (1,-1,0) => (1,0,0)]
face #2:
  face boundary:
    halfedge cycle [(1,-1,0) => (0,-1,0), (0,-1,0) => (-1,0,0), (-1,0,0) => (0,1,0), (0,1,0) => (1,0,0), (1,0,0) => (1,-1,0)]
 

Еще одна проблема в вашем коде — вы не должны использовать double числа с Simple_cartesian ядром, потому что это ядро требует рациональных чисел, поддерживающих точные вычисления. Кроме того, если у вас есть нецелочисленные входные данные, вам нужно преобразовать их в рациональные числа:

 SPoint p(1, 0, CGAL::Exact_rational(119, 31)); 
 

Если вместо этого вы используете 119/31 here , эта z-координата будет округлена до 3 .

ОБНОВЛЕНИЕ. Вы можете перемещаться только по той части DCEL, которая важна для вас — например, по всем ребрам созданного вами многоугольника. Один из способов сделать это приведен ниже:

 int main()
{
  const SPoint p1(1, 0, 0),  p2(0, 0, 1), p3(1, -1, 0);
  const std::vector<SSegment> v{{p1, p2}, {p2, p3}, {p3, p1}};
  const SPolygon p(v.cbegin(), v.cend());
  // ------ locate the starting vertex of this polygon and output its halfedge cycle
  const auto h = p.locate(p1);
  SPolygon::SVertex_const_handle vh = nullptr;
  if (CGAL::assign(vh, h))
  {
    cout << halfEdgeCycleAsString(vh->out_sedge()) << endl;
  }
  else
  {
    cout << "Vertex is not found" << endl;
  }
}