Generating a 3D object (e.g. via Mayavi) and exporting it as 3D image stack (e.g. tiff)

Question

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)

Answer 1

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)