Image alignment by point set registration

It is often the case that images are acquired with different modalities, pehaps different spatial resolutions on the same specimen and we want to register them together to have a consistent description across all available images. In this recipe we will see how to achieve that using point set registration between different images.

Get the complete Python source code: pointset_registration.py

To compute the affine transform between an image and the reference image, at least 3 pair of points must be identified. The matching pairs are usually obtained by selecting unique features in both images and measuring their coordinates. Here we will illustrate the method compute_affine_transform and then apply it on a real example.

Let’s consider 5 points given by their coordinates \((x, y)\) in the reference 2D space.

random.seed(13)
n = 5
ref_points = np.empty((n, 2))
for i in range(n):
    ref_points[i, 0] = random.randint(0, 10)
    ref_points[i, 1] = random.randint(0, 10)

For the sake of the example, let’s transform the 5 points by an isotropic scaling of \(s=1.4\), a rotation of \(\theta=22\) degrees and a translation of \((t_x, t_y)=(5, 3)\). The transformation is can be achieved by multiplying the augmented coordinates vector \((x, y, 1)^T\) by the affine transform matrix \(\mathbf{A}=\mathbf{T}.\mathbf{S}.\mathbf{R}\) obtained by the composition of the rotation, scaling and translation (the order of the translation is important, here it is applied after the rotation and scaling).

For this example we have chosen the following transformation:

\[\begin{split}\begin{pmatrix} x' \\ y' \\ 1 \end{pmatrix} =\mathbf{T}.\mathbf{S}.\mathbf{R}= \begin{pmatrix} 1 & 0 & t_x \\ 0 & 1 & t_y \\ 0 & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} s & 0 & 0 \\ 0 & s & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \\ \end{pmatrix} \begin{pmatrix} x \\ y \\ 1 \end{pmatrix}\end{split}\]

# compute the affine transform by composing R, S and T
s = 1.4
angle = 22 * pi / 180.
tx, ty = 5, 3
S = np.array([[s, 0.0, 0.0],
              [0.0, s, 0.0],
              [0.0, 0.0, 1.0]])
R = np.array([[cos(angle), -sin(angle), 0.0],
              [sin(angle), cos(angle), 0.0],
              [0.0, 0.0, 1.0]])
T = np.array([[1.0, 0.0, tx],
              [0.0, 1.0, ty],
              [0.0, 0.0, 1.0]])
A = np.dot(T, np.dot(S, R))
print('full affine transform:\n{:s}'.format(A))

Which is also equivalent to:

\[\begin{split}\begin{matrix} x' = t_x + x s \cos\theta - y s \sin\theta \\ y' = t_y + x s \sin\theta + y s \cos\theta \\ \end{matrix}\end{split}\]

In this case, one can verify that the full transformation is:

\[\begin{split}\mathbf{A} = \mathbf{T}.\mathbf{S}.\mathbf{R}= \begin{pmatrix} 1.298 & -0.524 & 5. \\ 0.524 & 1.298 & 3. \\ 0. & 0. & 1. \\ \end{pmatrix}\end{split}\]

With either solution, we can transform the points to obtain the new coordinates:

# transform the points
tsr_points = np.empty_like(ref_points)
for i in range(n):
    p = np.dot(A, [ref_points[i, 0], ref_points[i, 1], 1])
    tsr_points[i, 0] = p[0]
    tsr_points[i, 1] = p[1]

With the two lists (or arrays) of points, one can compute the affine transform:

# compute the affine transform from the point set
translation, transformation = compute_affine_transform(ref_points, tsr_points)
invt = np.linalg.inv(transformation)
offset = -np.dot(invt, translation)

And register the points using the inverse affine transform and the offset computed from the centroids.

ref_centroid = np.mean(ref_points, axis=0)
tsr_centroid = np.mean(tsr_points, axis=0)
new_points = np.empty_like(ref_points)
for i in range(n):
    new_points[i] = ref_centroid + np.dot(transformation, tsr_points[i] - tsr_centroid)

The 5 transformed points registered by applying the inverse affine transform.