Using ANTS for matching to the template : A story of some of the things that can go wrong¶
In this Notebook, we show some issue that arise with spatial normalization, using ANTS. http://stnava.github.io/ANTs/
Note that most of these issues are not specific to ANTS : any spatial normalisation program will have simlar hurdles.
It is important to start by stating that we don't know how to normalize one brain on another, and while we don't discuss the fundamental proble here, some of the results shown examplify this. In particular, we dont know in general how to set the regularisation parameter that will balance the amount of deformation (applied to the brain to be normalized) and the similarity to the template.
Later in this notebook, we show results with different regularizations.
%pylab inline
import numpy as np
import matplotlib.pylab as plt
import os
import nibabel as nib
# assume we have the current directory the "mni_icbm152_nlin_asym_09a" directory that contains the MNI template
# if not, this can be downloaded from : http://www.bic.mni.mcgill.ca/ServicesAtlases/ICBM152NLin2009 -
# choose the ICBM 2009a Nonlinear Asymmetric 1×1x1mm template
MNI_PATH = 'mni_icbm152_nlin_asym_09a'
os.listdir(MNI_PATH)
Show the template :
mni_fname = os.path.join(MNI_PATH, 'mni_icbm152_t1_tal_nlin_asym_09a.nii')
mni_img = nib.load(mni_fname)
plt.imshow(mni_img.get_data()[:, :, 30], cmap="gray")
print mni_img.get_affine()
print mni_img.shape
# where is our image to register:
DATA_PATH = os.path.join(os.path.expanduser('~'), 'data', 'ds105', 'sub001', 'anatomy')
os.listdir(DATA_PATH)
anat_fname = os.path.join(DATA_PATH, 'highres001.nii.gz')
print anat_fname
anat_img = nib.load(anat_fname)
plt.imshow(anat_img.get_data()[:, 70, :], cmap="gray")
print anat_img.get_affine()
print anat_img.shape
ANTS is picky about the affine information in the image. The subject image doesn't have a proper affine in the header: look at the srow_{x,y,z} fields
anat_hdr = anat_img.get_header()
print(anat_hdr)
Load our new fixnifti library
#- Make sure we have pna-utils on the path
!echo $PYTHONPATH
#- import the fixing tool
import fixnifti
For ANTS to work properly, we will write a new copy with a proper affine, by default prefixed by "f". In fact nibabel will guess a good-enough affine for us; we use that.
fixnifti.fixup_nifti_file(anat_fname)
fixed_fname = os.path.join(DATA_PATH, 'fhighres001.nii.gz')
# check that the new header contains information for the srow_{x,y,z}
fixed_anat = nib.load(fixed_fname)
fixed_anat_hdr = fixed_anat.get_header()
print(fixed_anat_hdr)
#save the fixed image in the current directory, call it 'my_fixed_anat.nii'
nib.save(fixed_anat, 'my_fixed_anat.nii')
# make a nibabel image
PNA_PATH = '/home/jb/code/pna-notebooks'
fixed_fname = os.path.join(PNA_PATH, 'my_fixed_anat.nii')
fixed_img = nib.load(fixed_fname)
First we check ANTS is on the path. There are many ways of usign ANTS, because there exist a number of commandes (shell scripts invoking eventually ANTS) with different options. We will show the use of 2 of these scripts.
!ants.sh
We set ANTSPATH in the environment if it isn't
os.environ['ANTSPATH'] = '/home/jb/bin/ants/'
The ANTS program itself and its help can be lauched with :
ANTS --help¶
But we are going to use the script ants.sh first with its default parameters. ANTS takes some time to run ... on my laptop it takes between an hour and an hour and a half.
# RUN THIS IF YOU HAVE TIME ...
# Output images will be prefixed with nrmlzd_"
# !ants.sh 3 my_fixed_anat.nii mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii nrmlzd_
Now comes the critical step. Look at the results of the registration:¶
# we use viz from nipy:
from nipy.labs import viz
normed_img_fname = os.path.join(PNA_PATH, 'nrmlzd_deformed.nii.gz')
normed_img = nib.load(normed_img_fname)
slicer = viz.plot_anat(normed_img.get_data(), normed_img.get_affine(), dim=.2, black_bg=True)
slicer.edge_map(mni_img.get_data(), mni_img.get_affine())
Look at the normalized image : you see above that it didn't register the two brains too well...¶
To help compare two or more images, let's do a little helper function:
def show_imgs(imgs, pos, titles=None, col="gray", interp="nearest"):
"""
"""
Nimg = len(imgs)
if (titles is None):
titles = ['No title']*Nimg
if not isinstance(pos, list):
pos = [pos]*Nimg
fig, ax = plt.subplots(1,Nimg)
for ix,im in enumerate(imgs):
ax[ix].imshow(im.get_data()[pos[ix]], cmap=col)
ax[ix].get_xaxis().set_ticks([])
ax[ix].get_yaxis().set_ticks([])
ax[ix].set_title(titles[ix])
pos = (slice(None), slice(None), 64)
show_imgs([mni_img, normed_img], [pos,pos], ['mni_img', 'normed_img'])
print normed_img.get_affine(), normed_img.shape
print mni_img.get_affine(), mni_img.shape
Note that only the viz plot can tell us about the registration: the show_imgs above doesn't take into account the affine. But since the results are ugly on the viz plot let's try something else.
We see that the global position is not right. We can try:
- The image has a lot of the neck, the template has not. This can trick our registration.
- The 0,0,0 position of the template and the image can be far off. making the origin close to AC in the image could help.
Below, we first try the first step, and register using a rigid transform to see if we can at least get this.
Cut his neck (a good French tradition)¶
#PNA_PATH = '/home/jb/code/pna-notebooks'
#fixed_fname = os.path.join(PNA_PATH, 'my_fixed_anat.nii')
#fixed_img = nib.load(fixed_fname)
data = fixed_img.get_data()
aff = fixed_img.get_affine().copy()
# cut from z=40 to z=225
lz, hz = 40, 225
data_cut_neck = data[:,:,lz:hz]
aff[2,3] = -(hz-lz)/2. # put the origin in the middle of the z
print aff
cut_neck_fname = os.path.join(PNA_PATH,"neck_cut.nii.gz")
cut_neck_img = nib.Nifti1Pair(data_cut_neck, aff)
nib.save(cut_neck_img, cut_neck_fname)
# re-read and display
new_img = nib.load(cut_neck_fname)
plt.imshow(new_img.get_data()[:, 80, :], cmap="gray")
print new_img.get_affine()
print new_img.shape
Let's use another ants script (antsIntroduction.sh) to do a first rigid normalisation to the template:¶
# get help with : !antsIntroduction.sh -h
# launch it with option : -t RI (for Rigid) and prefix the results with "rigid_"
!antsIntroduction.sh -d 3 -i neck_cut.nii.gz -r mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii -o rigid_ -t RI
# !ANTS 3 -i fixed_anat.nii -r mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii -o rigid_ -t RI
#!ANTS 3 -i neck_cut.nii.gz -r mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii -o nrmlzd_ -s MI -t RI
#--number-of-affine-iterations 10000x10000x10000
# look at the results:
rigid_fname = os.path.join(PNA_PATH, 'rigid_deformed.nii.gz')
rigid_img = nib.load(rigid_fname)
slicer = viz.plot_anat(rigid_img.get_data(), rigid_img.get_affine(), dim=.2, black_bg=True)
slicer.edge_map(mni_img.get_data(), mni_img.get_affine())
This is just a rigid transform, but it has worked. Let's start our non linear transform from here!¶
pos = (slice(None), slice(None), 64)
show_imgs([mni_img, rigid_img], [pos,pos], ['mni_img', 'rigid_img'])
# re-try with the rigid transform as the input
# Use antsIntroduction that gives us more flexibility : the -m 10x30x5 states the number of iterations
# by default 30x90x20
# !antsIntroduction.sh -d 3 -i rigid_deformed.nii.gz -r mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii \
# -o from_rigid_ -t SY -m 10x30x5
renormed_fname = os.path.join(PNA_PATH, 'from_rigid_deformed.nii.gz')
renormed_img = nib.load(renormed_fname)
slicer = viz.plot_anat(renormed_img.get_data(), renormed_img.get_affine(), dim=.2, black_bg=True)
slicer.edge_map(mni_img.get_data(), mni_img.get_affine())
This sounds better. But ... is it not too good ?¶
pos = (slice(None),40,slice(None))
show_imgs([mni_img, renormed_img], pos, ['mni_img', 'renormed_img'])
What we see here is that the elasticity of the deformation was too large - it has done some unatural deformation on our brain. Let's have a look at the deformations fields:¶
# load the deformation maps
# use slice objects
warp_fname = os.path.join(PNA_PATH, 'from_rigid_Warp.nii.gz')
warp_img = nib.load(warp_fname)
print warp_img.shape
print renormed_img.shape
print mni_img.shape
# The warp image is 5 dimensional : TO DO : EXPLAIN ...
pos = (slice(None),50,slice(None))
pos = [pos, pos + (0, 0), pos]
show_imgs([renormed_img, warp_img, mni_img], pos, ['renormed', 'warp', 'mni'])
The field of deformation is not smooth ... it has done crazy deformations to our brain...¶
# increase regularisation using the following command:
#
# Compare the results with the previous normalized file
renormed_25_fname = os.path.join(PNA_PATH, 'from_rigid_25_deformed.nii.gz')
renormed_25_img = nib.load(renormed_25_fname)
pos = (slice(None),50,slice(None))
show_imgs([renormed_img, renormed_25_img, rigid_img], [pos, pos, pos], ['full elastic', '.25 regul', 'rigid'])
Rigid and .25 look similar - try less regularization: 0.1¶
renormed_01_fname = os.path.join(PNA_PATH, 'from_rigid_01_deformed.nii.gz')
renormed_01_img = nib.load(renormed_01_fname)
pos = (slice(None),50,slice(None))
show_imgs([mni_img, renormed_img, renormed_01_img, rigid_img],
[pos, pos, pos, pos],
['MNI', '0 Regul', '.01 Regul', 'Rigid'])
To conclude: we should investigate something between the 0 regularization and the .01 regularization, or other ANTS options. At the moment,.01 seems to be very close to .25, and to rigid.¶