Herve Lombaert

Research Projects Algorithms School Presentations Schedule

Manual Registration with Thin Plates

This page tries to highlight the use of thin plate splines in a registration framework. Given two images, the goal is to deform an image so it matches the second one. Here a manual registration is done, that is the user manually provides a few matching points on the two images. First, thin plate splines will be explained: what they are, how they can be used to deform an image. At last, they are used in a typical registration framework where a few matching points are provided.

Thin Plates

There are quite a few techniques out there that provide a smooth interpolation between a set of control points. Thin Plate Splines is one such technique. It interpolates a surface that passes through each control point. A set of 3 points will thus generate a flat plane. It is easy to think of the control points as position constraints on a bending surface. The ideal surface is one that bends the least. Here is an example of such a surface with 7 control points. The surface is forced to pass through all these 7 control points.

Figure 1: A Thin Plate Splines that passes through a set of control points (Image from the Jarno Elonen demo)

This least bent surface is given by the following equation:

$$a_1 + a_2 x + a_3 y + \sum_{i=1}^{n}{w_i U\left(\left\vert P_i - (x,y) \right\vert\right)}$$ f(x,y) = (1)

The first three terms correspond to the linear part which defines a flat plane that best matches all control points (this can be seen as a least square fitting). The last term corresponds to the bending forces provided by n control points. There is a coefficient w_i for each control point. Also, $\left\vert P_i- (x,y)\right\vert$ is the distance between the control point P_i and a position (x,y). This distance is used in the function U defined by U(r) = r² log r .

So far, what is unknown are the coefficients a₁, a₂, a₃, and w_i for every control point. All w_i forms the vector W. These unknowns are defined by:

$$\begin{eqnarray*} L^{-1} V & = & \left(W \vert a_1 a_2 a_3\right)^T \end{eqnarray*}$$

What is known are the control point positions (x_i, y_i) of all points P_i and their heights v_i. Let us simplify the writing with:

$$\begin{eqnarray*} P & = & \left[ \begin{array}{ccc}1 & x_1 & y_1\\ 1 & x_2 & y_2... ...ray}\right], \textrm{the control point heights, padded with 0s.} \end{eqnarray*}$$

Let us use a matrix K that evaluates the function U(r_ij) where r_ij is the distance between two control points, $r_{ij} = \left\vert P_i-P_j\right\vert$. Matrix K is:

$$\begin{eqnarray*} K & = & \left[ \begin{array}{ccc}U(r_{11}) & U(r_{12}) & ...\\... ...21}) & U(r_{22}) & ...\\ ... & ... & U(r_{nn})\end{array}\right] \end{eqnarray*}$$

The matrix L mentioned earlier is composed with the matrix K on its top-left corner, with the matrix P on its right side, with the matrix P^T on its bottom side, and with zeros on its bottom-right corner:

$$\begin{eqnarray*} L & = & \left[ \begin{array}{cc}K & P\\ P^T & 0\end{array}\right] \end{eqnarray*}$$

To find the unknown coefficients, that is the matrix (W | a₁ a₂ a₃), we can either find the inverse L^-1, or solve L (W | a₁ a₂ a₃) = Y which is a linear system. Since L is by definition symmetric, the last system can be easily solved with for instance a LU decomposition.

Once (W | a₁ a₂ a₃) is computed, we can go back to the thin plate equation 1 and find the height v = f(x,y) for any point (x,y).

So far, interpolation of a surface height has been done at a 2D point, but the math can be extended to N-D by adding more terms a_i and by using a larger vector P. The result would still be the interpolation of the height of the hypersurface, that means a 1D variable.

It is interesting to note that the number of terms in the thin plate equation depends on the number of control points. The more control points we have, the longer it will take to interpolate the surface height.

Deformation

To deform an image, the displacement of every pixel is needed. That is, we need to know where a pixel (x,y) on image A is going to be on image B at position (x+dx, y+dy). For instance, in the image below, given the position of the smile corner (x,y), we need to know its displacement (dx, dy) on the next image. However, our interpolation method only gives one variable, the surface height.

Figure 2: Thin plates can be used to deform image (a) to (b)

(a) before	(b) after

It is safe enough to consider the coordinate x independent of y. One surface can then be used to interpolate the displacement dx and another surface can be used to interpolate the displacement dy. More hypersurfaces can be used as the dimension of the image increases. In our smile image we define 6 points: 4 in the image corners, 1 in an eye, and 1 in the smile corner. At the corners, there is no displacement, the surfaces of dx and dy have no height. The smile corner moved upward, so the surface dy should have a significant height while the surface dx should have close to no height. The eye position moved left, so the surface dx should have a negative height while the height of surface dy should be close to 0. This gives the following two surfaces dx and dy:

Figure 3: Interpolation of displacement in the (a) x direction and in the (b) y direction

(a) $\Delta x$	(b) $\Delta y$

With these two surfaces dx and dy it is possible to interpolate the displacement (dx, dy) at any point. For example, the height at the nose position of surface dx is fairly 0 while the height of surface dy is important. That means the nose should have no horizontal displacement but should be moved upward. Interpolating for every pixel of image B gives the following image:

Figure 4: The displacement of the previous control points creates a smile in the original image

The two surfaces dx and dy are the keys that deform image A into image B. Seeking the solution to the registration of image A and image B means finding these two surfaces.

Matching Points

We know now the key to the registration process is two thin plate surfaces. To find these surfaces, we need some control points. And to find these control points we need the displacement of a point from one image to another. That is to say we need to find matching points in the two images.

There are also a lot of method to find matching points. We have chosen the fairly simple least square method. With an area W around a point (x₀, y₀) on image A, we try to find a very similar area on image B.

$$\begin{eqnarray*} p(x,y) & = & \exp\left(- \sum_{i,j \in W} \left(A(x_0+i,y_0+j) - B(x+i,y+j)\right)^2 / {2 \sigma} \right) \end{eqnarray*}$$

For that we slide a window on image B and do a squared difference between the small area on image A and the sliding window on image B. The least squared difference will give the matching position on image B.

Figure 5: (a) area around a control point in the first image and (b) sliding window in the second image

(a) image A	(b) image B

For instance, on the left image below, a point is selected on the top of a building. A sliding window gives the potential matching score p(x,y) of the current window on image B. Image B is in the middle, and the score map is on the right. Top of buildings usually have high score, but the maximum score gives where our building has moved.

Figure 6: Image A on the left with the original control point position in red. Image B in the middle with the matched position in green. Probability map of possible matches in image B.

Manual Registration

As we only seek to highlight the use of thin plates, the registration framework is here very simple. The user defines a couple of points on an image, and matching points on the second image provides control points for the thin plate surfaces.

For example, with the deformed image of New York City below, we would like to know the deformation required to match the correct image of the city's skyline.

Figure 7: (a) image to deform so it matches the target image (b)

(a) image to deform	(b) target image

To find the deformation, we select a couple of points that we find interesting on image A on the left. Matching points are found in the image in the middle. And image A is deformed with the thin plate surface. The result is shown on the right.

Figure 8: Image A on the left, image B in the middle, deformation of image A on the right

Here is the deformation of image A so it looks like image B.

Figure 9: Deformation of the original image so it matches the target image

To better appreciate the registered image A, below are the difference of the original image A and image B, and the difference with the registered image A. Initially where the deformation where strong in the middle, differences are high, and with the registered image A, there is much less difference in the same area.

Figure 10: Difference (registration error) of the target image with (a) the original image and with (b) the deformed image

(a) before	(b) after

Thin Plate Splines or Surfaces provides a very handy interpolator for image deformation. Recovering a deformation is also very handy with thin plates provided that we have the key control points. An automatic scheme for registration will need to automatically find the meaningful control points of the deformation.

Bibliography

1

Principal Warps: Thin-Plate Splines and the Decomposition of Deformations, F.L. Bookstein, PAMI 1989 (Vol. 11, No. 6) pp. 567-585

2

Thin Plate Spline editor, Jarno Elonen, http://elonen.iki.fi/code/tpsdemo