form solution with Euler angles parameterization is not
feasible. Luo and Hancock [6], [13] find the rotation matrix
through singular value decomposition (SVD). They ignore
some terms of the objective function, which leads to only an
approximate solution. We shall derive the exact closed form
solution (M-step) for the rigid point set registration and
discuss its difference from the related methods in Section 4.
2.2 Nonrigid Point Set Registration Methods
Earlier works on nonrigid point set registration include
Hinton et al. [19], [20], who used the probabilistic GMM
formulation. The GMM centroids are uniformly positioned
along the contour (modeled using splines), which allows for
nonrigid transformations. In practice, the method is applic-
able only to contour-like point sets. One of the most popular
nonrigid point set registration method is by Chui and
Rangarajan [8]. They proposed using Thin Plate Spline
(TPS) [21], [22] parameterization of the transformation,
following RPM, which results in the TPS-RPM method.
Similarly to the rigid case, they use deterministic annealing
and alternate updates for soft assignment and TPS para-
meters estimation. They also showed that TPS-RPM is
equivalent (with several modifications) to EM for GMM [9].
Tsin and Kan ade [ 23] proposed a correlation-ba sed
approach to point set registration, which was later
improved by Jian and Vemuri [24]. The method considers
the registration as an alignment between two distributions,
where each of the point sets represents the GMM centroids.
One of the point sets is parameterized by rigid/affine
parameters (in the rigid/affine case) or TPS (in the nonrigid
case). The transformation parameters are estimated to
minimize the L
2
norm between the distributions. These
methods all use explicit TPS parameterization, which is
equivalent to a regularization of second order derivatives of
the transformation [21], [22]. The TPS parameterization
does not exist when the dimension of points is higher than
three, which limits the applicability of such methods.
Huang et al. [25] proposed implicitly embedding the
shape (or point sets in our case) into distance transform
space and aligning them similarly to nonrigid image
registration algorithms. The authors use sum-of-squared-
differences similarity measure between the embedded
spaces and incremental free form deformation (FFD) to
parameterize the transformation. The method performs
well on relatively simple data sets.
Finally, in our previous work, we presented the Coherent
Point Drift (CPD) algorithm [26] for nonrigid point set
registration. The algorithm regularizes the displacement
(velocity) field between the point sets following the motion
coherence theory (MCT) [27], [28]. Using variational
calculus, we obtained that the optimal displacement field
has an elegant kernel form in multiple dimensions. In this
paper, we shall elaborate and analyze the CPD algorithm.
We also extend the general nonrigid registration framework
and show that CPD and TPS-RPM are its special cases.
Among other contributions, we estimate the width of GMM
components without the time-consuming deterministic
annealing and show a fast CPD implementation to reduce
the computational complexity to linear. We shall discuss
and compare o ur method to the works of Chui and
Rangarajan [8] and Jian and Vemuri [24] in Section 5.
3GENERAL METHODOLOGY
We consider the alignment of two point sets as a probability
density estimation problem, where one point set represents
the GMM centroids and the other one represents the data
points. At the optimum, two point sets become aligned and
the correspondence is obtained using the maximum of the
GMM posterior probability for a given data point. Core to
our method is to force GMM centroids to move coherently
as a group to preserve the topological structure of the point
sets. Throughout the paper, we use the following notations:
. D—dimension of the point sets,
. N;M—number of points in the point sets,
. X
ND
¼ðx
1
; ...; x
N
Þ
T
—the first point set (the data
points),
. Y
MD
¼ðy
1
; ...; y
M
Þ
T
—the second point set (the
GMM centroids),
. TðY;Þ—Transformation T applied to Y, where is
a set of the transformation parameters,
. I—identity matrix,
. 1—column vector of all ones,
. dðaÞ—diagonal matrix formed from the vector a.
We consider the points in Y as the GMM centroids and the
points in X as the data points generated by the GMM. The
GMM probability density function is
pðxÞ¼
X
Mþ1
m¼1
PðmÞpðxjmÞ; ð1Þ
where pðxjmÞ¼
1
ð2
2
Þ
D=2
exp
xy
m
kk2
2
2
. We also added an
additional uniform distribution pðxjM þ 1Þ¼
1
N
to the
mixture model to account for noise and outliers. We use
equal isotropic covariances
2
and equal membership
probabilities P ðmÞ¼
1
M
for all GMM components
(m ¼ 1; ...;M). Denoting the weight of the uniform
distribution as w, 0 w 1, the mixture model takes
the form
pðxÞ¼w
1
N
þð1 wÞ
X
M
m¼1
1
M
pðxjmÞ: ð2Þ
We reparameterize the GMM centroid locations by a set of
parameters and estimate them by maximizing the
likelihood or, equivalently, by minimizing the negative
log-likelihood function
Eð;
2
Þ¼
X
N
n¼1
log
X
Mþ1
m¼1
PðmÞpðx
n
jmÞ; ð3Þ
where we make the i.i.d. data assumption. We define the
correspondence probability between two points y
m
and x
n
as the posterior probability of the GMM centroid given the
data point: P ðmjx
n
Þ¼PðmÞpðx
n
jmÞ=pðx
n
Þ.
We use the EM algorithm [29], [30] to find and
2
. The
idea of EM is first to guess the values of parameters (“old”
parameter values) and then use the Bayes’ theorem to
compute a posteriori probability distributions P
old
ðmjx
n
Þ of
mixture components, which is the expectation or E-step of
the algorithm. The “new” parameter values are then found
2264 IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 32, NO. 12, DECEMBER 2010