In [1]:
# load the module
import irtk
# read an image
img = irtk.imread("/home/kevin/Imperial/PhD/DATASETS/Originals/5066_4.nii")
print img[0:2,0:2,0:2] # we crop before printing to save space on paper...
In [2]:
# we previously read the image in its native type (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("/home/kevin/Imperial/PhD/DATASETS/Originals/5066_4.nii",
dtype='float32')
print "Maximum with slope/intercept", float(img.max())
In [3]:
# show us a quick view of the image with a little bit of saturation
irtk.imshow(img.saturate(0.01,0.99))
Out[3]:
In [4]:
# now a quick segmentation
irtk.imshow(img.saturate(0.01,0.99), img > 2000, opacity=0.4)
Out[4]:
In [5]:
# for a better visualisation, launch rview or display:
irtk.rview( img, img > 2000 )
irtk.display( img, img > 2000 )
# 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
irtk.imshow( img, irtk.slic(img, 1000, 10) )
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,
interpolation='linear' )
irtk.imshow(new_img)
Out[6]:
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
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)
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("/home/kevin/Imperial/PhD/DATASETS/Originals/890_*.nii")
print "Number of files:", len(filenames)
img1 = irtk.imread(filenames[0], dtype='float32')
img1 = img1.rescale()
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()
t = img2.register(img1) # rigid registration
img2 = img2.transform(t,target=img1) # linear interpolation by default
img1 += img2
weights += (img2 > 0).astype('int')
stop = time()
img1 /= weights.astype('float32') + 0.000001 # division by zero is bad
# write image to disk
irtk.imwrite( "average.nii", img1)
print "Elapsed time:", (stop - start) / 60 , "minutes"
irtk.imshow(img1)
Out[8]:
In [9]:
# it is still a work in progress, with a focus so far on rigid transformation
# here are the functions currently available:
print irtk.__all__