Chester VonWinchester Chester VonWinchester - 10 months ago 99
Python Question

Delaunay Triangulation of points from 2D surface in 3D with python?

I have a collection of 3D points. These points are sampled at constant levels (z=0,1,...,7). An image should make it clear:

point collection

These points are in a numpy ndarray of shape

(N, 3)
called
X
. The above plot is created using:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

X = load('points.npy')
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(X[:,0], X[:,1], X[:,2])
ax.scatter(X[:,0], X[:,1], X[:,2])
plt.draw()


I'd like to instead triangulate only the surface of this object, and plot the surface. I do not want the convex hull of this object, however, because this loses subtle shape information I'd like to be able to inspect.

I have tried
ax.plot_trisurf(X[:,0], X[:,1], X[:,2])
, but this results in the following mess:

mess

Any help?

Example data



Here's a snippet to generate 3D data that is representative of the problem:

import numpy as np
X = []
for i in range(8):
t = np.linspace(0,2*np.pi,np.random.randint(30,50))
for j in range(t.shape[0]):
# random circular objects...
X.append([
(-0.05*(i-3.5)**2+1)*np.cos(t[j])+0.1*np.random.rand()-0.05,
(-0.05*(i-3.5)**2+1)*np.sin(t[j])+0.1*np.random.rand()-0.05,
i
])
X = np.array(X)


Example data from original image



Here's a pastebin to the original data:

http://pastebin.com/YBZhJcsV

Here are the slices along constant z:

enter no image description here

Answer Source

Ok, here's the solution I came up with. It depends heavily on my data being roughly spherical and sampled at uniformly in z I think. Some of the other comments provide more information about more robust solutions. Since my data is roughly spherical I triangulate the azimuth and zenith angles from the spherical coordinate transform of my data points.

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

X = np.load('./mydatars.npy')
# My data points are strictly positive. This doesn't work if I don't center about the origin.
X -= X.mean(axis=0)

rad = np.linalg.norm(X, axis=1)
zen = np.arccos(X[:,-1] / rad)
azi = np.arctan2(X[:,1], X[:,0])

tris = mtri.Triangulation(zen, azi)

fig = plt.figure()
ax  = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(X[:,0], X[:,1], X[:,2], triangles=tris.triangles, cmap=plt.cm.bone)
plt.show()

Using the sample data from the pastebin above, this yields:

no thanks


update

I should update this for future people who come across it. While the above method works most of the time, it does assume (via spherical coordinate transformation) that no two points lie along the same ray. If you notice the artifact in the middle-left of the image above, this is why.

A better approach is to do "surgery" on the surface. Thinking about the surface like a orange peel, you cut it down one side, and then splay it open, stretching it. You then have a 2D plane which you can triangulate and interpolate. You simply have to keep track of how to get back to the corresponding place in 3D. To implement this idea takes quite a bit of work, and the implementation also requires some special care unique to how my data is represented.

Anyway, this is just to give an indication of how one might approach this more robustly.

update 2

I now do this as follows:

  1. I use the fact that the paths in each z-slice are closed and simple and use matplotlib.path to determine points inside and outside of the contour. Using this idea, I convert the contours in each slice to a boolean-valued volume.
  2. Next, I use skimage's marching_cubes method to obtain a triangulation of the surface for visualization.

Here's an example of the method. I think the data is slightly different, but you can definitely see that the results are much cleaner.