Main menu

Post navigation

NumPy to VTK: Converting your NumPy arrays to VTK arrays and files

In this post I will show how to ‘convert’ NumPy arrays to VTK arrays and files by means of the vtk.util.numpy_support module and the little-known PyEVTK package respectively.

Intro: The Conundrum

I wrote, or rather ranted, in my previous post about the value of VTK. Now lets say you were convinced (ha!) and decided to start including VTK in your scripts for visualization and processing.

Well, as if there weren’t enough deterrents in employing VTK, you will quickly realize that using your precious data – which let’s face it – will be stored in NumPy ndarray objects, with VTK ain’t all that straightforward. And why would it be? VTK was made in C++ and C++ isn’t about ease-of-use and concise programing. C++ is about putting hair on your chest :).

The traditional/ugly way, is creating new VTK objects, setting a bunch of properties like dimensions etc, and looping over your NumPy data to copy and populate your new objects. Since, looping in Python must be avoided like the black plague I will be focusing on the two ways I prefer.

The first way is using the vtk.util.numpy_support module that comes with VTK and allows you to ‘easily’ convert your data. The second way is by means of exporting your data into VTK-readable files using the PyEVTK package, a way which as you’ll see is great if you want to process and/or visualize that data in VTK-based applications.

Using the numpy_support module

So, given the popularity of Python and the fact that VTK is exposed in its near entirety to Python, the VTK folk decided to create the numpy_support module which resides under vtk.util. Of course, given the near-absence of documentation and/or examples, using it is as convoluted as doing anything in VTK. However, I’m here to try and elucidate their usage.

Usage

The functions of interest to us are numpy_to_vtk and vtk_to_numpy. Let us first inspect the docstring of the first function which can be accessed as follows, assuming you have VTK installed in your Python distro:

from vtk.util import numpy_support
help(numpy_support.numpy_to_vtk)

with the result being

numpy_to_vtk(num_array, deep=0, array_type=None)
Converts a contiguous real numpy Array to a VTK array object.
This function only works for real arrays that are contiguous.
Complex arrays are NOT handled. It also works for multi-component
arrays. However, only 1, and 2 dimensional arrays are supported.
This function is very efficient, so large arrays should not be a
problem.
If the second argument is set to 1, the array is deep-copied from
from numpy. This is not as efficient as the default behavior
(shallow copy) and uses more memory but detaches the two arrays
such that the numpy array can be released.
WARNING: You must maintain a reference to the passed numpy array, if
the numpy data is gc'd and VTK will point to garbage which will in
the best case give you a segfault.
Parameters
----------
- num_array : a contiguous 1D or 2D, real numpy array.

Upon first inspection, one might think that 3D NumPy arrays weren’t possible to convert. At least that’s what I thought (yeah, yeah, I suck). However, all that one needs to do is create a 1D representation of the array using ‘numpy.ndarray’ methods such as flatten or ravel. So here’s how to use this function assuming you mean to export an numpy.ndarray object named NumPy_data of type float32 :

As you can see we use ravel to flatten NumPy_data using the default C, or row-major, ordering (especially important on 2D and 3D arrays). In addition, we specify that we do want the data to be deep-copied by setting deep=True, while we also define the data type by array_type=vtk.VTK_FLOAT. Note, that we keep a copy of the shape of NumPy_data which we will use later to reshape the result of vtk_to_numpy. Converting back is much easier and can be done as such:

CAUTION: You may choose to allow for shallow-copies by setting deep=False but be warned: If for any reason, the array you pass is garbage-collected then the link will break and your nice VTK array will be useless. Should that occur, if you end up using vtk_to_numpy to back-convert, you will simply get memory trash. That is especially the case if you use flatten instead of ravel as the former always returns a ‘flattened’ copy of the array (check docs here) which is bound to get gc’ed if you use it directly in the call. ravel is slightly safer as it typically returns a new ‘view’ to the array and only copies when it has to (check docs here) but you don’t want to depend on that. Bottom line: disable deep-copying at your own risk and only if you know what you’re doing.

Using the PyEVTK package

Some time ago, I was struggling with the numpy_support module discussed above, mostly cause I sucked and didn’t think to ravel the array, so I started googling and came across PyEVTK a great little package by a Paulo Herrera. While back then, the code relied on C routines which refused to work on Windows platforms, since v0.9 – v1.0.0 being the current version – the package is built on pure-Python code and works like a charm.

What PyEVTK does, is allow you to save NumPy arrays straight to different types of VTK XML-based file formats (not the old legacy ones), without even needing to have VTK installed in your system, making NumPy the only dependency. This way, you can easily save your NumPy arrays into files that can be visualized and processed with any of the flagship VTK applications such as ParaView, VisIt, Mayavi, while obviously you can easily load said files with VTK and use all the goodies the toolkit offers.

However, while I’m supremely grateful to Mr. Herrera for creating PyEVTK, the package wasn’t hosted on PyPI and could only be used by checking out the code from the original repository in BitBucket, and using distutils to build/install it via the good ol’ fashioned python setup.py install.

Thus, possessed by the noble spirit of plagiarism, I took it upon myself to rip off, i.e., fork, the repository, re-package it, and upload it on PyPI. You can now find my fork on BitBucket here, while the package is hosted on PyPI under PyEVTKhere. Therefore, you can install with pip with the following command:

pip install pyevtk

As you’ll see on PyPI, I claim no ownership, authorship, or entitlement of PyEVTK. Paulo Herrera wrote the thing and the only reason I ‘stole’ it was cause I use it a lot and wanted to make it widely available.

Usage

I’m gonna give a quick example here to show you how to save a NumPy array as a rectilinear-grid file with a .vtr extension. In the interest of consistently misappropriating other people’s code, I modified the code for creating a Julia set from Ted Burke’s post here. I just thought I’d use a pretty dataset for my demo 🙂

The first important part of this code is from pyevtk.hl import gridToVTK which uses the high-level functions of the pyevtk package and which reside under the pyevtk.hl module. For the most-part, this module has all you need but the package also provides a pyevtk.vtk module for low-level calls. As we want to export data residing on a rectilinear grid we use the gridToVTK method.

The next 5 lines aren’t important but what we do is use the numpy.dstack function to stack 5 copies of the 2D julia array and create a 3D array to export. Note the numpy.arange calls through which we calculate ‘fake’ axes for the array to be exported. An important point to make is that since we’re defining a grid we need axes with N+1 coordinates which define the grid edges (as the data resides in the centers of that grid). If you want to export point-data, use the pyevtk.hl.pointsToVTK function.

As you guessed the magic happens with this one line:

gridToVTK("./julia", x, y, z, cellData = {'julia': juliaStacked})

As you can see, we pass the path to the file without an extension which will be automatically defined by PyEVTK depending on which function you used in order to create the appropriate VTK file. Subsequently, we pass the three coordinate axes, and provide a Python dictionary with a name for the array and the array itself. Doesn’t get easier than that :). The result of the above is .vtr file which you can find here while the IPython Notebook of the entire code can be found here.

Resources

A quick note here: while Paulo’s repo was named PyEVTK, he had originally named the package evtk but in the interest of consistency I renamed it to pyevtk so his examples won’t work directly if you installed PyEVTK through pip (see above section). If you want to see more examples of its usage, I refactored them a tad and can be seen in my repo here.

Well that’s it for today. Thanks for reading once more and I hope you gained something from this. Happy Python-ing 😀

Many thanks for that. I needed exactly what you want to do (tridimensional NumPy array -> VTK rectilinear grid, with the goal of doing volume rendering using ParaView, or other tools), and I managed to create the file and load it into ParaView, but the moment I select a real rendering method (other than Outline) paraview hangs with the following error printing many times :

The same problem happens with the julia.vtr file you provide. But not with other files coming from other sources. I don’t know if it is a bug in pyevtk, ParaView or a bit of both, but it looks like a problem

Hey Alex, that is indeed weird. What type of data are we talking about? Some of the volume-rendering algorithms in VTK only work with given data types. If you check my volume-rendering post (https://pyscience.wordpress.com/2014/11/16/volume-rendering-with-python-and-vtk/) you’ll see that classes like vtkVolumeRayCastMapper only work with unsigned char and unsigned short. You wouldn’t happen to have floats/signed ints or something would you 🙂

Hehe I was wondering when would someone call me on this one. Great question Darya and its indeed not an intuitive process. I’ve been meaning to write a full post of SimpleITKNumPyVTK conversions but I haven’t the time to touch my blog in months. I’m at work right now so I can’t give you the code but I’ll do it when I get off

You won’t find squat I’m afraid. The secret is in creating a new vtkImageData object, setting the size and spacing based on the NumPy array, and then using the setScalars method of the vtkImageData object to feed it the vtkDataArray you got from the initial conversion.

You can prob figure it out but I’ll write a snippet when I get the chance

At the link below you can find an IPython Notebook with the process to convert a dataset from SimpleITK to NumPy and then VTK (the process we described). Its the one I started preparing for that post I mentioned but its rather half-assed as I never got the time to complete it but it has the process you need.

For some reason, I can’t access the link you sent… For now I tried to write my code avoiding image conversions and I have my VTK pipeline. It isn’t working, because it runs out of memory (Windows 10 – 8 GB). I thought I was linking the filters properly by means of SetInputConnection(…GetOutputPort())…

So, I tried it on Linux and it is working (adding some corrections to ScalarTypes), but still nothing on Windows. Does it happen often?

I found that the Windows/OSX VTK is a little more temperamental than the Linux version where development and testing happens. If you’re running out of memory there’s not much you can do. You can verify that, however, by testing it on a smaller portion of your data.

I successfully converted my VTK image data no numpy array and back to VTK image data!
I only have one issue: the final image data is not of the same ScalarType as the original one (instead of short I get float). I guess the problem is with the functions SetScalarType() and SetNumberOfScalarComponents(). Both of them require 2 arguments: an ‘int’ and a ‘vtkInformation’ object. Probably I messed up with that, but it works 🙂

I’m glad it worked :). The ‘SetScalarType’ is the obvious suspect (the ‘SetNumberOfScalarComponents’ should only affect the ‘depth’ and not the type). If you can’t figure out what’s off just cast it down to an integer/short with the ‘vtkImageCast’ class (e.g., http://www.vtk.org/Wiki/VTK/Examples/Cxx/Images/Cast).

First of all amazing post ! I am also interested in what Darya wants to achieve . Specifically I want to feed a 3d numpy array (which is a grid of of isovalues basicly ) to vtkMarchingCubes , and I noticed that the input to the latter has to be of type vtkImageData . Since you already solved this I would be interested in looking at that script ! Have a good one!

Hello again , I was referring to Darya’s code . I already read the notebook you provided , and had a small issue but solved my problem already ! Thanks again!
Keep up the good work file mou!
Darya if you are reading this I would still be interested to read your code !

For comparison, I generated a mesh using the scikit-image marching cubes, and the marching cubes in PyMCubes. BTW, these two produce very similar results, which are worse than the mesh output from ITK-Snap, which I assume uses the VTK marching cubes algorithm.

My input volume is a pre-segmented and thresholded array, stored in a .nii (nifti) file format. It contains an artery, segmented from a contrast enhanced CT scan of a head. There is a large aneurysm visible in this vascular tree.

I load this volume, which has voxel values of either 0 or 1, then apply the various marching cubes algorithms to it. The results are saved as a .obj file (PyMCubes) or as a .stl file (VTK).

Both output files were loaded into MeshLab for visualization (and also visualized using mlab in the notebook). There is a screen capture showing the MeshLab renderings of these meshes – both the correct PyMCubes output (right panel) and the skewed VTK output (left panel).

I assume I am doing something wrong in in this block of code (from my notebook):

Solved, you didn’t set the Origin correctly mate. The original origin of the image is (0.0, 0.0, -93.79782104492188) and you set it to all 0s. Hence the squashing. I’ve updated your notebook and re-run the process to get an updated STL. You can find the thing under https://db.tt/WnBhtjQg. The new STL is potatoes.stl (don’t ask why). Lemme know if you’re having any more trouble 🙂

Hey, nice tutorial 🙂 But I have some problems. For my project i need the coordinate combinations for the ranges x:10-310, y:10-310, z: 0-65. I thought to use numpy.mgrid for this: coords = np.mgrid[xmin:xmax, ymin:ymax, zmin:zmax] and convert it to vtk_data_array like : vtk_data_array = numpy_support.numpy_to_vtk(num_array=coords.ravel(),deep=True,array_type=vtk.VTK_FLOAT). And after this to points: points = vtk.vtkPoints()
points.SetData(vtk_data_array). But this just crashes my python. Do you have an idea how to manage this?

Hi Varlor. Off the top of my head that sounds like you haven’t initialised the object properly. I can’t recall exactly but you might need to set the dimension, data-type, and allocate the point arrays prior to plonking an array in the `vtkPoints` object. Try looping over and inserting the points through the `InsertNextPoint` method. Also perhaps `vtkPoints` isn’t the right class for it and you might want to consider using the `vtkStructuredGrid` class instead.

Hello! I’m really glad you created such a tutorial! I’d be extremely grateful if you could help me with one issue that I cannot overcome. To cut a long story short, I’ve got a 3D numpy array filled after some segmentation procedures (e.g. value 255 is for the artery, whereas 0 is just the empty space) and I’d like to obtain a model of that vessel. Basing on your tutorial I managed to create a nice STL file, but I’d like to have a volumetric object (I guess it’s the .vtp or .vtk file) that can be used in VMTK program. I haven’t got the foggiest idea how to do that… Let’s assume STL creation firstly:

So that works perfectly and I can observe the surface of the artery. However, when I try to leave all ‘voxels’ as some objects and obtain a volumetric model composed of voxel-blocks (and save it as .vtk or .vtp file), I don’t have a clue what to do next…

Considering the pyEVTK module, I managed to export a .vtr file, but I didn’t manage to export .vtk file (don’t know how to do it – gridToVTK saves as .vtr, while pointsToVTK doesn’t work). Additionally, I guess this module saves all 0-elements also as some ‘objects’.

I would be extremely grateful if you could help me in creating .VTP or .VTK from 3D numpy array – it’s a huge must-be for my work. Thank you in advance!

Oh man I’m so sorry I was on vaca when it came through and must’ve missed it at the time. Have you managed to solve it?

Btw STLs are surface-deep but what you’re talking for is a volumetric mesh that can either consist of voxels or tetrahedra. I was under the impression that VMTK works with surface models (I mean why not, its just a vessel) but I’ll have to do some more looking into it.

Unfortunately I didn’t manage to solve that problem. Indeed, VMTK works with surface models (even after loading the vtp file – so PolyData, e.g. voxels – it applies marching cubes algorithm producing a surface). Hence, I guess I can operate on the STL files that I’m able to generate thanks to your tutorial. However, after applying all the code lines provided in VMTK tutorial (to tetrahedralize the mesh), I obtain only the surface mesh again, contrary to the tutorial guys. I haven’t got the foggiest idea whether VMTK somehow ‘remembers’ that it operated on the PolyData, and automatically seeks for the creation of the 3D volume mesh through Tetgen – initial file of every single tutorial is .vtp (volumetric) not .stl (surface). I’ve asked how to obtain a volume mesh from the pure STL file in the vmtk-users group, however, my question remains unanswered and unsolved. That’s why I wanted to test everything using .vtp files.
P.S. I think I’ll remind of myself in that group 🙂 Hopefully someone will provide me just a one line of code I was lacking.

Unfortunately my experience with VMTK is extremely old and limited but out of curiosity, is the vessel surface model closed or do you actually have open surfaces at the vessel interfaces? I’m not sure how VMTK meshes surfaces but I’d expect the code to operate on closed surfaces

In terms of loading the STL it has to be a waterproof surface (although I’m not 100% sure if it really has to be closed). Then you have to perform clipping to open inlet/outlets surfaces – there is a user interface and you perform clipping manually (as perdpendicularly to the channel as possible). The next algorithm that is triggered automatically closes all the previously cut surfaces and as a result there is a smooth plane for the inlet/outlets with displayed IDs. You can also approximate inlet/outlets with a circular shape and protrude it, so as to obtain pure circle cross-section for inlet/outlets. Afterwards you call a meshgenerator function that should create a volumetric mesh and then you can tetrahedralize it.

Hello!
I’ve got a problem with displaying the VTR files. I convert my DICOM (numpy arrays) into the VTR file and then I try to display it, I can see a rectangular cuboid of constant colour, instead of the most outer layers of the image set (e.g. part of the aorta, some bones, soft tissue, etc). Here are my scripts (maybe I’m doing something wrong):

To be honest, it was just a test of conversion function (from numpy to VTR) and I guess it failed. What’s the most probable reason behind that? My target objective is to create a VTK rendering window where I could see three 2D projections that cross one another (something like 3 separate matplotlib imshows, each one in each anatomical plane). I’d be extremely grateful if you could give me a little tip on how to do that since even after looking in the internet, I haven’t got the foggiest idea how to do that… I was thinking about generating VTI file instead of VTR, but it failed as well – I could have observed a red rectangular cuboid instead of normal Hounsfield-voxels…

Hey, a quick skim over the code didn’t show anything wrong with the conversion but I may have missed something. Did you try opening the VTR file in Paraview? If it doesn’t work there then at least you know the issue is indeed in the conversion

Thank you for your fast response. After reinstalling the ParaView, both VTR and VTI seem to be correct. When I select ‘aorta’ instead of the Solid Color, the displayed image was correct (otherwise I could have observed only a grey rectangular cuboid). However, when I load those files and try to display them in VTK rendering window (embedded in my application), I can observe only the grey surface of the rectangular cuboid. Do you know by any chance how to make vtk rendering window display true voxel intensities?