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 :

NumPy_data_shape = NumPy_data.shape
VTK_data = numpy_support.numpy_to_vtk(num_array=NumPy_data.ravel(), deep=True, array_type=vtk.VTK_FLOAT)

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:

NumPy_data = numpy_support.vtk_to_numpy(VTK_data)
NumPy_data = NumPy_data.reshape(NumPy_data_shape)

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 PyEVTK here. 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 🙂

Here, however, I’m only going to show the PyEVTK part of the code while you can see the full thing under this IPython Notebook. The result of the Julia code is the following pretty 2D array which was visualized with Plotly’s Heatmap method (see this previous post on how to use Plotly):

A Julia set visualized with Plotly

So, here’s the PyEVTK part of the notebook:

from pyevtk.hl import gridToVTK

noSlices = 5
juliaStacked = numpy.dstack([julia]*noSlices)

x = numpy.arange(0, w+1)
y = numpy.arange(0, h+1)
z = numpy.arange(0, noSlices+1)

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

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 😀

Advertisements

28 thoughts on “NumPy to VTK: Converting your NumPy arrays to VTK arrays and files

  1. Pingback: Surface Extraction: Creating a mesh from pixel-data using Python and VTK | PyScience

  2. Pingback: Image Segmentation with Python and SimpleITK | PyScience

  3. Pingback: Multi-Modal Image Segmentation with Python & SimpleITK | PyScience

  4. Hi,

    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 :

    ERROR: In /build/paraview-B2WpY_/paraview-4.1.0+dfsg+1/VTK/Rendering/OpenGL/vtkOpenGLPainterDeviceAdapter.cxx, line 340
    vtkOpenGLPainterDeviceAdapter (0x1d98360): Unsupported type for vertices: 8

    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

    Like

  5. Hi!

    Thank you so much for this post. I only have one question: the ‘numpy_support.numpy_to_vtk’ function returns a vtkDataArray, but I need to convert it to vtkImageData. How can I do this?

    Best regards,
    Darya

    Like

    • 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

      Like

      • That would be an amazing post, since I couldn’t find much information on this topic 🙂 I guess I’ll play around with vtkImageImport and vtkImageData. Thank you again and good luck at your work!

        Like

      • 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

        Like

      • Btw be extra careful with the setScalarType or whatever it’s called. You need to use the proper data type or VTK will crash wo explanation due to memory overflows

        Like

      • Hey Darya,

        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.

        https://dl.dropboxusercontent.com/u/3091507/ConvertSimpleITKtoVTK.ipynb

        Take a looksie and lemme know if you run into trouble 🙂

        Cheers,
        Adam

        Like

      • Hi Adam!

        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?

        Thank you!

        Best regards,
        Darya

        Like

      • Ok that’s weird. Maybe its cause Dropbox gets mixed up with the format. I’ve zipped and uploaded the file on both my DropBox and GoogleDrive. One of the following links should work:

        https://dl.dropboxusercontent.com/u/3091507/ConvertSimpleITKtoVTK.ipynb.zip

        https://drive.google.com/file/d/0B_IXVGLixuWMVkw4c1JsOGhqWnc/view?usp=sharing

        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.

        Anyway, I hope that one of the above links will do the trick.

        Lemme know how you go

        Cheers,
        Adam

        Like

      • Thank you so much, your script was very helpful!

        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 🙂

        Like

    • Hello!

      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!

      Like

      • 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 !

        Like

      • Hi Michael!

        Sorry for not answering in time to solve your problem! I’m glad that everything worked out. In my code I use something like this (just as Adamos explained):


        CT_Data = reader.GetOutput()

        … convert to numpy array and do some stuff …

        NP_data = numpy_support.numpy_to_vtk(np_image_data.ravel(), deep=True, array_type=vtk.VTK_SHORT)

        imageVTK = vtk.vtkImageData()
        imageVTK.SetSpacing(CT_Data.GetSpacing())
        imageVTK.SetOrigin(CT_Data.GetOrigin())
        imageVTK.SetDimensions(CT_Data.GetDimensions())
        imageVTK.GetPointData().SetScalars(NP_data)

        And then I just feed the imageVTK to the rest of the pipeline (by the way, it is also vtkDiscreteMarchingCubes). I hope this helps!

        Thank you all!

        Liked by 1 person

  6. Hi:

    Thanks for the great post and follow up comments. I am having an issue with the mesh produced by vtk – it appears ‘skewed’.

    I have uploaded a zip that contains an ipython notebook, input file, and output files:

    https://www.dropbox.com/s/0jy3e5idk61tr49/vtk_mesh.zip?dl=0

    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):

    —-
    # voxel_array is a numpy array obtained from the .nii file. The array has type np.int16
    # and dimensions of (256, 256, 255)
    #
    vtk_data = numpy_support.numpy_to_vtk(num_array=voxel_array.ravel(), deep=True, array_type=vtk.VTK_SHORT)

    # nii_spacing = (0.33868799, 0.33868799, 0.33868408)
    # nii_shape = (256, 256, 255)

    imageVTK = vtk.vtkImageData()
    imageVTK.SetSpacing(nii_spacing)
    imageVTK.SetOrigin([0.0, 0.0, 0.0])
    imageVTK.SetDimensions(nii_shape)
    imageVTK.GetPointData().SetScalars(vtk_data)

    dmc = vtk.vtkDiscreteMarchingCubes()
    dmc.SetInputData(imageVTK)
    dmc.GenerateValues(1, 1, 1)
    dmc.Update()
    —-

    Any insight appreciated!

    Sincerely,

    Ross Mitchell

    Like

  7. Really cool tutorial. Didn’t manage to get it working for .vtp (vector_scatter, or quiver – just plots one component as the arrow length…), but very useful otherwise!

    Thanks for your time.

    Like

  8. 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?

    Like

    • 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.

      Like

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s