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.
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:
Below is a simple comparison between itk, SimpleITK and irtk:
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 );
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" )
reader.Update()
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" )
writer.Update()
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" )
In [2]:
import irtk
# read image
img = irtk.imread( "input.nii" )
# do some processing
print type(img)
# write image
irtk.imwrite( "output.nii", img )
Out[2]:
Installation¶
The code is available for download as part of IRTK: https://github.com/sk1712/IRTK/tree/master/wrapping/cython
Get the full path of
git clone https://github.com/sk1712/IRTK.git
cd IRTK
mkdir build
cd build
cmake -D BUILD_WITH_PNG=ON -D WRAP_CYTHON=ON ..
make -j 3Get the full path of
irtk/build/lib and for bash, add to your ~/.bashrcexport PYTHONPATH=full_path/irtk/build/lib:$PYTHONPATH~/.cshrcsetenv PYTHONPATH full_path/irtk/build/lib:$PYTHONPATHExample 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
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())
In [3]:
# show us a view of the image with some saturation
irtk.imshow( img.saturate(0.01,0.99), filename="quickview.png" )
Out[3]:
In [4]:
# now a segmentation
irtk.imshow( img.saturate(0.01,0.99), # input image
img > 2000, # segmentation labels
opacity=0.4,
filename="segmentation.png" )
Out[4]:
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" )
Out[5]:
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[1]-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,
target=img,
interpolation='linear' )
irtk.imshow( new_img, filename="rotation.png" )
Out[6]:
In [7]:
# which can be compared to the un-centered rotation:
new_img2 = img.transform( rotation,
target=img,
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
center_of_rotation[shape[0]/2-5:shape[0]/2+5,
shape[1]/2-5:shape[1]/2+5,
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']
center_of_rotation[:,
world_center[1]-5:world_center[1]+5,
world_center[0]-5:world_center[0]+5] = 2 # will be green
irtk.imshow( diff,
center_of_rotation,
opacity=1.0,
filename="diff.png" )
Out[7]:
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" )
Out[8]:
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,
img1)
for f in filenames[1:] )
)
stop = time()
print "Elapsed time:", (stop - start) / 60 , "minutes"
print transformations[-1]
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.
References¶
Blogging with the IPython notebook, by Fernando Perez.
cyitk, A simple ITK-to-Python bridge written in Cython, by Luca Antiga
SimpleITK wiki
cyitk, A simple ITK-to-Python bridge written in Cython, by Luca Antiga
SimpleITK wiki
No comments:
Post a Comment