CGAL Как использовать глобальные индексы на вершинах в трехмерной триангуляции?

Эта проблема

Используя 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 выпуклые оболочки и триангуляции Делоне пакет для вычисления триангуляции Делоне для этого облака точек в 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);
}
...

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

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 классы, которые можно использовать для добавления любой дополнительной информации в вершину. Как показано Вот(вопрос stackoverflow).

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

Поиск

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

• Магистр Тезис Реализация четырехмерных триангуляций в 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-измерениях.