Point-Set Registration: Coherent Point Drift
aa r X i v : . [ c s . C V ] M a y Point Set Registration: Coherent Point Drift
Andriy Myronenko and Xubo Song
Abstract —Point set registration is a key component in many computer vision tasks. The goal of point set registration is to assigncorrespondences between two sets of points and to recover the transformation that maps one point set to the other. Multiple factors,including an unknown non-rigid spatial transformation, large dimensionality of point set, noise and outliers, make the point setregistration a challenging problem. We introduce a probabilistic method, called the Coherent Point Drift (CPD) algorithm, for bothrigid and non-rigid point set registration. We consider the alignment of two point sets as a probability density estimation problem. Wefit the GMM centroids (representing the first point set) to the data (the second point set) by maximizing the likelihood. We force theGMM centroids to move coherently as a group to preserve the topological structure of the point sets. In the rigid case, we impose thecoherence constraint by re-parametrization of GMM centroid locations with rigid parameters and derive a closed form solution of themaximization step of the EM algorithm in arbitrary dimensions. In the non-rigid case, we impose the coherence constraint by regularizingthe displacement field and using the variational calculus to derive the optimal transformation. We also introduce a fast algorithm thatreduces the method computation complexity to linear. We test the CPD algorithm for both rigid and non-rigid transformations in thepresence of noise, outliers and missing points, where CPD shows accurate results and outperforms current state-of-the-art methods.
Index Terms —Registration, correspondence, matching, alignment, rigid, non-rigid, point sets, Coherent Point Drift (CPD), Gaussianmixture model (GMM), coherence, regularization, EM algorithm. F NTRODUCTION R EGISTRATION of point sets is a key component inmany computer vision tasks including stereo match-ing, content-based image retrieval, image registrationand shape recognition. The goal of point set registrationis to assign correspondences between two sets of pointsand/or to recover the transformation that maps onepoint set to the other. For example, in stereo matching,in order to recover depth and infer structure from apair of stereo images, it is necessary to first define aset of points in each image and find the correspondencebetween them. An example of point set registrationproblem is shown in Fig. 1. The “points” in a point setare often features extracted from an image, such as thelocations of corners, boundary points or salient regions.The points can represent both geometric and intensitycharacteristics of an image.Practical point set registration algorithms should haveseveral desirable properties: (1) Ability to accuratelymodel the transformation required to align the point setswith tractable computational complexity; (2) Ability tohandle possibly high dimensionality of the point sets;(3) Robustness to degradations such as noise, outliersand missing points that occur due to imperfect imageacquisition and feature extraction.The transformation usually falls into two categories:rigid or non-rigid. A rigid transformation allows only fortranslation, rotation and scaling. The simplest non-rigid • A. Myronenko and X. Song are with the Department of Science andEngineering, School of Medicine, Oregon Health and Science University,Portland, OR, 97201.E-mail: [email protected], [email protected] transformation is affine, which also allows anisotropicscaling and skews. Non-rigid transformation occurs inmany real-world problems including deformable motiontracking, shape recognition and medical image regis-tration. The true underlying non-rigid transformationmodel is often unknown and challenging to model.Simplistic approximations of the true non-rigid trans-formation, including piece-wise affine and polynomialmodels, are often inadequate for correct alignment andcan produce erroneous correspondences. Due to the usu-ally large number of transformation parameters, the non-rigid point sets registration methods tend to be sensitiveto noise and outliers and are likely to converge into localminima. They also tend to have a high computationalcomplexity. A practical non-rigid point set registrationmethod should be able to accurately model the non-rigidtransformation with tractable computational complexity.Multidimensional point sets are common in many realworld problems. Most current rigid and non-rigid pointsets registration algorithm are well suited for 2D and3D cases, but their generalization to higher dimensionsare not always trivial. Furthermore, degradations suchas noise, outliers and missing points significantly com-plicates the problem. Outliers are the points that areincorrectly extracted from the image; outliers have nocorrespondences in the other point set. Missing pointsare the features that are not found in the image due toocclusion or inaccurate feature extraction. A point setregistration method should be robust to these degrada-tions.We present a robust probabilistic multidimensionalpoint sets registration algorithm for both rigid and non-rigid transforms. We consider the alignment of twopoint sets as a probability density estimation problem,
X: −0.1183Y: 2.111 ? Fig. 1. The point set registration problem: Given two setsof points, assign the correspondences and the transfor-mation that maps one point set to the other. where one point set represents the Gaussian MixtureModel (GMM) centroids, and the other one representsthe data points. We fit the GMM centroids to the databy maximizing the likelihood. At the optimum, thepoint sets become aligned and the correspondence isobtained using the posterior probabilities of the GMMcomponents. Core to our method is to force GMM cen-troids to move coherently as a group, which preservesthe topological structure of the point sets. We imposethe coherence constraint by explicit re-parametrizationof GMM centroid locations (for rigid and affine transfor-mations) or by regularization of the displacement field(for smooth non-rigid transformation). We also showhow the computational complexity of the method canbe reduced to linear, which makes it applicable forlarge data sets. The rest of the paper is organized asfollows. In Section 2, we overview the current rigid andnon-rigid point set registration methods and state ourcontributions. In Section 3, we formulate a probabilisticpoint set registration framework. In Sections 4 and 5, wedescribe our algorithms for rigid, affine and non-rigidregistration cases, and relate them to existing works. InSection 6, we discuss the computational complexity ofthe method and introduce its fast implementation. InSection 7, we evaluate the performance of our algorithm.Section 8 concludes with some discussions.
REVIOUS W ORK
Many algorithms exist for rigid and for non-rigid pointset registration. They aim to recover the correspondenceor the transformation required to align the point setsor both. Many algorithms involve a dual-step update,iteratively alternating between the correspondence andthe transformation estimation. Here, we briefly overviewthe rigid and non-rigid point set registration methodsand state our contributions.
Iterative Closest Point (ICP) algorithm, introduced byBesl and McKay [1] and Zhang [2], is the most popularmethod for rigid point set registration due to its simplic-ity and low computational complexity. ICP iterativelyassigns correspondences based on a closest distancecriterion and finds the least-squares rigid transformation relating the two point sets. The algorithm then redeter-mines the correspondences and continues until it reachesthe local minimum. Many variants of ICP have beenproposed that affect all phases of the algorithm fromthe selection and matching of points to the minimizationstrategy [3], [4]. ICP requires that the initial position ofthe two point sets be adequately close.To overcome the ICP limitations, many probabilisticmethods were developed [5], [6]. These methods usesoft-assignment of correspondences that establishes cor-respondences between all combinations of points accord-ing to some probability, which generalizes the binaryassignment of correspondences in ICP. Among thesemethods are Robust Point Matching (RPM) algorithmintroduced by Gold et al. [7], and its later variants [5],[8], [9]. In [10] it was shown that in RPM alternatingsoft-assignment of correspondences and transformationis equivalent to the Expectation Maximization (EM) al-gorithm for GMM, where one point sets is treated asGMM centroids with equal isotropic covariances andthe other point set is treated as data points. In fact,several rigid point set methods, including Joshi andLee [11], Wells [12], Cross and Hancock [13], Luo andHancock [6], [14], McNeill and Vijayakumar [15] andSofka et al. [16], explicitly formulate the point setsregistration as a maximum likelihood (ML) estimationproblem, to fit the GMM centroids to the data points.These methods re-parameterize GMM centroids by aset of rigid transformation parameters (translation androtation). The EM algorithm used for optimization of thelikelihood function consists of two steps: E-step to com-pute the probabilities, M-step to update the transforma-tion. Common to such probabilistic methods is the inclu-sion of an extra distribution term to account for outliers(e.g. large Gaussian [5] or uniform distribution [12]) anddeterministic annealing on the Gaussian width to avoidpoor local minima. These probabilistic methods performbetter than conventional ICP, especially in presence ofnoise and outliers.Another class of rigid point sets registration methodsare the spectral methods. Scott and Longuet-Higgins [17]introduced a non-iterative algorithm to associate pointsof two arbitrary patterns, exploiting some properties ofGaussian proximity matrix (Gram matrix) of point sets.The algorithm works well with translation, shearing andscaling deformations, but performs poorly with non-rigid transformations. Li and Hartley showed that corre-spondence and transformation are two factors of Grammatrices, and can be found iteratively using Newton-Schulz factorization [18]. This method performs wellfor moderate linear transformations. In spite of its el-egance, the large computational effort required for spec-tral methods prohibits its wide applicability. There area few other nonspectral methods worth mentioning. Hoet al. [19] proposed an elegant non-iterative algorithmfor 2D affine registration by searching for the roots ofthe associated polynomials. Unfortunately this methoddoes not generalize to higher dimensions. Belongie et al. [20] introduced the “shape context” descriptor, whichincorporates the neighborhood structure of the point setand thus helps to recover the correspondence betweenthe point sets.Our approach to the rigid point sets registration isprobabilistic and most closely related to the works ofRangarajan et al. [5], Wells [12] and Luo and Han-cock [14]. Despite extensive work in rigid probabilisticregistration, none of the methods, to our best knowledge,provides a closed form solution to the maximizationstep (M-step) of the EM algorithm for a general mul-tidimensional case. The fact that the rotation matrixshould be constrained to be orthogonal and to have apositive determinant further complicates its estimation.Rangarajan and collaborators [5] showed the solution for2D case only, where rotation is parametrized by a singleangle. In higher dimensions the closed form solutionwith Euler angles parametrization is not feasible. Luoand Hancock [6], [14] find the rotation matrix throughsingular value decomposition (SVD). They ignored someterms of the objective function, which leads to only anapproximate solution. We shall derive the exact closedform solution (M-step) for the rigid point set registrationand discuss its difference from the related methods inSection 4.
Earlier works on non-rigid point set registration in-clude Hinton et al. [21], [22], who used the probabilisticGMM formulation. The GMM centroids are uniformlypositioned along the contour (modeled using splines),which allows for non-rigid transformations. In practice,the method is applicable only to contour-like pointsets. One of the most popular non-rigid point set reg-istration method is by Chui and Rangarajan [8], [9].They proposed to use Thin Plate Spline (TPS) [23], [24]parametrization of the transformation, following RPM,which results into the TPS-RPM method. Similar to therigid case, they use deterministic annealing and alternateupdates for soft-assignment and TPS parameters estima-tion. They also showed that TPS-RPM is equivalent (withseveral modifications) to EM for GMM [10]. Tsin andKanade [25] proposed a correlation-based approach topoint set registration, which was later improved by Jianand Vemuri [26]. The method considers the registrationas an alignment between two distributions, where eachof the point sets represents the GMM centroids. One ofthe point sets is parametrized by rigid/affine parameters(in rigid/affine case) or TPS (in non-rigid case). Thetransformation parameters are estimated to minimize the L norm between the distributions. These methods alluse explicit TPS parametrization, which is equivalentto a regularization of second order derivatives of thetransformation [23], [24]. The TPS parametrization doesnot exist when the dimension of points is higher thanthree, which limits the applicability of such methods.Huang et al. [27] proposed to implicitly embed theshape (or point sets in our case) into distance transform space, and align them similar to non-rigid image reg-istration algorithms. The authors use sum-of-squared-differences similarity measure between the embeddedspaces and incremental free form deformation (FFD) toparameterize the transformation. The method performswell on relatively simple data sets.Finally, in our previous work we presented the Coher-ent Point Drift (CPD) algorithm [28] for non-rigid pointsets registration. The algorithm regularizes the displace-ment (velocity) field between the point sets following themotion coherence theory (MCT) [29], [30]. Using varia-tional calculus, we obtained that the optimal displace-ment field has an elegant kernel form in multiple di-menstions. In this paper, we shall elaborate and analyzethe CPD algorithm. We also extend the general non-rigidregistration framework, and show that CPD and TPS-RPM are its special cases. Among other contributions,we estimate the width of GMM components withoutthe time consuming deterministic annealing and showa fast CPD implementation to reduce the computationalcomplexity to linear. We shall discuss and compare ourmethod to the works of Chui and Rangarajan [9], andJian and Vemuri [26] in Section 5. ENERAL M ETHODOLOGY
We consider the alignment of two point sets as a proba-bility density estimation problem, where one point setrepresents the Gaussian mixture model (GMM) cen-troids, and the other one represents the data points.At the optimum, two point sets become aligned andthe correspondence is obtained using the maximum ofthe GMM posterior probability for a given data point.Core to our method is to force GMM centroids to movecoherently as a group to preserve the topological struc-ture of the point sets. Throughout the paper we use thefollowing notations: • D - dimension of the point sets, • N, M - number of points in the point sets, • X N × D = ( x , . . . , x N ) T - the first point set (the datapoints), • Y M × D = ( y , . . . , y M ) T - the second point set (theGMM centroids), • T ( Y , θ ) - Transformation T applied to Y , where θ is a set of the transformation parameters, • I - identity matrix, • - 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 theGMM. The GMM probability density function is p ( x ) = M +1 X m =1 P ( m ) p ( x | m ) (1)where p ( x | m ) = πσ ) D/ exp − k x − y m k σ . We also addedan additional uniform distribution p ( x | M + 1) = N tothe mixture model to account for noise and outliers. We use equal isotropic covariances σ and equal member-ship probabilities P ( m ) = M for all GMM components( m = 1 , . . . , M ). Denoting the weight of the uniformdistribution as w , ≤ w ≤ , the mixture model takesthe form p ( x ) = w N + (1 − w ) M X m =1 M p ( x | m ) (2)We re-parameterize the GMM centroid locations by a setof parameters θ and estimate them by maximizing thelikelihood, or equivalently by minimizing the negativelog-likelihood function E ( θ, σ ) = − N X n =1 log M +1 X m =1 P ( m ) p ( x | m ) (3)where we make the i.i.d. data assumption. We definethe correspondence probability between two points y m and x n as the posterior probability of the GMM centroidgiven the data point: P ( m | x n ) = P ( m ) p ( x n | m ) /p ( x n ) .We use Expectation Maximization (EM) algorithm [31],[32] to find θ and σ . The idea of EM is first to guess thevalues of parameters (“old” parameter values) and thenuse the Bayes’ theorem to compute a posteriori proba-bility distributions P old ( m | x n ) of mixture components,which is the expectation or E-step of the algorithm. The“new” parameter values are then found by minimizingthe expectation of the complete negative log-likelihoodfunction [32] Q = − N X n =1 M +1 X m =1 P old ( m | x n ) log( P new ( m ) p new ( x n | m )) (4)with respect to the “new” parameters, which is calledthe maximization or M-step of the algorithm. The Q function, which we call the objective function , is also anupper bound of the negative log-likelihood function in(3). The EM algorithm proceeds by alternating betweenE- and M-steps until convergence. Ignoring the constantsindependent of θ and σ , we rewrite (4) as Q ( θ, σ ) = 12 σ N X n =1 M X m =1 P old ( m | x n ) k x n − T ( y m , θ ) k + N P D σ (5)where N P = P Nn =1 P Mm =1 P old ( m | x n ) ≤ N (with N = N P only if w = 0 ) and P old denotes the posteriorprobabilities of GMM components calculated using theprevious parameter values: P old ( m | x n ) = exp − ‚‚‚‚ x n −T ( y m,θold ) σold ‚‚‚‚ P Mk =1 exp − ‚‚‚‚ x n −T ( y k,θold ) σold ‚‚‚‚ + c (6)where c = (2 πσ ) D/ w − w MN . Minimizing the function Q , we necessarily decrease the negative log-likelihoodfunction E , unless it is already at a local minimum. To proceed, we specify the transformation T for rigid, affineand non-rigid point set registration cases separately. IGID & A
FFINE P OINT SET R EGISTRATION
For rigid point set registration, we define the transforma-tion of the GMM centroid locations as T ( y m ; R , t , s ) = s Ry m + t , where R D × D is a rotation matrix, t D × isa translation vector and s is a scaling parameter. Theobjective function (5) takes the form: Q ( R , t , s, σ ) = 12 σ M,N X m,n =1 P old ( m | x n ) k x n − s Ry m − t k + N P D σ , s.t. R T R = I , det( R ) = 1 (7)Note that the first term is similar to the one in theabsolute orientation problem [33], [34], which is definedas min P Nn =1 k x n − ( s Ry n + t ) k in our notations. Equa-tion (7) can be seen as a generalized weighted absoluteorientation problem, because it includes weighted dif-ferences between all combinations of points. The exactminimization solution of the objective function (7) iscomplicated due to the constraints on R . To obtain theclosed form solution we shall use Lemma 1 [35]: Lemma 1:
Let R D × D be an unknown rotation ma-trix and A D × D be a known real square matrix. Let U SS V T be a Singular Value Decomposition (SVD) of A , where UU T = VV T = I and SS = d( s i ) with s ≥ s ≥ , . . . , ≥ s D ≥ . Then the optimal rotationmatrix R that maximizes tr ( A T R ) is R = UCV T , where C = d(1 , , . . . , , det( UV T )) .To apply this lemma, we need to simplify the Q functionto a form equivalent to tr ( A T R ) . First, we eliminatetranslation t from Q . Taking partial derivative of Q withrespect to t and equate it to zero, we obtain: t = 1 N P X T P T − s R N P Y T P1 = µ x − s R µ y , where the matrix P has elements p mn = P old ( m | x n ) in(6) and the mean vectors µ x and µ y are defined as: µ x = E ( X ) = 1 N X T P T , µ y = E ( Y ) = 1 N Y T P1 . Substituting t back into the objective function andrewritting it in matrix form, we obtain Q = 12 σ [tr( ˆ X T d( P T ) ˆ X ) − s tr( ˆ X T P T ˆ YR T )+ s tr( ˆ Y T d( P1 ) ˆ Y )] + N P D σ (8)where ˆ X = X − µ T x and ˆ Y = Y − µ T y are the centeredpoint set matrices. We use the fact that trace is invariantunder cyclic matrix permutations and R is orthogonal.We can rewrite (8) as Q = − c tr(( ˆ X T P T ˆ Y ) T R ) + c ,where c , c are constants independent of R and c > .Thus minimization of Q with respect to R is equivalentto maximization of max tr( A T R ) , A = ˆ X T P T ˆ Y , s.t. R T R = I , det( R ) = 1 . Rigid point set registration algorithm: • Initialization: R = I , t = 0 , s = 1 , ≤ w ≤ σ = DNM P Nn =1 P Mm =1 k x n − y m k • EM optimization, repeat until convergence: • E-step: Compute P , p mn = exp − σ k x n − ( s Ry m + t ) k P Mk =1 exp − σ k x n − ( s Ry k + t ) k +(2 πσ ) D/ w − w MN • M-step: Solve for R , s, t , σ : · N P = T P1 , µ x = N P X T P T , µ y = N P Y T P1 , · ˆ X = X − µ T x , ˆ Y = Y − µ T y , · A = ˆ X T P T ˆ Y , compute SVD of A = U SS V T , · R = UCV T , where C = d(1 , .., , det( UV T )) , · s = tr( A T R )tr( ˆ Y T d( P1 ) ˆ Y ) , · t = µ x − s R µ y , · σ = N P D (tr( ˆ X T d( P T ) ˆ X ) − s tr( A T R )) . • The aligned point set is T ( Y ) = s YR T + 1 t T , • The probability of correspondence is given by P . Fig. 2. Rigid point set registration algorithm.
Now we are ready to use Lemma 1, and the optimal R is in the form R = UCV T , where U SS V T = svd ( ˆ X T P T ˆ Y ) (9)and C = d(1 , .., , det( UV T )) . To solve for s and σ , weequate the corresponding partial derivative of (8) to zero.Solving for R , s, t , σ is the M-step of the EM algorithm.We summarize the rigid point sets registration algorithmin Fig. 2.The algorithm has one free paramater, w ( ≤ w ≤ ),which reflects our assumption on the amount of noisein the point sets. The solution for the rotation matrix isgeneral D -dimensional. Affine point set registration:
Affine registration caseis simpler compared to the rigid case, because theoptimization is unconstrained. Affine transformation isdefined as T ( y m ; R , t , s ) = By m + t , where B D × D is anaffine transformation matrix, t D × is translation vector.The objective function takes the form: Q ( B , t , σ ) = 12 σ M,N X m,n =1 P old ( m | x n ) k x n − ( By m + t ) k + N P D σ (10)We can directly take the partial derivatives of Q , equatethem to zero, and solve the resulting linear system ofequations. The solution is straightforward and similarto the rigid case. We summarize the affine point setregistration algorithm in Fig. 3. Here, we discuss the probabilistic rigid point set regis-tration methods most closely related to ours. Rangarajanet al. [5] presented the RPM method for rigid point setregistration. The method is shown for 2D case, where
Affine point set registration algorithm: • Initialization: B = I , t = 0 , ≤ w ≤ σ = DNM P Nn =1 P Mm =1 k x n − y m k • EM optimization, repeat until convergence: • E-step: Compute P , p mn = exp − σ k x n − ( By m + t ) k P Mk =1 exp − σ k x n − ( By k + t ) k +(2 πσ ) D/ w − w MN • M-step: Solve for B , t , σ : · N P = T P1 , µ x = N P X T P T , µ y = N P Y T P1 , · ˆ X = X − µ T x , ˆ Y = Y − µ T y , · B = ( ˆ X T P T ˆ Y )( ˆ Y T d( P1 ) ˆ Y ) − , · t = µ x − B µ y , · σ = N P D (tr( ˆ X T d( P T ) ˆ X ) − tr( ˆ X T P T ˆ YB T )) . • The aligned point set is T ( Y ) = YB T + 1 t T , • The probability of correspondence is given by P . Fig. 3. Affine point set registration algorithm. rotation matrix is parametrized by a single rotationangle, which allows to find an explicit update. Such Eu-ler’s angles approach is not feasible in multidimensionalcases. RPM also uses deterministic annealing on σ ,which requires to set the starting and stopping criteriaas well as the annealing rate. The EM iterations has to berepeated at each annealing step, which can be slow. Weargue that it is preferable to estimate σ instead of us-ing deterministic annealing. The deterministic annealinghelps to overcome poor local minima, but for the rigidpoint set registration problem the rigid parametrizationis a strong constraint that alleviates the advantages ofthe deterministic annealing.Luo and Hancock [14], [36] introduced the rigid pointsets registration algorithm that is the most similar toours. The authors optimize the objective function ratherintuitively than rigorously, which leads to several as-sumptions and approximate minimization. They ignore afew terms of the objective function (see Eqs.10,11 in [36]),where the last term does depend on transformationparameters, and must not be ignored. If such optimiza-tion converge, the M-step of the EM algorithm is onlyapproximate. Among other differences, we want to men-tion that the authors use structural editing , a techniqueto remove some undesirable points, instead of usingan additional uniform distribution to account for thesepoints. Some other authors [15] also follow the rigidparameters estimation steps of Luo and Hancock [36]. ON -R IGID P OINT SET R EGISTRATION
Non-rigid point set registration remains a challengingproblem in computer vision. The transformation thataligns the point sets is assumed to be unknown andnon-rigid, which is generally broad class of transforma-tions that can lead to an ill-posed problem. In order todeal with the problem we use Tikhonov regularizationframework [37]–[39]. We define the transformation as the initial position plus a displacement function v : T ( Y , v ) = Y + v ( Y ) , (11)We regularize the norm of v to enforce the smoothnessof the function [38]. Such approach is also supported bythe Motion Coherence Theory (MCT) [29], [30], whichstates that points close to one another tend to movecoherently, and thus, the displacement function betweenthe point sets should be smooth. This is mathematicallyformulated as a regularization on the displacement (alsocalled velocity) function.Additing a regularization term to the negative log-likelihood function we obtain f ( v, σ ) = E ( v, σ ) + λ φ ( v ) (12)where E is the negative log-likelihood function (3), φ ( v ) is a regularization term and λ is a trade-off parame-ter. Such regularization is well formulated in Bayesianapproach, where the regularization term comes from aprior on displacement field: p ( v ) = exp − λ φ ( v ) .We estimate the displacement function v using varia-tional calculus. We shall define the regularization term φ ( v ) in different but equivalent forms and show thatthe optimal functional form of v is a linear combinationof particular kernel functions. A particular choice ofthe regularization will lead to our non-rigid point setregistration method. A norm of v in the Hilbert space H m is defined as: k v k H m = Z R m X k =0 (cid:13)(cid:13)(cid:13)(cid:13) ∂ k v∂x k (cid:13)(cid:13)(cid:13)(cid:13) dx. (13)Alternatively, we can define the norm in the Reproduc-ing Kernel Hilbert Space (RKHS) [38], [40] as k v k H m = Z R D | ˜ v ( s ) | ˜ G ( s ) d s (14)where G is a unique kernel function associated withthe RKHS, and ˜ G is its Fourier transform. Function ˜ v indicates the Fourier transform of the function v and s isa frequency domain variable. The Fourier domain normdefinition has been used in the Regularization Theory(RT) [40] to regularize the smoothness of a function.Regularization theory defines smoothness as a measureof the “oscillatory” behavior of a function. Within theclass of differentiable functions, one function is said tobe smoother than another if it oscillates less; in otherwords, if it has less energy at high frequency. The highfrequency content of a function can be measured byfirst high-pass filtering the function, and then measuringthe resulting power. This can be represented by (14),where ˜ G represents a symmetric positive definite low-pass filter, which approaches zero as k s k → ∞ . Forconvenience, we shall write the regularization term as φ ( v ) = k v k H m = k P v k (15) where an operator P “extracts” a part of the function forregularization, in our case, the high frequency contentpart [38], [39]. We find the functional form of v using calculus of varia-tion. Minimization of regularized negative log-likelihoodfunction in (12) boils down to minimization of the fol-lowing objective function (M-step): Q ( v, σ ) = 12 σ M,N X m,n =1 P old ( m | x n ) k x n − ( y m + v ( y m )) k + N P D σ + λ k P v k (16)A function v that minimizes (16) must satisfy the Euler-Lagrange differential equation σ λ N X n =1 M X m =1 P old ( m | x n )( x n − ( y m + v ( y m ))) δ ( z − y m )= ˆ P P v ( z ) (17)for all vectors z , where ˆ P is the adjoint operator to P .The solution to such partial differential equation can bewritten as the integral transformation of its left side witha Green’s function G of the self-adjoint operator ˆ P P . v ( z ) = 1 σ λ M,N X m,n =1 P old ( m | x n )( x n − ( y m + v ( y m ))) G ( z , y m )= M X m =1 w m G ( z , y m ) (18)where w m = σ λ P Nn =1 P old ( m | x n )( x n − ( y m + v ( y m ))) .Note that this solution is incomplete. In general, thesolution also includes the term ψ ( z ) that lies in the nullspace of P [40], [41]. Thus, we achieve Lemma 2. Lemma 2:
The optimal displacement function thatminimizes (16) is given by linear combination of theparticular kernel functions centered at the points Y plusthe term ψ ( z ) in the null space of P : v ( z ) = M X m =1 w m G ( z , y m ) + ψ ( z ) (19)where the kernel function is a Green’s function of theself-adjoint operator ˆ P P . We choose the regularization term according to (14): φ ( v ) = Z R D | ˜ v ( s ) | ˜ G ( s ) d s (20)where G is a Gaussian (note it is not related to theGaussian form of the distribution chosen for the mix-ture model). There are several motivations for such aGaussian choice: First, the Green’s function (the kernel) corresponding to the regularization term in (20) is alsoa Gaussian (and remains a Gaussian for an arbitrary di-mensional case); the Gaussian kernel is positive definiteand the null space term ψ ( z ) = 0 [40]. Second, by choos-ing an appropriately sized Gaussian function we havethe flexibility to control the range of filtered frequenciesand thus the amount of spatial smoothness. Third, thechoice of the Gaussian makes our regularization termequivalent to the one in the Motion Coherence Theory(MCT) [30]: φ MCT ( v ) = Z R d ∞ X l =0 β l l !2 l (cid:13)(cid:13) D l v ( x ) (cid:13)(cid:13) d x (21)where D is a derivative operator so that D l v = ∇ l v and D l +1 v = ∇ ( ∇ l v ) , where ∇ is the gradient operator and ∇ is the Laplacian operator. Lemma 3:
The regularization term in (20) with aGaussian choice of low-pass filter G is equivalent to thethe regularization term in (21). Both terms represent thenorm of the function v , after applying the operator P ,and the corresponding Green’s function is a Gaussian inboth cases [38].The equivalence of our regularization term with thatof the Motion Coherence Theory implies that we areimposing motion coherence among the points and thuswe call the non-rigid point set registration method theCoherent Point Drift (CPD) algorithm.We can obtain the coefficients w m by evaluating (19)at y m points ( G + λσ d ( P1 ) − ) W = d ( P1 ) − PX − Y (22)where W M × D = ( w , . . . , w M ) T is a matrix of coeffi-cients, G M × M is a kernel matrix with elements g ij = G ( y i , y j ) = e − ‚‚‚ y i − y jβ ‚‚‚ and d − ( · ) is the inverse diag-onal matrix. The transformed position of y m are foundaccording to (11) as T = T ( Y , W ) = Y + GW . We obtain σ by equating the corresponding derivative of Q to zero σ = 1 N P D N X n =1 M X m =1 k x n − T ( y m , W ) k =1 N P D (tr( X T d( P T ) X ) − PX ) T T )+tr( T T d( P1 ) T )) (23)We summarize the CPD non-rigid point set registrationalgorithm in Fig. 4. Analysis:
The CPD algorithm includes three free pa-rameters: w, λ and β . Parameter w ( ≤ w ≤ ) reflectsour assumption on the amount of noise in the pointsets. Parameters λ and β both reflect the amount ofsmoothness regularization. A discussion on the differ-ence between λ and β can be found in [29], [30]. Brieflyspeaking, parameter β defines the model of the smooth-ness regularizer (width of smoothing Gaussian filter in(20)). Parameter λ represents the trade-off between thegoodness of maximum likelihood fit and regularization. Non-rigid point set registration algorithm: • Initialization: W = 0 , σ = 1 DN M
M,N X m,n =1 k x n − y m k • Initialize w ( ≤ w ≤ ), β > , λ > , • Construct G : g ij = exp − β k y i − y j k , • EM optimization, repeat until convergence: • E-step: Compute P , p mn = exp − σ k x n − ( y m + G ( m, · ) W ) k P Mk =1 exp − σ k x n − ( y k + G ( k, · ) W ) k + w − w (2 πσ D/ MN • M-step: · Solve ( G + λσ d ( P1 ) − ) W = d ( P1 ) − PX − Y · N P = T P1 , T = Y + GW , · σ = N P D (tr( X T d( P T ) X ) − PX ) T T )+tr( T T d( P1 ) T )) , • The aligned point set is T = T ( Y , W ) = Y + GW , • The probability of correspondence is given by P . Fig. 4. The Coherent Point Drift algorithm for non-rigidpoint set registration.
We note that solution of (22) gives the exact minimumof Q (16), if σ is assumed fixed. As far as we areestimating σ , (22) and (23) should be solved simulta-neously. The non-linear dependency of σ on W andvice-verse does not allow for simultaneous analyticalsolution. Iterative exact solution can be obtained byperforming a few cyclic iterations on (22) and (23) withina single EM step. Practically, a single iteration, given by(22) and (23), decrease the Q function almost to the exactminimum. Such an iterative procedure, which decreasesthe Q function but not to exact minimum, is often calledthe generalized EM algorithm [31], [42]. The CPD algorithm follows our previous work [28] onnon-rigid point sets registration. However, previouslywe have used deterministic annealing on σ , whereashere, we estimate the Gaussian width σ within MLframework. This allows us to significantly speed upthe algorithm, alleviating the repeated EM-iterationsfor every single annealing step. We have not observedany decrease in accuracy of the method related to thischange. In [28], we used a slightly different notation forthe GMM centroid locations: we called Y the initialcentroids position (which we call Y here), and Y forthe finall GMM centroid position (which we call T ( Y ) here).The most relevant non-rigid point sets registrationalgorithm to ours is TPS-RPM, more precisely itsGMM formulation [10]. TPS-RPM uses Thin Plate Spline(TPS) [23], [24] parametrization of the transformation,which can be obtained by adding the regularization termthat penalizes second order derivatives of the transfor- mation. For instance, in 2D such regularization term is k L T k = Z Z [( ∂ T ∂x ) + 2( ∂ T ∂x∂y ) + ( ∂ T ∂y ) ] dxdy (24)This term can be equivalently formulated in the Fourierspace as: k L T k = Z R k s k | ˜ T ( s ) | d s (25)which is a special case of the Duchon splines [43]. Thenull space of such regularization includes affine trans-formations. Using the variational approach we can showthat the optimal transformation T for such regularizationis in the form T ( Y ) = YA + KC , where A is matrix ofaffine transformation cooefficients, C is a matrix of non-rigid cooefficients. For 2D case, matrix K M × M is kernelmatrix with elements k ij = k y i − y j k log k y i − y j k . For3D case, matrix K has elements k ij = k y i − y j k . For4D or higher dimensions the TPS kernel solution doesnot exist [44]. Finally, to link such regularization toour non-rigid registration framework, we note that theregularization of the displacement field v , instead of thetransformation itself, is exactly the same, because, (24)is invariant under affine transformations, in other words k L T ( z ) k = k L ( z + v ( z )) k = k Lv ( z ) k . This means thatboth CPD and TPS-RPM regularizes the displacementfunction, but using different regularization terms.The advantage of CPD regularization (as given by (20)or (21)) comparing to TPS ((24) or (25)), is that it easilygeneralizes to N dimensions. Also we can control thelocality of spatial smoothness by changing the Gaussianfilter width β , whereas TPS does not have such flexibility.This, however, introduces one extra parameter to themethod, but TPS-RPM has to uses one extra parameterto regularize affine matrix after all. Among other differ-ences, TPS-RPM approximates the M-step solution of theEM algorithm [10] for simplicity and use deterministicannealing on σ .Finally, Jian and Vemuri [26] consider the registrationas an alignment between the distributions of two pointsets, where a separate GMMs are used to model thedistribution for the point sets. One of the point setsis parametrized by TPS. The transformation parametersare estimated to minimize the L norm between thedistributions. In our case, the CPD method maximizesthe likelihood function, which is equivalent to KL diver-gence minimization between two mixture distributions:GMM and mixture of delta functions. KL divergence ismore appropriate similarity measure for the densitiesthan L norm, because it weights the error accordingto its probability. AST I MPLEMENTATION
Here we show that CPD computational complexity canbe reduced to linear up to a multiplicative constant.We use the fast Gauss transform (FGT) [45] to computethe matrix-vector products P1 , P T , PX , which are thebottlenecks for both rigid and non-rigid cases. We use Compute P T , P1 and PX : • Compute K T (using FGT), • a = 1 ./ ( K T + c ) , • P T = − c a , • P1 = Ka (using FGT), • PX = K ( a . ∗ X ) (using FGT), Fig. 5. Matrix-vector products computation through FGT. low-rank matrix approximation to speed-up the solutionof the linear system of equations (22) for the non-rigidcase.
The fast Gauss transform:
Greengard and Strain [45]introduced the fast Gauss transform (FGT) for fast com-putation of the sum of exponentials: f ( y m ) = N X n =1 z n exp − σ k x n − y m k , ∀ y m , m = 1 , . . . , M. The naive approach takes O ( M N ) operations, while FGTtakes only O ( M + N ) . The basic idea of FGT is to expandthe Gaussians in terms of truncated Hermit expansion,and approximate (6) up to the predefined accuracy.Rewriting (6) in matrix form, we obtain f = Kz , where z is some vector and K M × N is a Gaussian affinity matrixwith elements: k mn = exp − σ k x n −T ( y m ) k , which wehave already used in our notations. We simplify thematrix-vector products P1 , P T and PX , to the formof Kz and apply FGT. Matrix P (6) can be partitionedinto P = K d( a ) , a = 1 ./ ( K T + c ) (26)where d( a ) is diagonal matrix with a vector a along thediagonal. Here, we use Matlab element-wise division ( ./ )and element-wise ( . ∗ ) multiplication notations. We showthe algorithm to compute the bottleneck matrix-vectorproducts P1 , P T and PX using FGT in Fig. 5. We notethat for dimensions higher than three, we can use the im-proved fast Gauss transform (IFGT) method [46], whichis a faster alternative to FGT for higher dimensions.During the finall EM iterations, the width of the Gaus-sians σ becomes small. The Hermitian expansion thusrequires many terms to approximate highly multimodalGaussian distribution for a given precision. At the finaliterations, the Gaussian becomes very narrow, and wecan switch to the truncated Gaussian approximation (setzeros outside a predefined box). Low-rank matrix approximation:
In the non-rigid case,we need to solve the linear system (22), which is O ( M ) using direct matrix inversion. We note that the left handside matrix of (22) is symmetric and positive definite.We use low-rank matrix approximation of G , where G is a Gaussian affinity matrix with elements g ij =exp − β k y i − y j k . We approximate the matrix G as ˆ G = Q Λ Q T (27)where Λ K × K is a diagonal matrix with K largest eigen-values and the matrix Q M × K is formed from the cor- a)b)c) Initialization Iteration 10 Iteration 30 Iteration 40 Result (iteration 50) Fig. 6. Fish data set, rigid registration examples. We align Y (blue circles) onto X (red stars). The columns showthe iterative alignment progress. a) Registration of the point sets with missing non-overlapping parts ( w = 0 . ); b)Registration of the point sets corrupted by random outliers ( w = 0 . ); c) A challenging rigid registration example,where both point sets are corrupted by outliers and biased to different sides of the point sets. We have also deletedsome parts from both point sets. We set w = 0 . and fix scaling s = 1 . CPD registration is robust and accurate in allexperiments. responding eigenvectors. ˆ G is the closest K -rank ma-trix approximation to G both in L and Frobeniusnorms [47]. To solve the linear system in (22) we usethe Woodbury identity and invert the first term as ( Q Λ Q T + λσ d( P1 ) − ) − = 1 λσ d( P1 ) − λσ ) d( P1 ) Q (Λ − + 1 λσ Q T d( P1 ) Q ) − Q T d( P1 ) (28)The inside matrix inversion is of O ( K ) , where K ≪ M . For instance choosing K = M / largest eigenvalues,we reduce the computational complexity to linear. Wecan pre-compute K largest eigenvalues and eigenvectorsof G using deflation techniques [48]. It requires severaliterations with the matrix-vector product Gz , which canbe implemented explicitly or through FGT.The low-rank matrix approximation intuitively con-straints the space of the non-rigid transformations, andcan be even desirable to further constrain the non-rigidtransformation. If the number of points is large and wellclustered, then an extremely small percent of eigenvalueswill be sufficient for an accurate approximation. ESULTS
We implemented the algorithm in Matlab, and tested iton a Pentium4 CPU 3GHz with 4GB RAM. We imple-mented the matrix-vector products in C as a Matlab mex functions to avoid the storage of P . The code is avail-able at .We shall refer to our method as Coherent Point Drift(CPD) both for rigid and non-rigid point sets registrationmethods presented in this paper. We have also imple-mented the matrix-vector products through FGT usingthe Matlab FGT implementation by Sebastien Paris [49].We consider rigid and non-rigid experiments sepa-rately below. We shall always pre-align both point setsto zero mean and unit variance before the registration. We show the CPD rigid registration on several examples,test the fast CPD implementation and evaluate its per-formance in comparison with LM-ICP [3], which is oneof the most popular robust rigid point set registrationmethods.
Rigid fish point set registration:
Fig. 6 shows severalrigid regsitration tests on 2D fish point sets. In Fig. 6awe deleted non-overlapping parts in both point setsand set w = 0 . , where w is a weight of the uniformditribution that accounts for noise and outliers. In Fig. 6bwe corrupted the point sets by outliers. We generateoutliers randomly from a normal zero-mean distribution.CPD demonstrates robust and accurate performance inall examples. Fig. 6c demonstrates a challenging exam-ple, where both point sets have missing points and are a)b) Initialization Iteration 10 Iteration 20 Iteration 30 Result (iteration 50) Fig. 7. 3D bunny point set rigid registration examples. We align Y (blue circles) onto X (red dots). The columns showthe iterative alignment progress. We initialized one of the point sets with 50 degree rotation and scaling equal . a)Registration of the point sets with missing points ( w = 0 . ); b) A challenging example of CPD rigid registration withmissing points, outliers and noise. CPD shows robust and accurate registration result in all experiments. E rr o r CPDLM−ICP0 0.05 0.1 0.15 0.200.020.040.060.080.10.120.140.160.18 Noise E rr o r CPDLM−ICP
Performance 0.04 Noise STD 0.12 Noise STD 0.2 Noise STD
Fig. 8. A comparison of CPD and LM-ICP rigid registration performances with respect to noise in the X (first row) andthe Y point sets (second row). We align Y (blue circles) onto X (red dots). The columns 2,3 and 4 show the examplesof initial point sets for different random noise stds added to the point set positions. The first column shows the error inestimating the rotation matrix for CPD (blue) and TPS-RPM (red). CPD outperforms LM-ICP in all cases. corrupted by outliers. The most challenging here is thatwe biased the outliers to the different sides of fish pointsets. We were able to register such point sets only byfixing the scaling to be constant (estimating rotation andtranslation only). CPD demonstrates accurate and robustregistration performance. Rigid bunny point set registration:
We test 3D rigidpoint sets registration on the Stanford “bunny” dataset [50]. We use a subsampled bunny version of × points. In Fig. 7a, we have deleted the front and backparts of the bunny point sets. In Fig. 7b, we have addedrandom outliers to different sides of the point sets. Weset w = 0 . . CPD registration is accurate and robust inall examples.We compare the CPD rigid algorithm to the LM-ICPmethod [3], a robust version of ICP. Fig. 8 shows theperformance of CPD and LM-ICP with respect to noisein the point sets. We align the Y point set (blue circles)onto the X point set (red dots). We set w = 0 . . Theknown initial rotation discrepancy between the point sets is degrees. The first and second rows shows thealignment performance when a random noise is addedto the X and Y point set positions respectively. We usea norm of the difference between the true and estimatedrotation matrix as an error measure. A few initial pointsets examples with different noise std are shown in thecolumns 2, 3 and 4 of Fig. 8. For each level of the noisestds we made independent runs. The first columnplots the error values (mean and standard deviation)in the estimated rotation matrix as a function of noiselevels. The CPD rigid algorithm outperforms the robustLM-ICP method, especially when the noise is present inthe X point set.Fig. 9 shows the performance of CPD and LM-ICPwith respect to the outliers in the point sets. We adddifferent number of outliers (irrelevant random points)to the point sets. An examples of such initial point setsare shown in columns 2, 3 and 4 of Fig. 9 for 600, 1800and 3000 outlier points added respectively. The first andsecond row show the cases of outliers present in the E rr o r CPDLM−ICP0 500 1000 1500 2000 2500 3000−0.0200.020.040.060.080.10.120.140.160.18 Outliers E rr o r CPDLM−ICP
Performance Outliers 600 Outliers 1800 Outliers 3000
Fig. 9. A comparison of CPD and LM-ICP rigid registration performances with respect to outliers in the X (first row)and the Y (second row) point sets. We align Y (blue circles) onto X (red dots). The columns 2,3 and 4 show theexamples of initial point sets with different number of outliers added. The first column show the error in estimating therotation matrix. CPD outperforms LM-ICP. a)b)c) Initialization Iteration 10 Iteration 20 Iteration 40 Result (iteration 50) Fig. 10. Non-rigid CPD registration of 2D fish point sets. a) Noiseless fish point sets registration ( × points, w = 0 );b) Registration of 2D fish point set with missing points ( w = 0 . ); c) Registration of 2D fish point set in presence ofoutliers ( w = 0 . ). CPD registration is robust and accurate in all experiments. X and Y point sets respectively. CPD performs wellin all experiments, whereas LM-ICP performance is lessaccurate. Fast rigid CPD implementation:
We also test the CPDperformance with FGT implementation of the matrix-vector products. We use four Stanford bunny sets ofsizes: × , × , × and × . For eachof the cases we add a small amount of noise and outliersto both point sets, initialized them with degreesrotation and set w = 0 . . For the FGT parameters,we used “ratio of far field”= , “number of centers”= , “order of truncation”= . Table. 1 shows the registrationtime with and without FGT. The FGT implementationis significantly faster. We note that there are severaldownsides of using the FGT: a) FGT requires its ownparameter initialization; b) CPD (with FGT) aligns thepoint sets to 0.1 degree error rotation and then startsbeing unstable. This is because σ becomes small andthe FGT approximation error becomes significant. Atthis point one can either stop (the alignment already isreasonably accurate) or proceed with ICP or truncatedGaussian CPD. Deformation level E rr o r RPMCPD
Noise level E rr o r RPMCPD
Outliers to data ratio E rr o r RPMCPD (a) (b) (c)
Fig. 11. A comparison of CPD and TPS-RPM on the 2D fish point sets with respect to a) Deformation level; b) Noiselevel; c) Outliers. CPD shows more accurate registration performance compared to TPS-RPM, especially in presenceof outliers and complex non-rigid deformations.
N, M
Naive FGT × ×
11s 3s ×
4m 10s × TABLE 1The rigid CPD registration time for naive (no FGT) andFGT implementations. The FGT-based implementation issignificantly faster.
We show CPD non-rigid registration on several exam-ples, test the fast CPD implementation and evaluate CPDperformance in comparison to TPS-RPM [9], which is oneof the best performing non-rigid point set registrationmethods. We set λ = 2 , β = 2 . Non-rigid fish point set registration:
Fig. 10a showsnon-rigid CPD registration of two fish point sets withclean data. Fig. 10b is with missing points ( w = 0 . ).Fig. 10c is with both point sets are corrupted by outliers( w = 0 . ). The non-rigid CPD registration results areaccurate in all experiments.We test CPD against TPS-RPM [9] on synthenticallygenerated 2D fish non-rigid examples with respect toa) level of non-rigid deformation, b) amount of noisein the point sets locations c) number of outliers. Weset w = 0 . in all experiments. Since we know thetrue correspondences, we use the mean squared distancebetween the corresponding points after the registrationas an error measure. For each set of parameters wehave conducted runs. Fig. 11a shows the methodsperformances with respect to the level of initial non-rigid deformation between the point sets. To generatethe non-rigid transformation, we parameterize the pointsets domain by a mesh of control points, perturb thepoints and use splines to interpolate the deformation.The higher level of mesh point perturbations producethe higher deformation. CPD shows accurate registrationperformance and outperforms the TPS-RPM. Fig. 11bshows the methods performances with respect to theamount of noise. We add a zero-mean white noise withincreasing levels of stds to the point sets. Both CPD and N, M × D Naive FGT Low-rank FGT & Low-rank ×
2s 2.3s 1.7s 1.8s × × × – – 40m 10m TABLE 2Registration time required for non-rigid registration of 3Dbunny point sets. The time is shown when using onlyFGT of vector-matrix products, only low-rank matrixapproximation of Gaussian kernel matrix or both.
TPS-RPM show accurate performances. We note that,due to deterministic annealing used by TPS-RPM, itsconvergence takes significantly more iterations and time.Fig. 11a shows the methods performances with respectto the number of outliers. We add random outliers to thepoint sets and plot the registration error with respect tothe ration of number of outliers to the number of datapoints. CPD shows robust registration performance andoutperforms the TPS-RPM.
Non-rigid 3D face registration:
We show the CPD per-formance on 3D face point sets. Fig. 12a shows two 3Dface point sets related through non-rigid deformation.Fig. 12b shows two 3D face point sets point sets withadded outliers and non-rigid deformation. Non-rigidCPD registration is accurate in all experiments.
Non-rigid 3D LV point set registration:
Finally, wedemonstrate the CPD performance on non-rigid a 3D leftventricle (LV) contours segmented from 3D ultrasoundimages, using active contour based segmentation [51].Fig. 13 shows (a) two LV point sets at different timeinstances, (b) the registration result, (d) the displacementfield required for CPD alignment. That the registrationresult is accurate.
Fast non-rigid CPD implementation:
We test the com-putational time of the fast CPD non-rigid implementa-tion on several subsampled 3D Stanford bunny pointsets. We use FGT of the matrix-vector products, thelow-rank matrix approximations of the kernel matrix, orboth. We applied a moderate non-rigid deformation tothe bunny point sets. The registration time of the non-rigid CPD is shown in Table 2. We were unable to run a)b) Initialization Result Fig. 12. Non-rigid registration of 3D face point sets. a)Registration of clean point sets b) Registration of pointset with outliers. CPD shows accurate alignment. (a) Initialization (b) Result (c) Displacement
Fig. 13. Non-rigid registration of 3D left ventricle (LV)point sets. (a) two LV point sets at different time instances,(b) the registration result, (c) displacement field betweenthe corresponding points. −20 −15 −10 −5 Fig. 14. The log-plot of the eigenspectrum of the kernelmatrix G for the bunny point sets of size × . the test without the low-rank matrix approximation forthe largest bunny set ( × ), because of the largeRAM requirements to construct the kernel matrix G . Weused only leading eigenvalues and eigenvectors inall cases. Table 2 shows that the main computational bot-tleneck is in solving the linear system of equations (22),because the low-rank matrix approximation alone can re-duce the computational time significantly. Both FGT andlow-rank approximations provide further speed-up withonly moderate loss of accuracy. We note that almost of the time required to complete the CPD registrationusing the low-rank matrix approximation were requiredto pre-compute the eigenvalues and eigenvectors of the kernel matrix G .We also show the eigenvalues for a particular exampleof the bunny point set of size × in Fig. 14.Eigenvalues drops quickly below − only after 10largest eigenvalues, and drops below − after about100 eigenvalues. The approximation error of using a lowrank approximate matrix (constructed from leadingeigenvectors and eigenvalues), is only − in terms ofits Frobenius norm. ISCUSSION AND C ONCLUSION
We introduce a probabilistic method for rigid and non-rigid point set registration, called the Coherent PointDrift algorithm. We consider the alignment of two pointsets as a probability density estimation, where one pointset represents the Gaussian Mixture Model centroids,and the other represents the data points. We iterativelyfit the GMM centroids by maximizing the likelihood andfind the posterior probabilities of centroids, which pro-vide the correspondence probability. Core to our methodis to force the GMM centroids to move coherently as agroup, which preserves the topological structure of thepoint sets.Our contribution includes the following aspects. Forthe non-rigid point set registration, we formulate themotion coherence constraint and derive a solution ofthe regularized ML estimation through the variationalapproach, which leads to an elegant kernel form. CPDsimultaneously finds both the transformation and thecorrespondence between two point sets without makingany prior assumption of the non-rigid transformationmodel except that of motion coherence. For the rigidcase, we derived a closed form multidimensional solu-tion (of the M-step of the EM algorithm), which has notbeen derived exactly before. Finally, we introduced thefast CPD implementation using fast Gauss transform andlow-rank matrix approximation to reduce the computa-tional complexity of the method to as low as linear. Ontop of the computational advantage, the low-rank kernelapproximation provides more stable solutions in caseswhere the matrix G is poorly conditioned. To our bestknowledge, CPD is the only method capable of non-rigidregistration of large data sets. Both rigid and non-rigidCPD registration methods are multidimensional and canbe applied to arbitrary dimensional data sets.We estimate the GMM width, σ , within the MLformulation. We have not observed any decrease inperformance compared to the deterministic annealingapproach. Estimation σ allows to reduce the numberof free parameters and, most importantly, to significantlyreduce the number of iterations and the processing time.We have used an addition uniform distribution toaccount for noise and outliers. The weight, w , of thisdistribution provides a flexible control over the methodrobustness and allows accurate CPD performance, espe-cially in presence of severe outliers and missing points.We have tested CPD on various synthetic and realexamples and comare it to LM-ICP (in rigid case) and TPS-RPM (in non-rigid case). CPD shows robust andaccurate performance with respect to noise, outliers andmissing points. Our method is of general interest withnumerous computer vision applications. We provide theMatlab code of the CPD algorithm free for academicresearch. R EFERENCES [1] P. J. Besl and N. D. McKay, “A method for registration of 3-Dshapes,”
IEEE PAMI , vol. 14, no. 2, pp. 239–256, Feb. 1992.[2] Z. Zhang, “Iterative point matching for registration of free-formcurves and surfaces,”
IJCV , vol. 13, no. 2, pp. 119–152, Oct. 1994.[3] A. W. Fitzgibbon, “Robust registration of 2D and 3D point sets,”
Image and Vision Computing , vol. 21, pp. 1145–1153, 2003.[4] S. Rusinkiewicz and M. Levoy, “Efficient variants of the ICPalgorithm,” in
International Conference on 3D Digital Imaging andModeling (3DIM) , 2001, pp. 145–152.[5] A. Rangarajan, H. Chui, E. Mjolsness, L. Davachi, P. S. Goldman-Rakic, and J. S. Duncan, “A robust point matching algorithm forautoradiograph alignment,”
MIA , vol. 1, no. 4, pp. 379–398, 1997.[6] B. Luo and E. R. Hancock, “Structural graph matching using theem algorithm and singular value decomposition,”
IEEE Trans.Pattern Anal. Mach. Intell. , vol. 23, no. 10, pp. 1120–1136, 2001.[7] S. Gold, C. P. Lu, A. Rangarajan, S. Pappu, and E. Mjolsness,“New algorithms for 2D and 3D point matching: Pose estimationand corresp.” in
NIPS , vol. 7. The MIT Press, 1995, pp. 957–964.[8] H. Chui and A. Rangarajan, “A new algorithm for non-rigid pointmatching,” in
CVPR , vol. 2. IEEE Press, Jun. 2000, pp. 44–51.[9] ——, “A new point matching algorithm for non-rigid registra-tion,” in
CVIU , vol. 89, no. 2-3. Elsevier, Feb. 2003, pp. 114–141.[10] ——, “A feature registration framework using mixture models,”in
IEEE Workshop on MMBIA , Jun. 2000, pp. 190–197.[11] A. Joshi and C.-H. Lee, “On the problem of correspondence inrange data and some inelastic uses for elastic nets,”
IEEE Trans.on Neural Networks , vol. 6, no. 3, pp. 716–723, May 1995.[12] W. M. Wells, “Statistical approaches to feature-based object recog-nition,”
IJCV , vol. 22, no. 1-2, pp. 63–98, Jan. 1997.[13] A. D. Cross and E. R. Hancock, “Graph matching with dual stepem algorithm,”
IEEE PAMI , vol. 20, no. 1, pp. 1236–1253, 1998.[14] B. Luo and E. R. Hancock, “A unified framework for alignmentand correspondence,”
CVIU , vol. 92, no. 1, pp. 26–55, 2003.[15] G. McNeill and S. Vijayakumar, “A probabilistic approach torobust shape matching,” in
IEEE ICIP , 2006, pp. 937–940.[16] M. Sofka, G. Yang, and C. V. Stewart, “Simultaneous covariancedriven correspondence (CDC) and transformation estimation inthe expectation maximization,” in
CVPR , Jun. 2007.[17] G. L. Scott and C. Longuet-Higgins, “An algorithm for associatingthe features of two images,”
Proceedings of the Royal Society London:Biological Sciences , vol. 244, no. 1309, pp. 21–26, Apr. 1991.[18] H. Li and R. Hartley, “A new and compact algorithm for simul-taneously matching and estimation,” in
IEEE ICASSP , 2004.[19] J. Ho, M.-H. Yang, A. Rangarajan, and B. Vemuri, “A newaffine registration algorithm for matching 2D point sets,” in
IEEEWorkshop on ACV , 2007, p. 25.[20] S. Belongie, J. Malik, and J. Puzicha, “Shape matching and objectrecognition using shape contexts,”
IEEE Transaction on PatternAnalysis Machine Intelligence , vol. 24, no. 4, pp. 509–522, 2002.[21] G. E. Hinton, C. K. I. Williams, and M. D. Revow, “Adaptiveelastic models for hand-printed character recognition,” in
NIPS ,vol. 4, 1992, pp. 512–519.[22] M. Revow, C. K. I. Williams, and G. E. Hinton, “Using generativemodels for handwritten digit recognition,”
IEEE PAMI , vol. 18,no. 6, pp. 592–606, 1996.[23] G. Wahba,
Spline Models for Observational Data . Philadelphia:SIAM, 1990.[24] F. L. Bookstein, “Principal warps: Thin-plate splines and thedecomposition of deformations,”
IEEE PAMI , vol. 11, no. 6, pp.567–585, Jun. 1989.[25] Y. Tsin and T. Kanade, “A correlation-based approach to robustpoint set registration,” in
ECCV , vol. 3, 2004, pp. 558–569.[26] B. Jian and B. C. Vemuri, “A robust algorithm for point setregistration using mixture of gaussians,” in
IEEE ICCV , vol. 2,Oct. 2005, pp. 1246–1251. [27] X. Huang, N. Paragios, and D. N. Metaxas, “Shape registration inimplicit spaces using information theory and free form deforma-tions,”
IEEE PAMI , vol. 28, no. 8, pp. 1303–1318, 2006.[28] A. Myronenko, X. Song, and M. ´A. Carreira-Perpi˜n´an, “Non-rigidpoint set registration: Coherent Point Drift,” in
NIPS , 2007, pp.1009–1016.[29] A. L. Yuille and N. M. Grzywacz, “The motion coherence theory,”in
ICCV , vol. 3, 1988, pp. 344–353.[30] ——, “A mathematical analysis of the motion coherence theory,”
IJCV , vol. 3, no. 2, pp. 155–175, Jun. 1989.[31] A. Dempster, N. Laird, and D. Rubin, “Maximum likelihood fromincomplete data via the EM algorithm,”
Journal of Royal StatisticalSociety. Series B (Methodological) , vol. 39, no. 1, pp. 1–38, 1977.[32] C. M. Bishop,
Neural Networks for Pattern Recognition . OxfordUniversity Press, Inc., 1995.[33] K. Arun, T. S. Huang, and S. D. Blostein, “Least-squares fitting oftwo 3-D point sets,”
IEEE PAMI , vol. 9, no. 5, pp. 698–700, 1987.[34] S. Umeyama, “Least-squares estimation of transformation param-eters between two point patterns,”
IEEE PAMI , vol. 13, no. 4, pp.376–380, Apr. 1991.[35] A. Myronenko and X. Song, “On the closed-form solution of therotation matrix arising in computer vision problems,” OregonHealth and Science Univ., Tech. Rep. arXiv:0904.1613v1, 2009.[36] B. Luo and E. R. Hancock, “Iterative procrustes alignment withthe em algorithm,”
Image and Vision Computing , vol. 20, no. 5-6,pp. 377–396, Apr. 2002.[37] A. N. Tikhonov and V. I. Arsenin,
Solutions of Ill-Posed Problems .Washington, D.C: Winston and Sons, 1977.[38] Z. Chen and S. Haykin, “On different facets of regularizationtheory,”
Neural Computation , vol. 14, no. 12, pp. 2791–2846, 2002.[39] B. Sch¨olkopf and A. J. Smola,
Learning with Kernels . Cambridge,MA: The MIT Press, 2002.[40] F. Girosi, M. Jones, and T. Poggio, “Regularization theory andneural networks architectures,”
Neural Computation , vol. 7, no. 2,pp. 219–269, 1995.[41] G. Kimeldorf and G. Wahba, “Some results on tchebycheffianspline functions,”
Journal of Mathematical Analysis and Applications ,vol. 33, no. 1, pp. 82–95, 1971.[42] R. M. Neal and G. E. Hinton, “A view of the EM algorithm thatjustifies incremental, sparse, and other variants,” in
Learning inGraphical Models , M. I. Jordan, Ed. Kluwer, 1998.[43] J. Duchon, “Spline minimizing rotation-invariant semi-norms insobolev spaces,” in
Constructive Theory of Functions of SeveralVariables , vol. 571, 1977, pp. 85–100.[44] R. Sibson and G. Stone, “Computation of thin-plate splines,”
SIAM JSSC , vol. 12, no. 6, pp. 1304–1313, Nov. 1991.[45] L. Greengard and J. Strain, “The fast gauss transform,”
SIAMJSSC , vol. 12, no. 1, pp. 79–94, 1991.[46] C. Yang, R. Duraiswami, N. A. Gumerov, and L. Davis, “Improvedfgt and efficient kernel density estimation,” in
ICCV , 2003, p. 464.[47] G. H. Golub and C. F. V. Loan,
Matrix Computations , 2nd ed.Baltimore, MD: Johns Hopkins University Press, 1989.[48] J. K. Cullum and R. A. Willoughby,
Lanczos Algorithms for LargeSymmetric Eigenvalue Computations