2
votes

I was reading about Delaunay (scipy) and came across the code:

import numpy as np
points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])

from scipy.spatial import Delaunay
tri = Delaunay(points)

import matplotlib.pyplot as plt
plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.plot(points[:,0], points[:,1], 'o')
plt.show()

As far as I understand, a simplex is the generalization of a triangle to higher dimensions.

I don't understand the meaning of the code below and would like help understanding it:

# Point indices and coordinates for the two triangles forming the triangulation:

tri.simplices
array([[3, 2, 0],
       [3, 1, 0]], dtype=int32)

points[tri.simplices]
array([[[ 1. ,  1. ],
        [ 1. ,  0. ],
        [ 0. ,  0. ]],
       [[ 1. ,  1. ],
        [ 0. ,  1.1],
        [ 0. ,  0. ]]])

Triangle 0 is the only neighbor of triangle 1, and it’s opposite to vertex 1 of triangle 1:


tri.neighbors[1]
# array([-1,  0, -1], dtype=int32)

points[tri.simplices[1,1]]
array([ 0. ,  1.1])

Thanks!

1

1 Answers

2
votes

This code is creating a Delaunay triangulation from four vertices containing two triangles. The triangulation looks like this:

triangulation

The code begins by defining the four vertices in an array:

points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])

Next, scipy constructs a Delaunay triangulation for these points:

from scipy.spatial import Delaunay
tri = Delaunay(points)

Now, tri.simplices contains a list of triangles (in this 2D case) in the Delaunay triangulation. Each triangle is represented as three integers: each value represents a index in to the original points array.

tri.simplices
array([[3, 2, 0],
       [3, 1, 0]], dtype=int32)

So [3,2,0] is the triangle between vertex 3 (1,1), vertex 2 (1,0) and vertex 0 (0,0). The next code connects the points and tri data structures to compute the coordinates of the vertices of each triangle, eliminating the indirection:

points[tri.simplices]
array([[[ 1. ,  1. ],
        [ 1. ,  0. ],
        [ 0. ,  0. ]],
       [[ 1. ,  1. ],
        [ 0. ,  1.1],
        [ 0. ,  0. ]]])

The tri.neighbors array contains information about which triangles are adjacent to each other.

tri.neighbors[1]
# array([-1,  0, -1], dtype=int32)

Recall that triangle 1 (tri.simplices in position 1) has vertices [3,1,0]. Triangle 0 is neighboring Triangle 1 opposite vertex 1, which is why the result has value 0 in the second element (corresponding to 1 in [3,1,0]). There are no triangles opposite vertex 3 (i.e., connected along the edge between vertex 0 and 1) or opposite vertex 0, so the neighbors array contains -1 in those positions.

Finally, is this code.

points[tri.simplices[1,1]]
array([ 0. ,  1.1])

Recalling the tri.simplices data structure above, there is a value 1 contained in position 1 in simplex 1 (i.e., [3,1,0]. This line is just looking up the coordinates of vertex 1.

As a final note, the order of vertices in the simplices returned is not required to match this original example and can vary from version to version. Here is a recent run matching an observation in the comment below which is inconsistent with the original vertex order (provided in the original documentation):

enter image description here