SciPy Delaunay triangulation changes multiple points of simplex for tiny change in parameters

Question

I'm using a Delaunay triangulation to interpolate in some function values evaluated at a set of parameters on a regular 4-dimensional grid. Sometimes, when a parameter value changes by a small amount that takes it to a new simplex, more than one point in the simplex changes. I expect that as I vary one of the parameter continuously, I'd move from simplex to simplex by changing just one point in the simplex at a time (and that's usually the case in my code, too). Instead, consider this script:

import numpy as np
from scipy.spatial import Delaunay

# hideous construction to get the desired 4d grid of points
# with points at [-1, -0.5, 0, 0.5, 1] along each axis
X = np.vstack(list(map(np.ravel, np.meshgrid(*[np.linspace(-1, 1, 5) for i in range(4)])))).T

tri = Delaunay(X)

delta = 1e-10

print(np.sort(tri.vertices[tri.find_simplex([-0.25, -0.05, 0.5+delta, 0.1])]))
print(np.sort(tri.vertices[tri.find_simplex([-0.25, -0.05, 0.5-delta, 0.1])]))

which produces

[192 292 317 318 322]
[167 292 293 313 317]

Note that these two simplices differ by 3 points, where I expect one, and I haven't devised a 2- or 3-D example where more than one vertex would change.

I'm 99% sure this is because my points are on a regular grid but I can't find a detailed answer of why or how to avoid the problem. I know that the triangulation isn't unique but that isn't fundamentally a problem. Various tricks appear to change where I encounter this issue but I haven't yet found a "fix" that prevents the issue from appearing anywhere.

Edit

I've managed to find an example in 3D, which allows me to visualise the problem.

import numpy as np
from scipy.spatial import Delaunay

X = np.vstack(list(map(np.ravel, np.meshgrid(*[np.linspace(-1, 1, 5) for i in range(3)])))).T
tri = Delaunay(X)

delta = 1e-6
x = np.array([-0.25, 0, 0.07])

fig = pl.figure()
ax = fig.add_subplot(projection='3d')

ax.scatter(*x, c='k')

x[1] = delta
s = X[tri.vertices[tri.find_simplex(x)]]

for i, si in enumerate(s):
    for j, sj in enumerate(s[i:]):
        ax.plot3D(*np.vstack([si, sj]).T, c='C0')

x[1] = -delta
s = X[tri.vertices[tri.find_simplex(x)]]
for i, si in enumerate(s):
    for j, sj in enumerate(s[i:]):
        ax.plot3D(*np.vstack([si, sj]).T, c='C1')

ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')

and here are two aspects on the output.

The blue and orange simplices are contain the interpolation point before and after it crosses y=0 (from positive to negative). My assumption was that both simplices would have the same triangular face along the y=0 plane but apparently this is incorrect. Is this fundamental to the Delaunay triangulation in this degenerate case or something about the implementation? And is there a way to avoid it? The QHull option Qx (which is in SciPy's default options for D>4) seems to help for this example but I'm not sure about globally.

Answer 1

Your question is not really about the implementation of triangulation in scipy.spatial. and it is more about the math of the Delaunay triangulation as a mathematical object.

Delaunay triangulations in dimension D are very well defined, ... when the points are "general position". That means that no D+2 points from the input points are on a common sphere. If that happens, one says the Delaunay triangulation is "degenerate". When the triangulation is degenerate, the Delaunay triangulation is not well defined, and there exists multiple ways to triangulate the convex hull of the points while preserving the Delaunay property.

What is what you observe: your points are on a regular grid, and that is a very degenerated point set (for the Delaunay property). Any slight modification of the coordinates can trigger the flip of multiple simplices, to restore the Delaunay property.

Maybe you can understand that behavior by having a look at the dual object of the Delaunay triangulation: its Voronoi diagram. For points sets close to a regular grid, the diagram is degenerated: it has Voronoi edges that are zero-length, or with a length close to zero. And any small modification of the coordinates of the point can change the topology of the Voronoi diagram (and thus of the Delaunay triangulation as well).