Tuesday, 17 December 2013

IRTK Python

Python for Medical Imaging

There are several reasons for choosing Python for Medical Imaging. First of all, the main biomedical imaging library out there, ITK, offers both exhaustive SWIG wrappers and a simplified interface called SimpleITK, which can be directly installed from PyPI. Other major libraries, such as VTK for data visualization and OpenCV for Computer Vision also provide excellent Python wrappers.

Then comes the Python eco-system itself, whose ressources for Machine Learning and Image Processing have been growing in the last years with the main libraries: scipy.ndimage, scikit-image and scikit-learn. There are also Python modules for input/output of medical images, such as pydicom and nibabel.

Finally, enough is never enough, there is cython to easily mix python and C/C++ code.

But there's more to Python than a set of existing libraries. It is simple to write, yet powerful. It has a wide community, much wider than the limited scope of scientific computing. It is Open Source, which is a key element for research. It offers several possibilities for optimizing code, for instance making an efficient use of numpy, using joblib for parallelization or cython for C++ optimization.

Why another Medical Imaging library in Python?

IRTK is a Medical Imaging library developped within the BioMedIA group at Imperial College London, initially built around Daniel Rueckert's non-rigid registration using free-form deformations. While ITK deals with ND-images, can read most image formats and has a wide range of applications ranging from cell microscopy to satellite images, IRTK is much more focused, its core component being image registration. It can read NIFTI files, and write NIFTI and PNG files, and all images in IRTK are 4D (X, Y, Z, T), with flat dimensions where necessary. The code is only templated other VoxelType.

The main reason I chose to write a Python interface for IRTK instead of using what already existed is that it is much easier to collaborate within the lab if we all use a same code base. Moreover, I wanted medical images to be as much pythonic as possible, namely a subclass of numpy.ndarray which would automatically update its spatial coordinates when cropped.

As I wanted the C++ code to be free of Python syntax, so that others can use it without Python, and I wanted the Python code to contain pure Python objects, I chose the following code organization:
Python <--> Cython <--> C++

template <class dtype>
void irtk2py( irtkGenericImage<dtype>& irtk_image,
              dtype* img,
              double* pixelSize,
              double* xAxis,
              double* yAxis,
              double* zAxis,
              double* origin,
              int* dim );

Considering that a typical image registration in C++ is not instantaneous, I conveniently copy and reshape arrays quite often in the background instead of sharing buffers (which I would not know how to do properly) or maintaining a C++ object accessible through Python (which would not be a numpy.ndarray anymore).

Below is a simple comparison between itk, SimpleITK and irtk:
In [2]:
import itk

# read an image
pixelType = itk.F
imageType = itk.Image[ pixelType, 3 ]
readerType = itk.ImageFileReader[ imageType ]
reader = readerType.New()
reader.SetFileName( "input.nii" )

itk_image = reader.GetOutput()

# get an numpy array
itk2np = itk.PyBuffer[ imageType ]
data = itk2np.GetArrayFromImage( itk_image )

# do some processing
print type( itk_image )
print type( data )

# create a new itk image, preserving spatial information
new_itk_image = itk2np.GetImageFromArray( data )
new_itk_image.SetSpacing( itk_image.GetSpacing() )
new_itk_image.SetOrigin( itk_image.GetOrigin() )
new_itk_image.SetDirection( itk_image.GetDirection() )

# write to a file
writerType = itk.ImageFileWriter[ imageType ]
writer = writerType.New()
writer.SetInput( new_itk_image )
writer.SetFileName( "output.nii" )
<class 'itkImage.itkImageF3_PointerPtr'>
<type 'numpy.ndarray'>

In [1]:
import SimpleITK as sitk

# read image and get numpy array
img = sitk.ReadImage( "input.nii" )
data = sitk.GetArrayFromImage( img )

# do some processing
print type( img )
print type( data )

# create a new sitk image, preserving spatial information
output = sitk.GetImageFromArray( data )
output.SetSpacing( img.GetSpacing() )
output.SetOrigin( img.GetOrigin() )
output.SetDirection( img.GetDirection() )

# write to a file
sitk.WriteImage( output, "output.nii" )
<class 'SimpleITK.Image'>
<type 'numpy.ndarray'>

In [2]:
import irtk

# read image
img = irtk.imread( "input.nii" )

# do some processing
print type(img)

# write image
irtk.imwrite( "output.nii", img )
<class 'irtk.image.Image'>



The code is available for download as part of IRTK: https://github.com/sk1712/IRTK/tree/master/wrapping/cython
git clone https://github.com/sk1712/IRTK.git
mkdir build
cd build
make -j 3
Lastly, you need to update your PYTHONPATH:
Get the full path of irtk/build/lib and for bash, add to your ~/.bashrc
export PYTHONPATH=full_path/irtk/build/lib:$PYTHONPATH
for tcsh, add to your ~/.cshrc
setenv PYTHONPATH full_path/irtk/build/lib:$PYTHONPATH
The documentation for IRTK Python is available online.

Example usage

In [1]:
# load the module
import irtk

# read an image
img = irtk.imread( "input.nii" )

print img[0:2,0:2,0:2] # we crop before printing to save space on the page
Image data:
array([[[0, 0],
        [0, 0]],

       [[0, 0],
        [0, 0]]], dtype=int16)

Image shape: (2, 2, 2)

{   'dim': array([2, 2, 2, 1], dtype=int32),
    'orientation': array([[-0.97086735,  0.10419369,  0.2157784 ],
       [ 0.23358087,  0.2106676 ,  0.94924135],
       [-0.05344743, -0.9719891 ,  0.22886794]]),
    'origin': array([ 118.94322308,   90.36445636, -187.99950024,    0.        ]),
    'pixelSize': array([ 1.18055558,  1.18055558,  1.265625  ,  1.        ])}

In [2]:
# we previously read the image in its type on disk (short == int16)
# we thus had the warning on stderr:
# irtkGenericImage<short>::Read: Ignore slope and intercept, 
# use irtkGenericImage<float> or irtkGenericImage<double> instead
print "Maximum without slope/intercept", float(img.max())

# let's read it again requesting for float ( float == float32 )
img = irtk.imread( "input.nii", dtype='float32' )
print "Maximum with slope/intercept", float(img.max())
Maximum without slope/intercept 414.0
Maximum with slope/intercept 5461.76513672

In [3]:
# show us a view of the image with some saturation
irtk.imshow( img.saturate(0.01,0.99), filename="quickview.png" )
In [4]:
# now a segmentation
irtk.imshow( img.saturate(0.01,0.99), # input image
            img > 2000,               # segmentation labels
            filename="segmentation.png" )
In [5]:
# Modules from irtk.ext need to be built separately
# by running make in irtk/wrapping/cython/ext

# segmentation colormaps are automatically generated
# the first 10 colors are predefined, the remaining ones are random.
# this works for rview, display and the quick imshow tool
# (the latter works only for ipython notebook or qtconsole, 
# unless you specify a filename for writing to disk)
# let's try a more comlex segmentation like SLIC supervoxels
from irtk.ext.slic import slic
irtk.imshow( img,
             slic( img.gaussianBlurring(2.0), 2000, 10), # segmentation labels
             filename="slic.png" )
In [6]:
# the advantage of a this new interface over the old SWIG wrapper is
# to offer a pythonesque access to IRTK.
# a rule of thumb is that the Python code should be shorter,
# and easier to read than C++

# for instance, if we want to rotate an image around its center (pixel coordinates):

# get the translation
tx,ty,tz = img.ImageToWorld( [(img.shape[2]-1)/2,
                              (img.shape[0]-1)/2] )
translation = irtk.RigidTransformation( tx=-tx, ty=-ty, tz=-tz )

# form a rotation
rotation = irtk.RigidTransformation( ry=60 )

# apply the transformations
new_img = img.transform( translation.invert()*rotation*translation,
                         interpolation='linear' )

irtk.imshow( new_img, filename="rotation.png" )
In [7]:
# which can be compared to the un-centered rotation:
new_img2 = img.transform( rotation,
                          interpolation='linear' )

import numpy as np
diff = new_img-new_img2

center_of_rotation = irtk.zeros(img.get_header(), dtype='int')
shape = center_of_rotation.shape
                   shape[2]/2-5:shape[2]/2+5] = 1 # will be red

world_center = img.WorldToImage([0,0,0])
# it happens to be outside of the image, so we represent it as a "cylinder"
print "World center:", world_center
print "Image dimensions (XYZT):", img.header['dim']
                   world_center[0]-5:world_center[0]+5] = 2 # will be green
irtk.imshow( diff, 
             filename="diff.png" )
World center: [125 112 109]
Image dimensions (XYZT): [288 288  80   1]

In [8]:
# here is a second example:
# we read a bunch of scans, register them to the first one,
# and compute an average volume in the coordinate system of the first file
from glob import glob
from time import time

filenames = glob("890_*.nii")

print "Number of files:", len(filenames)

img1 = irtk.imread(filenames[0], dtype='float32')
img1 = img1.rescale().resample(2)
weights = (img1 > 0).astype('int')

# while the registration is running, you can see its usual outpout on stdout
start = time()
for f in filenames[1:]:
    img2 = irtk.imread( f, dtype='float32')
    img2 = img2.rescale().resample(2)
    t = img2.register(img1) # rigid registration
    img2 = img2.transform(t,target=img1, interpolation='nearest') 
    img1 += img2
    weights += (img2 > 0).astype('int')
stop = time()

img1 /= weights.astype('float32') + 0.000001 # dividing by zero is bad

# write image to disk
irtk.imwrite( "average.nii", img1)

print "Elapsed time:", (stop - start) / 60 , "minutes"
irtk.imshow( img1, filename="average.png" )
Number of files: 10
Elapsed time: 3.58047141631 minutes

In [9]:
# finally, an example where joblib is used to parallelise the registrations
import irtk
from joblib import Parallel, delayed
from glob import glob
from time import time

def register( f, img1):
    img2 = irtk.imread( f, dtype='float32')
    img2 = img2.rescale().resample(2)
    t = img2.register(img1)
    return t

filenames = glob( "890_*.nii")
img1 = irtk.imread(filenames[0], dtype='float32')
img1 = img1.rescale().resample(2)

transformations = [irtk.RigidTransformation()]

start = time()
transformations.extend( Parallel(n_jobs=2)(delayed(register)( f,
                         for f in filenames[1:] )
stop = time()
print "Elapsed time:", (stop - start) / 60 , "minutes"
print transformations[-1]
Elapsed time: 3.21967080037 minutes
tx: 3.11568652093 (mm)
ty: 4.90160108265 (mm)
tz: -8.44972466072 (mm)
rx: -2.60710439831(degrees)
ry: 3.20516000688(degrees)
rz: -7.51380234957 (degrees)

This Python interface is still far from complete compared to what is available in IRTK, but I am only taking it as far as I need for my PhD. Do let me know if you use it, if you have ideas for improvement, and you can directly modify it as you like on github.


No comments:

Post a Comment