CODEWITHSUNDEEP
pythonscipycloudpointsdelaunay
s.t.e.a.l.t.h
s.t.e.a.l.t.h

Reputation: 520

Return surface triangle of 3D scipy.spatial.Delaunay

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

Upvotes: 13

Views: 31454

Answers (3)

gmagno
gmagno

Reputation: 1990

Following Jaime's answer, but elaborating a bit more with an example:

import matplotlib as mpl
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3
import numpy as np
import scipy as sp
from scipy import spatial as sp_spatial


def icosahedron():
    h = 0.5*(1+np.sqrt(5))
    p1 = np.array([[0, 1, h], [0, 1, -h], [0, -1, h], [0, -1, -h]])
    p2 = p1[:, [1, 2, 0]]
    p3 = p1[:, [2, 0, 1]]
    return np.vstack((p1, p2, p3))


def cube():
    points = np.array([
        [0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1],
        [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1],
    ])
    return points


points = icosahedron()
# points = cube()

hull = sp_spatial.ConvexHull(points)
indices = hull.simplices
faces = points[indices]

print('area: ', hull.area)
print('volume: ', hull.volume)


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.dist = 30
ax.azim = -140
ax.set_xlim([0, 2])
ax.set_ylim([0, 2])
ax.set_zlim([0, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

for f in faces:
    face = a3.art3d.Poly3DCollection([f])
    face.set_color(mpl.colors.rgb2hex(sp.rand(3)))
    face.set_edgecolor('k')
    face.set_alpha(0.5)
    ax.add_collection3d(face)

plt.show()

Which should depict the following figure:

enter image description here

Upvotes: 3

Juha
Juha

Reputation: 2115

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()

Pyramid what the above example plots

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], ]);

Upvotes: 11

Jaime
Jaime

Reputation: 67427

It looks like you want to compute the convex hull of your point cloud. I think this is what you want to do:

from scipy.spatial import ConvexHull

hull = ConvexHull(points)
indices = hull.simplices
vertices = points[indices]

Upvotes: 4

Related Questions