0

My current task is to generate a 3D image space where there are 3D objects (of iso-surfaces) that I designed and export it as an image stack (numpy or tiff).

I came down to using Mayavi to generate 3D iso-surfaces. I know Mayavi is originally designed to provide 3D visualizations on its own, but I would like to find a way that I can export a 3D object to a 3D image stack as in a numpy form of (z,y,x). My initial idea was to iteratively take snapshots of the sliced volume from the Mayavi mlab object along z-axis but I am not sure if there's any option to save a sliced image of an iso-surface as a snapshot.

The best case scenario would be to export a 3D image stack (tiff) exactly from what I see from a Mayavi window. Otherwise, I will take any suggestions to carry out this task in general.

Here's an example code.

import numpy as np
from mayavi import mlab

# Produce some nice data.
n_mer, n_long = 6, 11
pi = np.pi
dphi = pi/1000.0
phi = np.arange(0.0, 2*pi + 0.5*dphi, dphi, 'd')
mu = phi*n_mer
x = np.cos(mu)*(1+np.cos(n_long*mu/n_mer)*0.5)
y = np.sin(mu)*(1+np.cos(n_long*mu/n_mer)*0.5)
z = np.sin(n_long*mu/n_mer)*0.5

# Init plot
source = mlab.points3d(x, y, z)

Example of a 3D object generated by Mayavi for visualization

peterpark828
  • 177
  • 1
  • 1
  • 9

1 Answers1

2

You might go for the vtk class vtkImplicitModeller. E.g.:

import numpy as np
from vedo import Points, Volume

n_mer, n_long = 6, 11
dphi = np.pi/1000.0
phi = np.arange(0.0, 2*np.pi + 0.5*dphi, dphi, 'd')
mu = phi*n_mer
x = np.cos(mu)*(1+np.cos(n_long*mu/n_mer)*0.5)
y = np.sin(mu)*(1+np.cos(n_long*mu/n_mer)*0.5)
z = np.sin(n_long*mu/n_mer)*0.5

source = Points([x, y, z], r=4)
modl = source.implicitModeller(
         distance=0.15,
         res=(60,60,30),
         bounds=(-1.8,1.8,-1.8,1.8,-0.7,0.7),
)
modl.smoothLaplacian().computeNormals()
modl.c("blue9").lw(1).lighting("metallic").show(axes=1)

#######################################################
import vtk
imp = vtk.vtkImplicitModeller()
imp.SetInputData(source.polydata())
imp.SetSampleDimensions(50,50,30)
imp.SetModelBounds(-1.8,1.8,-1.8,1.8,-0.7,0.7)
imp.Update()

vol = Volume(imp.GetOutput())
arr = np.clip(vol.getDataArray(), 0, 1.2)
print(arr.shape)

enter image description here

mmusy
  • 1,292
  • 9
  • 10
  • it works for me! I just have one more question. It seems the clipping at the end is necessary to preserve the shading effect (otherwise it's all filled to the brim). Is the clipping threshold set empirically or is there a way to calculate an effective threshold for clipping? – peterpark828 Jul 02 '21 at 07:02
  • the `np.clip` is due to the fact that `vtkImplicitModeller` generates very large values (not sure why this happens, it should not..). The bounds are set a bit empirically, they should roughly be augmented by the desired distance value when isosurfacing the volume. – mmusy Jul 02 '21 at 11:07