3

The Problem

Using CGAL, I am creating a random point cloud in d-Dimensional space like so: 

#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;
    }
...

Then I'm using the dD Convex Hulls and Delaunay Triangulations package to compute a Delaunay triangulation of that point cloud in d-Dimensions:

...

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);
}
...

Then I would like to write the triangulation to a .txt file in a format like this:

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

Where v1, ..., vd+1 are the global indices of the vertices in one simplex. Like it is described in this question.

Now, for 2D and 3D, there are the CGAL::Triangulation_vertex_base_with_info_2 and CGAL::Triangulation_vertex_base_with_info_3 classes which can be used to add any additional info to a vertex. Like shown here(stackoverflow question).

However, I could not find a clear explanation how to do the same in d-Dimensions.

The Search

Searching Google gave me the following results:

  • A Master's thesis Implementing four-dimensional triangulations in CGAL. Here, a user defined simple geometric kernel is implemented, which provides the necessary objects for working in R^4 (Point_4, Orientation_4 etc.). An implementation of a class Triangulation_vertex_base_4 is mentioned, but I couldn't find the code to see how it's done.

  • A paragraph in the CGAL manual mentioning the flexibility of the design of the Vertex and Face classes used in 2D and 3D triangulations. There is an example there of extending the Vertex Base class with custom data, but if I understood correctly it only works with the already defined Triangulation_vertex_base_with_info for 2D and 3D.

The Tried Solutions

  • I tried modifying the following classes:

    • Triangulation_ds_vertex_base_3.h
    • Triangulation_vertex_base_3.h
    • Triangulation_vertex_base_with_info_3.h
    • Dummy_tds_3.h
      Basically, replacing the 3 with d but I'm not sure I'm doing the right thing. Then there is the Triangulation_data_structure_3.hwhich looks pretty challenging and I wouldn't go down the road of modifying it for d-Dimensions. There must be an easier way?

  • I tried a simple solution of looking up the point for each vertex in the point cloud container that I created in the beginning. The idea is to look at the distance from the point given by the vertex (Vertex_handle::point()) to all of the points in the std::vector<Point> container and choose the index of the closest point in that container. Here is some code:

    ...
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;
}

The problem is that when I try to create that vector, in order to calculate the distance, CGAL throws an unhandled exception (Assertion_exception). I am guessing that this occurs when the two points (from the vhandle and the point cloud) are actually the same. But I might be wrong. Anyways, I can't compute the distance, and Google hasn't been of any help for this either. (I have also tried with creating a squared_distance_d_object() but it does the same thing).

The Question

Is there an easy way to plug additional information to vertices in d-Dimensional triangulations in CGAL?

Is there another way to get the global indices of the vertices?

What am I doing wrong with the euclidean distance?

Any help would be greatly appreciated.

Edit 1

I solved the problem for the euclidean distance approach. The point that I was getting with Point pt = vh->point() was not compatible with the Cartesian_d::Point_d class. This created problems when creating the vector. I solved the problem with using Delaunay_d::point_at_simplex(si, i) method, which I somehow overlooked before. Now the search for the indices works fine.

It would be nice, however, to find out how to add data to the vertices in d-dimensions.

Community
  • 1
  • 1
Stefan Bojarovski
  • 31
  • 1
  • 4

0 Answers0