How to display a volume with non-cubic voxels correctly in mayavi

Question

I'm using mayavi (3.3.2) to display volume isosurfaces.

Generally, my volumes do not have cubic voxels; for example, the sampling grid might be 1mm x 1mm in X and Y, but 1.4mm in the Z direction.

How can I get such volumes to display with the correct spatial proportions using mayavi's mlab.contour3d or mlab.pipeline.iso_surface ? I'd really prefer to not resample the volumes to a cubic grid.

Another way of stating the problem: what can I do to get the below code to display a sphere instead of a flattened elipsoid (taking the volume with it's intended 1:1:2 aspect-ratio voxels as a given, and without regenerating or resampling the volume).

import numpy as np
from enthought.mayavi import mlab

def sqr(x): return x*x

s=64
x,y,z = np.ogrid[0:s,0:s,0:s/2]

volume = np.sqrt(sqr(x-s/2)+sqr(y-s/2)+sqr(2*z-s/2))

isos = mlab.contour3d(volume,contours=[5,15,25],transparent=True)
mlab.show()

I'm guessing there ought to be some way of getting at the underlying VTK graphics pipeline (its transforms etc) and inserting the appropriate anisotropic scaling (if there isn't some way of doing it more directly through the mlab API).

Answer 1

For this, it's easiest to explicitly create a scalar_field object from the input data.

I actually do this quite frequently, as we like to put things in depth (where positive is downwards) in geology. That means that you need a negative increment in the z-direction. It would be nice if it was just an argument to the various mlab functions, but its still not too hard to do.

from mayavi import mlab
import numpy as np

s=64
x,y,z = np.ogrid[0:s,0:s,0:s/2]

data = np.sqrt((x-s/2)**2 + (y-s/2)**2 + (2*z-s/2)**2)

grid = mlab.pipeline.scalar_field(data)
grid.spacing = [1.0, 1.0, 2.0]

contours = mlab.pipeline.contour_surface(grid, 
                         contours=[5,15,25], transparent=True)
mlab.show()