I Introduction
The registration of point sets is an essential methodology in computer vision, computer graphics, robotics, and medical image analysis. The vast majority of existing techniques solve the pairwise (two sets) registration problem, e.g.,
[1, 2, 3, 4, 5, 6], while the multiple-set registration problem has comparatively received less attention, e.g., [7, 8, 9]. Solutions to this problem are often approximated by repeatedly solving for pairwise registration, either sequentially [10, 11, 12], or via a one-versus-all strategy [13, 14, 15].Independently of the particular two-set registration algorithm that is used, the above mentioned approximate solutions have their own limitations. On the one hand, sequential strategies suffer from the well known drift accumulation owing to the chain-based optimization, i.e., sequential registration between pairs of point sets. On the other hand, one-versus-all strategies lead to a biased estimator since the registration is governed by a single reference set. In addition, both strategies lack closed-loop information and one needs to further consider this constraint. Therefore, an unbiased solution that treats all the point sets on an equal footing and that implicitly enforces a loop constraint is particularly desirable.
Such an unbiased solution is targeted by motion averaging approaches that build on pairwise registration schemes and aim to evenly distribute the total error across the network of point sets, either as a post-processing step [16] or as an over-successive registration between pairs of point sets [9]. We rather aim to jointly register all the point sets and not re-distribute the error from a pairwise registration. To this end, we propose a generative approach to the joint registration of multiple point sets. An arbitrary number of point sets, observed from different sensor locations, are assumed to be generated from a singleGaussian mixture model (GMM). The problem is cast into a data clustering problem which, in turn, is solved via maximum likelihood and leads to an expectation maximization (EM) algorithm, whereby both the mixture and registration parameters are optimally estimated. We present batch and incremental EM algorithms: both can deal with point sets of different cardinalities and contaminated by noise and outliers.
Pairwise probabilistic registration methods constrain the GMM means to coincide with the points of one set, e.g., [4, 5]. Note that such a coincidence is inherently problematic, as long as both point sets are noisy and may include outliers. Even if one includes a uniform component in the mixture to deal with outliers [17], one of the sets is supposed to be “perfect”. Instead, the means of the proposed formulation are not tight to a particular set: they result from fitting a mixture model to the data sets that are appropriately rotated and translated. In addition to registration, this also achieves scene reconstruction, since the cluster means may be viewed as the scene model. The proposed formulation implicitly enforces a closed-loop constraint. In other words, the proposed model assumes a star network topology, while the pairwise registration schemes assume a ring topology or a fully connected network
This article is an extended version of [18]. Several aspects of the proposed model are discussed into more detail, namely initialization, behavior, complexity and advantages over existing methods. In addition to the batch EM described in [18], we introduce an incremental version of EM, which solves the parameter estimation problem more efficiently at the cost of less accurate results. Experiments with novel datasets and benchmarks with several recent methods are included as well.
The remainder of this paper is organized as follows: Section II discusses the related work. Section III formulates the problem in a generative probabilistic framework. Section IV describes the batch EM together with an algorithm analysis while the incremental version is presented in Section V. Section VI describes in detail various initialization procedures. Section VII presents the experimental results and Section VIII concludes the paper.^{1}^{1}1Matlab code, datasets and videos are available at https://team.inria.fr/perception/research/jrmpc/.
Ii Related work
The two-set registration problem is usually solved by ICP [1, 19] or by one of its variants [20, 2, 3, 21, 22]
. While ICP alternates between hard assignments and transformation estimation, more sophisticated registration approaches replace the binary assignments with probabilities
[23, 24, 25, 4, 5]. Nevertheless, whether based on ICP or on soft assignments, these methods consider one set as the “model” and the other set as the “data”, thus leading to solutions that are biased as long as both sets may contain noise and outliers. Alternatively, [6, 26] consider two Gaussian mixtures, one per point set, and the rigid transformation is applied to one of these mixtures. This leads to a non-linear optimization problem.Multiple point set registration is often addressed by a sequential pairwise registration strategy [10, 11, 19], in particular when an online solution is required. Whenever an additional set is available, the model set is updated using either an ICP-like or a probabilistic scheme. Apart from the drawbacks associated with pairwise registration, this incremental mode of operation is subject to error propagation, while it fails to close any existing loop. As for offline applications, several approaches have been proposed, being mostly based on the underlying network, a.k.a. viewgraph, defined by the sets (represented as nodes) and their relative overlap (represented as edges). The majority of these methods initialize the poses via a pairwise registration.
The first solution for the problem in question was proposed in [13], where the sets are organized in a star-shaped network with one of the sets in the center, and such that any two sets are linked via two edges, hence by combining two rigid transformations. An algorithm computes the transformations incrementally based on a point-to-plane ICP algorithm [19]. [27] proposed to accelerate this algorithm by allowing incremental updates once pairwise registration within the loop has been performed. [15] starts with pairwise registrations to build the set graph, while a global registration step eliminates inconsistent matches and leads to the model graph whereby poses are provided. All these methods, however, consider in practice one set as reference, thus favoring a biased and non-symmetric solution.
Alternatively, [28]
proposes a method to register multiple range images based on shape modeling: a point-to-surface distance is defined, the signed distance field, and the algorithm alternates between alignment and registration. Measurement errors and wrong correspondences are handled by a robust loss function. This method is well suited for dense range data since a surface representation is necessary.
Other methods consider known and fixed correspondences across multiple sets, thus updating only the transformations to balance the global error over the viewgraph [7, 14, 31, 32, 16, 33]. The main principle of these methods is that transformations along a network cycle ideally compose to the identity transformation. The cycles may refer to either minor loops between two adjacent sets or a larger cycle over the network.^{2}^{2}2When a spanning tree is used, an unused edge is added to obtain a cycle. Provided an approximate alignment, the goal is to minimize the on-cycle accumulated error from registering pairs of relevant (nearby) views. However, when data are ignored, a low inconsistency between coordinate frames does not necessarily mean better surface registration, in particular when good initialization is not available. As a consequence, these methods just “spread” any existing bias across the network without any correspondence refinement.
An alternative approach consists of considering a dense sequence of depth images and of estimating slight misalignments between these images. If the images are linearly correlated, the image alignment can be obtained via low-rank decomposition of a large matrix which has as columns the misaligned images. This formulation has been successfully applied to 2D image alignment [29] and extended to align images gathered with a RGB-D sensor [30]. The method is however limited to small camera motions such as to preserve the necessary condition that the images are linearly correlated. Our method addresses a different scenario, because it can handle large camera displacements and it does not necessitate dense RGB-D data. We conclude that our method and [30] are complementary.
Several recent methods built on the motion averaging principle introduced in [34] and based on rotation averaging [35]. Provided the view network, [9] suggests a motion averaged ICP algorithm. This algorithm alternates between the correspondence step and a double motion update. Any edge of the network implies an ICP run that updates a relative motion, and the redundancy information from all the relative motions in turn lead to a new global motion (one transformation per set) through the Lie-algebraic motion averaging principle. Then, the global motion information is back propagated in order to re-update the relative motions in a globally consistent manner. Again, the main assumption behind averaging is that traversing a cycle on the view network implies no motion. However, point correspondences are also updated here. [36] adopts the same technique but it employs trimmed-ICP [37] to compute pairwise motions. Note that an existing closed-loop may need to be pre-defined or pre-detected.
Probabilistic methods have been also proposed. As in [6], [8] represents each point set as a GMM and the non-rigid transformations are applied to cluster centers rather than to raw points. The model parameters are estimated by minimizing the Jensen-Shannon divergence of multiple densities and a probabilistic mean shape is built (as a by-product) from the convex combination of the aligned sets. This method vitally depends on each set’s clusters, thus requiring highly and well structured point sets with no outliers. More closely to our method, [38] proposed an EM algorithm that alternates between the reconstruction of the object’s mean shape and the registration between the sets and this shape. Despite the same principle, i.e., an emerging mean shape generates the sample sets, [38] considers given correspondences as well as several simplifications. KinectFusion [12] would roughly fall into this category owing to its model-to-frame registration strategy. Unlike these approaches, [39] and [40] build on pairwise registrations. The former generalizes [4] to align multiple super-resolved depth images by jointly optimizing many pairwise alignments along with compensating for pixel-dependent systematic bias. The latter extends the objective function of [7] and [32] by considering correspondences as missing data that are inferred along with the pairwise transformations in an EM fashion. Recently, [41] proposed an extension of [18] that integrates RGB information which enables better initial matches. In a large-scale outdoor context, [42] exploits positioning and map data to a pre-detect closed loop, while it proposes a multiple point set extension of [21] to simultaneously refine the intra-loop poses of the range sensor.
Iii Problem formulation
Let be data points that belong to point set and let be the number of point sets. We denote with the union of all these sets. A rigid transformation
, i.e. a rotation matrix and a translation vector, maps
from a set-centered frame to a model-centered frame, such that all the points form all the sets are expressed in the same coordinate frame. The objective is to estimate the data-set-to-model-set transformations under the assumption that the observed points are generated from the same mixture model(1) |
where is a rotation matrix and is a translation vector such that , are the mixing coefficients with , and are the mean vectors and covariance matrices respectively, and
is the uniform distribution parameterized by the volume
of the 3D convex hull encompassing the data [5]. We now define as the ratio between outliers and inliers(2) |
This allows to balance the outlier/inlier proportion via a judicious choice of . To summarize, the model parameters are
(3) |
This problem can be solved using an EM algorithm. We define hidden variables , such that means that observation is assigned to the -th component of the mixture, and we seek to estimate the parameters by maximizing the expected complete-data log-likelihood given the observed data
(4) |
Iv Batch Registration
Assuming that the observed data are independent and identically distributed, it is straightforward to write (III) as
(5) |
where are the posteriors. By replacing the standard expressions of the likelihoods [43] and by ignoring constant terms, (5) can be written as an objective function of the form
(6) |
where denotes the determinant and . The model is farther restricted to isotropic covariances, i.e. , since this leads to closed-form maximization solutions for all the model parameters (3), while non-isotropic covariances lead to a more complex convex optimization problem with no significant gain in accuracy [5]. Particular care must be given to the estimation of the rotation matrices, namely a constrained optimization problem
(7) |
which can be solved via EM. Notice that the standard M steps for Gaussian mixtures are augmented with a step that estimates the rigid transformation parameters. We will refer to this algorithm as joint registration of multiple point clouds (JRMPC). The batch version will be referred to as JRMPC-B and it is outlined in Algorithm 1. This leads to a conditional maximization procedure [44]. Each M-step first estimates the transformation parameters, given the current responsibilities and Gaussian mixture parameters, and then estimates the new mixture parameters, given the new transformation parameters. It is of course possible to adopt a reverse order, in particular when rough rigid transformations are available. However, the proposed order does not assume such prior information.
Iv-a E-step
The posterior probability of point
to be associated with cluster , e.g. an inlier, is(8) |
where accounts for the uniform component in the mixture, and with the notation:
(9) |
Therefore, the posterior probability of being an outlier is simply given by
.
As shown in Algorithm 1, the posterior probability at the -th iteration, , is computed from (8) using the parameter set .
Iv-B M-rigid-step
This step estimates the rotations and translations that maximize , given current values for , , , and . Notice that this estimation can be carried out independently for each set . By setting the GMM parameters to their current values, we reformulate the problem of estimating the rotations and translations. The rigid transformation parameters that maximize can be estimated from the following constrained minimization
(10) |
where is a diagonal matrix with entries , , is a vector of ones, denotes the Frobenius norm, and , where is the weighted average of the -th point set assigned to the -th mixture component
(11) |
The minimization (10) can be solved in closed-form and is a weighted version of the solution [45]. The optimal rotation matrices are
(12) |
where and
are the left and right matrices respectively, obtained from the singular value decomposition of matrix
, with(13) |
is a projection matrix and . Once the optimal rotation matrices are estimated, the optimal translation vectors are easily computed with
(14) |
Note that each rigid transform aligns the GMM means with K virtual points (one virtual point per component). Therefore, the proposed method can deal with point sets of different cardinalities and the number of components in the mixture, , can be chosen independently of these cardinalities. This is an important advantage over pairwise registration methods that assume that the cardinalities of the two point sets must be similar.
Iv-C M-GMM-step
Given rigid transformation estimates and posterior probabilities, one can use standard optimization techniques to compute the optimal means and covariances:
(15) | ||||
(16) |
where is a small scalar to avoid singularities. As for the priors, we introduce a Lagrange multiplier to take into account the constraint . This leads to the following dual function
(17) |
and its optimization yields
(18) | ||||
(19) |
with and . Note that if , which means that there is no uniform component in the mixture, then , which is in agreement with [43].
Iv-D Algorithm Analysis
The leading complexity of JRMPC-B is owing to E-step and equation (9). If is the average cardinality of a point set, the complexity can be written as . Typically, owing to underlying clustering while could be close to or even greater than when many non-ovelapping sets cover a large volume.
The proposed algorithm has a number of advantages over pairwise registration methods. Such methods, e.g. [7, 9, 36], are intrinsically more time-consuming than joint registration because one has to consider all point-set-to-point-set combinations. Either using EM or ICP when registering each pair of sets, the evaluation of all point-to-point distances is needed. Such a strategy requires operations in principle. This complexity can be decreased by structuring the data, e.g. using KD-trees, at the cost of data structure building. Approximate solutions, e.g., sequential or one-versus-all approaches, consider pairs of sets, thus requiring operations at the expense of performance.
Another important difference between joint and pairwise registrations is that the former puts all the point sets on an equal footing and registration is truly cast into clustering, i.e. (1), whereas the latter performs an unbalanced treatment of the point sets, i.e. one set constitutes the data and the other set constitutes the model. More precisely, when EM is used for registering pairs [4, 5], the generative model in (1) is replaced with , where belongs to point set (the data), belongs to point set (the model) and is the rigid alignment that maps onto . Hence, in the joint case, the mixture is modeled by a set of free parameters, while in the pairwise case, the means directly depend on the rigid parameters. The immediate consequence is that noisy points or outliers that may be present in the data will propagate via this dependency and will give rise to bad means in the mixture. When ICP is used, again one of the sets is identified with the model.
The minimal configuration required by the proposed method consists of two sets with at least three overlapping points. The algorithm can be applied to a large number of point sets. However in this case, the computation time increases linearly with on the premise that does not depend on the number of point sets to be aligned. For this reason, it is interesting to provide an incremental version of the algorithm, on the following ground: once point sets () are aligned with the JRMPC-B algorithm described above, new sets can be added incrementally (one at a time) and aligned with the current model, at a lower computational cost than the batch algorithm, using update formulae for mixture parameters. The incremental version of the algorithm is described in detail in the next section.
V Incremental Registration
The incremental version of the proposed registration method, referred to as JRMPC-I, is outlined in Algorithm 2. This algorithm considers the new -th point set to be aligned with already registered sets. The latter and the corresponding transformations are not updated and therefore are not used as input and output arguments. The JRMPC-I algorithm starts with computing the responsibilities , it then estimates the rigid transformation that aligns the set with the already aligned sets, and , and finally updates the mixture parameters. This process can be optionally repeated for a few iterations (). While the mixture parameters are initialized with those previously calculated, the integration of the new set requires initialization of and , referred to as and in Algorithm 2. One possible strategy that has been successfully used with a moving RGB-D camera, e.g., [12], is to initialize the rigid transformation with and . In the more general case, one can use the initialization strategy discussed in Section VI.
We denote with the GMM parameters estimated with JRMPC-B, where the notation denotes the sets from to . The incremental registration algorithm proceeds iteratively. The E-step computes the responsibilities associated with the -th set:
(20) |
The M-rigid-step uses equations (11)-(14) to calculate and in closed-form. This rigid transformation aligns the -th set with the GMM means that explain the previously aligned sets, hence the joint alignment of all the sets. The M-GMM-step updates the means, covariances and priors:
(21) | ||||
(22) | ||||
(23) |
with
The number of iterations that JRMPC-I needs to converge depends on its initialization. It was noticed that a small number of iterations are sufficient when data are gathered from a smoothly moving camera. Once a few sets have been integrated with JRMPC-I, it may be useful to run JRMPC-B in order to obtain a globally optimal alignment and to reject the outliers. Also, it is worthwhile to remark that JRMPC-I is not meant to grow the model, i.e. it is not designed to increase the number of components of the Gaussian mixture as point sets, possibly with no overlap, are incrementally added. JRMPC-I should be merely used when an efficient algorithm is needed. In the particular case of a large number of sets, e.g., depth sequences, a temporally hierarchical scheme that benefits from both versions is recommended to cope with the large memory requirements.
Vi Initialization
It is well known that initialization plays a crucial role in EM procedures. Therefore, we discuss here initialization options well suited for point set registration. We assume no prior information about the position and orientation of the camera(s) with respect to the scene. However, information such as the calibration parameters of a network of static cameras, or transformations between pairs of point sets, could be used if available. We also assume that there is sufficient overlap between pairs of point sets. The sensitivity of our method to the amount of overlap is tested an analyzed in Sec. VII.
When the point sets have a sufficient joint overlap, the translation vectors can be initialized by centroid differences, i.e., , where is the centroid of the cluster centers and the centroid of the th set. If the point sets suffer from strong artifacts, e.g., flying pixels, the difference of medians can be used instead. Rotation matrices can be simply initialized with . Instead, when many non-overlapping pairs exist, a pairwise registration is preferred, i.e., the minimum number of pairs that leads to a rough global alignment can be registered beforehand.
Several strategies may be adopted for initializing the mixture parameters. One way to do it is to initialize the means with the points of one set. Another way is to distribute the means on the surface of a sphere that encompasses the convex hull of the point sets already centered at . Concerning the variance, we found that starting with a high value yields very good results and that the variances quickly converge to the final values. Our algorithm converges much faster than EM algorithms that adopt a deterministic annealing behavior, i.e. the variance is decreased according to an annealing schedule. Finally the priors are initialized with , where we remind that is the number of Gaussian components. While update formulae for the priors are provided with both our algorithms, in practice it was found that keeping the priors constant affect neither the convergence nor the quality of the registrations. Notice that any rough pre-alignment of the sets results in a very good initilazation of the mixture parameters, that is, the means can be intitialized from re-sampling the registered set while the variances can be intiialized such that each cluster encompasses a sufficient number of points.
In order to choose the number of components, we propose the following empirical strategy. If the cardinalities of the point sets are similar, one can use (recall that is the average number of points in a set). However, may be chosen to be smaller than if the sets highly overlap, or larger if there are many non-overlapping sets. One should notice that the number of components in the mixture merely depends on the data and on the application at hand. Experimentally we found that yields excellent alignment results in the presence of dense depth data, as provided by depth sensors.
(a) | (b) | (c) |
(a) noise | (b) noise+outliers | (c) noise | (d) noise+outliers |
Vii Experiments
In this section, we test and benchmark the proposed algorithms with widely used and publicly available 3D data, as well as with time-of-flight (TOF) and structured-light (Kinect) data. First, we compare the proposed algorithm with pairwise registration methods, which illustrates the behavior of the algorithm in comparison with other algorithms and in particular its robustness to noise and to outliers. Second, we compare our algorithm with recently proposed joint-registration algorithms. Third, we test and evaluate the best performing algorithms on challenging TOF data and on data captured with a moving Kinect sensor.
Vii-a Simulated Data
Vii-A1 Comparison with pairwise registration algorithms
We use 3D models from the Stanford 3D scanning repository^{3}^{3}3https://graphics.stanford.edu/data/3Dscanrep/, i.e., “Bunny”, “Lucy” and “Armadillo”, and we proceed as follows in order to synthesize multiple point sets from different viewpoints. The model is shifted around the origin, the points are downsampled and then rotated in the -plane; points with negative coordinates are rejected. This way, only a part of the object is viewed in each set, the point sets do not fully overlap, and the extent of the overlap depends on the rotation angle, as in real scenarios. It is important to note that downsampling differs over the sets, such that different points are present in each set as well as different cardinalities (from the range ) are obtained. We add Gaussian noise to point coordinates based on a predefined signal-to-noise ratio (SNR), and more importantly, we add outliers to each set which are uniformly distributed around five randomly chosen points of the set. A tractable case of registering four point sets () is considered here, the angle between the first set and the other sets being , and respectively. We include JRMPC-I in the latter experiments where more point-sets are registered.
Bunny | Lucy | Armadillo | |||||||
---|---|---|---|---|---|---|---|---|---|
ICP [1] | |||||||||
GMMReg [6] | |||||||||
CPD [4], ECMPR [5] | |||||||||
SimReg [7] | |||||||||
JRMPC-B |
respectively, while the third column shows the standard deviation of the two errors (
, outliers).For comparison, we consider the following baselines that follow the one-vs-all approach: ICP [1], CPD [4], ECMPR [5], GMMReg [6]. In addition, we include a sequential version of ICP (seqICP) and a modification of [7], abbreviated here as SimReg. Unlike the original version, the latter allows updating the matches at each iteration. Recall that CPD is exactly equivalent to ECMPR when it comes to rigid registration.^{4}^{4}4CPD considers common variance for all components, while each component has its own variance with ECMPR; that latter case is considered here. As showed in [6], Levenberg-Marquardt ICP [2] performs similarly with GMMReg, while [8] shows that GMMReg is superior to Kernel Correlation [3]. As a consequence, we implicitly assume a variety of baselines. All the competitors employ registrations between the first and rest sets, while SimReg considers all the pairs of (overlapping) sets.
To evaluate the performance. we use the root of the mean squared error (RMSE) of the rotation parameters averaged by the number of sets. For all algorithms, we implicitly initialise the translations by transferring the centroids of the point clouds into the same point, while identity matrices initialize the rotations. GMMReg and SimReg are kind of favored in the comparison, since the former benefits from a two-level optimization (the first level initializes the second one) while the latter starts from the point where the pairwise ICP ends. Notice that the proposed method provides a transformation for every point set, while ground rotations are typically expressed in terms of the first set. Hence, the product of estimations is compared with the ground rotation , i.e., the error for the -th set is .
JRMPC starts here from a completely unknown GMM where the initial means are distributed on a sphere that spans the convex hull of the sets. The variances are here initialized with the median distance between and all the points in . We found that updating the priors does not drastically improve the registration. We therefore keep them constant and equal to during EM, while is chosen to be the volume of a sphere whose radius is ; the latter is not an arbitrary choice because the point coordinates are normalized by the maximum distance between points of the convex hull of . The number of the components, , is here equal to of the mean cardinality. We use iterations for all algorithms while GMMReg performs and function evaluations for the first and second optimization levels respectively. However, the current authors’ implementation allows to extract the parameters after the latest evaluation.
(a) | (b) | (c) |
(d) | (e) | (f) |
(g) | (h) | (i) |
Fig. 2 shows the final log-RMSE, averaged over realisations and all views, as a function of outlier percentage for each 3D model. Apparently, ICP and SimReg are more affected by the presence of outliers owing to one-to-one correspondences. CPD and GMMReg are affected in the sense that the former assigns outliers to any of the GMM components, while the latter may merge outliers into clusters. The proposed method is more robust to outliers and the registration is successful even with densely present outliers. The behavior of the proposed algorithm in terms of the outliers is discussed in detail below and showed on Fig. 4. To visualize the convergence rate of the algorithms, we show curves for a challenging setting ( and outliers). Regarding GMMReg, we just plot a line that shows the error in steady state, since the author’s implementation allow to extract the final parameters only. There is a performance variation as the model’s surface changes. “Lucy” is more asymmetric than “Bunny” and “Armadillo”, thus a lower floor is achieved. Unlike the competitors, JRMPC-B may show a minor perturbation in the first iterations owing to the joint solution and the initialization of the means and the variances.
It is also important to show the estimation error between sets whose geometric relation is not directly estimated. This also shows how biased each algorithm is. Based on the above experiment (SNR=, outliers), Table I reports the average rotation error for the pairs and , as well as the standard deviation of these two errors as a measure of bias. All but seqICP do not estimate these individual mappings alone. The proposed scheme, not only provides the lowest error, but it also offers the most symmetric solution.
A second experiment evaluates the robustness of the algorithms in terms of the rotation angle between two point sets, that is, the extent of their overlap. This also allows us to show how the proposed algorithm deals with the simple case of two point sets. Recall that JRMPC-B does not reduce to CPD/ECMPR in the two-set case, but it still computes the poses of the two sets with respect to the “central” GMM. Fig. 3 plots the average RMSE over realizations of ”Lucy“ and “Armadillo”, when the relative rotation angle varies from to . As for an acceptable registration error, the proposed scheme achieves the widest and shallowest basin for “Lucy”, and competes GMMReg for “Armadillo”. Since “Armadillo” consists of smooth and concave surfaces, the performance of the proposed scheme is better with multiple point sets than the two-set case, hence the difference with GMMReg. The wide basin of GMMReg is also due to its sophisticated initialization.
As mentioned, a by-product of the proposed method is the reconstruction of an outlier-free model. In addition, we are able to detect the majority of the outlying points based on the variance of the component they most likely belong to. To show this effect, we use the results of one realization of the first experiment with outliers. Fig. 4 shows in (a) and (b) two out of four point sets, thereby one verifies the distortion of the point sets, as well as how different the sets may be, e.g., the right hand is missing in the first set. The progress of estimation is shown in (d-f). Apparently, the algorithm starts by reconstructing the scene model (observe the presence of the right hand). Notice the size increment of the hull of the points , during the progress. This is because the posteriors in the first iteration are very low and make the means shrink into a very small cell. While the two point sets are around the points and , we build the scene model around the point . The distribution of the final deviations is shown in (c). We get the same distribution with any model and any outlier percentage, as well as when registering real data. Although one can fit a pdf, e.g., Rayleigh, it is convenient enough here to split the components using the threshold , where . Accordingly, we build the scene model and we visualize the binary classification of points . Apparently, whenever components attract outliers, even not far from the object surface, they tend to spread their hull by increasing their scale. Based on the above thresholding, we can detect such components and reject points that are assigned with high probability to them, as shown in (g). Despite the introduction of the uniform component that prevents the algorithm from building clusters away from the object surface, locally dense outliers are likely to create components outside the surface. In this example, most of the point sets contain outliers above the shoulders, and the algorithm builds components with outliers only, that are post-detected by their variance. The integrated surface is shown in (h) and (i) when “bad” points are automatically removed. Of course, the surface can be post-processed, e.g., smoothing, for a more accurate representation, but this is beyond of our goal.
(a) MAICP | (b) MATrICP | (c) JRMPC-B | (d) JRMPC-I |
Vii-A2 Comparison with joint registration algorithms
We here compare our method with the joint registration algorithms of [9] and [36].^{5}^{5}5The code was kindly provided by the authors. Recall that both rely on the motion averaging strategy using the ICP and the trimmed-ICP algorithm, respectively, hence abbreviated as MAICP and MATrICP. According to the literature, MAICP and MATrICP seem to outperform the methods of [31, 27, 7, 13]. The method of [7] is also included here as a baseline that considers fixed matches between the sets, and is referred to as multi-view ICP (MV-ICP). As mentioned above, MV-ICP considers all the pairs of overlapping views. While [8] generalizes GMMReg [6] for multiple point-sets, the authors provide the code for two-set case only.
Raw-data | Initialization | MV-ICP [7] | MAICP [9] | MATrICP [36] | JRMPC-B | JRMPC-I | |
---|---|---|---|---|---|---|---|
Bunny | 3.45 | 2.10 | 1.54 | 0.95 | 0.27 | 0.37 | 0.69 |
Dragon | 7.28 | 4.37 | 3.75 | 1.95 | 0.62 | 0.47 | 0.73 |
Happy Buddha | 10.77 | 3.18 | 2.45 | 0.64 | 0.43 | 0.36 | 0.77 |
For consistency reasons, the experimental setup of [9] is adopted, i.e., the point-sets have been roughly pre-aligned using a standard pair-wise ICP scheme. The error metric is the angle (in degrees) obtained from the composition of true and inverse estimation averaged over all point-sets, that is, it should ideally vanish. As with [9], we use the “Bunny”, “Dragon” and “Happy Buddha” models from Stanford scanning repository owing to the availability of the ground truth motions. While “Bunny” is asymmetrically captured from viewpoints, the last two sets contain scans from evenly spaced view angles (every degrees). To get the point-sets, true transformations first apply to the sets and then, we deform each set by a random yet known transformation. Finally, we down-sample the point-sets so that the cardinalities vary from to points. Unlike [9] and [36], we also evaluate the registration performance, when the point-sets have been further perturbed by noise. We deliberately avoid adding outliers since any mis-registration in the initialization step would make the motion averaging methods completely fail.
cross-section area | pairwise ICP | MAICP | MATrICP | JRMPC-B | JRMPC-I |
JRMPC-B | JRMPC-I |
Raw-data | Initialization | MV-ICP [7] | MAICP [9] | MATrICP [36] | JRMPC-B | JRMPC-I | |
---|---|---|---|---|---|---|---|
Bunny | 3.45 | 2.42 | 2.66 | 2.87 | 2.37 | 1.07 | 1.41 |
Dragon | 7.28 | 7.34 | 7.37 | 3.28 | 1.55 | 0.64 | 0.89 |
Happy Buddha | 10.77 | 6.88 | 6.86 | 4.13 | 1.92 | 1.18 | 1.69 |
Both JRMPC-B and JRMPC-I consider the same number of components () while the initial centers are randomly selected points from roughly aligned sets. For JRMPC-I, when a new set appears, components are rejected and re-initialized with points from the new set. This is to enforce the displacement of some GMM means towards the new data, as long as model growing is not considered here. Several conditions may apply to this rejection stage. Here, we first reject degenerate clusters ( (if any) and we randomly select old components to replace. One iteration of integration step and refinement iterations with JRMPC-B are allowed owing to the different viewpoints, while we let the algorithm run cycles to register the two first sets. Note that the current implementations of MV-ICP and MAICP consider the closed-loop known, that is, pairing the last with the first set, while MAICP also considers the scan boundaries known and rejects such points for potential matching. Instead, both versions of our algorithm as well as MATrICP make no use of any prior knowledge about the loop and the overlap.
seqICP | MVICP | MAICP | |||
MATrICP | JRMPC-B | JRMPC-B after outlier removal |
Table II shows the registration error of the methods. As expected, MV-ICP fails to provide accurate registration owing to fixed matches. The proposed algorithm along with MATrICP achieve the most accurate registration, while JRMPC-I provides results of sufficient quality. Indeed, as claimed in [36], it seems that motion averaging benefits from more robust versions of ICP. The corresponding integrated models of the best performing algorithms are shown in Fig. 5. Likewise, MAICP is less accurate while JRMPC schemes and MATrICP provide very good reconstructions.
Fig. 6 shows cross-sections of the reconstructions obtained by the proposed and motion-averaging algorithms (best viewed on-screen). The more “clean” and solid the sketch, the more accurate the alignment. The algorithms achieve to correct the initial sketch of the pair-wise ICP method. A detailed look verifies the superiority of the proposed batch method and the potential of the incremental version. Note that down-sampling makes short lines intersect in the cross-sections, even when using the ground truth motions.
Despite its incremental nature, JRMPC-I achieves comparable reconstructions and closes the loop successfully. However, the components are not distributed in the same way. Fig. 7 shows the distribution of means after running both versions of JRMPC. Despite the refinement step and the rejection stage, the means seem to remain a little biased towards initial sets, which might be problematic with long data sequences. In such a scenario, one should enforce a constraint so that new components that replace the rejected ones entirely belong to new scene surface. Note that detecting the points that may belong to the new part of the scene/object when the depth sensor is moving is easy with today hybrid sensors that deliver visual and inertial data.
Table III provides a quantitative comparison between the methods when the point-sets are further perturbed with noise of SNR=dB. As seen, the motion averaging methods seem to be more sensitive than the proposed ones. This is mainly because the GMM means get cleaned over time and the registration module in JRMPC is more robust to noise. As a consequence, even JRMPC-I outperforms the motion averaging methods. The presence of noise make the illustration of cross-sections and integrated models meaningless.
Remarkably, we experimentally found that fixing the variance for the initial iterations make JRMPC-B converge at a lower level. When the sets are roughly aligned, a fixed and reasonable value of the variance (that make each cluster include a few points) leads to better distributed means in terms of the object skeleton, which in turn lead to more accurate transformations. This is because the skeleton carries more informative points than the surface itself. Then, the update of the variance leads to better reconstruction of the object and to “safe” refinements of the rotations. From a mathematical point of view, this strategy helps avoiding local minima in the variance-rotations subspace.
Vii-B Real Data
In [18], we tested JRMPC-B along with pairwise strategies on EXBI dataset, that contains depth data captured with a time-of-flight camera rigidly attached to two color cameras. Once calibrated [46, 47], this TOF-stereo sensor provides RGB-D data. The EXBI data consist of ten point clouds gathered by manually moving the TOF-stereo sensor in front of a scene, e.g., Fig. 8. Each point cloud contains approximatively points. While JRMPC-B only uses the depth data, color information is used the final assessment and also shows the potential for fusing RGB-D data.
The comparison in [18] showed that, unlike JRMPC, all the pairwise strategies suffer from misalignments and need further processing, e.g., motion averaging.^{6}^{6}6we also refer the reader to the supplementary material of [18] Therefore, we test the performance of MAICP, MATrICP and MVICP on EXBI data-set and compare with JRMPC. SeqICP is used to roughly initialize the transformations of the point clouds.
Fig. 8 shows the front and top view of the integrated sets obtained by seqICP (initialization) as well as by MVICP, MAICP, MATrICP and JRMPC algorithms. Both versions of the proposed algorithm provide visually similar results. As verified, the motion averaging method cannot fully compensate for the misalignments of the initialization. This is shown even in front views, e.g., on the dummy head area. Again, MATrICP is more robust than MAICP, while MVICP clearly underperforms. The proposed scheme, however, achieves to register the point clouds accurately. Despite the large number of outliers, we are also able to get an outlier-free reconstruction of the scene based on the above thresholding principle. Of particular note, finally, is that JRMPC obtains these results with only components, a fact that further validates its potential.
Finally, we evaluate the performance JRMPC-I with a large number of point clouds collected with a moving sensor, namely the TUM dataset [48]. In particular, we converted the depth sequences fr1desk and fr2desk from this dataset into two sequences of and point sets, respectively. The first sequence includes several sweeps (local loops) over four desks in a typical office environment while the second sequence includes a full loop around a desk. The sequences are captured with a Kinect and ground-truth camera poses are provided with the help of a motion-capture system. As in the previous experiment, color data are not used by the registration algorithm.
To enable the algorithm to deal with a large number of sets, we considered a hierarchical scheme with two modules, a front-end and a back-end module. The front-end registers groups of successive point-sets and provides an outlier-free GMM, whose means are referred to as the mean set. The back-end module uses JRMPC-I on a temporal window of mean sets, that is, a new mean set is integrated at every times tamps and the local model instance of the window is refined with the batch method. We noticed that applying JRMPC-B on a small number of temporally non-overlapping integrated sets, e.g. one downsampled registered set per point sets, further improves the joint alignment.
All the initial point sets have been downsampled by a factor of before running the algorithm. We used , while and for fr1desk and fr2desk, respectively owing to differences in motion patterns. The overlap between successive groups in the front-end module is one set. The number of components is and for the back-end and front-end modules, respectively. The batch method refines the window model for iterations owing to its relatively small window. An optimized implementation of this procedure may lead to real-time performance.
ORB-SLAM2 (RGB-D) [49] | Elastic fusion [51] | RGBD SLAM [50] | JRMPC-I | |
---|---|---|---|---|
fr1/desk | 0.016 | 0.020 | 0.026 | 0.047 |
fr2/desk | 0.009 | 0.071 | 0.057 | 0.034 |
Table IV shows the performance of the proposed algorithm based on the protocol of [48]. Typically, the RMSE of the translation is used for evaluating SLAM methods. The error of JRMPC is computed for all frames while the other algorithms use only keyframes. We also provide the error of state-of-the-art RGB-D SLAM methods [49, 50, 51] as a reference (a direct comparison is not fair), as reported in [49]. Although SLAM methods use both modalities (RGB and depth) and invoke several modules to achieve accurate camera localization, the performance of the proposed algorithm is quite close to theirs.
Fig. 9 shows the camera trajectories obtained with the proposed algorithm as well as from RGBD-SLAM [50]. SLAM methods generally provide smooth trajectories owing to their internal tracking module and pose graph optimization. Instead, our algorithm simply registers depth data in a model-to-frame manner and one may observe local perturbations. Fig. 10 shows the final alignment for the fr1/desk sequence obtained with JRMPC-I and the corresponding ground truth. Interestingly, the proposed scheme delivers promising reconstructions despite the fact that it uses only depth data without pose graph optimization.
fr1/desk | fr2/desk |
JRMPC-I | Ground-truth |
Viii Conclusions
We presented a probabilistic generative model and its associated algorithm to jointly register multiple point sets. The vast majority of state-of-the-art techniques select one of the sets as the model and attempt to align the other sets onto this model. Instead, the proposed method treats all the point sets on an equal footing: any point is considered as realization of a single GMM and the registration is cast into a clustering problem. We formally derived an expectation-maximization algorithm that estimates the GMM parameters as well as the rotations and translations between each individual set and the initially unknown GMM means. An incremental version of the algorithm that efficiently integrates new point sets into the registration pipeline was also derived. We thoroughly validated the proposed method on challenging data sets gathered with depth cameras, we compared it with several state-of-the-art methods, and we showed its potential for effectively fusing depth data. In the future we plan to investigate the use of more efficient representations of generative models, e.g., [52] and an incremental registration method allowing the number of clusters to grow.
References
- [1] P. J. Besl and N. D. McKay, “A method for registration of 3-D shapes,” IEEE-TPAMI, vol. 14, pp. 239–256, 1992.
- [2] A. W. Fitzgibbon, “Robust registration of 2D and 3D point sets,” IVC, vol. 21, no. 12, pp. 1145–1153, 2001.
- [3] Y. Tsin and T. Kanade, “A correlation-based approach to robust point set registration,” in ECCV, 2004.
- [4] A. Myronenko and X. Song, “Point-set registration: Coherent point drift,” IEEE-TPAMI, vol. 32, no. 12, pp. 2262–2275, 2010.
- [5] R. Horaud, F. Forbes, M. Yguel, G. Dewaele, and J. Zhang, “Rigid and articulated point registration with expectation conditional maximization,” IEEE-TPAMI, vol. 33, no. 3, pp. 587–602, 2011.
- [6] B. Jian and B. C. Vemuri, “Robust point set registration using gaussian mixture models,” IEEE-TPAMI, vol. 33, no. 8, pp. 1633–1645, 2011.
- [7] J. Williams and M. Bennamoun, “Simultaneous registration of multiple corresponding point sets,” CVIU, vol. 81, no. 1, pp. 117–142, 2001.
- [8] F. Wang, B. C. Vemuri, A. Rangarajan, and S. J. Eisenschenk, “Simultaneous nonrigid registration of multiple point sets and atlas construction,” IEEE-TPAMI, vol. 30, no. 11, pp. 2011–2022, 2008.
- [9] V. M. Govindu and A. Pooja, “On averaging multiview relations for 3d scan registration,” IEEE-TIP, vol. 23, no. 3, pp. 1289–1302, 2014.
- [10] G. Blais and M. D. Levine, “Registering multiview range data to create 3d computer objects,” IEEE-TPAMI, vol. 17, no. 8, pp. 820–824, 1995.
- [11] T. Masuda and N. Yokoya, “A robust method for registration and segmentation of multiple range images,” CVIU, vol. 61, no. 3, pp. 295–307, 1995.
- [12] S. Izadi, D. Kim, O. Hilliges, D. Molyneaux, R. Newcombe, P. Kohli, J. Shotton, S. Hodges, D. Freeman, A. Davison, and A. Fitzgibbon, “Kinectfusion: Real-time 3d reconstruction and interaction using a moving depth camera,” in ACM Symposium on UIST, 2011.
- [13] R. Bergevin, M. Soucy, H. Gagnon, and D. Laurendeau, “Towards a general multi-view registration technique,” IEEE-TPAMI, vol. 18, no. 5, pp. 540–547, 1996.
- [14] U. Castellani, A. Fusiello, and V. Murino, “Registration of multiple acoustic range views for underwater scene reconstruction,” CVIU, vol. 87, no. 1-3, pp. 78–89, 2002.
- [15] D. F. Huber and M. Hebert, “Fully automatic registration of multiple 3d data sets,” IVC, vol. 21, no. 7, pp. 637–650, 2003.
- [16] S.-W. Shih, Y.-T. Chuang, and T.-Y. Yu, “An efficient and accurate method for the relaxation of multiview registration error,” IEEE-TIP, vol. 17, no. 6, pp. 968–981, 2008.
- [17] J. D. Banfield and A. E. Raftery, “Model-based Gaussian and non-Gaussian clustering,” Biometrics, vol. 49, no. 3, pp. 803–821, 1993.
- [18] G. Evangelidis, D. Kounades-Bastian, R. Horaud, and E. Psarakis, “A generative model for the joint registration of multiple point sets,” in ECCV, 2014.
- [19] Y. Chen and G. Medioni, “Object modelling by registration of multiple range images,” IVC, vol. 10, no. 3, pp. 145–155, 1992.
- [20] D. Chetverikov, D. Svirko, D. Stepanov, and P. Krsek, “The trimmed iterative closest point algorithm,” in ICPR, 2002.
- [21] A. V. Segal, D. Haehnel, and S. Thrun, “Generalized-icp,” in Robotics: Science and Systems, 2009.
- [22] J. Yang, H. Li, and Y. Jia, “Go-icp: Solving 3d registration efficiently and globally optimally,” in ICCV, 2013.
- [23] W. Wells III, “Statistical approaches to feature-based object recognition,” IJCV, vol. 28, no. 1/2, pp. 63–98, 1997.
- [24] S. Granger and X. Pennec, “Multi-scale EM-ICP: A fast and robust approach for surface registration,” in ECCV, 2002.
- [25] H. Chui and A. Rangarajan, “A new point matching algorithm for non-rigid registration,” CVIU, vol. 89, no. 2-3, pp. 114–141, 2003.
- [26] J. Hermans, D. Smeets, D. Vandermeulen, and P. Suetens, “Robust point set registration using em-icp with information-theoretically optimal outlier handling,” in CVPR, 2011.
- [27] R. Benjemaa and F. Schmitt, “A solution for the registration of multiple 3d point sets using unit quaternions,” in 5th European Conference on Computer Vision, 1998, pp. 34–50.
- [28] T. Masuda, “Registration and integration of multiple range images by matching signed distance fields for object shape modeling,” CVIU, vol. 87, no. 1, pp. 51–65, 2002.
- [29] Y. Peng, A. Ganesh, J. Wright, W. Xu, and Y. Ma, “Rasl: Robust alignment by sparse and low-rank decomposition for linearly correlated images,” IEEE TPAMI, vol. 34, no. 11, pp. 2233–2246, 2012.
- [30] D. Thomas, Y. Matsushita, and A. Sugimoto, “Robust simultaneous 3D registration via rank minimization,” in IEEE 3DIMPVT, 2012, pp. 33–40.
- [31] G. C. Sharp, S. W. Lee, and D. K. Wehe, “Multiview registration of 3d scenes by minimizing error between coordinate frames,” IEEE T PAMI, vol. 26, no. 8, pp. 1037–1050, August 2004.
- [32] S. Krishnan, P. Y. Lee, and J. B. Moore, “Optimisation-on-a-manifold for global registration of multiple 3d point sets,” Int. J. Intelligent Systems Technologies and Applications, vol. 3, no. 3/4, pp. 319–340, 2007.
- [33] A. Torsello, E. Rodol , and A. A, “Multiview registration via graph diffusion of dual quaternions,” in CVPR, 2011.
- [34] V. M. Govindu, “Combining two-view constraints for motion estimation,” in CVPR, 2001.
- [35] R. Hartley, J. Trumpf, Y. Dai, and H. Li, “Rotation averaging,” IJCV, vol. 103, no. 3, pp. 267–305, 2013.
- [36] Z. Li, J. Zhu, K. Lan, C. Li, and C. Fang, “Improved techniques for multi-view registration with motion averaging,” in 3DV, 2014.
- [37] D. Chetverikov, D. Stepanov, and P. Krsek, “Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm,” IVC, vol. 23, no. 3, pp. 299–309, 2005.
- [38] G. Jacob, “Registration of multiple point sets using the em algorithm,” in ICCV, 1999.
- [39] Y. Cui, S. Schuon, S. Thrun, D. Stricker, and C. Theobalt, “Algorithms for 3d shape scanning with a depth camera,” TPAMI, vol. 35, no. 5, pp. 1039–1050, 2013.
- [40] X. Mateo, X. Orriols, and X. Binefa, “Bayesian perspective for the registration of multiple 3d views,” CVIU, vol. 118, pp. 84 – 96, 2014.
- [41] M. Danelljan, G. Meneghetti, F. Khan, and M. Felsberg, “A probabilistic framework for color-based point set registration,” in CVPR, 2016.
- [42] T. Shiratori, J. Berclaz, M. Harville, C. Shah, T. Li, Y. Matsushita, and S. Shiller, “Efficient large-scale point cloud registration using loop closures,” in 3DV, 2015.
- [43] C. Bishop, Pattern Recognition and Machine Learning. Springer, 2006.
- [44] X.-L. Meng and D. B. Rubin, “Maximum likelihood estimation via the ECM algorithm: a general framework,” Biometrika, vol. 80, pp. 267–278, 1993.
- [45] S. Umeyama, “Least-squares estimation of transformation parameters between two point patterns,” IEEE-TPAMI, vol. 13, no. 4, pp. 376–380, 1991.
- [46] M. Hansard, R. Horaud, M. Amat, and G. Evangelidis, “Automatic detection of calibration grids in time-of-flight images,” CVIU, vol. 121, pp. 108–118, 2014.
- [47] M. Hansard, G. Evangelidis, Q. Pelorson, and R. Horaud, “Cross-Calibration of Time-of-flight and Colour Cameras,” CVIU, vol. 134, pp. 105–115, 2015.
- [48] J. Sturm, N. Engelhard, F. Endres, W. Burgard, and D. Cremers, “A benchmark for the evaluation of rgb-d slam systems,” in IROS, 2012.
- [49] J. D. T. Raul Mur-Artal, “ORB-SLAM2: An open-source SLAM system for monocular, stereo and RGB-D cameras,” arXiv, 2016.
- [50] F. Endres, J. Hess, J. Sturm, D. Cremers, and W. Burgard, “3-d mapping with an rgb-d camera,,” IEEE TR, vol. 30, no. 1, pp. 177–187, 2014.
- [51] T. Whelan, R. F. Salas-Moreno, B. Glocker, A. J. Davison, and S. Leutenegger, “Elasticfusion: Real-time dense slam and light source estimation,” IJRR, vol. 35, no. 14, pp. 1697–1716, 2016.
- [52] B. Eckart, K. Kim, A. Troccoli, and A. Kelly, “Accelerated generative models,” in CVPR, 2016.