Matching as an Optimization Problem

The following method of registration of point sets is part of many publications so only a short summery is given here. The complete algorithm was published 1992 first and can be found, e.g., in[2]. The method is called Iterative Closest Points (ICP) algorithm.

Given two independently derived set of 3D points, $M$ (model set, $\vert M\vert = N_m$) and $D$ (data set, $\vert D\vert = N_d$), which correspond to a single shape, we want to find the transformation consisting of a rotation $\M R$ and a translation $\V t$ which minimizes the following cost function:

$$E(\M R, \V t) = \sum_{i=1}^{N_m}\sum_{j=1}^{N_d}w_{i,j}\norm {\V m_{i}-(\M R \V d_j+\V t)}^2.$$ (1)

The value of 1 is assigned to $w_{i,j}$ if the $i$-th point of set $M$ describes the same point in space as the $j$-th point of set $D$. Otherwise $w_{i,j}$ is set to 0. Two things have to be calculated: First the corresponding points and second the transformation $\M R$ and $\V t$ that minimizes $E(\M R, \V t)$ using the point correspondeces.

The ICP calculates iteratively the point correspondences. In each iteration step the algorithm selects the closest points as correspondences and calculates the transformation $(\M R, \V t)$ for Eq. (1). It is shown that the iteration terminates in a (local) minima[2]. The assumption is that in the last iteration step the point correspondences are correct.

In each iteration the transformation is calculated by the quaternion based method of Horn[5]. A unit quaternion is a 4 vector $\dot q = (q_0,q_x,q_y,_z)^T$, where $q_0 \geq 0, q_0^2+q_x^2+q_y^2+q_z^2 = 1$. It describes a rotation axis and an angle to rotate around that axis. A $3 \times 3$ rotation matrix $\M R$ is calculated from the unit quaternion according the the following scheme:

$$\begin{eqnarray*} \M R = \left( \begin{array}{ccc} (q_0^2 + q_x^2 - q_y^2 - q_z^... ... + q_xq_0) & (q_0^2 - q_x^2 - q_y^2 + q_z^2) \end{array}\right). \end{eqnarray*}$$

To determine the transformation, the mean values (centroid vectors) $\V c_m$ and $\V c_d$ of the points that contribute to the matching are subtracted from all points in $M$ and $D$ respectively, resulting in sets $M'$ and $D'$. The rotation expressed as quaternion that minimizes equation (1) is the largest eigenvalue of the cross-covariance matrix

$$\begin{eqnarray*} \begin{scriptsize} \M N = \left( \begin{array}{cccc} (S_{xx} +... ...(- S_{xx} - S_{yy} + S_{zz}) \end{array}\right), \end{scriptsize}\end{eqnarray*}$$

with $S_{xx} = \sum_{i=1}^{N_m}\sum_{j=1}^{N_d} w_{i,j} \m'_{ix} d'_{jx}, \S_{xy} = \sum_{i=1}^{N_m}\sum_{j=1}^{N_d} w_{i,j} \ m'_{ix} d'_{jy}, ldots$. After the calculation of the rotation $\M R$ the translation is $\V t = \V c_m - \M R \V c_d$ [5]. Fig. 1 shows three steps of the ICP algorithm¹, The corresponding surface meshes are given in Fig. 2.

Figure 1: Registration of two 3D reliefs with the ICP algorithms. Left: Initial alignment. Middle: Alignment after 4 iterations. Right: Final alignment after 85 iterations.