# Images and affines

## Contents

# 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)
```