3D Contour plot from data using Mayavi / Python

Question

I would like to do a 3D contour plot using Mayavi in exactly the same way as the third figure on this page (a hydrogen electron cloud model) :

http://www.sethanil.com/python-for-reseach/5

I have a set of data points which I created using my own model which I would like to use. The data points are stored in a multi-dimensional numpy array like so:

XYZV = [[1, 2, 3, 4],
 [6, 7, 8, 9],
 ...
 [4, 5, 6, 7]]

The data points are not uniformly spread in XYZ space and not stored in any particular order. I think the example uses a meshgrid to generate the data points - I have looked this up but totally don't understand it. Any help would be much appreciated?

H
_{(source: sethanil.com)}

• Don Question: Show us, what you tried so far. We'll help from there onwards.
• David Z: Just for future reference, questions like this are a great fit for Computational Science.

Answer 1

The trick is to interpolate over a grid before you plot - I'd use scipy for this. Below R is a (500,3) array of XYZ values and V is the "magnitude" at each XYZ point.

from scipy.interpolate import griddata
import numpy as np
# Create some test data, 3D gaussian, 200 points
dx, pts = 2, 100j
N = 500
R = np.random.random((N,3))*2*dx - dx
V = np.exp(-( (R**2).sum(axis=1)) )
# Create the grid to interpolate on
X,Y,Z = np.mgrid[-dx:dx:pts, -dx:dx:pts, -dx:dx:pts]
# Interpolate the data
F = griddata(R, V, (X,Y,Z))

From here it's a snap to display our data:

from mayavi.mlab import *
contour3d(F,contours=8,opacity=.2 )

This gives a nice (lumpy) Gaussian.

enter image description here

Take a look at the docs for griddata, note that you can change the interpolation method. If you have more points (both on the interpolated grid, and on the data set), the interpolation gets better and better represents the underlying function you're trying to illustrate. Here is the above example at 10K points and a finer grid:

enter image description here

• joshlk: Thanks so much. It works like a charm! Just one question: if I wanted to double the amount of points on the fitting grid what would I change in the line 'dx, pts = 2, 100j'?
• Hooked: @Josh you would change it to dx, pts = 2, 200j, however this would double the number of points in each dimension so you would have 2^3=8 times as many points to interpolate over. dx controls the extent of the plot grid. For finer control just mgrid for each linear dimension.
• joshlk: I find the griddata function can take a very long time to compute. Do you know any tips for speeding the process up?
• Hooked: @Josh you could try changing the interpolation to nearest rather than linear, I haven't tried a timing test so I'm not sure how much of a difference this will make. In addition, you can split the plot into sections - high density regions can have more points/voxel while low density regions could be more rarefied. You'll have to have separate griddata commands for that though.

Answer 2

You can use delaunay3d filter to create cells from points. Then you can create an iso_surface() for the output UnstructuredGrid of delaunay3d. If you want ImageData, you can use image_data_probe filter.

import numpy as np
from tvtk.api import tvtk
from mayavi import mlab
points = np.random.normal(0, 1, (1000, 3))
ug = tvtk.UnstructuredGrid(points=points)
ug.point_data.scalars = np.sqrt(np.sum(points**2, axis=1))
ug.point_data.scalars.name = "value"
ds = mlab.pipeline.add_dataset(ug)
delaunay = mlab.pipeline.delaunay3d(ds)
iso = mlab.pipeline.iso_surface(delaunay)
iso.actor.property.opacity = 0.1
iso.contour.number_of_contours = 10
mlab.show()

enter image description here

• joshlk: Hi when I use this method with my own data I keep getting loads of errors popping up in a new windows with the last saying: "Warning: In C:\pisi\tmp\VTK-5.6.0-2\work\VTK\Graphics\vtkDelaunay3D.cxx, line 488 vtkDelaunay3D (03AD2250): 27 degenerate triangles encountered, mesh quality suspect". Do you know what's happening?