Linear transformations in Numpy

geometry geometric-transformations python numpy matplotlib

A linear transformation of the plane $\mathbb R^2$ is a geometric transformation of the form

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}a&b\\c&d\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

where $a$, $b$, $c$ and $d$ are real constants.

Linear transformations leave the origin fixed and preserve parallelism. Scaling, shearing, rotation and reflexion of a plane are examples of linear transformations.

Applying a geometric transformation to a given matrix in Numpy requires applying the inverse of the transformation to the coordinates of the matrix, create a new matrix of indices from the coordinates and map the matrix to the new indices. Since this can be tricky, let's start with a simple example involving a matrix that represents the indices itself.

A simple example. Indices transformation

>>> import numpy as np
>>> M, N = 3, 4
>>> matrix = np.arange(M*N).reshape((M,N))
>>> matrix
array([[ 0,  1,  2,  3],
[ 4,  5,  6,  7],
[ 8,  9, 10, 11]])

Then, we need to obtain the indices pairs of the matrix in a matrix form. The new indices of the matrix will result from the product of the inverse of the transformation matrix and this matrix, therefore the indices pairs in this case need to be a 2x12 matrix as

$\mathbf p = \begin{pmatrix} 0 & 0 & 0 & 1 & 1 & 1 & 2 & 2 & 2 & 3 & 3 & 3 \\ 0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2 & 0 & 1 & 2 \end{pmatrix},$

then:

>>> points = np.mgrid[0:N, 0:M].reshape((2, M*N))
>>> points
array([[0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3],
[0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2]])

An alternative way to get the indices pairs, longer but maybe more obvious:

>>> x, y = np.mgrid[0:N, 0:M]
>>> points = np.vstack([x.ravel(), y.ravel()])

Now, apply the transformation to the indices pairs. The new indices pairs need to be integers to map the given matrix to the indices. In this example the result will be just casted to integers for simplicity, then it will be easy to spot the result of the transformation since it involve points halfway between two integers. Other methods require interpolations of the given matrix from the indices. As an example, the transformation matrix will be

$\mathbf A = \begin{pmatrix}2&0\\0&1\end{pmatrix},$

which corresponds with scaling the plane in the $x$-axis by a factor of $2$. Hence, the new indices pairs are

$\mathbf p' = \lfloor \mathbf A \mathbf p \rfloor = \begin{pmatrix} 0&0&0&0&0&0&1&1&1&1&1&1\\ 0&1&2&0&1&2&0&1&2&0&1&2 \end{pmatrix},$

then:

>>> a = np.array([[2, 0],
[0, 1]])
>>> new_points = np.linalg.inv(a).dot(points).astype(int)
>>> new_points
array([[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
[0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2]])

To finish this example, convert the indices pairs to a matrix of indices which in this example corresponds with the resulting matrix.

The $x$ components is

$\mathbf x = \begin{pmatrix} 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{pmatrix},$

the $y$ component is

$\mathbf y = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \\ 2 & 2 & 2 & 2 \end{pmatrix},$

and the resulting matrix from the $x$ and $y$ components is

$\mathbf x + N \mathbf y = \begin{pmatrix} 0 & 0 & 1 & 1 \\ 4 & 4 & 5 & 5 \\ 8 & 8 & 9 & 9 \end{pmatrix}.$

Here it is easy to spot the scaling of the values along the $x$-axis in the matrix. These final steps correspond to:

>>> x, y = new_points.reshape((2, M, N), order='F')
>>> x + N*y
array([[0, 0, 1, 1],
[4, 4, 5, 5],
[8, 8, 9, 9]])

A more visual example. Matrix transformation

In the following example we will use a bigger matrix, represented as an image for visual support. Once we calculate the new indices matrix we will map the original matrix to the new indices, wrapping the out-of-bounds indices to obtain a continuous plane using numpy.take with mode='wrap'.

import matplotlib as mpl
import matplotlib.pyplot as plt

Some visual settings:

mpl.rcParams.update({'image.cmap': 'Accent',
'image.interpolation': 'none',
'xtick.major.width': 0,
'xtick.labelsize': 0,
'ytick.major.width': 0,
'ytick.labelsize': 0,
'axes.linewidth': 0})

Create a 200x200 matrix for this example:

aux = np.ones((100, 100), dtype=int)
src = np.vstack([np.c_[aux, 2*aux], np.c_[3*aux, 4*aux]])
plt.imshow(src)
plt.show()

The linear transformation function, which includes the operations of the previous examples but rounding the new indices pairs and mapping the source matrix to the new indices might be written as follows:

def linear_transformation(src, a):
M, N = src.shape
points = np.mgrid[0:N, 0:M].reshape((2, M*N))
new_points = np.linalg.inv(a).dot(points).round().astype(int)
x, y = new_points.reshape((2, M, N), order='F')
indices = x + N*y
return np.take(src, indices, mode='wrap')

Let's see some linear transformations we can do.

Scaling the plane in the $x$-axis by a factor of 1.5:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1.5 & 0\\0 & 0\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

a = np.array([[1.5, 0],
[0, 1]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()

Dilating the plane by a factor of 1.8:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1.8 & 0\\0 & 1.8\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

a = 1.8*np.eye(2)
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()

Dilating the plane by a factor of 0.5:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}0.5 & 0\\0 & 0.5\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

a = .5*np.eye(2)
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()

Scaling the plane in the $y$-axis by a factor of 0.5:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1 & 0\\0 & 0.5\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

a = np.array([[1, 0],
[0, .5]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()

Shearing about the $y$-axis with a vertical displacement of $+x/2$:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix}1 & 0\\{1 \over 2} & 0\end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

a = np.array([[1, 0],
[.5, 1]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()

Rotation through $45^{\circ}$ about the origin:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix} \cos{\pi \over 4} & -\sin{\pi \over 4} \\ \sin{\pi \over 4} & \cos{\pi \over 4} \end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

alpha = np.pi/4
a = np.array([[np.cos(alpha), -np.sin(alpha)],
[np.sin(alpha), np.cos(alpha)]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()

Reflexion in a line with inclination of $45^{\circ}$ through the origin:

$f \begin{pmatrix}x\\y\end{pmatrix} = \begin{pmatrix} \cos{\pi \over 2} & \sin{\pi \over 2} \\ \sin{\pi \over 2} & -\cos{\pi \over 2} \end{pmatrix} \begin{pmatrix}x\\y\end{pmatrix},$

in Numpy:

alpha = np.pi/4
a = np.array([[np.cos(2*alpha), np.sin(2*alpha)],
[np.sin(2*alpha), -np.cos(2*alpha)]])
dst = linear_transformation(src, a)
plt.imshow(dst)
plt.show()