How to include all points into error-less triangulation mesh with scipy.spatial.Delaunay?

Question

I am testing scipy.spatial.Delaunay and not able to solve two issues:

1. the mesh has errors
2. the mesh doesn't include all points

Code and image of plot:

import numpy as np
from scipy.spatial import Delaunay,delaunay_plot_2d
import matplotlib.pyplot as plt    
#input_xyz.txt contains 1000 pts in "X Y Z" (float numbers) format
points = np.loadtxt("input_xyz.txt", delimiter=" ", usecols=(0, 1))
tri = Delaunay(points)
delaunay_plot_2d(tri)
plt.plot(points[:,0], points[:,1], 'o')
plt.show()

As mentioned under scipy.spatial.Delaunay:

"Unless you pass in the Qhull option “QJ”, Qhull does not guarantee that each input point appears as a vertex in the Delaunay triangulation."

But if I use QJ:

tri = Delaunay(points, qhull_options = "QJ")

I get Qhull error and if I use QJn (n=some high number):

tri = Delaunay(points, qhull_options = "QJ200")

to overcome that error the generated mesh looks terrible - triangles all over the place crossing each other.

How to include all points into error-less triangulation mesh with scipy.spatial.Delaunay?

Answer 1

The problem is that your data set is not centered. Qhull (used to do the Delaunay triangulation) does not center the data set for you under the default options, so it runs to rounding errors far away from origin.

You can center it yourself before triangulation

points -= points.mean(axis=0)
tri = Delaunay(points)