Exporting a 3D numpy to a VTK file for viewing in Paraview/Mayavi

Question

For those that want to export a simple 3D numpy array (along with axes) to a .vtk (or .vtr) file for post-processing and display in Paraview or Mayavi there's a little module called PyEVTK that does exactly that. The module supports structured and unstructured data etc.. Unfortunately, even though the code works fine in unix-based systems I couldn't make it work (keeps crashing) on any windows installation which simply makes things complicated. Ive contacted the developer but his suggestions did not work

Therefore my question is: How can one use the from vtk.util import numpy_support function to export a 3D array (the function itself doesn't support 3D arrays) to a .vtk file? Is there a simple way to do it without creating vtkDatasets etc etc?

Thanks a lot!

Wirer

Mayavi can display 3d numpy array and export it to couple of 3d formats. I'm not sure what are you after – Zellner

Paraview doesn't work with numpy arrays and especially if you want to throw a rectilinear 3D array with its axes at it. Such data need to be converter to a format like .vtr first. This conversion is what Im after. Any ideas? Also what did you mean about Mayavi supporting numpy arrays? The idea was to create a file from a numpy array created outside Mayavi and import it into such software as a file – Wirer

ParaView does work with numpy (read the manual about "Python Programmable Filter") and Mayavi can read numpy array of course, and export it to various 3d formats (which are BTW supported by ParaView) – Zellner

Answer 1

It's been forever and I had entirely forgotten asking this question but I ended up figuring it out. I've written a post about it in my blog (PyScience) providing a tutorial on how to convert between NumPy and VTK. Do take a look if interested:

pyscience.wordpress.com/2014/09/06/numpy-to-vtk-converting-your-numpy-arrays-to-vtk-arrays-and-files/

Wirer

I'd like to talk to you about this @somada141, I am trying this now in python 3.5 and am finding errors-- are you still working with these methodologies? – Gregggreggory

@Gregggreggory last I checked Python 3 support for many of those libraries was a tad iffy but feel free to reach out. As SO supports no PM mail me at <my username> at Gmail if you need help w sth. – Wirer

Answer 2

It's not a direct answer to your question, but if you have tvtk (if you have mayavi, you should have it), you can use it to write your data to vtk format. (See: http://code.enthought.com/projects/files/ETS3_API/enthought.tvtk.misc.html )

It doesn't use PyEVTK, and it supports a broad range of data sources (more than just structured and unstructured grids), so it will probably work where other things aren't.

As a quick example (Mayavi's mlab interface can make this much less verbose, especially if you're already using it.):

import numpy as np
from enthought.tvtk.api import tvtk, write_data

data = np.random.random((10,10,10))

grid = tvtk.ImageData(spacing=(10, 5, -10), origin=(100, 350, 200), 
                      dimensions=data.shape)
grid.point_data.scalars = np.ravel(order='F')
grid.point_data.scalars.name = 'Test Data'

# Writes legacy ".vtk" format if filename ends with "vtk", otherwise
# this will write data using the newer xml-based format.
write_data(grid, 'test.vtk')

And a portion of the output file:

# vtk DataFile Version 3.0
vtk output
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS 10 10 10
SPACING 10 5 -10
ORIGIN 100 350 200
POINT_DATA 1000
SCALARS Test%20Data double
LOOKUP_TABLE default
0.598189 0.228948 0.346975 0.948916 0.0109774 0.30281 0.643976 0.17398 0.374673 
0.295613 0.664072 0.307974 0.802966 0.836823 0.827732 0.895217 0.104437 0.292796 
0.604939 0.96141 0.0837524 0.498616 0.608173 0.446545 0.364019 0.222914 0.514992 
...
...

Reckford

Thanks a lot for the response Joe! That seems to be doing what I wanted more or less. I just need to use a format for non-uniform grids but Im sure i'll pull that one off :) – Wirer

line grid.point_data.scalars = np.ravel(order='F') should be grid.point_data.scalars = data.ravel(order='F') – Nonperishable

Answer 3

TVTK of Mayavi has a beautiful way of writing vtk files. Here is a test example I have written for myself following @Joe and tvtk documentation. The advantage it has over evtk, is the support for both ascii and html.Hope it will help other people.

from tvtk.api import tvtk, write_data
import numpy as np

#data = np.random.random((3, 3, 3))
#
#i = tvtk.ImageData(spacing=(1, 1, 1), origin=(0, 0, 0))
#i.point_data.scalars = data.ravel()
#i.point_data.scalars.name = 'scalars'
#i.dimensions = data.shape
#
#w = tvtk.XMLImageDataWriter(input=i, file_name='spoints3d.vti')
#w.write()

points = np.array([[0,0,0], [1,0,0], [1,1,0], [0,1,0]], 'f')
(n1, n2)  = points.shape
poly_edge = np.array([[0,1,2,3]])

print n1, n2
## Scalar Data
#temperature = np.array([10., 20., 30., 40.])
#pressure = np.random.rand(n1)
#
## Vector Data
#velocity = np.random.rand(n1,n2)
#force     = np.random.rand(n1,n2)
#
##Tensor Data with 
comp = 5
stress = np.random.rand(n1,comp)
#
#print stress.shape
## The TVTK dataset.
mesh = tvtk.PolyData(points=points, polys=poly_edge)
#
## Data 0 # scalar data
#mesh.point_data.scalars = temperature
#mesh.point_data.scalars.name = 'Temperature'
#
## Data 1 # additional scalar data
#mesh.point_data.add_array(pressure)
#mesh.point_data.get_array(1).name = 'Pressure'
#mesh.update()
#
## Data 2 # Vector data
#mesh.point_data.vectors = velocity
#mesh.point_data.vectors.name = 'Velocity'
#mesh.update()
#
## Data 3 additional vector data
#mesh.point_data.add_array( force)
#mesh.point_data.get_array(3).name = 'Force'
#mesh.update()

mesh.point_data.tensors = stress
mesh.point_data.tensors.name = 'Stress'

# Data 4 additional tensor Data
#mesh.point_data.add_array(stress)
#mesh.point_data.get_array(4).name = 'Stress'
#mesh.update()

write_data(mesh, 'polydata.vtk')

# XML format 
# Method 1
#write_data(mesh, 'polydata')

# Method 2
#w = tvtk.XMLPolyDataWriter(input=mesh, file_name='polydata.vtk')
#w.write()

Scleroprotein

Answer 4

I know it is a bit late and I do love your tutorials @somada141. This should work too.

def numpy2VTK(img, spacing=[1.0, 1.0, 1.0]):
 # evolved from code from Stou S.,
 # on http://www.siafoo.net/snippet/314
 # This function, as the name suggests, converts numpy array to VTK
 importer = vtk.vtkImageImport()

 img_data = img.astype('uint8')
 img_string = img_data.tostring()  # type short
 dim = img.shape

 importer.CopyImportVoidPointer(img_string, len(img_string))
 importer.SetDataScalarType(VTK_UNSIGNED_CHAR)
 importer.SetNumberOfScalarComponents(1)

 extent = importer.GetDataExtent()
 importer.SetDataExtent(extent[0], extent[0] + dim[2] - 1,
                       extent[2], extent[2] + dim[1] - 1,
                       extent[4], extent[4] + dim[0] - 1)
 importer.SetWholeExtent(extent[0], extent[0] + dim[2] - 1,
                        extent[2], extent[2] + dim[1] - 1,
                        extent[4], extent[4] + dim[0] - 1)

 importer.SetDataSpacing(spacing[0], spacing[1], spacing[2])
 importer.SetDataOrigin(0, 0, 0)


 return importer

Hope it helps!

Already

hey cheers Rick! I haven't quite seen this approach before. What VTK version were you using for this? – Wirer

hey @Wirer I am currently using VTK 7. Its Sidh though ;) – Already

@RickM. VTK_UNSIGNED_CHAR is not defined in your example. Any pointers? – Nonperishable

@Nonperishable You could use 5 instead of this or call, importer.SetDataScalarTypeToUnsignedChar – Already

Answer 5

Here's a SimpleITK version with the function load_itk taken from here:

import SimpleITK as sitk
import numpy as np


if len(sys.argv)<3:
    print('Wrong number of arguments.', file=sys.stderr)
    print('Usage: ' + __file__ + ' input_sitk_file' + ' output_sitk_file', file=sys.stderr)
    sys.exit(1)


def quick_read(filename):
    # Read image information without reading the bulk data.
    file_reader = sitk.ImageFileReader()
    file_reader.SetFileName(filename)
    file_reader.ReadImageInformation()
    print('image size: {0}\nimage spacing: {1}'.format(file_reader.GetSize(), file_reader.GetSpacing()))
    # Some files have a rich meta-data dictionary (e.g. DICOM)
    for key in file_reader.GetMetaDataKeys():
        print(key + ': ' + file_reader.GetMetaData(key))



def load_itk(filename):
    # Reads the image using SimpleITK
    itkimage = sitk.ReadImage(filename)
    # Convert the image to a  numpy array first and then shuffle the dimensions to get axis in the order z,y,x
    data = sitk.GetArrayFromImage(itkimage)
    # Read the origin of the ct_scan, will be used to convert the coordinates from world to voxel and vice versa.
    origin = np.array(list(reversed(itkimage.GetOrigin())))
    # Read the spacing along each dimension
    spacing = np.array(list(reversed(itkimage.GetSpacing())))
    return data, origin, spacing


def convert(data, output_filename):
    image = sitk.GetImageFromArray(data)
    writer = sitk.ImageFileWriter()
    writer.SetFileName(output_filename)
    writer.Execute(image)


def wait():
    print('Press Enter to load & convert or exit using Ctrl+C')
    input()


quick_read(sys.argv[1])
print('-'*20)
wait()
data, origin, spacing = load_itk(sys.argv[1])
convert(sys.argv[2])

Nonperishable