As I had some problems in figuring out how to use VTK to render data contained in a three dimensional numpy array, I have decided to share my code. This code is based on the otherwise excellent documentation for VTK and the now outdated vtkImageImportFromArray-class created by David Gobbi found at http://public.kitware.com/cgi-bin/cvsweb.cgi/vtk/python/?cvsroot=vtk

The example is very simple, for more advanced functionality: read the documentation.

   1 import vtk
   2 from numpy import *
   3 
   4 # We begin by creating the data we want to render.
   5 # For this tutorial, we create a 3D-image containing three overlaping cubes. 
   6 # This data can of course easily be replaced by data from a medical CT-scan or anything else three dimensional.
   7 # The only limit is that the data must be reduced to unsigned 8 bit or 16 bit integers.
   8 data_matrix = zeros([75, 75, 75], dtype=uint8)
   9 data_matrix[0:35, 0:35, 0:35] = 50
  10 data_matrix[25:55, 25:55, 25:55] = 100
  11 data_matrix[45:74, 45:74, 45:74] = 150
  12 
  13 # For VTK to be able to use the data, it must be stored as a VTK-image. This can be done by the vtkImageImport-class which
  14 # imports raw data and stores it. 
  15 dataImporter = vtk.vtkImageImport()
  16 # The preaviusly created array is converted to a string of chars and imported.
  17 data_string = data_matrix.tostring()
  18 dataImporter.CopyImportVoidPointer(data_string, len(data_string))
  19 # The type of the newly imported data is set to unsigned char (uint8)
  20 dataImporter.SetDataScalarTypeToUnsignedChar()
  21 # Because the data that is imported only contains an intensity value (it isnt RGB-coded or someting similar), the importer
  22 # must be told this is the case.
  23 dataImporter.SetNumberOfScalarComponents(1)
  24 # The following two functions describe how the data is stored and the dimensions of the array it is stored in. For this
  25 # simple case, all axes are of length 75 and begins with the first element. For other data, this is probably not the case.
  26 # I have to admit however, that I honestly dont know the difference between SetDataExtent() and SetWholeExtent() although
  27 # VTK complains if not both are used.
  28 dataImporter.SetDataExtent(0, 74, 0, 74, 0, 74)
  29 dataImporter.SetWholeExtent(0, 74, 0, 74, 0, 74)
  30 
  31 # The following class is used to store transparencyv-values for later retrival. In our case, we want the value 0 to be
  32 # completly opaque whereas the three different cubes are given different transperancy-values to show how it works.
  33 alphaChannelFunc = vtk.vtkPiecewiseFunction()
  34 alphaChannelFunc.AddPoint(0, 0.0)
  35 alphaChannelFunc.AddPoint(50, 0.05)
  36 alphaChannelFunc.AddPoint(100, 0.1)
  37 alphaChannelFunc.AddPoint(150, 0.2)
  38 
  39 # This class stores color data and can create color tables from a few color points. For this demo, we want the three cubes
  40 # to be of the colors red green and blue.
  41 colorFunc = vtk.vtkColorTransferFunction()
  42 colorFunc.AddRGBPoint(50, 1.0, 0.0, 0.0)
  43 colorFunc.AddRGBPoint(100, 0.0, 1.0, 0.0)
  44 colorFunc.AddRGBPoint(150, 0.0, 0.0, 1.0)
  45 
  46 # The preavius two classes stored properties. Because we want to apply these properties to the volume we want to render,
  47 # we have to store them in a class that stores volume prpoperties.
  48 volumeProperty = vtk.vtkVolumeProperty()
  49 volumeProperty.SetColor(colorFunc)
  50 volumeProperty.SetScalarOpacity(alphaChannelFunc)
  51 
  52 # This class describes how the volume is rendered (through ray tracing).
  53 compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
  54 # We can finally create our volume. We also have to specify the data for it, as well as how the data will be rendered.
  55 volumeMapper = vtk.vtkVolumeRayCastMapper()
  56 volumeMapper.SetVolumeRayCastFunction(compositeFunction)
  57 volumeMapper.SetInputConnection(dataImporter.GetOutputPort())
  58 
  59 # The class vtkVolume is used to pair the preaviusly declared volume as well as the properties to be used when rendering that volume.
  60 volume = vtk.vtkVolume()
  61 volume.SetMapper(volumeMapper)
  62 volume.SetProperty(volumeProperty)
  63 
  64 # With almost everything else ready, its time to initialize the renderer and window, as well as creating a method for exiting the application
  65 renderer = vtk.vtkRenderer()
  66 renderWin = vtk.vtkRenderWindow()
  67 renderWin.AddRenderer(renderer)
  68 renderInteractor = vtk.vtkRenderWindowInteractor()
  69 renderInteractor.SetRenderWindow(renderWin)
  70 
  71 # We add the volume to the renderer ...
  72 renderer.AddVolume(volume)
  73 # ... set background color to white ...
  74 renderer.SetBackground(1, 1, 1)
  75 # ... and set window size.
  76 renderWin.SetSize(400, 400)
  77 
  78 # A simple function to be called when the user decides to quit the application.
  79 def exitCheck(obj, event):
  80     if obj.GetEventPending() != 0:
  81         obj.SetAbortRender(1)
  82 
  83 # Tell the application to use the function as an exit check.
  84 renderWin.AddObserver("AbortCheckEvent", exitCheck)
  85 
  86 renderInteractor.Initialize()
  87 # Because nothing will be rendered without any input, we order the first render manually before control is handed over to the main-loop.
  88 renderWin.Render()
  89 renderInteractor.Start()

To exit the application, simply press q.

In my opinion, the volume renderer creates extremely ugly images if not the following option is used:

   1 volumeProperty.ShadeOn()


CategoryCookbook

Cookbook/vtkVolumeRendering (last edited 2010-09-08 08:25:22 by PauliVirtanen)