1. pna2015
  2. day11
In [1]:
# - compatibility with Python 3
from __future__ import print_function  # print('me') instead of print 'me'
from __future__ import division  # 1/2 == 0.5, not 0
In [2]:
# - show figures inside the notebook
%matplotlib inline
In [3]:
# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
In [4]:
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
In [5]:
# - import numpy.linalg with a shorter name
import numpy.linalg as npl
In [6]:
# - make numpy print out values in less verbose form
np.set_printoptions(suppress=True, precision=4)

Applying a deformation field

We should have already calculated the deformation field for the image ds107_sub012_highres.nii.

We did this using the SPM12 Normalise: Estimate option from the GUI.

This should have left an image called y_ds107_sub012_highres.nii in the directory containing this notebook.

If not:

1) Ask for help; 2) If you get stuck, you can download the image from https://nipy.bic.berkeley.edu/pna/results/y_ds107_sub012_highres.nii

We make an image object for this deformations image by loading with nibabel.

In [7]:
# - Load y_ds107_sub012_highres.nii with nibabel
# Get the image array data
import nibabel as nib
deformations_img = nib.load('y_ds107_sub012_highres.nii')
deformations_data = deformations_img.get_data()
deformations_data.shape
Out[7]:
(121, 145, 121, 1, 3)

We are going to work out how to apply this deformations image to reslice our original image ds107_sub012_highres.nii.

Oddly - this is a 5 dimensional image, where the 4th dimension is length 1.

The length-1 4th dimension is an artefact of the NIfTI image format - so let's get rid of this dimension with np.squeeze:

In [8]:
# - Remove length-1 4th dimension from deformation data
deformations_data = np.squeeze(deformations_data)
deformations_data.shape
Out[8]:
(121, 145, 121, 3)

The data is now a 4-dimensional image, containing 3 volumes. These volumes are:

  • x coordinates;
  • y coordinates;
  • z coordinates

respectively.

Put another way, the vector deformations_data[i, j, k, :] gives the (x, y, z) coordinates for the voxel [i, j, k]. More on this later.

If you were looking carefully on the SPM interface, SPM has calculated the distortions necessary to go from a template of tissue probability maps (called TPM.nii) to the ds107_sub012_highres.nii image.

We can get the original 3D shape and affine of TPM.nii because SPM stored them in y_ds107_sub012_highres.nii:

In [9]:
# - get original TPM.nii 3D shape and affine
tpm_shape = deformations_data.shape[:3]
tpm_affine = deformations_img.affine
tpm_affine
Out[9]:
array([[  -1.5,    0. ,    0. ,   90. ],
       [   0. ,    1.5,    0. , -126. ],
       [   0. ,    0. ,    1.5,  -72. ],
       [   0. ,    0. ,    0. ,    1. ]])

First we look at the images before the normalization has been applied.

To do that, we will make a new copy of the MNI template, with the same shape as the TPM image.

The MNI template image we will use is mni_icbm152_t1_tal_nlin_asym_09a.nii

In [10]:
# - load the template image we will resample from
template_img = nib.load('mni_icbm152_t1_tal_nlin_asym_09a.nii')
template_data = template_img.get_data()
template_data.shape
Out[10]:
(197, 233, 189)

To do this, we need the mapping from voxels in TPM.nii to voxels in mni_icbm152_t1_tal_nlin_asym_09a.nii.

We can break this down into two transforms:

  • from voxels in TPM.nii to mm;
  • from mm to voxels in mni_icbm152_t1_tal_nlin_asym_09a.nii.

We have these already. So the full transform is:

In [11]:
# - voxels in TPM.nii to voxels in mni_icbm152_t1_tal_nlin_asym_09a.nii
# Matrix multiplication is right to left
vox2vox = npl.inv(template_img.affine).dot(tpm_affine)
vox2vox
Out[11]:
array([[  -1.5,    0. ,    0. ,  188. ],
       [   0. ,    1.5,    0. ,    8. ],
       [   0. ,    0. ,    1.5,    0. ],
       [   0. ,    0. ,    0. ,    1. ]])

We break this down into 3x3 mat and length 3 vec components:

In [12]:
# - to mat and vec
mat, vec = nib.affines.to_matvec(vox2vox)
mat, vec
Out[12]:
(array([[-1.5,  0. ,  0. ],
        [ 0. ,  1.5,  0. ],
        [ 0. ,  0. ,  1.5]]), array([ 188.,    8.,    0.]))

Then we resample from the MNI template, into the voxel grid of the TPM.nii:

In [13]:
# - resample MNI template onto TPM grid
from scipy.ndimage import affine_transform
template_into_tpm = affine_transform(template_data, mat, vec, output_shape=tpm_shape)
template_into_tpm.shape
Out[13]:
(121, 145, 121)

Here is the new version of the template image:

In [14]:
# - plot the template image resampled onto the TPM grid
plt.imshow(template_into_tpm[:, :, 60])
Out[14]:
<matplotlib.image.AxesImage at 0x1092fdb90>

Now, what do with with the SPM distortion field in deformations_data?

By checking in the SPM source code, it is possible work out that deformations_data contains, for every voxel in TPM, the corresponding mm coordinate in the mm space of the subject image.

That is, deformations_data[i, j, k] is a length 3 vector [x, y, z] where [x, y, z] are the mm coordinates of voxel [i, j, k] when mapped into millimeters for the subject image.

Here is the subject data, and the image, which contains the affine:

In [15]:
# - load the subject data that we will resample from
subject_img = nib.load('ds107_sub012_highres.nii')
subject_data = subject_img.get_data()
subject_data.shape
Out[15]:
(256, 208, 192)

With this information, and with the scipy.ndimage.map_coordinates function, you should be able to:

  • get the mapping from voxels in TPM to voxels in the subject image and
  • resample the subject image into the grid of the TPM image using this mapping.

Hint: remember that map_coordinates expects the 3-length coordinate dimension to be first, but deformations_data - at the moment - has the 3-length coordinate dimension last.

In [16]:
# * get the mapping from voxels in TPM to voxels in the subject image and
# * resample the subject image into the grid of the TPM image using this mapping.
from scipy.ndimage import map_coordinates
vox2vox_mapping = nib.affines.apply_affine(npl.inv(subject_img.affine), deformations_data)
for_map_coords = vox2vox_mapping.transpose(3, 0, 1, 2)
subject_into_tpm = map_coordinates(subject_data, for_map_coords)
subject_into_tpm.shape
Out[16]:
(121, 145, 121)

Show an example slice from the template resampled into TPM and the subject resampled into TPM

In [17]:
# - show an example slice from the resampled template and resampled subject
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
axes[0].imshow(template_into_tpm[:, :, 60])
axes[1].imshow(subject_into_tpm[:, :, 60])
Out[17]:
<matplotlib.image.AxesImage at 0x10a0c6a50>