Triangulation__vertex__base__with__info__2 for Delaunay_triangulation_on_sphere_2

[AlekseyMuraviev]

Do I understand correctly that now there is no way to use something like CGAL::Triangulation_on_sphere_vertex_base_2 for triangulation on a sphere, which can be used when constructing triangulation on a plane?

[MaelRL]

The triangulation on the sphere is templated by the same template parameters (geometric traits + triangulation data structure) as other CGAL triangulations. You can use a "with_info" vertex base in the stack :

...
using Vbb = Triangulation_vertex_base_with_info_2<Your_info, Gt>;
using Vb = Triangulation_on_sphere_vertex_base_2<Gt, Vbb>;
using Tds = Triangulation_data_structure_2<Vb, Triangulation_on_sphere_face_base_2<Gt> > >;
using Dt = Delaunay_triangulation_on_sphere_2<Gt, Tds>;
...

There are no insert() overloads with zip iterator, but you can just extract code from the range insertion function (https://github.com/CGAL/cgal/blob/master/Triangulation_on_sphere_2/include/CGAL/Delaunay_triangulation_on_sphere_2.h#L684=) and add your info after insertion:

...
Vertex_handle v = insert(p, hint);
if(v != Vertex_handle()) // could happen if the point is not on the sphere
{
  hint = v->face(); 
  v->info() = ...;
}
...

[AlekseyMuraviev]

Thank you for quick answer!

But when I try to construct a triangulation like this:

using K = CGAL::Exact_predicates_inexact_constructions_kernel;                
using Gt = CGAL::Delaunay_triangulation_on_sphere_traits_2<K>;
using Vbb = CGAL::Triangulation_vertex_base_with_info_2<size_t, Gt>;
using Vb = CGAL::Triangulation_on_sphere_vertex_base_2<Gt, Vbb>;
using Tds = CGAL::Triangulation_data_structure_2<Vb, CGAL::Triangulation_on_sphere_face_base_2<Gt>>;
using Dt = CGAL::Delaunay_triangulation_on_sphere_2<Gt, Tds>;

I get an error:

Error C2385 ambiguous access of 'Point_2' D:\cgal-5.3\include\Triangulation_2\include\CGAL\Triangulation_vertex_base_2.h 33
I'm using cgal-5.3 and Microsoft Visual Studio Professional 2019 Version 16.10.4.

[MaelRL]

Sorry my vertex stack was wrong: I thought Triangulation_vertex_base_with_info_2 inherited Triangulation_ds_vertex_base_2 by default, but it is Triangulation_vertex_base_2 actually.

Use

using Vbb = CGAL::Triangulation_on_sphere_vertex_base_2<Gt>;
using Vb = CGAL::Triangulation_vertex_base_with_info_2<std::size_t, Gt, Vbb>;

instead of

using Vbb = CGAL::Triangulation_vertex_base_with_info_2<std::size_t, Gt>;
using Vb = CGAL::Triangulation_on_sphere_vertex_base_2<Gt, Vbb>;

(so, just a switch)

[AlekseyMuraviev]

Thanks, everything compiled!

But I want to store the index of the point in the original array, i.e. even before spatial_sort_on_sphere(). The code below doesn't work the way I want:

using Dt = CGAL::Delaunay_triangulation_on_sphere_2<Gt, Tds>;
using Point = Dt::Point;

std::vector<OtherPoint> original_points = get_points();

std::vector<Point> points = some_transformation(original_points);

Gt traits(Point(0, 0, 0), 1);
Dt dt(traits);

CGAL::spatial_sort_on_sphere(points.begin(), points.end(), traits, CGAL::square(traits.radius()), traits.center());

Dt::Face_handle hint;
for (size_t i = 0; i < points.size(); i++)
{
    Dt::Vertex_handle v = dt.insert(points[i], hint);
    if (v != Dt::Vertex_handle())
    {
        hint = v->face();
        v->info() = i;
    }
}

I tried using CGAL::Spatial_sort_traits_adapter_3 but it doesn't seem to work.

Probably I could use something like Point_3_with_iterator. But it seems that in this case I will have to make an extra copy of the original data, and I would like to avoid this.

And do I understand correctly that random_shuffle is not needed here? After all, it is already in spatial_sort_on_sphere().

[MaelRL]

An easy and not-so-costly fix would be to keep a map that associates the geometric position of the point to your info and then it's just v->info() = my_info[v->point].

It should also work with the adapted spatial sort, have a look at the Delaunay_triangulation_3 vector for inspiration, and you can send me a gist of your code if you are still stuck.

It looks like that random_shuffle is redundant indeed.

[AlekseyMuraviev]

This code works, but maybe there is a simpler solution?

I will then make the resulting mesh from the originPoints and triangles.

std::vector<Eigen::Vector3d> generateRandomPoints(size_t numPoints)
{
    using namespace std;
    random_device rd;
    mt19937 gen(rd());
    normal_distribution<double> dirDis;
    uniform_real_distribution<double> distDis(100, 150);
    vector<Eigen::Vector3d> points(numPoints);
    for (size_t i = 0; i < numPoints; i++)
    {
        auto& p = points[i];
        p.x() = dirDis(gen);
        p.y() = dirDis(gen);
        p.z() = dirDis(gen);
        p.normalize();
        p *= distDis(gen);
    }
    return points;
}

int main()
{
    using K = CGAL::Exact_predicates_inexact_constructions_kernel;
    using Gt = CGAL::Delaunay_triangulation_on_sphere_traits_2<K>;
    using Vbb = CGAL::Triangulation_on_sphere_vertex_base_2<Gt>;
    using Vb = CGAL::Triangulation_vertex_base_with_info_2<size_t, Gt, Vbb>;
    using Tds = CGAL::Triangulation_data_structure_2<Vb, CGAL::Triangulation_on_sphere_face_base_2<Gt>>;
    using Dt = CGAL::Delaunay_triangulation_on_sphere_2<Gt, Tds>;
    using Point = Dt::Point;

    struct PointWithIndex : public Point
    {
        PointWithIndex(const Point& p, size_t index)
            : Point(p)
            , index(index)
        {}
        size_t index;
    };

    const std::vector<Eigen::Vector3d> originPoints = generateRandomPoints(1000);

    std::vector<PointWithIndex> points;
    points.reserve(originPoints.size());
    for (size_t i = 0; i < originPoints.size(); i++)
    {
        const auto& p = originPoints[i];
        Eigen::Vector3d q(p);
        q.normalize();
        points.emplace_back(Point(q.x(), q.y(), q.z()), i);
    }

    Gt traits(Point(0, 0, 0), 1);
    Dt dt(traits);

    CGAL::spatial_sort_on_sphere(points.begin(), points.end(), traits, CGAL::square(traits.radius()), traits.center());

    Dt::Face_handle hint;
    for (size_t i = 0; i < points.size(); i++)
    {
        Dt::Vertex_handle v = dt.insert(*(Point*)(&points[i]), hint);
        if (v != Dt::Vertex_handle())
        {
            hint = v->face();
            v->info() = points[i].index;
        }
    }

    std::vector<Eigen::Vector3i> triangles;
    triangles.reserve(dt.number_of_solid_faces());
    for (const auto& fh : dt.all_face_handles())
    {
        if (!fh->is_ghost())
        {
            triangles.emplace_back(fh->vertex(0)->info(), fh->vertex(1)->info(), fh->vertex(2)->info());
        }
    }
}