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?

_{(source: sethanil.com)}

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.

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:

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