Hi !
I aim trying to voxelize in python a 3D mesh of tetrahedron elements stored as a vtk unstructured grid (produced by tetgen) where I would like that each elements be graduatly incremented (numbered from 1 to n_cells +1)
I was therefore thinking that this might be possible with the vtkResampleToImage class but I struggle making it work as nothing seems to be written in the output file. I did the following :
from vtk import *
# The source file
file_name = "./patients/01/grid.vtk"
# Read the source file.
reader = vtkUnstructuredGridReader()
reader.SetFileName(file_name)
reader.Update()
output_vtk = reader.GetOutput()
# the number of tetrahedrons
n_tetrahedrons = output_vtk.GetNumberOfCells()
resample = vtkResampleToImage()
#resample.SetSamplingDimensions(700, 700, 700) ## i guess I should define the size of grid but that leads to a crash of kernel
resample.AddInputDataObject(output_vtk)
resample.Update()
writer = vtk.vtkMetaImageWriter()
writer.SetInputData(resample.GetOutput())
writer.SetFileName('grid_voxelized.mha')
writer.SetCompression(True)
writer.Write()
Hi @paul-bd, it’s not exactly clear to me what you are trying to achieve and why resampling to an image is needed - could you be a bit more specific?
You don’t use this in the code above… what would you like to do with it?
From what I’ve taken so far, you want to add an array on the cell data of the grid that has values from 1 to n_cells+1 (basically matching the cell IDs but off by one). You could do this with PyVista by:
import pyvista as pv
import numpy as np
output_vtk = pv.read("./patients/01/grid.vtk")
n_tetrahedrons = output_vtk.n_cells
output_vtk["my array"] = np.arange(1, n_tetrahedrons+1)
output_vtk.save("./patients/01/grid_with_array.vtk")
If you really want to resample to an image, you could do this with the vtkResampleWithDataSet algorithm and build the UniformGrid first. Note that it needs to have the same spatial extent as the input mesh. We can do this pretty easily in PyVista leverage some helper methods:
In pure-VTK this gets a bit verbose, but it look something like:
from vtk import *
import numpy as np
... # read the mesh
# Create a uniform grid of the same spatial extent
bounds = np.array(mesh.GetBounds())
dimensions = np.array([700,700,700], dtype=int)
image = vtkImageData()
image.SetDimensions(dimensions)
dims = (dimensions - 1)
imageSetSpacing((bounds[1::2] - bounds[:-1:2]) / dims)
image.SetOrigin(bounds[::2])
# Run the sampling
alg = vtkResampleWithDataSet()
alg.SetInputData(mesh) # Set the Input data (actually the source i.e. where to sample from)
alg.SetSourceData(image) # Set the Source data (actually the target, i.e. where to sample to)
alg.SetPassCellArrays(True)
alg.SetPassPointArrays(True)
alg.Update() # Perfrom the resampling
result = alg.GetOutput()
# ... Do something with the resampled image data
Thank you very much for your answer and your advices. It is exactly what I wanted to do.
I therefore added a scalar to the mesh to identify each cell as followed :
I want to do the resampling part so that I can extract the numpy array of shape = dimensions and were the value at x,y, z is the integer value of mesh[“my_array”] corresponding to x,y,z (so the tetrahedron number + 1).
Is there a method in pyvista or vtk to do so or shall I loop within image.points and image[“my_array”] to fill that numpy array.