1. pna2015
  2. day9
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'

Rotations and rotation matrices

See wikipedia on rotation matrices

In two dimensions, rotating a vector $\theta$ around the origin can be expressed as a 2 by 2 transformation matrix:

$$ R(\beta) = \begin{bmatrix} \cos \beta & -\sin \beta \\ \sin \beta & \cos \beta \\ \end{bmatrix} $$

This rotates column vectors by matrix multiplication:

$$ \begin{bmatrix} x' \\ y' \\ \end{bmatrix} = \begin{bmatrix} \cos \beta & -\sin \beta \\ \sin \beta & \cos \beta \\ \end{bmatrix}\begin{bmatrix} x \\ y \\ \end{bmatrix} $$

So the coordinates $(x',y')$ of the point $(x,y)$ after rotation are

$$ x' = x \cos \beta - y \sin \beta \\ y' = x \sin \beta + y \cos \beta $$

You might be able to see why this is true by looking at this diagram, and remembering:

$$ \cos(\alpha) = x \\ \sin(\alpha) = y $$

(Image copyright Stackexchange user Blue licensed with cc-by-sa)

Rotations in three dimensions simply extend from 2D. For example, in 3D, the rotation above would be a rotation around the z axis (z stays constant, (x, y) change):

$$ \begin{alignat}{1} R_x(\theta) &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta \\[3pt] 0 & \sin \theta & \cos \theta \\[3pt] \end{bmatrix} \\[6pt] R_y(\theta) &= \begin{bmatrix} \cos \theta & 0 & \sin \theta \\[3pt] 0 & 1 & 0 \\[3pt] -\sin \theta & 0 & \cos \theta \\ \end{bmatrix} \\[6pt] R_z(\theta) &= \begin{bmatrix} \cos \theta & -\sin \theta & 0 \\[3pt] \sin \theta & \cos \theta & 0\\[3pt] 0 & 0 & 1\\ \end{bmatrix} \end{alignat} $$

We can combine rotations with matrix multiplication. For example, here is an rotation around the x axis of $\gamma$ radians:

$$ \begin{bmatrix} x'\\ y'\\ z'\\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\gamma) & 0 & -\sin(\gamma) \\ 0 & \sin(\gamma) & 0 & \cos(\gamma) \\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} $$

We could then apply a rotation around the y axis of $\phi$ radians:

$$ \begin{bmatrix} x''\\ y''\\ z''\\ \end{bmatrix} = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \begin{bmatrix} x'\\ j'\\ k'\\ \end{bmatrix} $$

But we could also write this as:

$$ \begin{bmatrix} x''\\ y''\\ z''\\ \end{bmatrix} = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\gamma) & 0 & -\sin(\gamma) \\ 0 & \sin(\gamma) & 0 & \cos(\gamma) \\ \end{bmatrix} \begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} $$

And, because matrix multiplication is associative:

$$ \mathbf{Q} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\gamma) & 0 & -\sin(\gamma) \\ 0 & \sin(\gamma) & 0 & \cos(\gamma) \\ \end{bmatrix} $$$$ \mathbf{P} = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \\ \end{bmatrix} $$$$ \mathbf{M} = \mathbf{P} \cdot \mathbf{Q} $$$$ \begin{bmatrix} x''\\ y''\\ z''\\ \end{bmatrix} = \mathbf{M} \begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} $$

Here $\mathbf{M}$ is our rotation matrix that encodes a rotation by $\gamma$ radians around the x axis followed by (right to left matrix multiplication) a rotation by $\phi$ radians around the y axis.

Have a look at the file rotations.py in pna_code. It has routines that will make these 3 by 3 matrices for given angles of rotation around the x, y, and z axes.

Reampling, pull, push

Let's say I have two images $\mathbf{I}$ and $\mathbf{J}$.

There is some spatial transform between them, such as a translation, or rotation.

We could either think of the transformation that maps $\mathbf{I} \to \mathbf{J}$ or $\mathbf{J} \to \mathbf{I}$.

The transformation maps voxel coordinates in one image to coordinates in the other.

For example, write a coordinate in $\mathbf{I}$ as $(x_i, y_i, z_i)$, and a coordinate in $\mathbf{J}$ as $(x_j, y_j, z_j)$.

The mapping $\mathbf{I} \to \mathbf{J}$ maps $(x_i, y_i, z_i) \to (x_j, y_j, z_j)$.

Now let's say that I want to move image $\mathbf{J}$ to match image $\mathbf{I}$.

To do this moving, I need to resample $\mathbf{J}$ onto the same voxel grid as $\mathbf{I}$.

Specifically, I am going to do the following:

  • make a new empty image $\mathbf{K}$ that has the same voxel grid as $\mathbf{I}$;
  • for each coordinate in $\mathbf{I} : (x_i, y_i, z_i)$ I will apply the transform $\mathbf{I} \to \mathbf{J}$ to get $(x_j, y_j, z_j)$;
  • I will probably need to estimate a value $v$ from $\mathbf{J}$ at $(x_j, y_j, z_j)$ because the coordinate values $(x_j, y_j, z_j)$ will probably not be integers, and so there is no exactly matching value in $\mathbf{J}$;
  • I place $v$ into $\mathbf{K}$ at coordinate $(x_i, y_i, z_i)$.

Notice that, in order to move $\mathbf{J}$ to match $\mathbf{I}$ I needed the opposite transform - that is $\mathbf{I} \to \mathbf{J}$. This is called pull resampling.

I will call the $\mathbf{I} \to \mathbf{J}$ transform the resampling transform.

Scipy ndimage and affine_transform

Scipy has a function for doing reampling with transformations, called scipy.ndimage.affine_transform.

In [5]:
from scipy.ndimage import affine_transform

It does all the heavy work of resampling for us.

For example, lets say I have an image volume $\mathbf{I}$:

In [6]:
import nibabel as nib
img = nib.load('ds107_sub012_t1r2.nii')
data = img.get_data()
I = data[..., 0]  # I is the first volume

And another image volume $\mathbf{J}$:

In [7]:
J = data[..., 1]  # I is the second volume

Let's say I know that the resampling transformation $\mathbf{I} \to \mathbf{J}$ is given by:

  • A rotation by 0.2 radians about the x axis
  • A translation of [1, 2, 3] voxels.
In [8]:
from rotations import x_rotmat  # rotations from `pna_code`
M = x_rotmat(0.2)  # rotation matrix for rotation of 0.2 radians around x axis
translation = [1, 2, 3]  # Translation from I to J
M, translation
Out[8]:
(array([[ 1.        ,  0.        ,  0.        ],
        [ 0.        ,  0.98006658, -0.19866933],
        [ 0.        ,  0.19866933,  0.98006658]]), [1, 2, 3])

affine_transform now does the work for me:

In [9]:
K = affine_transform(J, M, translation, order=1)  # order=1 for linear interpolation
K.shape
Out[9]:
(64, 64, 35)

Here affine_transform

  • makes the new empty volume K, assuming it will be the same shape as J;
  • for each coordinate $(x_i, y_i, z_i)$ implied by the volume K:
    • applies the transformations implied by M and translation to $(x_i, y_i, z_i)$ to get the corresponding point in J : $(x_j, y_j, z_j)$;
    • resamples J at $(x_j, y_j, z_j)$ to get $v$;
    • places $v$ at $(x_i, y_i, z_i)$ in K