CGAL Как использовать глобальные индексы для вершин в dDimensional Triangulation?

Эта проблема

Используя CGAL, я создаю случайное облако точек в d-мерном пространстве следующим образом:

#include <CGAL/Cartesian_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Delaunay_d.h>
#include <CGAL/basic.h>

typedef CGAL::Cartesian_d<double> K;    //Kernel
typedef CGAL::Point_d<K> Point;
typedef std::vector<Point> PointSet;

typedef CGAL::Delaunay_d<K> Delaunayd;
typedef Delaunayd::Simplex_handle S_handle;
typedef Delaunayd::Simplex_iterator S_iterator;
typedef Delaunayd::Vertex_handle V_handle;
typedef Delaunayd::Vertex_iterator V_iterator;

int main()
{
    int mDim = 3;    // will be higher, I use it now for simplicity
    int mNum = 10;
    PointSet mPs;
    CGAL::Random_points_in_cube_d<Point> rpg(d, range);

    for(int i = 0; i < mNum; ++i)
    {
        mPs.push_back(*rpg);
        ++rpg;
    }
...

Затем я использую пакет dD Convex Hulls and Delaunay Triangulations для вычисления триангуляции Делоне этого облака точек в d-измерениях:

...

Delaunayd * mT = new Delaunayd(mDim);

PointSet::iterator pt_it = mPs.begin();
PointSet::iterator pt_it_end = mPs.end();
for (; pt_it != pt_it_end; ++pt_it)
{
    V_handle v = mT->insert(*pt_it);
}
...

Затем я хотел бы записать триангуляцию в файл .txt в таком формате:

d        //number of dimensions
nv ns        // nv - number of vertices(points), ns - number of simplices
// followed by nv lines
coord_0_1  coord_0_2  coord_0_3  ... coord_0_d
...
coord_nv_1 coord_nv_2 coord_nv_3 ... coord_nv_d
//followed by ns lines
v1 v2 v3 ... vd+1
...
v1 v2 v3 ... vd+1

Где v1, ..., vd+1 — глобальные индексы вершин в одном симплексе. Как это описано в этом вопросе.

Теперь для 2D и 3D существуют классы CGAL::Triangulation_vertex_base_with_info_2 и CGAL::Triangulation_vertex_base_with_info_3, которые можно использовать для добавления любой дополнительной информации к вершине. Как показано здесь (вопрос о стеке).

Однако я не смог найти внятного объяснения, как сделать то же самое в d-Dimensions.

Поиск

Поиск в Google дал мне следующие результаты:

диссертация < em>Реализация четырехмерных триангуляций в CGAL. Здесь реализовано определяемое пользователем простое геометрическое ядро, которое предоставляет необходимые объекты для работы в R^4 (Point_4, Orientation_4 и т.д.). Упоминается реализация класса Triangulation_vertex_base_4, но я не смог найти код, чтобы увидеть, как это делается.
параграф в руководстве CGAL, в котором упоминается гибкость дизайна Vertex и Face классы, используемые в 2D и 3D триангуляции. Там есть пример расширения класса Vertex Base пользовательскими данными, но, если я правильно понял, он работает только с уже определенным Triangulation_vertex_base_with_info для 2D и 3D.

Проверенные решения

Я попытался изменить следующие классы:
- Triangulation_ds_vertex_base_3.h
- Triangulation_vertex_base_3.h
- Triangulation_vertex_base_with_info_3.h
- Dummy_tds_3.h
  По сути, замена 3 на d, но я не уверен, что поступаю правильно. Затем есть Triangulation_data_structure_3.h, который выглядит довольно сложно, и я бы не стал модифицировать его для d-Dimensions. Должен же быть более простой способ?

Я попробовал простое решение поиска точки для каждой вершины в контейнере облака точек, который я создал вначале. Идея состоит в том, чтобы посмотреть на расстояние от точки, заданной вершиной (Vertex_handle::point()), до всех точек в контейнере std::vector<Point> и выбрать индекс ближайшей точки в этом контейнере. Вот код:

...
int ns = T->all_simplices().size();
vertIndices [ns][mDim + 1];  // number of simplices and d+1 vertices per simplex
S_iterator s_it = mT->simplices_begin();
S_iterator s_it_end = mT->simplices_end();
int currSimp = 0;
int numVert = mDim + 1; // number of vertices per simplex
for (; s_it != s_it_end; ++s_it)
{
    for(int i = 0; i < numVert; ++i)
    {
        V_handle vh = mT->vertex_of_simplex(s_it, i);
        int idx = findIndexOf(vh, mPs);
        vertIndices[currSimp][i] = idx;
    }
    currSimp++;
}
...
int findIndexOf(V_handle vh, PointSet mPs)
{
    Point pt = vh->point();
    int idx = -1;
    int count = 0;
    PointSet::iterator pt_it = mPs.begin();
    PointSet::iterator pt_it_end = mPs.end();
    for (; pt_it != pt_it_end; ++pt_it)
    {
        Point opt = *pt_it;
        K::Vector_d v = pt - opt;         // The problem is here!
        double dst = v.squared_length();
        if (dst < 0.01)
        {
            idx = count;
        }
        count++;
    }
    return idx;
}

Проблема в том, что когда я пытаюсь создать этот вектор, чтобы вычислить расстояние, CGAL выдает необработанное исключение (Assertion_exception). Я предполагаю, что это происходит, когда две точки (из vhandle и облака точек) на самом деле одинаковы. Но я могу ошибаться. В любом случае, я не могу вычислить расстояние, и Google тоже не помог. (Я также пытался создать squared_distance_d_object(), но он делает то же самое).

Вопрос

Есть ли простой способ добавить дополнительную информацию к вершинам в d-мерных триангуляциях в CGAL?

Есть ли другой способ получить глобальные индексы вершин?

Что я делаю не так с евклидовым расстоянием?

Любая помощь будет принята с благодарностью.

Изменить 1

Я решил проблему для подхода евклидова расстояния. То, что я получил с помощью Point pt = vh->point(), было несовместимо с классом Cartesian_d::Point_d. Это создало проблемы при создании вектора. Решил проблему методом Delaunay_d::point_at_simplex(si, i), который раньше как-то упускал из виду. Теперь поиск по индексам работает нормально.

Однако было бы неплохо узнать, как добавлять данные к вершинам в d-измерениях.

Stefan Bojarovski 25.06.2014 источник