\(\newcommand{L}[1]{\| #1 \|}\newcommand{VL}[1]{\L{ \vec{#1} }}\newcommand{R}[1]{\operatorname{Re}\,(#1)}\newcommand{I}[1]{\operatorname{Im}\, (#1)}\)

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. We are using the rotations.py module to make the 3D rotation matrices. Remember – matrix multiplication works right to left.

>>> # get functions to make 3D rotation matrices
>>> from rotations import x_rotmat, y_rotmat, z_rotmat

Here is a rotation matrix (3 x 3) for a rotation of -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 – 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.    ]])

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):

\[\mathbf{E^{-1}} \mathbf{E} = \mathbf{I}\]

and:

\[\begin{split}\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}\end{split}\]

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)  # doctest: +SKIP
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.    ]])
>>> np.allclose(E_inv.dot(combined), first_affine)
True

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 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}\):

\[\begin{split}\mathbf{LA} = \mathbf{I}\\ \mathbf{AR} = \mathbf{I}\\\end{split}\]

then:

\[\begin{split}\mathbf{LAR} = \mathbf{LAR}\\ \mathbf{L(AR)} = \mathbf{(LA)R}\\ \mathbf{L} = \mathbf{R}\end{split}\]

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

\[\begin{split}\mathbf{F} \triangleq \mathbf{C} \mathbf{B} \\ \mathbf{D} = \mathbf{F} \mathbf{A} \\ \mathbf{D} \mathbf{A^{-1}} = \mathbf{F} \mathbf{A} \mathbf{A^{-1}} \\ \mathbf{D} \mathbf{A^{-1}} = \mathbf{F}\end{split}\]

For our actual affines:

>>> third_with_second = combined.dot(npl.inv(first_affine))
>>> third_with_second  # doctest: +SKIP
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.dot(second_affine)
>>> F  # doctest: +SKIP
array([[  0.9211,   0.    ,   0.3894,  10.    ],
       [  0.    ,   1.    ,   0.    ,  20.    ],
       [ -0.3894,   0.    ,   0.9211,  30.    ],
       [  0.    ,   0.    ,   0.    ,   1.    ]])
>>> np.allclose(third_with_second, F)
True