Robust Multi-view Registration of Point Sets with Laplacian Mixture Model

Point set registration is an essential step in many computer vision applications, such as 3D reconstruction and SLAM. Although there exist many registration algorithms for different purposes, however, this topic is still challenging due to the increasing complexity of various real-world scenarios, such as heavy noise and outlier contamination. In this paper, we propose a novel probabilistic generative method to simultaneously align multiple point sets based on the heavy-tailed Laplacian distribution. The proposed method assumes each data point is generated by a Laplacian Mixture Model (LMM), where its centers are determined by the corresponding points in other point sets. Different from the previous Gaussian Mixture Model (GMM) based method, which minimizes the quadratic distance between points and centers of Gaussian probability density, LMM minimizes the sparsity-induced L1 distance, thereby it is more robust against noise and outliers. We adopt Expectation-Maximization (EM) framework to solve LMM parameters and rigid transformations. We approximate the L1 optimization as a linear programming problem by exponential mapping in Lie algebra, which can be effectively solved through the interior point method. To improve efficiency, we also solve the L1 optimization by Alternating Direction Multiplier Method (ADMM). We demonstrate the advantages of our method by comparing it with representative state-of-the-art approaches on benchmark challenging data sets, in terms of robustness and accuracy.


Introduction
Point set registration is a fundamental problem that has wide applications in computer vision [26], robotics [25], computer graphics [12], medical image analysis [19] and so on. With the advent of sensors such as LiDAR (Light Detection and Ranging) and depth cameras, it becomes relatively easy to capture realworld 3D scene data. However, one usually attains partial point clouds at once due to the 3D nature of the objects. To accurately reconstruct the 3D model, it is necessary to align multiple point clouds acquired from multi views of the object into a unified coordinate system.
Point set registration has been extensively studied in literature, and many efforts are devoted to the pair-wise (two sets) registration problem. Among these methods, the most classical one is the Iterative Closest Point (ICP) algorithm [18], in which the registration is decomposed as the alternative implementation of point correspondence and transformation estimation. Given a set of data points, the rigid transformation of ICP is solved by minimizing the summation of the squared distance of the closest point pairs. However, ICP tends to get trapped in a local optimum because of the hard assignment scheme, improper initialization, occlusions, and the interference of noise and outliers.
In contrast to the hard assignment scheme of ICP, previous works also explore the use of statistical models for registration, such as GMM [13], which replaces the previous binary assignment with probability. These methods formulate point set registration as a probability density estimation, where the mean of each component is initialized as the point location. In principle, probability registration can provide better estimation for convergence and geometric matching, thereby they are further extended for multi-view registration [8]. Most existing probability methods rely on GMM for registration, nevertheless, Gaussian distribution minimizes the quadratic distance between the data points and their means, which makes GMM susceptible to noise with heavy tail and sensitive to outliers [10].
In this paper, to address the aforementioned problems, we propose a novel and robust multi-view registration method based on the LMM composed by the Laplacian distribution. Due to the heavy-tail property and the sparsity-induced L 1 norm, LMM is more robust against outliers. We formulate point set registration as a likelihood estimation problem, which can be solved by EM framework. To handle the L 1 optimization, we customize the rotation estimation to a linear programming problem by exponential mapping and name it as linear programming approximation (LPA). Moreover, inspired by [4], we also use ADMM to solve L 1 optimization, which has higher efficiency without significant accuracy decreasing. We test and compare the proposed methods (including LPA and ADMM) with representative state-of-the-art approaches in terms of accuracy and robustness, and the results indicate that our methods outperform previous ones with higher accuracy and are more robust against noise and outliers.
In a nutshell, the main contributions of this paper are twofold as follows: -We propose a novel multi-view registration method for point clouds based on the sparsity-induced Laplacian mixture model. Due to the L 1 norm of Laplacian distribution, the proposed method is more robust against noise and outliers; -To handle L 1 optimization, we customize it as a linear programming problem via exponential mapping, which can be effectively solved by the interior point method. For efficiency, we further deduce it into the ADMM framework, thereby it can also be solved by the ADMM algorithm.

Related Work
We briefly review the related work of point set registration from the perspective of pair-wise registration and multi-view registration.

Pair-wise Registration
The objective of pair-wise registration is to align two point sets by solving a transformation matrix. The most popular algorithm of pair-wise registration is ICP [18], which composes two iterative steps: 1) correspondence step in which the point correspondence between two point sets is established, and 2) transformation step in which the transformation matrix based on the current correspondence is updated. ICP takes the least-squares estimator as the objective function, thereby the closed-form solution in each step can be attained through singular value decomposition (SVD). However, as the Gauss-Markov theorem [21] pointed, the least-squares estimator is sensitive to outliers. Moreover, ICP requires good initialization and fails to handle non-overlapping point sets. To improve performance, many variants of ICP have been proposed. Fitzgibbon [9] adopts M-estimator for registration error minimization and solves the problem by non-linear Levenberg-Marquardt optimization, thereby proper initialization is required. Chetverikov et al. [6] propose a trim scheme named TrICP to automatically remove non-overlapping regions for accurate registration. For efficiency, Rusu et al. [22] introduce Fast Point Feature Histograms (FPFH) to describe the local geometry around a point to reduce the computational complexity. Lei et al. [15] compute eigenvalues and normals from multiple scales and take them as local descriptors to speed up matching. However, descriptor-based methods [15,11] are susceptible to point sets with noise and low overlapping.
Recently, deep learning-based methods such as [1] are also proposed to estimate the transformation matrix. Nevertheless, the lack of comprehensive registration datasets results in a hard time for the learning methods to grasp all shape variations. Alternatively, probability registration methods adopt GMM to represent data points, such as Robust Point Matching [7], Coherent Point Drift [17], and GMM-Reg [13]. These methods either simultaneously model the source and the target point sets by GMM, and the transformation matrix is solved by minimizing the discrepancy between the two GMMs, or singly model the source point set as GMM, and then evoke maximum likelihood estimation to fit the target point set. However, GMM fails to attain accurate registration results under the contamination of noise with a heavy tail. Moreover, since GMM minimizes the quadratic distance between the points and its means, it suffers from severe outliers. In contrast, we introduce Laplacian mixture models to model data distribution, which adopt the L 1 norm for error evaluation, thereby it is more robust against outliers for registration.

Multi-view registration
Multi-view registration aims to simultaneously register multiple point sets from different views, which is often solved through sequential pair-wise registration. Transformation parameters are sequentially updated by ICP or probability methods if a new point set is added. Except for the drawbacks from pair-wise registration, these methods also suffer from error accumulation and propagation. Bergevin et al. [3] propose a star-network and sequentially put one point set in the center of it, then pair-wise registration by ICP is implemented to align the central point set with the other ones. In contrast, Williams et al. [24] simultaneously compute the correspondence between all point sets, which is timeconsuming. Mateo et al. [16] introduce the Bayesian perspective to assign different weights for different correspondences, with the detection of false correspondences. Although the registration accuracy is improved, it needs to compute many variables.
Additionally, information theoretic measures are also customized for multiview registration task. Wang et al. [23] first represent point sets by cumulative distribution function (CDF), and then minimize the Jensen-Shannon divergence between CDFs for multi-view registration. Later, Chen et al. [5] use the Havrda-Charvat divergence to evaluate the differences between CDFs, compared with previous Jensen-Shannon divergence, this method is more efficient. Nevertheless, information theory based approaches are still with low efficiency.
Recently, GMM is generalized for multi-view registration by [8] named JRMPC, which assumes data points are generated by a central GMM, and then casts the registration task as a clustering process. By this, global information of point sets is combined to avoid the error accumulation. Zhu et al. [27] propose a method named EMPMR, which assumes that each data point is generated from a GMM whose Gaussian centroids are composed of corresponding points from other point sets. Nevertheless, GMM suffers from heavy-tail noise, and is sensitive to severe outliers due to the L 2 norm. Based on JMRPC, TMM [20] achieves better robustness by replacing the Gaussian distribution with a t-distribution, however, its time consumption is relatively considerable. To address these drawbacks, we propose a novel and robust multi-view registration method based on the sparsityinduced L 1 norm, which is effectively solved by LPA or ADMM.

Multivariate Laplacian Distribution
Suppose x ∈ R d is a random variable following the multivariate Laplacian distribution, with the probability density function as where µ and b represent the mean and the scale parameter of the d-dimensional Laplacian distribution, respectively, and · 1 denotes the sparsity-induced L 1 norm, summarizing the absolute values of all elements of a vector. Previous work has shown that the Gaussian distribution is sensitive to heavy outliers because of its short tails [10,2], while the Laplacian distribution has heavier tails.

Laplacian Mixture Model
The task of multi-view registration is to align the multiple point sets in X to the same center frame. To this end, we solve the rotation matrix R i and the translation vector t i for each X i , and represent the transformed coordinates asx il = R i x il + t i . Due to the influence of noise, it is hard to hope that corresponding points in different frames will have the same coordinate after alignment. In contrast, we assume the multi-view point sets will constitute multiple clusters.
Suppose each data pointx il is generated from a unique LMM, where the Laplacian centers are composed of the corresponding points in other frames after alignment. However, it is difficult to directly get the accurate correspondence relationship between different point sets. Here we adopt the nearest neighbor search based on the kd-tree to approximate the corresponding points. In specific, for pointx of the center frame, we denote c j (x) as the nearest neighbor ofx from the j-th point set, which can be defined as: For simplicity, we use isotropic covariance b and equal membership probability for all LMM components. Thereby, the LMM of x ∈ X i is defined as Note that different from previous GMM based methods, in which a uniform distribution has to be added to account for noise and outliers, there is no need for our LMM based method, since Laplacian distribution has a heavy tail and is sparse enough to accommodate outliers. We adopt the maximum likelihood estimation (MLE) to solve unknown parameters. The log-likelihood function of observed data points is , b represents all parameters. It is intractable to directly solve the maximum likelihood solution for such a complex model. Instead, we adopt the effective EM framework by introducing latent variables for solving, as presented in the following.

Registration by EM Algorithm
To estimate the parameters by EM algorithm, we first define a set of latent variables where Z il = j means thatx il is generated from the j-th component of the LMM. Given all observed data X , model parameters can be estimated by maximizing the expectation of the loglikelihood function: Since we regard each Laplacian component equal, the prior probability P (Z; Θ) is a constant term. In addition, P (X | Z; Θ) can be derived from Eq. (3), namely, After ignoring constant terms, the objective function of Eq. (5) is reformulated as where α ilj = P (Z il = j | x il ; Θ) denotes the posterior and d denotes the data dimension (d = 3 in our case). Therefore, the multi-view registration problem is cast into a constrained optimization problem as follows: In order to maximize f (Θ) by EM algorithm, we alternatively perform E-step and M-step as follows, after which, the transformation parameters of multi-view registration are attained.

E-Step
Given data point set X and currently estimated parameters Θ (k) , E-step calculates the distribution of latent variable Z, which can be divided into two steps. The first step is to establish point correspondence. For point x il in i-th point set, it is initially transformed into the center frame with new coordinatê x i . Then the corresponding point in j-th frame can be approximated by the nearest neighbour ofx The second step updates α (k+1) ilj = P (Z il = j | x il ; Θ (k) ), which represents the posterior probability of point x il generated from j-th component of LMM. According to Bayesian formula, we have where β ilj denotes the probability density of Laplacian distribution:

M-Step
M-step estimates the parameters Θ (k+1) by maximizing the expectation f (Θ) with current values α , b}, we first solve transformation parameters {R i , t i } and then solve the variance scale b. Transformation parameters are estimated by The above equation is a weighted least absolute value (WLAV) problem with SO(3) constraint. Although there is lack of a closed-form solution, we propose two methods to iteratively solve it in Section 4.3. After updating transformation parameters, we take the partial derivative of f (Θ) with respect to b and equate it to zero, then the update of b (k+1) is where N i denotes the cardinality of the i-th point set. As can be seen, b presents the weighted average of the deviations.

WLAV subproblem
Previous GMM based approaches usually need to solve a weighted least square problem, and the closed-form solution can be attained through SVD. However, due to the L 1 norm optimization, WLSV problem requires iterative solving. In this work, we propose two methods to solve this kind of problem, namely, LPA based and ADMM based methods, as presented in the following.

LPA Method
We first transform the sub-problem (11) into a canonical form: where P = [p 1 , . . . , p n ] ∈ R 3×n and Q = [q 1 , . . . , q n ] ∈ R 3×n represent the source point set and the target point set respectively, while w i represents the weight. By exponential mapping, the rotation matrix R can be written as: where [r] × represents the skew-symmetric matrix of r ∈ R 3 , expressed as To eliminate the SO(3) constraint, we linearly approximate the rotation matrix R by neglecting the higher-order terms of (14). Then the deviation can be approximated as Then the objective function can be reformulated as: where x denotes [r, t] T . Let u = |Ax − b|, the original WLAS problem is transformed into a linear programming problem: which can be efficiently solved by interior point methods [14]. Moreover, we deduce another method to solve L 1 optimization based on the ADMM in the following. We name the two optimizations as Ours-LPA and Ours-ADMM, respectively, and compare them with previous approaches in Section 5. (13) can be reformulated as:

ADMM Method Considering the constraint
Then the augmented Lagrangian function of ADMM is where we replace w i (Rp i + t − q i ) by s i to simplify the formula. Due to the space limit, we directly present the iterative steps as where sub-problem (21) can be solved efficiently by the following shrinkage operator: while sub-problem (22) can be solved by SVD.

Experiments
In this section, we compare the performance of the proposed method with three representative state-of-the-art approaches for multi-view registration, namely, JRMPC [8], TMM [20], and EMPMR [27]. The implementation of all compared methods is publicly available. All experiments in the following are performed on a laptop with a 6-core 2.2GHz Intel CPU and 16GB RAM.

Data Sets and Evaluation Measure
We use six 3D data sets from the Stanford 3D Scanning Repository 5 (Bunny, Buddha, Dragon, and Armadillo) and the AIM@SHAPE Repository 6 (Bimba and Olivier hand) for test. To quantitatively evaluate the registration performance of different methods, we compute the error by where denote the ground truth and the estimated rigid transformation, respectively. Table 1. Statistics of registration errors of all compared methods, where the red and blue fonts indicate the best and the second-best performance for each metric. The proposed method Ours-LPA attains the overall best performance.

Methods
Armadillo

Comparison Results
We first downsample each data set to 4,000 points, and crop them along the xyplane to generate the missing overlapping, then we rotate them around the x, y and z-axes with the rotation angle uniformly distributed between [−20 • , 20 • ]. Moreover, we add Gaussian noise to each data point with the signal-to-noise ratio (SNR) equal to 70dB and 30% outliers following the uniform distribution.
The comparison results are reported in Table 1. We color the best in red and the second-best in blue for each metric. As observed, compared with initial errors, TMM has larger rotation deviations for all data sets, JRMPC and EMPMR also suffer from rotation estimation, such as the Bunny and the Armadillo data sets. In contrast, our proposed methods including Ours-ADMM and Ours-LPA achieve fewer rotation errors than initial cases for all data sets. Moreover, the proposed Ours-LPA attains the highest accuracy for rotation estimation than all competitors. EMPMR has the overall lowest translation error, but its rotation errors are relatively large. Our proposed method Ours-LPA attains the second highest accuracy for translation estimation, and the deviations are small enough. Thereby, the proposed method Ours-LPA achieves the overall best performance. We present several test examples in Fig. 1.
We further report the time consumption of different methods in Table 2. As observed, due to the closed-form solution of EMPMR based on the L 2 norm, it consumes relatively less time. Our proposed method Ours-ADMM has the second fastest speed, and its time consumption is quite close to EMPMR. Moreover, from Table 1, we find that Ours-ADMM has fewer rotation errors than EMPMR, meanwhile with acceptable translation deviations. Ours-LPA consumes relatively more time than Ours-ADMM, however, it is still faster than TMM and JRMPC. Therefore, to register point sets with high efficiency and comparable accuracy, we suggest Ours-ADMM for use.

Robustness against Outliers
Subsequently, we evaluate the robustness of the proposed method against outliers. We adopt the bunny data set for this purpose. We add 1% − 80% outliers to the point cloud. Besides, we contaminate the point cloud with 70dB Gaussian noise. Several examples are illustrated in Fig. 2. The test results are reported in Table 3. As can be seen, with outlier increasing, the proposed method Ours-LPA exhibits higher robustness than compared ones. EMPMR has fewer deviations at low outlier contamination, but it suffers from severe outliers (more than 30%). JRMPC also shows certain robustness against outliers, whereas its performance is unstable. In contrast, Ours-ADMM attains the quite stable performance for all outlier tests, and it has the a very similar registration accuracy to the first two winners. TMM has the largest deviations than the others, indicating its weakness in handling outliers.

Robustness against Noise
We also test the performance of the proposed methods in term of noise. To this end, we first contaminate the point cloud with 30% outliers, and then increase the noise intensity by decreasing the signal-to-noise ratio. The test results are reported in Table 4. As observed, the proposed method Ours-LPA achieves the overall best performance, and its registration errors are even 1,000 times fewer than TMM. JRMPC and Ours-ADMM attain similar results, thereby they have the second-best performance. EMPMR shows shortcomings for noise, especially for low SNR noise.

Conclusion
We have presented a novel and robust multi-view registration method for point clouds based on the Laplacian mixture models (LMM). We adopt Laplacian distribution to represent each data point, and then cast the multi-view registration task as a density estimation problem, which can be efficiently solved through the expectation-maximization framework. Due to the heavy tail and the sparsity-induced L 1 norm, LMM is more robust against GMM. To solve the L 1 problem, we deduce our objective function into two optimization paradigms, namely, linear programming and ADMM. We test the proposed methods on challenging data sets with contamination of noise and outliers, and compare it with three representative state-of-the-art approaches, and results demonstrate the salient advantages of the proposed method: more robust against noise and outliers, as well as higher accuracy.