[SciPy-User] creating a 3D surface plot from collected data

[Joe Kington]

Hi there,

 I don't understand this, how can my list of the x-positions of the vertices
> be a
> 2D array?
>

mlab's mesh expects 2d grids of x,y, and z values so that it can infer how
the points are connected.

Because you just have lists of x,y,z verticies, there's no pre-defined way
to connect them.  Therefore, you need some way to define how they connect
into a surface.  (The program you exported them from probably stored them as
triangle-strip mesh.  You have the verticies, but don't have the definition
of how they link into individual triangles)

The surface reconstruction filter Robert mentioned is one way (and probably
the best way) of doing this.  (Credit where it's due... I had no idea vtk
had that particular filter before Robert mentioned it!)

It's slightly confusing, as it's not quite as simple as just replacing the
delaunay2d in Gael's earlier example with a call to the
SurfaceReconstructionFilter.

The surface reconstruction filter returns a volume that you need to extract
an isosurface at a value of 0 from.  (It took me awhile to figure that
out... Guess I should have read the vtk docs a bit closer!)

If this helps, here's a very artificial example using the surface
reconstruction filter.  It will produce some artifacts, but it should (?)
work for your data.... I think...  The first half just builds a set of
points on the outside of a sphere.  (capital) X,Y,Z would be representative
of your data.

from enthought.mayavi import mlab
import numpy as np

# First off, let's make 2d grids of phi and theta
phi = np.linspace(-np.pi, np.pi, 100)
theta = np.linspace(0, np.pi*2, 100)
phi, theta = np.meshgrid(phi, theta)

# Now lets make points on the outside of a sphere and
# put things back into cartesian coordinates
r = 1000
x = r * np.sin(phi) * np.cos(theta)
y = r * np.cos(phi)
z = r * np.sin(phi) * np.sin(theta)

# Since x,y,z are all regularly spaced 2d arrays, there's
# a defined way to link up each x,y,z coordinate into a
# surface (if that makes sense).  We can plot this using mesh
mlab.figure()
mlab.mesh(x,y,z)
mlab.title('The ideal surface')

# Now, lets "lose" the connectivity information by flattening the
# arrays.  This is more like what you have... A bunch of x,y,z
# coordinates of points on a surface with no clear definition
# of how to connect them.
X, Y, Z = [item.flatten() for item in [x,y,z]]

# Now we have to reconstruct how they connect into a surface
mlab.figure()
# Make a point source from the x,y,z coordinates
pts = mlab.pipeline.scalar_scatter(X,Y,Z,Z)

# This creates a volume of the distance to the reconstructed surface
vol = mlab.pipeline.user_defined(pts, filter='SurfaceReconstructionFilter')

# We can then extract an iso_surface at 0 that should be close to the
# original surface
mlab.pipeline.iso_surface(vol, contours=[0])
mlab.title('The reconstructed surface\nNotice the artifacts')
mlab.show()

Hope that helps,
-Joe

On Wed, Feb 17, 2010 at 8:54 AM, URI
<zallen at urologicresearchinstitute.org>wrote:

> Joe Kington <jkington <at> wisc.edu> writes:
>
> >
> >
> > No worries about the bottom post. I actually think I'm breaking protocol
> by
> top-posting, but I use gmail, and it's just more natural to
> top-post.Anyway, I
> think your problem is coming from the fact that you're working with an
> isosurface that fully encloses a volume. I assumed that you had scattered
> data
> points (say, elevation measurements) that you needed to interpolate
> between.
> The interpolation example I posted implicitly assumes that z = f(x,y). In
> other
> words, that there is only one z value for any given x and y.  If you're
> working
> with a 3D kidney-shape, this is clearly not the case.  Therefore, you're
> getting
> a singluar matrix when you try to interpolate (more than one value of z for
> identical rows of x and y in the matrix).The good news is that this may
> make the
> problem much simpler. (Your data is already in some sort of triangle-strip
> format if it's been exported from an imaging program. No need to
> interpolate.)What happens when you just do:fig = plt.figure()ax =
> Axes3D(fig)ax.plot_surface(x,y,z)plt.show()Where x,y,z are your raw data?
> There's a reasonable chance that the program you used to export the x,y,z
> isosurface left the verticies in the order that plot_surface expects them
> in...On a side note, you might also look into Mayavi's mlab module for 3D
> plotting.  I prefer it to matplotlib's native 3D plotting. If you have it
> installed, you'd just do:from enthought.mayavi import mlabs = mlab.mesh(x,
> y,
> z)mlab.show()Hope that helps some,
> > -Joe
> > On Tue, Feb 16, 2010 at 6:00 AM, URI <zallen <at>
> urologicresearchinstitute.org> wrote:
> > >Joe Kington <jkington <at> wisc.edu> writes:
> >
> > >
>
> I tried using mlab.mesh from mayavi but apparently the x y z components
> need to
> be 2D arrays (see
>
> http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/auto/mlab_helper_functions.html#enthought.mayavi.mlab.mesh
> ).
>
> I don't understand this, how can my list of the x-positions of the vertices
> be a
> 2D array?
>
> I did try to use mlab.plot3d (see
>
> http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/auto/mlab_helper_functions.html#enthought.mayavi.mlab.plot3d
> )
> and it gave me a nice line drawing of my shape, but an actual surface would
> be
> much nicer.
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100218/91e1ae2f/attachment.html>