Return surface triangle of 3D scipy.spatial.Delaunay

Question

I have this problem. I try to triangulate points cloud by scipy.spatial.Delaunay. I used:

tri = Delaunay(points) # points: np.array() of 3d points 
indices = tri.simplices
vertices = points[indices]

But, this code return tetrahedron. How is it possible return triangle of surface only?

Thanks

Answer 1

To get it to work as in code form, you have to parametrize the surface to 2D. For example in the case of ball (r,theta, psi), radius is constant (drop it out) and points are given by (theta,psi) which is 2D.

Scipy Delaunay is N-dimensional triangulation, so if you give 3D points it returns 3D objects. Give it 2D points and it returns 2D objects.

Below is a script that I used to create polyhedra for openSCAD. U and V are my parametrization (x and y) and these are the coordinates that I give to Delaunay. Note that now the "Delaunay triangulation properties" apply only in u,v coordinates (angles are maximized in uv -space not xyz -space, etc).

The example is a modified copy from http://matplotlib.org/1.3.1/mpl_toolkits/mplot3d/tutorial.html which originally uses Triangulation function (maps to Delaunay eventually?)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri
from scipy.spatial import Delaunay

# u, v are parameterisation variables
u = np.array([0,0,0.5,1,1]) 
v = np.array([0,1,0.5,0,1]) 

x = u
y = v
z = np.array([0,0,1,0,0])

# Triangulate parameter space to determine the triangles
#tri = mtri.Triangulation(u, v)
tri = Delaunay(np.array([u,v]).T)

print 'polyhedron(faces = ['
#for vert in tri.triangles:
for vert in tri.simplices:
    print '[%d,%d,%d],' % (vert[0],vert[1],vert[2]),
print '], points = ['
for i in range(x.shape[0]):
    print '[%f,%f,%f],' % (x[i], y[i], z[i]),
print ']);'


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')

# The triangles in parameter space determine which x, y, z points are
# connected by an edge
#ax.plot_trisurf(x, y, z, triangles=tri.triangles, cmap=plt.cm.Spectral)
ax.plot_trisurf(x, y, z, triangles=tri.simplices, cmap=plt.cm.Spectral)


plt.show()

Below is the (slightly more structured) text output:

polyhedron(
    faces = [[2,1,0], [3,2,0], [4,2,3], [2,4,1], ], 

    points = [[0.000000,0.000000,0.000000], 
              [0.000000,1.000000,0.000000], 
              [0.500000,0.500000,1.000000], 
              [1.000000,0.000000,0.000000], 
              [1.000000,1.000000,0.000000], ]);