8

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

timday
  • 24,582
  • 12
  • 83
  • 135

1 Answers1

7

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

Spherical shell with non-cubic voxels

Joe Kington
  • 275,208
  • 71
  • 604
  • 463
  • Thanks that's exactly what I was looking for. I really need to put some effort into learning my way around the mlab API. – timday May 24 '11 at 16:42
  • 1
    @timday - This is a good place to start learning the "pipeline" side of things: http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/mlab_pipeline.html#controlling-the-pipeline-with-mlab-scripts Also, if you're familiar with VTK, and there's a filter that isn't wrapped directly by `mlab`, you can call it with `mlab.pipeline.user_defined`, which can be quite handy if you can find what you need in VTK but not directly in mayavi/mlab. – Joe Kington May 24 '11 at 16:53
  • 1
    Also, have at look at the different VTK data structures: http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/data.html#data-structures-used-by-mayavi – Joe Kington May 24 '11 at 16:54
  • I'm guessing that `data` is supposed to `== volume`? – foobarbecue Jul 23 '15 at 03:16
  • @foobarbecue - Yep, absolutely! Thanks for the catch! – Joe Kington Jul 23 '15 at 13:37