\(\newcommand{L}[1]{\| #1 \|}\newcommand{VL}[1]{\L{ \vec{#1} }}\newcommand{R}[1]{\operatorname{Re}\,(#1)}\newcommand{I}[1]{\operatorname{Im}\, (#1)}\)
Resampling with scipy.ndimage¶
>>> # import common modules
>>> import numpy as np
>>> np.set_printoptions(precision=4) # print arrays to 4DP
>>> import matplotlib.pyplot as plt
>>> # gray colormap and nearest neighbor interpolation by default
>>> plt.rcParams['image.cmap'] = 'gray'
>>> plt.rcParams['image.interpolation'] = 'nearest'
Reampling, pull, push¶
Let us say we 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 transformations map 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 us say that we want to move image \(\mathbf{J}\) to match image \(\mathbf{I}\).
To do this moving, we need to resample \(\mathbf{J}\) onto the same voxel grid as \(\mathbf{I}\).
Specifically, we are 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)\) we will apply the transform \(\mathbf{I} \to \mathbf{J}\) to get \((x_j, y_j, z_j)\);
- we 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}\);
- we 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}\) we needed the opposite transform – that is: \(\mathbf{I} \to \mathbf{J}\). This is called pull resampling.
We 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.
>>> from scipy.ndimage import affine_transform
It does all the heavy work of resampling for us.
For example, lets say we have an image volume \(\mathbf{I}\)
(ds107_sub012_t1r2.nii
):
>>> import nibabel as nib
>>> img = nib.load('ds107_sub012_t1r2.nii')
>>> data = img.get_data()
>>> I = data[..., 0] # I is the first volume
We have another image volume \(\mathbf{J}\):
>>> J = data[..., 1] # I is the second volume
Let’s say we 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.
See Rotations and rotation matrices for more on 2D and 3D rotation matrices.
Download the file rotations.py
. It has routines that will make
3 by 3 rotation matrices for rotations by given angles around the x, y, and z
axes.
Of course you will want to test these functions. Download
test_rotations.py
to the same directory as rotations.py
and
run the following from your terminal:
py.test test_rotations.py
We use the routines in rotation.py
to make the rotation matrix we need:
>>> from rotations import x_rotmat # from rotations.py
>>> # rotation matrix for rotation of 0.2 radians around x axis
>>> M = x_rotmat(0.2)
>>> M
array([[ 1. , 0. , 0. ],
[ 0. , 0.9801, -0.1987],
[ 0. , 0.1987, 0.9801]])
>>> translation = [1, 2, 3] # Translation from I to J
>>> translation
[1, 2, 3]
The affine_transform
function does the work of resampling. By default, it
implements the following algorithm:
makes the new empty volume
K
, assuming it will be the same shape asJ
;for each coordinate \((x_i, y_i, z_i)\) implied by the volume
K
:- apply the transformations implied by
M
andtranslation
to \((x_i, y_i, z_i)\) to get the corresponding point inJ
: \((x_j, y_j, z_j)\); - resample
J
at \((x_j, y_j, z_j)\) to get \(v\); - place \(v\) at \((x_i, y_i, z_i)\) in
K
- apply the transformations implied by
>>> # order=1 for linear interpolation
>>> K = affine_transform(J, M, translation, order=1)
>>> K.shape
(64, 64, 35)
>>> plt.imshow(K[:, :, 17])
<...>
Resampling with images of different shapes¶
Notice the assumption that affine_transform
makes above – that the
output image will be the same shape as the input image.
This need not be the case. In fact we can tell affine_transform
to start
with an empty volume K
with another shape. To do this, we use the
output_shape
parameter.
Now we can be more precise about the algorithm of affine_transform
.
affine_transform
accepts:
input
– an array to resample from. Say this array asn
dimensions (len(input.shape)
);matrix
– ann
byn
transformation matrix (the top leftmat
part of an affine);offset
– a optional lengthn
translation vector to be applied after the thematrix
transformation (thevec
part of an affine);output_shape
: an optional tuple giving the shape of the output array into which we will put the values resampled frominput
.output_shape
defaults toinput.shape
; this is the default we have been using above.
affine_transform
then generates all the voxel coordinates implied by the
output_shape
, and transforms them with the matrix
and offset
transforms to get a new set of coordinates C
. It then samples the
input
array at the coordinates given by C
to generate the output
array.
Here is the call to affine_transform
that we used above, where we allowed
the routine to assume that the output shape was the same as the input shape:
>>> K = affine_transform(J, M, translation, order=1)
>>> K.shape
(64, 64, 35)
The is the same as the following call, where we specify the shape explicitly:
>>> K = affine_transform(J, M, translation, output_shape=J.shape, order=1)
>>> K.shape
(64, 64, 35)
The output shape can be different from the input shape:
>>> K = affine_transform(J, M, translation,
... output_shape=(65, 65, 36), order=1)
>>> K.shape
(65, 65, 36)
Remember that the M
matrix and translation
vector apply to the
coordinates implied by the output shape.