Images and affines#
See: coordinate systems and affine transforms for background.
# import common modules
import numpy as np
# print arrays to 4DP
np.set_printoptions(precision=4, suppress=True)
import numpy.linalg as npl
import nibabel as nib
Affines, inverses#
We often have the situation where we compose an affine of several transformations. We do the composing using matrix multiplication. For example, the following code composes two rotations and a translation.
Remember: Matrix multiplication works right to left.
Here is a rotation matrix (3 x 3) for a rotation of -0.2 radians around the x axis:
def x_rotmat(theta):
""" Rotation matrix for rotation of `theta` radians around x axis
Parameters
----------
theta : scalar
Rotation angle in radians
Returns
-------
M : shape (3, 3) array
Rotation matrix
"""
cos_t = np.cos(theta)
sin_t = np.sin(theta)
return np.array([[1, 0, 0],
[0, cos_t, -sin_t],
[0, sin_t, cos_t]])
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:
def y_rotmat(theta):
""" Rotation matrix for rotation of `theta` radians around y axis
Parameters
----------
theta : scalar
Rotation angle in radians
Returns
-------
M : shape (3, 3) array
Rotation matrix
"""
cos_t = np.cos(theta)
sin_t = np.sin(theta)
return np.array([[cos_t, 0, sin_t],
[0, 1, 0],
[-sin_t, 0, cos_t]])
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 — matrix multiplication goes from right to left:
combined = third_affine @ second_affine @ 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. ]])
See The nibabel.affines module for a module with useful functions for working with affine matrices.
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):
and:
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 @ second_affine
E_inv = npl.inv(E)
E_inv @ 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. ]])
assert np.allclose(E_inv @ combined, first_affine)
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:
Luckily, it turns out that, for square matrices, if there is a left inverse \(\mathbf{E^{-1}}\) then this is also the right inverse:
It is a bit out of our way to prove that a matrix with a left inverse must also have a right inverse. If you accept that on faith for now, 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}\):
then:
So, in our case, where we want to find the transformations following the first affine, we can do this:
For our actual affines:
third_with_second = combined @ 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
F = third_affine @ second_affine
F
array([[ 0.9211, 0. , 0.3894, 10. ],
[ 0. , 1. , 0. , 20. ],
[-0.3894, 0. , 0.9211, 30. ],
[ 0. , 0. , 0. , 1. ]])
assert np.allclose(third_with_second, F)