This short tutorial will get you started using the Insight Toolkit (ITK) with Python in Ubuntu. We will explore basic file reading, image processing, and volume rendering. It is an introduction, the ITK Software Guide and Python Tutorial contain more in depth information.
Basics
This tutorial requires installation of the packages as described in my previous tutorial. Also, a few other packages are needed.
sudo apt-get install python-scipy insighttoolkit-examples
Insighttoolkit-examples contains images we are going to use in this demo. Uncompress the image first.
sudo gunzip /usr/share/doc/insighttoolkit-examples\
/examples/Data/BrainProtonDensity3Slices.raw.gz
I prefer to use IPython instead of python for the command line. IPython gives nice features like auto-completion, text coloring, and debugging. It’s great.
sudo apt-get install ipython
ipython
Alright, we are ready to start using python. All the following code samples are typed into the IPython command line. Begin by importing the ITK module.
import itk
ITK is a heavily templated system. As a result, everything requires an image type. For example, a three dimensional image of type unsigned char is defined by
image_type = itk.Image[itk.UC, 3]
A large combination of image types are available, in fact ITK puts no restriction of data types and image dimensions. Unfortunately, this template system doesn’t translate into python. I had to compromise and include a sub-set of dimensions and data-types. If you need more functionality, let me know and I can add it to subsequent versions of python-insighttoolkit. The following data-types for 2, 3, and 4 dimensional images are included
- unsigned char
- signed short
- unsigned short
- RGB unsigned short
- float
- complex float
- vector float
- covariant vector float
Reading and Writing Files
How do we read files? Just enter
file_name = '/usr/share/doc/insighttoolkit-examples/\
examples/Data/BrainProtonDensity3Slices.mha'
reader = itk.ImageFileReader[image_type].New()
reader.SetFileName( file_name )
reader.Update()
That’s it. See how the reader required an image_type. This is slightly restrictive because you need to know the dimensionality and data type of the image before you read it. But not so bad.
These four short lines of code will read TIFF, JPEG, PNG, BMP, DICOM, GIPL, Bio-Rad, LSM, Nifti, Analyze, SDT/SPR (Stimulate), Nearly Raw Raster Data (Nrrd), and VTK images. I don’t know what half of those file formats are, but someone will find them useful.
Simple Image Processing Example
Alright, an image is loaded into memory and ready to go. We’ll start with a simple example, a median filter.
median_filter = itk.MedianImageFilter[image_type, image_type].New()
Notice that image_type is used not once, but twice. The first and second image type refer to the input and output type, respectively. One image type definition was necessary for the reader above, because there is only an output image type. Filters require an input and output type.
Filters often expose parameters to adjust their function. For example, setting the radius to 1 specifies the median filter will operate on 3x3x3 neighborhoods.
median_filter.SetRadius( 1 )
Set the filter input to the output of the reader and update.
median_filter.SetInput( reader.GetOutput() )
median_filter.Update()
Remember earlier when we called reader.Update(). Well, this isn’t necessary. When median_filter.Update() is called, ITK will look at its input, and if necessary update it. If you put together a long string of filters, an update on the final filter will cascade all the way back to the beginning and update everything. Of course ITK is even smarter than this and will only rerun filters if its parameters or input has changed.
There are many more operations you can perform on your images, but we need to move on.
ITK and VTK
ITK is great for reading, writing, and processing images. However, ITK is not equipped for visualization. This is where the Visualization Toolkit (VTK) comes in. Fortunately, ITK and VTK were created by some of the same people, so the concepts are very similar. In addition, connecting ITK and VTK is trivial with the python-insighttoolkit-extras package.
itk_vtk_converter = itk.ImageToVTKImageFilter[image_type].New()
itk_vtk_converter.SetInput( median_filter.GetOutput() )
itk_vtk_converter.Update()
Now itk_vtk_converter.GetOutput() returns VTK data. Very handy, especially if you want to render the volume.
Volume Rendering with VTK
This section is a quick introduction to volume rendering with VTK. It is low on details because there are better introductions to VTK and python already available.
First we need to import VTK.
import vtk
A volume mapper determines how image data is rendered to the screen. In this case we use a high quality ray casting mapper. Depending on your video card, you can use the higher performance vtk.vtkVolumeTextureMapper3D() or vtk.vtkVolumeTextureMapper2D() mappers.
volume_mapper = vtk.vtkVolumeRayCastMapper()
volume_mapper.SetInput( itk_vtk_converter.GetOutput() )
The Ray Cast Mapper requires a composite function.
composite_function = vtk.vtkVolumeRayCastCompositeFunction()
volume_mapper.SetVolumeRayCastFunction( composite_function )
Our input image is 8 bit data with values from 0 to 255. A mapping is required between these values and the displayed color. For instance, the following code will map the image intensity to the blue channel. This is a linear function; 0 will map to black (0.0, 0.0, 0.0) and 255 will map to solid blue (0.0, 0.0, 1.0). Everything in between will be a linear interpolation between the two end points, i.e. 127 will map to dark blue (0.0, 0.0, 0.5).
color_transfer_func = vtk.vtkColorTransferFunction()
color_transfer_func.AddRGBPoint( 0, 0.0, 0.0, 0.0 )
color_transfer_func.AddRGBPoint( 255, 0.0, 0.0, 1.0 )
Each voxel in the input image must also map to opacity. The following maps image intensities of 0 to completely transparent (0.0), and 255 to completely opaque (1.0). Again, values between 0 and 1 are interpolated.
opacity_transfer_func = vtk.vtkPiecewiseFunction()
opacity_transfer_func.AddPoint( 0, 0.0 )
opacity_transfer_func.AddPoint( 255, 1.0 )
Encapsulate the above properties in a vtkVolumeProperty.
volume_properties = vtk.vtkVolumeProperty()
volume_properties.SetColor( color_transfer_func )
volume_properties.SetScalarOpacity( opacity_transfer_func )
If you are familiar with VTK, a VTK volume is similar to a VTK actor. It encapsulates the data, properties, and rendering method for one ‘object’ in the scene.
volume = vtk.vtkVolume()
volume.SetMapper( volume_mapper )
volume.SetProperty( volume_properties )
Now lets create a few objects used for rendering.
renderer = vtk.vtkRenderer()
render_window = vtk.vtkRenderWindow()
window_interactor = vtk.vtkRenderWindowInteractor()
The renderer handles the volume rendering, the render window is where the rendered volume will appear, and the window interactor handles user input to adjust the viewpoint (zoom, rotate, move, etc). Hook ‘em up, and add our volume.
render_window.AddRenderer( renderer )
window_interactor.SetRenderWindow( render_window )
renderer.AddVolume( volume )
Just render the scene and start the window interactor so you can rotate the volume with your mouse.
render_window.Render()
window_interactor.Start()
To stop the interactor hit ‘q’.
Numpy and ITK
Working with images with ITK is great, but when working in python, you want all the functionality of python. Fortunately, converting ITK images to python arrays is simple with the python-insighttoolkit-extras package.
itk_py_converter = itk.PyBuffer[image_type]
image_array = itk_py_converter.GetArrayFromImage( reader.GetOutput() )
How about converting a 10x10x10 python array to an ITK image.
import scipy
another_image_array = scipy.zeros( (10,10,10) )
itk_image = itk_py_converter.GetImageFromArray( another_image_array.tolist() )
Does it get any easier?
No Trackbacks
You can leave a trackback using this URL: http://paulnovo.us/wrapitktutorial/trackback
25 Comments
Thanks, Paul. These are great! My think my next big programming project will be in python now…
Thanks Paul. It is really a time saver, the Ubuntu repository.
I’m working on a thesis in Spain about MRI images processing. Just starting, I’d been thinking about working in C++ because I know it. But I think I should learn python, to code faster my experiments.
Will keep in touch.
Regards.
http://alexsavio.blogspot.com
terrific! his sure does justice to the excellent WrapITK project!
Hello and thanks for creating this tutorial.
I’m attempting this tutorial with pythonxy in windows and it seems to almost work, but I get the error listed below. Could you give some advice on how to fix this problem. I am quite new to using ITK. Any help is appreciated.
>>> median_filter = itk.MedianImageFilter[image_type, image_type].New()
KeyError: “itkTemplate : No template (<class ‘itkImage.itkImageUC3′>, <class ‘itkImage.itkImageUC3′ > ) for the itk::MedianImageFilter class”
I am not sure how pythonxy compiles ITK. If they compile ITK-WrapITK with the standard settings, as I suspect, it would not include support for Unsigned Chars. I have added support for extra data-types, and 4D data to my packages.
This is great Paul, I’m having more success in a few minutes with your package than a couple of days with the ‘raw’ itk src – thanks!!
I was wondering it your itk package is feature complete? There’s an image filter that I was looking forward to playing with: itkBayesianClassifierImageFilter which is in the cxx code but I can’t find it in the python lib. Am I missing something?
Thanks again, Dan
Glad to hear you are finding this useful. I have also found experimenting in python a powerful way to learn ITK.
You are right, some itk modules are not wrapped for python. Most of the time this is simple to fix, you just need to create a small file for each module that describes howto wrap the modules for python. ITK includes these for many of the modules, but not all. You can either: wait for ITK to expand their coverage of modules, convince me to figure this out and add a couple new modules, or learn how to do it and send me the changes and I’ll add them to the ubuntu packages.
-Paul
hi Paul,
Do you know how to Defining Origin and Spacing use ITK-WrapITK?
I use”itk_image.SetSpacing(1.5,1.5,2.0)” to define the space of a ITK image, but failed.
Can you give a help?
many thanks~
xh
@xhwang: The spacing is set a little different than that. Probably the easiest way to do it in Python is the following (for 1.5, 1.5, 2.0 spacing):
itk_image.GetSpacing().SetElement(0, 1.5)
itk_image.GetSpacing().SetElement(1, 1.5)
itk_image.GetSpacing().SetElement(2, 2.0)
To set the spacing of the output image using the input image, you can do such thing :
spacing = reader.GetOutput().GetSpacing()
output_image.SetSpacing(spacing)
However, I have some trouble to have access to the spacing values !
spacing[0] does not work for instance :
TypeError: ‘itkVectorD3Ptr’ object does not support indexing
Do you know how to do this ?
spacing isn’t a python array. To get the values, use GetElement:
spacing = reader.GetOutput().GetSpacing()
spacing_x = spacing.GetElement(0)
spacing_y = spacing.GetElement(1)
spacing_z = spacing.GetElement(2)
Hello Paul, I’m pretty new with ITK and Python, but not so new in programming and image analysis, I get an error in ipython when I write import itk: No module named itk. I did follow your first tutorial with no problems (I’m working in Ubuntu karmic). Do you know why this happen and can I fix it?
I’d like to use images of type met_short but if I write
image_type = itk.Image[itk.MS, 3]
it dosen’t work.
What’s the right command?
Thank you very much.
I am not sure what datatype met_short is? Are you referring to the datatype names used in Meta Images? I am pretty sure that a met_short corresponds to a normal signed short for meta images. So you would just need to use the following:
image_type = itk.Image[itk.SS, 3]
@paul:
That command doesn’t work!
If I write
image_type = itk.Image[itk.F, 3]
it works, but I don’t want use float type!
how can i access to the element of matrix type(e.g. image.GetDirection())? and do multiplication with the vector.
@gordon: For image direction, it is a little hairy. But you can get a scipy matrx with the following:
import scipy scipy_matrix = scipy.matrix(scipy.zeros((3,3))) vnl_matrix = i.GetDirection().GetVnlMatrix() for i in range(3): for j in range(3): scipy_matrix[i,j] = vnl_matrix.get(i,j)With the scipy matrix you can now do all sorts of matrix manipulations including vector multiplication. I would look for a scipy tutorial if you are not familiar with it.
@paul thanks for your work on maintaining this package! I noticed somewhere mentioned there is a optional patch to enable access by index to the structure such as itk.Size.
http://mima2.jouy.inra.fr/darcs/contrib-itk/WrapITK/patch/optional/
Is it possible to apply them to your build?
@gordon: Hmm, that looks like a great feature. I’ll look into it.
Thanks, this was really useful.
Has anyone tried using image registration with python? I’m trying to convert Imageregistration7 example and getting the error:
RuntimeError: /tmp/buildd/insighttoolkit-3.18.0/Code/Numerics/itkRegularStepGradientDescentBaseOptimizer.cxx:197:
itk::ERROR: RegularStepGradientDescentOptimizer(0xad25f0e8): The size of Scales is 6, but the NumberOfParameters for the CostFunction is 4.
Does anyone know how to increase the number of parameters in python?
Hi, I am wondering if itk.Mesh has a wrapper. I tried doing as follows
itk.Mesh[itk.F,3]
it say that there is no such template. Or I am not doing it correctly???
Thanks
@shan I don’t have floats wrapped for that type, but it looks like doubles are. If you use ipython you can figure out what types are wrapped and which aren’t. See my comment here.
I stumbled the problem that itk.Image.UL2/3 exist and is returned by filters such as itk.WatershedImageFilter. But practically no filter templates exist for this type. Especially frustrating is that neither itk.ImageFileWriter nor itk.RescaleIntensityImageFilter nor the most simple itk.CastImageFilter can handle the UL type.
Any suggestion how I could retrieve my data without having to write my own conversion filter?
Greetings.
@loli I am not compiling unsigned longs(UL) in these packages. It looks like the WatershedImageFilter has a hard-coded output of UL, so that is why it returns UL. Unfortunately, you would have to recompile the packages with UL turned on, so that all the other filters will also be able to use this data type.
@paul: That sadly poses a problem for me, as ITK won’t compile with the python bindings enabled on my machine.
I found two workarounds, which both have their problems:
I. Use itk.RelabelComponentImageFilter, as it supports UL input
The problem is here that when more than ushort (~65.000) labels are created, the ones out of range are concatenated
II. Cast to VTK and save as MetaIO image
This solves the problem of I, but one is restricted to the MetaIO format and the execution takes signif. longer
Any hint how I could recompile your pre-compiled package with UL support enabled?