How startinpy works

Original code written in rust

The library for calculating the Delaunay triangulation is originally written in Rust, it is called ‘startin’ and its source code is open.

Robust arithmetic for the geometric predicates is used (Shewchuk’s predicates, well the Rust port of the code), so startin/py is robust and shouldn’t crash (touch wood).

Insertion + deletion are possible

It uses an incremental algorithm for the construction of a Delaunay triangulation (constraints are not supported), that is each point is inserted one after another and the triangulation is updated between each insertion. The algorithm is based on flips.

The deletion of a vertex is also possible. The algorithm implemented is a modification of the one of Mostafavi, Gold, and Dakowicz (2003). The ears are filled by flipping, so it’s, in theory, more robust. I have also extended the algorithm to allow the deletion of vertices on the boundary of the convex hull. The algorithm is sub-optimal, but, in practice, the number of neighbours of a given vertex in a DT is only 6, so it doesn’t really matter.

The data structure

The data structure of the Rust code is a cheap implementation of the star-based structure defined in Blandford et al. (2003); cheap because the link of each vertex is stored a simple array and not in an optimised blob like they did. It results in a pretty fast library (comparison will come at some point), but it uses more space than the optimised one.

However, notice that the stars are not exposed in startinpy, to keep it a simple and higher-level library.

The data structure of startinpy is composed of 2 arrays:

  1. an array of Points, where each entry is an array of 3 floats (x-coordinate, y-coordinate, z-coordinate)

  2. an array of Triangles, where each Triangle is an array of 3 integers, the values of the indices of the 3 vertices (ordered counter-clockwise) in the array of Points (startinpy.DT.points(), which is 0-based, 0 being the infinite vertex).

A Vertex is an integer, it is the index in the array of points (startinpy.DT.points(), which is 0-based).

If you delete a vertex (with startinpy.DT.remove()), then the entry in the array of Points is not deleted (this would be slow because arrays are contiguous and a lot of copying would be necessary), instead the vertex/point is flagged as being removed and none of the Triangles will refer to it (startinpy.DT.is_vertex_removed() can be used to test).

Infinite vertex and triangles

_images/infinite.png

The implementation of startinpy has one infinite vertex and infinite triangles, this greatly simplifies the algorithms and ensures that one can insert new points outside the convex hull of a dataset (or even delete some vertices on the boundary of the convex hull). The CGAL library also does this, and the internal workings are well explained here.

The infinite vertex is the first vertex in the array of points (startinpy.DT.points()), and, thus, it has the index of 0 (zero). It has infinite coordinates ([inf inf inf]), those are of type numpy infinity.

An infinite triangle is a triangle having the infinite vertex as one of its vertices; a finite triangle doesn’t have the infinite vertex is part of the triangulation of the dataset.

In the figure, notice that there are 5 finite triangles (126, 236, 346, 456, 516), but the data structure actually stores 5 extra infinite triangles (102, 150, 540, 304, 203). Those are adjacent to the 5 edges on the boundary of the convex hull of the dataset. You can conceptualise the triangulation as being embedded on a sphere, and the infinite vertex is on the other side.

Some examples of the data structure and infinity

_images/tr.png

For instance, consider this 5-vertex Delaunay triangulation:

import startinpy
import numpy as np

np.set_printoptions(precision=10)

t = startinpy.DT()
t.insert_one_pt([0.5, 0.5, 1.0])
t.insert_one_pt([0.0, 0.0, 2.0])
t.insert_one_pt([1.0, 0.0, 3.0])
t.insert_one_pt([1.0, 1.0, 4.0])
t.insert_one_pt([0.0, 1.0, 5.0])

print(t.points)
print(t.triangles)

Which outputs this below. Notice first that there are a total of 6 vertices: the 5 we inserted plus the infinite vertex (at index-0 with infinity coordinates [inf inf inf]). Notice also that no finite triangle refers to the vertex 0.

[[inf inf inf]
 [0.5 0.5 1. ]
 [0.  0.  2. ]
 [1.  0.  3. ]
 [1.  1.  4. ]
 [0.  1.  5. ]]
[[1 2 3]
 [1 3 4]
 [1 4 5]
 [1 5 2]]

However, startinpy stores internally infinite triangles. For instance, if you retrieve the triangles incident to a given vertex on the convex hull:

re = t.incident_triangles_to_vertex(2)
for each in re:
    print(each)

[2 0 3]
[2 3 1]
[2 1 5]
[2 5 0]

you will notice that 2 triangles are infinite: [2 0 3] and [2 5 0].

Also, if you remove one vertex (eg the one in the middle of the square, vertex 1), observe that now its coordinates are “Not a Number (nan)”, and that no triangle in the DT refers to it anymore:

t.remove(1)
print(t.points)
print(t.triangles)
print(t.is_vertex_removed(1))

[[inf inf inf]
 [nan nan nan]
 [ 0.  0.  2.]
 [ 1.  0.  3.]
 [ 1.  1.  4.]
 [ 0.  1.  5.]]
[[2 3 4]
 [2 4 5]]
True

Finally, you can remove the unused/deleted vertices from the startinpy.DT.points() array by using startinpy.DT.collect_garbage(), which will assign a new ID to most vertices, and triangles will be updated too. Notice that now 5 vertices are in the array, and only 2 finite triangles are in the DT.

t.collect_garbage()
print(t.points)
print(t.triangles)

[[inf inf inf]
 [ 0.  0.  2.]
 [ 1.  0.  3.]
 [ 1.  1.  4.]
 [ 0.  1.  5.]]
[[1 2 3]
 [1 3 4]]