# - 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'
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):
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) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \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) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \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) & -\sin(\gamma) \\ 0 & \sin(\gamma) & \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.
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:
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 has a function for doing reampling with transformations, called scipy.ndimage.affine_transform
.
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}$:
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}$:
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:
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
(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:
K = affine_transform(J, M, translation, order=1) # order=1 for linear interpolation
K.shape
(64, 64, 35)
Here affine_transform
K
, assuming it will be the same shape as J
;K
:M
and translation
to $(x_i, y_i, z_i)$ to get the corresponding point in J
: $(x_j, y_j, z_j)$;J
at $(x_j, y_j, z_j)$ to get $v$;K