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
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
# 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)
['mni_icbm152_t1_tal_nlin_asym_09a_face_mask.nii', 'mni_icbm152_pd_tal_nlin_asym_09a.nii', 'mni_icbm152_gm_tal_nlin_asym_09a.nii', 'mni_icbm152_t1_tal_nlin_asym_09a.nii', 'mni_icbm152_csf_tal_nlin_asym_09a.nii', 'mni_icbm152_t1_tal_nlin_asym_09a_eye_mask.nii', 'mni_icbm152_t1_tal_nlin_asym_09a_InverseWarp.nii.gz', 'mni_icbm152_wm_tal_nlin_asym_09a.nii', 'mni_icbm152_t2_tal_nlin_asym_09a.nii', 'mni_icbm152_t1_tal_nlin_asym_09a_mask.nii']
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
[[ 1. 0. 0. -98.] [ 0. 1. 0. -134.] [ 0. 0. 1. -72.] [ 0. 0. 0. 1.]] (197, 233, 189)
# where is our image to register:
DATA_PATH = os.path.join(os.path.expanduser('~'), 'data', 'ds105', 'sub001', 'anatomy')
os.listdir(DATA_PATH)
['highres001_brain_mask.nii.gz', 'fhighres001.nii.gz', 'highres001.nii.gz', 'highres001_brain.nii.gz']
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
/home/jb/data/ds105/sub001/anatomy/highres001.nii.gz [[ -1.20000005 0. 0. 73.80000293] [ 0. 0.9375 0. -119.53125 ] [ 0. 0. 0.9375 -119.53125 ] [ 0. 0. 0. 1. ]] (124, 256, 256)
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)
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<' sizeof_hdr : 348 data_type : db_name : extents : 0 session_error : 0 regular : r dim_info : 0 dim : [ 3 124 256 256 1 1 1 1] intent_p1 : 0.0 intent_p2 : 0.0 intent_p3 : 0.0 intent_code : none datatype : int16 bitpix : 16 slice_start : 0 pixdim : [ 1. 1.20000005 0.9375 0.9375 1. 0. 0. 0. ] vox_offset : 352.0 scl_slope : 1.0 scl_inter : 0.0 slice_end : 0 slice_code : unknown xyzt_units : 10 cal_max : 0.0 cal_min : 0.0 slice_duration : 0.0 toffset : 0.0 glmax : 0 glmin : 0 descrip : FSL4.0 aux_file : qform_code : unknown sform_code : unknown quatern_b : 0.0 quatern_c : 0.0 quatern_d : 0.0 qoffset_x : 0.0 qoffset_y : 0.0 qoffset_z : 0.0 srow_x : [ 0. 0. 0. 0.] srow_y : [ 0. 0. 0. 0.] srow_z : [ 0. 0. 0. 0.] intent_name : magic : n+1
Load our new fixnifti
library
#- Make sure we have pna-utils on the path
!echo $PYTHONPATH
#- import the fixing tool
import fixnifti
/home/jb/code/pna-utils:/home/jb/.local/lib/python2.7/site-packages:/home/jb/code/nbconvert/nbconvert
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')
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<' sizeof_hdr : 348 data_type : db_name : extents : 0 session_error : 0 regular : r dim_info : 0 dim : [ 3 124 256 256 1 1 1 1] intent_p1 : 0.0 intent_p2 : 0.0 intent_p3 : 0.0 intent_code : none datatype : int16 bitpix : 16 slice_start : 0 pixdim : [-1. 1.20000005 0.9375 0.9375 1. 0. 0. 0. ] vox_offset : 352.0 scl_slope : 1.0 scl_inter : 0.0 slice_end : 0 slice_code : unknown xyzt_units : 10 cal_max : 0.0 cal_min : 0.0 slice_duration : 0.0 toffset : 0.0 glmax : 0 glmin : 0 descrip : FSL4.0 aux_file : qform_code : aligned sform_code : aligned quatern_b : 0.0 quatern_c : 1.0 quatern_d : 0.0 qoffset_x : 73.8000030518 qoffset_y : -119.53125 qoffset_z : -119.53125 srow_x : [ -1.20000005 0. 0. 73.80000305] srow_y : [ 0. 0.9375 0. -119.53125] srow_z : [ 0. 0. 0.9375 -119.53125] intent_name : magic : n+1
# 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
USAGE :: sh ants.sh ImageDimension fixed.ext moving.ext OPTIONAL-OUTPREFIX OPTIONAL-max-iterations OPTIONAL-Labels-In-Fixed-Image-Space-To-Deform-To-Moving-Image Option-DoANTSQC be sure to set ANTSPATH environment variable Max-Iterations in form : JxKxL where J = max iterations at coarsest resolution (here, reduce by power of 2^2) K = middle resolution iterations ( here, reduce by power of 2 ) L = fine resolution iterations ( here, full resolution ) -- this level takes much more time per iteration an extra Ix before JxKxL would add another level Default ierations is 30x90x20 -- you can often get away with fewer for many apps Other parameters are updates of the defaults used in the A. Klein evaluation paper in Neuroimage, 2009
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 :
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_
# 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())
/home/jb/.local/lib/python2.7/site-packages/nipy/labs/datasets/volumes/volume_img.py:284: FutureWarning: Numpy has detected that you (may be) writing to an array returned by numpy.diagonal or by selecting multiple fields in a record array. This code will likely break in the next numpy release -- see numpy.diagonal or arrays.indexing reference docs for details. The quick fix is to make an explicit copy (e.g., do arr.diagonal().copy() or arr[['f0','f1']].copy()). pixdim[0] = -pixdim[0]
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
[[ -1.20000005 -0. 0. 73.80000305] [ -0. 0.9375 -0. -119.53125 ] [ 0. 0. 0.9375 -119.53125 ] [ 0. 0. 0. 1. ]] (124, 256, 256) [[ 1. 0. 0. -98.] [ 0. 1. 0. -134.] [ 0. 0. 1. -72.] [ 0. 0. 0. 1.]] (197, 233, 189)
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:
Below, we first try the first step, and register using a rigid transform to see if we can at least get this.
#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)
[[ -1.20000005 0. 0. 73.80000305] [ 0. 0.9375 0. -119.53125 ] [ 0. 0. 0.9375 -92.5 ] [ 0. 0. 0. 1. ]]
# 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
[[ -1.20000005 0. 0. 73.80000305] [ 0. 0.9375 0. -119.53125 ] [ 0. 0. 0.9375 -92.5 ] [ 0. 0. 0. 1. ]] (124, 256, 185)
# 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
input files: checking -------------------------------------------------------------------------------------- Data integrity check failed -------------------------------------------------------------------------------------- There seems to be a problem with the header definitions of the input files. This script has tried to repair this issue automatically by invoking: /home/jb/bin/ants/ImageMath 3 rigid_repaired.nii.gz CompareHeadersAndImages mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii neck_cut.nii.gz The repaired image is: rigid_repaired.nii.gz and should be located in: /home/jb/code/pna-notebooks Try passing the original input image to antsaffine.sh. If the affine normalization succeeds, you can safely ignore this warning. Otherwise, you may need to reorient your date to correct the orientation incompatibility. -------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------- Mapping parameters -------------------------------------------------------------------------------------- ANTSPATH is /home/jb/bin/ants/ Dimensionality: 3 Fixed image: mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii Moving image: neck_cut.nii.gz Use label image: 0 N3BiasFieldCorrection: 0 DoANTS Quality Check: 0 Similarity Metric: PR Transformation: RI Regularization: MaxIterations: 30x90x20 Number Of MultiResolution Levels: 3 OutputName prefix: rigid_ -------------------------------------------------------------------------------------- -------------------------------------------------------------------------------------- ANTS command: /home/jb/bin/ants/ANTS 3 -m MI[mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii,neck_cut.nii.gz,1,32] -o rigid_.nii.gz -i 0 --use-Histogram-Matching --number-of-affine-iterations 10000x10000x10000x10000x10000 --MI-option 32x16000 --do-rigid: true -------------------------------------------------------------------------------------- WARNING: Unknown options --do-rigid: Run Reg values 1 Fixed image file: mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii Moving image file: neck_cut.nii.gz Metric 0: Not a Point-set Fixed image file: mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii Moving image file: neck_cut.nii.gz similarity metric weight: 1 Radius: [32, 32, 32] Radius: [32, 32, 32] Use identity affine transform as initial affine para. aff_init.IsNull()==1 Use identity affine transform as initial fixed affine para. fixed_aff_init.IsNull()==1 Continue affine registration from the input affine_opt.use_rotation_header = 0 affine_opt.ignore_void_orgin = 0 transform_initial: IsNotNull():0 OptAffine: metric_type=AffineWithMutualInformation MI_bins=32 MI_samples=16000 number_of_seeds=0 time_seed=1373763247 number_of_levels=5 number_of_iteration_list=[10000,10000,10000,10000,10000] graident_scales=[1,1,1,1,1,1,1,1,1,1,0.0001,0.0001,0.0001] is_rigid = 0 mask null: 1 maximum_step_length=0.1 relaxation_factor=0.5 minimum_step_length=0.0001 translation_scales=0.0001 opt.transform_initial.IsNull(): 1 opt.use_rotation_header: 0 opt.ignore_void_orgin: 0 input affine center: [-0.108839, 18.5778, 7.95606] input affine para: [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1.04609, -10.4327, -22.9379] level 0, iter 214, size: fix[32, 32, 32]-mov[32, 32, 32], affine para: [-0.136024, -0.021254, -0.0217099, 0.99024, 0.853847, 0.881971, 0.833728, -0.00904119, -0.00271425, -0.00119782, 0.599941, -9.45593, -13.6834] reach oscillation, current step: 9.76563e-05<0.0001 level 1, iter 81, size: fix[32, 32, 32]-mov[32, 32, 32], affine para: [-0.14322, -0.0233183, -0.0232153, 0.989144, 0.853777, 0.881911, 0.838739, -0.0104005, -0.00599301, -0.00712232, 0.711311, -9.29394, -13.8528] reach oscillation, current step: 9.76563e-05<0.0001 level 2, iter 84, size: fix[49, 58, 47]-mov[40, 64, 46], affine para: [-0.136036, -0.018551, -0.0225377, 0.990274, 0.858903, 0.887971, 0.843682, -0.0066924, -0.00698513, -0.00742868, 0.674914, -9.55677, -13.6831] reach oscillation, current step: 9.76563e-05<0.0001 level 3, iter 143, size: fix[99, 117, 95]-mov[79, 128, 93], affine para: [-0.135301, -0.0163787, -0.032014, 0.990152, 0.865197, 0.890771, 0.844465, -0.0131004, -0.00559151, -0.00524861, 0.895805, -9.38701, -13.4516] reach oscillation, current step: 9.76563e-05<0.0001 level 4, iter 137, size: fix[197, 233, 189]-mov[124, 256, 185], affine para: [-0.132035, -0.0195311, -0.0266981, 0.990693, 0.867563, 0.892625, 0.843917, -0.014863, -0.00534956, -0.012821, 0.794206, -9.40098, -13.5497] reach oscillation, current step: 9.76563e-05<0.0001 level 0, iter 251, size: fix[32, 32, 32]-mov[32, 32, 32], affine para: [0.0148373, -0.134694, 0.990391, -0.0276312, -1.1596, -1.12891, 1.16635, 0.0104492, 0.00448012, 0.030759, -0.0131428, 7.09968, 12.1662] reach oscillation, current step: 9.76563e-05<0.0001 level 1, iter 25, size: fix[32, 32, 32]-mov[32, 32, 32], affine para: [0.0170129, -0.130576, 0.990892, -0.0281714, -1.1597, -1.12912, 1.16746, 0.0108979, 0.00437984, 0.0305073, 0.0314324, 7.15781, 12.0045] reach oscillation, current step: 9.76563e-05<0.0001 level 2, iter 73, size: fix[40, 64, 46]-mov[49, 58, 47], affine para: [0.0185747, -0.135831, 0.990219, -0.0259118, -1.15781, -1.12365, 1.16641, 0.010386, 0.00248833, 0.0277585, -0.1085, 7.34215, 11.8598] reach oscillation, current step: 9.76563e-05<0.0001 level 3, iter 150, size: fix[79, 128, 93]-mov[99, 117, 95], affine para: [0.014477, -0.134876, 0.99052, -0.0216631, -1.15709, -1.1245, 1.1777, 0.00497012, 0.00276196, 0.0338989, -0.24611, 7.25086, 11.5775] reach oscillation, current step: 9.76563e-05<0.0001 level 4, iter 85, size: fix[124, 256, 185]-mov[197, 233, 189], affine para: [0.0165078, -0.137977, 0.989929, -0.0270101, -1.15299, -1.12417, 1.17417, 0.00349944, 0.00174478, 0.031315, -0.25929, 6.99053, 12.0323] reach oscillation, current step: 9.76563e-05<0.0001 v1 -0.61317 v2 -0.610286 final [-0.132035, -0.0195311, -0.0266981, 0.990693, 0.867563, 0.892625, 0.843917, -0.014863, -0.00534956, -0.012821, 0.794206, -9.40098, -13.5497] outputput affine center: [-0.108839, 18.5778, 7.95606] output affine para: [-0.132035, -0.0195311, -0.0266981, 0.990693, 0.867563, 0.892625, 0.843917, -0.014863, -0.00534956, -0.012821, 0.794206, -9.40098, -13.5497] initial measure value (MMI): rval = -0.312316 final measure value (MMI): rval = -0.620159 finish affine registeration... Requested Transformation Model: SyN : Using SyN diffeomorphic model for transformation. Grad Step 0.5 total-smoothing 0.5 gradient-smoothing 3 setting N-TimeSteps = 1 trunc 256 Registration Done begin writing rigid_.nii.gz writing rigid_ affine -------------------------------------------------------------------------------------- Applying rigid transformation to neck_cut.nii.gz -------------------------------------------------------------------------------------- AFFINE: rigid_Affine.txt moving_image_filename: neck_cut.nii.gz components 1 output_image_filename: rigid_deformed.nii.gz reference_image_filename: mni_icbm152_nlin_asym_09a/mni_icbm152_t1_tal_nlin_asym_09a.nii [0/1]: AFFINE: rigid_Affine.txt User Linear interpolation You are doing something more complex -- we wont check syntax in this case output origin: [98, 134, -72] output size: [197, 233, 189] output spacing: [1, 1, 1] output direction: -1 0 0 0 -1 0 0 0 1
# 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())
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())
pos = (slice(None),40,slice(None))
show_imgs([mni_img, renormed_img], pos, ['mni_img', 'renormed_img'])
# 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'])
(197, 233, 189, 1, 3) (197, 233, 189) (197, 233, 189)
# 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'])
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'])