Why does CGAL's Delaunay_triangulation_3 return a vertex handle that already exists when it calls insert()?

[liuxiaozhu01]

Issue Details

I am trying to build delaunay triangulation inserting points one by one. I've already knew that when inserting a point that already located in the target position will return a vertex handle to that old point. But however, there seem to be other situations where this behavior occurs as well. I do need different handles for different points in the coordinates, while new added point whose coordinate already exit would not be added. Any one can helps me?

Source Code

Here is some key code.

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_info_3<unsigned int, K> Vb;
typedef CGAL::Triangulation_data_structure_3<Vb> Tds;
typedef CGAL::Delaunay_triangulation_3<K, Tds> Triangulation;

typedef K::FT FT;
typedef CGAL::Parallel_if_available_tag Concurrency_tag;

typedef Triangulation::Cell_handle Cell_handle;
typedef Triangulation::Vertex_handle Vertex_handle;
typedef Triangulation::Locate_type Locate_type;
typedef Triangulation::Point Point;

class Builder {
public:
    Builder();
    ~Builder();

    void add_points(size_t num_points, float3 *points);

private:
    std::unordered_map<size_t, Vertex_handle> idx_to_vertex;
    Triangulation delaunay_tetra;
    unsigned int cur_size = 0;
};

And when i implement the function add_points() like this

void Builder::add_points(size_t num_points, float3 *points) {
    /* add  points sequentially */
    for (int i = 0; i < num_points; i++) {
        auto p = Point(points[i].x, points[i].y, points[i].z);
        auto vertex_handle = delaunay_tetra.insert(p);
        if (idx_to_vertex.find(vertex_handle->info()) != idx_to_vertex.end()) {
            std::cout << "\033[33m" << "Warning: There is a point already located in the target position (" << points[i].x << ", " << points[i].y << ", " << points[i].z << ")!" << std::endl;
            continue;
        }
        vertex_handle->info() = cur_size;
        idx_to_vertex[cur_size] = vertex_handle;
        cur_size++;
    }
}

some points would not be added because the vertex_handle returned by delaunay_tetra.insert(p) has already existed, but its corrdinate seems to be different from any points. But when i implemented it without checking like this

void TetrahedraBuilder::add_points(size_t num_points, float3 *points) {
    /* add  points sequentially */
    for (int i = 0; i < num_points; i++) {
        auto p = Point(points[i].x, points[i].y, points[i].z);
        auto vertex_handle = delaunay_tetra.insert(p);
        vertex_handle->info() = cur_size;
        idx_to_vertex[cur_size] = vertex_handle;
        cur_size++;
    }
}

everything goes well.

Environment

• Linux, 64 bits
• Compiler: g++ 8.4
• Release or debug mode: debug
• Specific flags used (if any): nope
• CGAL version: 5.5.2
• Boost version: libboost1.71-dev:amd64
• Other libraries versions if used (Eigen, TBB, etc.):

[liuxiaozhu01]

I made more attempts. A function that used to put vertices' coordinates into an array acoording to their indices( the unsigned int info as typedef CGAL::Triangulation_vertex_base_with_info_3<unsigned int, K> Vb;) is shown below

std::vector<float3> Builder::export_points_idx() {
    std::vector<float3> result(cur_size);
    float *result_float = reinterpret_cast<float*>(result.data());
    size_t cnt = 0;
    for (auto vertex_handle : delaunay_tetra.finite_vertex_handles()) {
        const auto point = vertex_handle->point();
        const unsigned int i = vertex_handle->info();
        result_float[i * 3 + 0] = point.x();
        result_float[i * 3 + 1] = point.y();
        result_float[i * 3 + 2] = point.z();
        cnt += 1;
    }
    std::cout<< "current points number: " << cnt << std::endl;
    return result;
}

I notice that the length of the list that delaunay_tetra.finite_vertex_handles() returns, not equal to cur_size used to record the number of vertice I added to the delaunay triangulation. I am totally confused now... Please help me!

[MaelRL]

When you write

        auto vertex_handle = delaunay_tetra.insert(p);
        if (idx_to_vertex.find(vertex_handle->info()) != idx_to_vertex.end()) {

if the vertex is new, the info field is uninitialized.

If you want to check if the triangulation has gained a new point, you can either split your insertion to call locate() first, or simply check the number of vertices of the triangulation after the insertion : if it did not increase, a vertex already existed at this position.

I notice that the length of the list that delaunay_tetra.finite_vertex_handles() returns, not equal to cur_size used to record the number of vertice I added to the delaunay triangulation.

With which version of add_points? Can you please post self-contained examples that can be easily run?

[liuxiaozhu01]

Thanks for your kindness! when the vertex is new, the info field is unintialized and i chose to initialize this field after checking the vertex handle returned by insert(p) is a new one. However, not all the time and that's what im confusing.

the tested pcd i used is download from here https://github.com/jkulhanek/tetra-nerf/blob/master/tests/assets/bottle.ply Here is my example code

triangulation.h

#include <CGAL/Delaunay_triangulation_3.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_3.h>
#include <CGAL/Triangulation_vertex_base_with_info_3.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/convex_hull_3_to_face_graph.h>
#include <CGAL/Surface_mesh/IO/PLY.h>

#include <CGAL/compute_average_spacing.h>

#include <vector>
#include <unordered_map>

#include "utils/vec_math.h" // include <sutil/vec_math.h>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_with_info_3<unsigned int, K> Vb;
typedef CGAL::Triangulation_data_structure_3<Vb> Tds;
typedef CGAL::Delaunay_triangulation_3<K, Tds> Triangulation;

typedef K::FT FT;
typedef CGAL::Parallel_if_available_tag Concurrency_tag;

typedef Triangulation::Cell_handle Cell_handle;
typedef Triangulation::Vertex_handle Vertex_handle;
typedef Triangulation::Locate_type Locate_type;
typedef Triangulation::Point Point;

typedef CGAL::Surface_mesh<Point> Surface_mesh;

// std::vector<uint4> triangulate(size_t num_points, float3 *points);
// float find_average_spacing(size_t num_points, float3 *points);

class Builder {
public:
    // Builder(size_t num_points, float3 *points);
    Builder();
    ~Builder();

    // std::vector<uint4> build();
    // float get_average_spacing();
    void add_points(size_t num_points, float3 *points);
    std::vector<uint4> build_tetrahedra();

    std::vector<float3> export_points_idx();

// #endif

private:
    std::unordered_map<size_t, Vertex_handle> idx_to_vertex;
    Triangulation delaunay_tetra;
    unsigned int cur_size = 0;
};

triangulation.cpp

#include "triangulation.h"

#include <cassert>
#include <fstream>
#include <iostream>
#include <list>
#include <string>
#include <vector>

Builder::Builder() {
    idx_to_vertex == std::unordered_map<size_t, Vertex_handle>();
    delaunay_tetra = Triangulation();
    cur_size = 0;
}

Builder::~Builder() {
}

void Builder::add_points(size_t num_points, float3 *points) {
    /* add  points sequentially */
    size_t cnt_strange = 0;
    for (int i = 0; i < num_points; i++) {
        auto p = Point(points[i].x, points[i].y, points[i].z);
        auto vertex_handle = delaunay_tetra.insert(p);
        auto iter = idx_to_vertex.find(vertex_handle->info());
        if (iter != idx_to_vertex.end()) {
            cnt_strange += 1;
            // continue;    // The difference between the different versions is whether or not this line is executed
        }
        vertex_handle->info() = cur_size;
        idx_to_vertex[cur_size] = vertex_handle;
        cur_size++;
    }
    std::cout << "There are " << cnt_strange << " vertex handle already exist!" << std::endl;
    std::cout << "There are " << delaunay_tetra.number_of_vertices() << " vertices in the triangulation!" << std::endl;

}

std::vector<uint4> Builder::build_tetrahedra() {
    std::vector<uint4> cells(delaunay_tetra.number_of_finite_cells());
    unsigned int *cells_uint = reinterpret_cast<unsigned int*>(cells.data());

    size_t i = 0;
    for (auto cell : delaunay_tetra.finite_cell_handles()) {
        for (int j = 0; j < 4; j++) {
            cells_uint[i * 4 + j] = cell->vertex(j)->info();
        }
        i++;
    }
    return cells;
}

std::vector<float3> Builder::export_points_idx() {
    std::vector<float3> result(cur_size);
    float *result_float = reinterpret_cast<float*>(result.data());
    // std::cout<<"hello in here 1!" << std::endl;

    // size_t i = 0;
    size_t cnt = 0;
    // for (auto vertex_handle : delaunay_tetra.all_vertex_handles()) {     /* the first element all_vertex_handles() returns donot contain any data! fuck! */
    for (auto vertex_handle : delaunay_tetra.finite_vertex_handles()) {

        const auto point = vertex_handle->point();
        // result_float[vertex_handle->info()] = make_float3(point.x(), point.y(), point.z());
        const unsigned int i = vertex_handle->info();
        // std::cout<<"hello in index " << i << std::endl;
        result_float[i * 3 + 0] = point.x();
        result_float[i * 3 + 1] = point.y();
        result_float[i * 3 + 2] = point.z();
        cnt += 1;
    }
    std::cout<< "current points number: " << cnt << std::endl;
    std::cout<< "Recorded cur_size: " << cur_size << std::endl;

    return result;
}

main.cpp

#include "triangulation.h"

#include <iostream>
#include <pcl/common/io.h>
#include <pcl/io/ply_io.h>
#include <pcl/point_types.h>

int main() {
    std::string filename = "./bottle.ply";
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);

    pcl::io::loadPLYFile<pcl::PointXYZ>(filename, *cloud);

    std::vector<float3> points;
    for (const auto& point : *cloud) {
        float x = point.x;
        float y = point.y;
        float z = point.z;
        points.push_back(make_float3(x, y, z));
    }

    Builder builder;
    builder.add_points(points.size(), points.data());
    auto cells = builder.build_tetrahedra();
    auto points_idx = builder.export_points_idx();

    return 0;
}

[liuxiaozhu01]

When you write

        auto vertex_handle = delaunay_tetra.insert(p);
        if (idx_to_vertex.find(vertex_handle->info()) != idx_to_vertex.end()) {

if the vertex is new, the info field is uninitialized.

If you want to check if the triangulation has gained a new point, you can either split your insertion to call locate() first, or simply check the number of vertices of the triangulation after the insertion : if it did not increase, a vertex already existed at this position.

I notice that the length of the list that delaunay_tetra.finite_vertex_handles() returns, not equal to cur_size used to record the number of vertice I added to the delaunay triangulation.

With which version of add_points? Can you please post self-contained examples that can be easily run?

Thank you for your help. I follow your second advice that check the number of vertices of the triangulation after the insertion. And there is indeed some points(the number of vertices of the triangulation increase) get an already existing vertex handle returned from insert() and i don't know whether it will influence the info field of vertex, resulting an messy squence which not what i expected.