# - 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

# - show figures inside the notebook
%matplotlib inline

# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package

# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'

# - import numpy.linalg with a shorter name
import numpy.linalg as npl

Making floating points numbers print nicely¶

By default, when numpy prints an array, it looks for very small or very large numbers. If it finds either, it uses exponents to show the numbers. This can be annoying:

np.pi

3.141592653589793

np.array([np.pi, 0.000001])

array([  3.14159265e+00,   1.00000000e-06])

In order to avoid this, you can tell numpy not to use exponential notation for small numbers:

np.set_printoptions(suppress=True)
np.array([np.pi, 0.000001])

array([ 3.14159265,  0.000001  ])

This setting stays in place until you change it:

np.array([np.pi, 0.000001])

array([ 3.14159265,  0.000001  ])

It can also be annoying to see many digits after the decimal point, if you know that these are not important. You can set the number of digits after the decimal point for numpy printing like this:

np.set_printoptions(precision=4)
np.array([np.pi, 0.000001])

array([ 3.1416,  0.    ])

This only affects printing, not calculations:

np.cos(np.pi)

-1.0

np.cos(3.1416)

-0.99999999997301514

The image affine¶

So far we have not payed much attention to the image header.

An image consists of:

the array data;
metadata (data about the array data).

The header contains the metadata for the image.

One piece of metadata, is the image affine.

import nibabel as nib
img = nib.load('ds107_sub012_highres.nii')
img.affine

array([[   1.    ,    0.    ,    0.    , -127.    ],
       [   0.    ,    1.    ,    0.    ,  -83.3253],
       [   0.    ,    0.    ,    1.    ,  -90.0533],
       [   0.    ,    0.    ,    0.    ,    1.    ]])

As you can imagine, nibabel is getting the affine from the header:

print(img.header)

<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 256 208 192   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.  1.  1.  0.  0.  0.  0.]
vox_offset      : 0.0
scl_slope       : nan
scl_inter       : nan
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      : scanner
sform_code      : scanner
quatern_b       : 0.0
quatern_c       : 0.0
quatern_d       : 0.0
qoffset_x       : -127.0
qoffset_y       : -83.325302124
qoffset_z       : -90.0532836914
srow_x          : [   1.    0.    0. -127.]
srow_y          : [  0.       1.       0.     -83.3253]
srow_z          : [  0.       0.       1.     -90.0533]
intent_name     : 
magic           : n+1

img.header.extensions

Nifti1Extensions()

Notice the srow_x, srow_y, srow_z fields in the header, that contain the affine for this image. It is not always this simple though - see http://nifti.nimh.nih.gov/nifti-1 for more details. In general, nibabel will take care of this for you, by extracting the affine from the header, and returning it via img.affine.

Nifti images can also be `.img, .hdr` pairs¶

So far, all the images we have seen have been NIfTI format images, stored in a single file with a .nii extension. The single file contains the header information, and the image array data.

The NIfTI format also allows the image to be stored as two files, one with extension .img storing the image array data, and another with extension .hdr storing the header. These are called NIfTI pair images.

For example, there are two files in the same directory as this notebook, called ds114_sub009_highres_moved.img and ds114_sub009_highres_moved.hdr. These two files together form one NIfTI image. You can load these with nibabel in the usual way:

pair_img = nib.load('ds114_sub009_highres_moved.img')
pair_img.affine

array([[   0.9416,   -0.4311,   -0.0586,  -98.8336],
       [   0.336 ,    1.1887,    0.2264, -164.1377],
       [  -0.0215,   -0.3028,    0.9723, -158.4178],
       [   0.    ,    0.    ,    0.    ,    1.    ]])

This form of the NIfTI image is getting less common, because it is inconvenient to have to keep the .img and .hdr files together, but you may still find them used. They have only one advantage, which is that, if some software wants to change only the header information, it only has to rewrite a small .hdr file, rather than the whole .nii file containing the image data and the header.

Identity matrices with `np.eye`¶

An identity matrix is a matrix with zeros off the diagonals, and ones on the leading diagonal. These are so common, than numpy (and MATLAB) have a function for making them. The parameter is the numbers of rows = number of columns:

np.eye(3)

array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

np.eye(6)

array([[ 1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1.]])

Affines, inverses¶

We often have the situation where we compose an affine of several transformations. We do the composing using matrix multiplication. For example, here we compose two rotations and a translation. Remember - matrix multiplication works right to left:

# - get our functions to make rotation matrices
from rotations import x_rotmat, y_rotmat, z_rotmat

A rotation matrix (3 x 3) for rotation -0.2 radians around the x axis:

first_rotation = x_rotmat(-0.2)
first_rotation

array([[ 1.    ,  0.    ,  0.    ],
       [ 0.    ,  0.9801,  0.1987],
       [ 0.    , -0.1987,  0.9801]])

We can make this rotation matrix into an affine transformation, by putting it into the top left of a 4 x 4 identity matrix:

first_affine = np.eye(4)  # The identity affine
first_affine[:3, :3] = first_rotation
first_affine

array([[ 1.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.9801,  0.1987,  0.    ],
       [ 0.    , -0.1987,  0.9801,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  1.    ]])

Now we made a second affine matrix for a rotation around y of 0.4 radians:

second_affine = np.eye(4)
second_affine[:3, :3] = y_rotmat(0.4)
second_affine

array([[ 0.9211,  0.    ,  0.3894,  0.    ],
       [ 0.    ,  1.    ,  0.    ,  0.    ],
       [-0.3894,  0.    ,  0.9211,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  1.    ]])

Finally we make a translation of 10 in x, 20 in y and 30 in z:

third_affine = np.eye(4)
third_affine[:3, 3] = [10, 20, 30]
third_affine

array([[  1.,   0.,   0.,  10.],
       [  0.,   1.,   0.,  20.],
       [  0.,   0.,   1.,  30.],
       [  0.,   0.,   0.,   1.]])

We compose these three affine matrices to give an affine implementing first a rotation of -0.2 around the x axis, then a rotation of 0.4 around the y axis, and finally a translation [10, 20, 30] in [x, y, z]. Note the order though - matrix multiplication goes from right to left:

combined = third_affine.dot(second_affine.dot(first_affine))
combined

array([[  0.9211,  -0.0774,   0.3817,  10.    ],
       [  0.    ,   0.9801,   0.1987,  20.    ],
       [ -0.3894,  -0.183 ,   0.9027,  30.    ],
       [  0.    ,   0.    ,   0.    ,   1.    ]])

nibabel.affines¶

In fact, nibabel has a short-cut routine to make a 4x4 affine matrix from a 3x3 matrix and an (optional) vector of translations:

# Affine from a 3x3 matrix (the 'mat' in 'matvec')
nib.affines.from_matvec(y_rotmat(0.4))

array([[ 0.9211,  0.    ,  0.3894,  0.    ],
       [ 0.    ,  1.    ,  0.    ,  0.    ],
       [-0.3894,  0.    ,  0.9211,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  1.    ]])

# Affine from a 3x3 matrix ('mat') and a translation vector ('vec')
nib.affines.from_matvec(y_rotmat(0.4), [10, 20, 30])

array([[  0.9211,   0.    ,   0.3894,  10.    ],
       [  0.    ,   1.    ,   0.    ,  20.    ],
       [ -0.3894,   0.    ,   0.9211,  30.    ],
       [  0.    ,   0.    ,   0.    ,   1.    ]])

Manipulating affines with inverses¶

Let us say we have an affine, like the one we just made:

combined

array([[  0.9211,  -0.0774,   0.3817,  10.    ],
       [  0.    ,   0.9801,   0.1987,  20.    ],
       [ -0.3894,  -0.183 ,   0.9027,  30.    ],
       [  0.    ,   0.    ,   0.    ,   1.    ]])

Imagine that we knew that this affine was composed of three affines, and we knew the last two, but not the first. How would we find what the first affine was?

Call our combined affine $\mathbf{D}$. We know that $\mathbf{D} = \mathbf{C} \cdot \mathbf{B} \cdot \mathbf{A}$. We know $\mathbf{C}$ and $\mathbf{B}$ but we want to find $\mathbf{A}$.

Above I've written matrix multiplication with a dot - as in $\mathbf{B} \cdot \mathbf{A}$, but in what follows I'll omit the dot, just writing $\mathbf{B} \mathbf{A}$ to mean matrix multiplication.

We find $\mathbf{A}$ using matrix inverses. Call $\mathbf{E} = \mathbf{C} \mathbf{B}$. Then $\mathbf{D} = \mathbf{E} \mathbf{A}$. If we can find the inverse of $\mathbf{E}$ (written as $\mathbf{E^{-1}}$) then (by the definition of the inverse):

$$ \mathbf{E^{-1}} \mathbf{E} = \mathbf{I} $$

and

$$ \mathbf{E^{-1}} \mathbf{D} = \mathbf{E^{-1}} \mathbf{E} \mathbf{A} \\ \mathbf{E^{-1}} \mathbf{D} = \mathbf{I} \mathbf{A} \\ \mathbf{E^{-1}} \mathbf{D} = \mathbf{A} $$

For reasons we do not have time to go into, our affine matrices are almost invariably invertible.

Let's see if we can reconstruct our first_affine from the combined affine, given we know the third_affine and second_affine:

E = third_affine.dot(second_affine)
E_inv = npl.inv(E)
E_inv.dot(combined)

array([[ 1.    ,  0.    , -0.    ,  0.    ],
       [ 0.    ,  0.9801,  0.1987,  0.    ],
       [ 0.    , -0.1987,  0.9801,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  1.    ]])

This is the same as our first affine:

first_affine

array([[ 1.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.9801,  0.1987,  0.    ],
       [ 0.    , -0.1987,  0.9801,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  1.    ]])

What about the situation where we know the first part of the affine, but we want to find the rest?

To solve this problem, we will need the right inverse.

The inverse we have used so far is the left inverse - so called because we apply it multiplying on the left of the original matrix:

$$ \mathbf{E^{-1}} \mathbf{E} = \mathbf{I} $$

Luckily, it turns out that, for square matrices, if there is a left inverse $\mathbf{E^{-1}}$ then this is also the right inverse:

$$ \mathbf{E^{-1}} \mathbf{E} = \mathbf{E} \mathbf{E^{-1}} = \mathbf{I} $$

It is a little tricky to prove that there must be a right inverse, but it is easy to prove that, if there is a right inverse, it must be the same as the left inverse. Call the left inverse $\mathbf{L}$ and the right inverse $\mathbf{R}$:

$$ \mathbf{LA} = \mathbf{I}\\ \mathbf{AR} = \mathbf{I}\\ $$

then:

$$ \mathbf{LAR} = \mathbf{LAR}\\ \mathbf{L(AR)} = \mathbf{(LA)R}\\ \mathbf{L} = \mathbf{R} $$

So, in our case, where we want to find the transformations following the first affine, we can do this:

$$ \mathbf{D} = \mathbf{C} \mathbf{B} \mathbf{A} \\ \mathbf{D} \mathbf{A^{-1}} = \mathbf{C} \mathbf{B} \mathbf{A} \mathbf{A^{-1}} \\ \mathbf{D} \mathbf{A^{-1}} = \mathbf{C} \mathbf{B} $$

For our actual affines:

third_with_second = combined.dot(npl.inv(first_affine))
third_with_second

array([[  0.9211,  -0.    ,   0.3894,  10.    ],
       [  0.    ,   1.    ,   0.    ,  20.    ],
       [ -0.3894,  -0.    ,   0.9211,  30.    ],
       [  0.    ,   0.    ,   0.    ,   1.    ]])

# This is the same as
third_affine.dot(second_affine)

array([[  0.9211,   0.    ,   0.3894,  10.    ],
       [  0.    ,   1.    ,   0.    ,  20.    ],
       [ -0.3894,   0.    ,   0.9211,  30.    ],
       [  0.    ,   0.    ,   0.    ,   1.    ]])