1 Introduction
Point cloud alignment is a fundamental problem for many applications in robotics [34, 22] and computer vision [43, 38, 51]. Finding the global transformation is generally hard: pointtopoint correspondences typically do not exist, the point clouds might only have partial overlap, and the underlying objects themselves are often nonconvex, leading to a potentially large number of alignment local minima. As such, popular local optimization techniques suffice only in circumstances with small true relative transformations and large overlap, such as in dense 3D incremental mapping [22, 38, 51]. Solving the alignment problem for large unknown relative transformations and small point cloud overlap calls for a global approach. Example applications are the loopclosure problem in SLAM [7] and the modelbased detection of objects in 3D scenes [28].
Motivated by the observation that surface normal distributions are translation invariant
[24] and straightforward to compute [37, 44], we develop a twostage branch and bound (BB) [30, 31] optimization algorithm for point cloud alignment. We model the surface normal distribution of each point cloud as a Dirichlet process (DP) [17, 47] vonMisesFisher (vMF) [19] mixture [45] (DPvMFMM). To find the optimal rotation, we minimize the distance between the distributions over the space of 3D rotations. We develop a novel refinable tessellation consisting of 4D tetrahedra (see Fig. LABEL:fig:600cell) which more uniformly approximates rotation space and is more efficient than the common axisangle tessellation [32, 21] during BB optimization. Given the optimal rotation and modeling the two point distributions as DP Gaussian mixtures [2, 10](DPGMM), we obtain the optimal translation similarly via BB over the space of 3D translations. The use of mixture models circumvents discretization artifacts, while still permitting efficient optimization. In addition to algorithmic developments, we provide corresponding theoretical bounds on the convergence of both BB stages, linking the quality of the derived rotation and translation estimates to the depth of the search tree and thus the computation time of the algorithm. Experiments on real data corroborate the theory, and demonstrate the accuracy and efficiency of BB as well as its robustness to realworld conditions, such as partial overlap, high noise, and large relative transformations.
2 Related Work
Local Methods There exists a variety of approaches for local point cloud alignment [9, 43]. Iterative closest point (ICP) [5], the most common of these, alternates between associating the points in both clouds and updating the relative transformation estimate under those associations. There are many variants of ICP [41] differing in their choice of cost function, how correspondences are established, and how the objective is optimized at each iteration. An alternative developed by Magnusson et al. [34] relies on the normal distribution transform (NDT) [6], which represents the density of the scans as a structured GMM. This approach has been shown to be more robust than ICP in certain cases [35]
. Approaches that use correlation of kernel density estimates (KDE) for alignment
[48] or GMMs [27] use a similar representation as the proposed approach. KDEbased methods scale poorly with the number of points. In contrast, we use mixture models inferred by nonparametric clustering algorithms (DPmeans [29] and DPvMFmeans [45]). This allows adaptive compression of the data, enabling the processing of large noisy point clouds (see Sec. 6 for experiments with more than k points). Straub et al. propose two local rotational alignment algorithms [45, 44] that, similarly to the proposed approach, utilize surface normal distributions modeled as vMF mixtures. Common to all local methods is the assumption of an initialization close to the true transformation and significant overlap between the two point clouds. If either of these assumptions are violated, local methods become unreliable as they tend to get stuck in suboptimal local minima [41, 43, 35].Global Methods Global point cloud alignment algorithms make no prior assumptions about the relative transformation or amount of overlap. For those reasons global algorithms, such as the proposed one, are often used to initialize local methods. 3Dsurfacefeaturebased algorithms [42, 20, 28, 1] involve extracting local features, obtaining matches between features in the two point clouds, and finally estimating the relative pose using RANSAC [18] or other robust estimators [25]. Though popular, featurebased algorithms are vulnerable to large fractions of incorrect feature matches, as well as repetitive scene elements and textures. A second class of approaches, including the proposed approach, rely on statistical properties of the two point clouds. Makadia et al. [36] separate rotational and translational alignment. Rotation is obtained by maximizing the convolution of the peaks of the extended Gaussian images (EGI) [24]
of the two surface normal sets. This search is performed using the spherical Fourier Transform
[16]. After rotational alignment, the translation is found similarly via the fast Fourier Transform. The use of histogrambased density estimates for the surface normal and point distributions introduces discretization artifacts. Additionally, the sole use of the peaks of the EGI makes the method vulnerable to noise in the data. For the alignment of 2D scans, Weiss et al. [50] and Bosse et al. [7] follow a similar convolutionbased approach. Early work by Li, Hartley and Kahl [32, 21] on BB for point cloud alignment used the axisangle (AA) representation of rotations. A drawback of this approach is that a uniform AA tessellation does not lead to a uniform tessellation in rotation space (see Sec. 4.1). As we show in Sec. 6, this leads to less efficient BB search. Parra et al. [39] propose improved bounds for rotational alignment by reasoning carefully about the geometry of the AA tessellation. GoICP [52] nests BB over translations inside BB over rotations and utilizes ICP internally to improve the BB bounds. GOGMA [8] uses a similar approach, but replaces the objective with a convolution of GMMs. Both GoICP and GOGMA involve BB over the joint 6dimensional rotation and translation space; since the complexity of BB is exponential in the dimension, these methods are relatively computationally expensive (see results Fig. LABEL:fig:apartment).3 The Point Cloud Alignment Problem
Our approach to point cloud alignment relies on the fact that surface normal distributions are invariant to translation [24] and easily computed [37, 44], allowing us to isolate the effects of rotation. Thus we decompose the task of finding the relative transformation into first finding the rotation using only the surface normal distribution, and then obtaining the translation given the optimal rotation.
Let a noisy sampling of a surface be described by the joint point and surface normal density , where and . A sensor observes two independent samples from this model: one from , and one from differing in an unknown rotation and translation . Given these samples, we model the marginal point densities , using the posterior of a Dirichlet process Gaussian mixture (DPGMM) [2], and model the marginal surface normal densities , using the posterior of a Dirichlet process von MisesFisher mixture (DPvMFMM) [4, 45]. Note that the formulation using DP mixture models admits arbitrarily accurate estimates of a large class of noisy surface densities (Theorem 2.2 in [14]). Given the density estimates, we formulate the problem of finding the relative transformation as
(1) 
where we represent rotations using unit quaternions in , the 4D sphere [23], and where denotes the rotation of a surface normal by a unit quaternion . Eq. (1) minimizes the metric via maximization of the convolution, which has been shown to be robust in practice [27]. This is a common approach for Gaussian MMs [48, 27, 8] but to our knowledge has not been explored for vMFMMs, nor for Bayesian nonparametric DP mixtures. In fact, the use of DP mixtures is critical, as it allows the automatic selection of a parsimonious, but accurate, representation of the point cloud data. This improves upon both kernel density estimates [48], which are highly flexible but make optimizing Eq. (1) intractable for large RGBD datasets, and fixedsized GMMs [27, 8]
, which require heuristic model selection and may not be rich enough to capture complex scene geometry. While exact posterior predictive DPMM densities cannot be computed tractably, excellent estimation algorithms are available, which we use in this work
[29, 45].Both optimization problems in Eq. (1) are nonconcave maximizations. Considering the geometry of the problem, we expect many local maxima, rendering typical gradientbased methods ineffective. This motivates the use of a global approach. We develop a twostep BB procedure [30, 31] that first searches over for the optimal rotation , and then over for the optimal translation . As BB may return multiple optimal rotations (e.g. if the scene has rotational symmetry) we estimate the optimal translation under each of those rotations, and return the joint transformation with the highest translational cost lower bound. Note that while is not necessarily the optimal transformation under rotation and translation jointly, the decoupling of rotation and translation we propose reduces the computational complexity of BB significantly. This is because the complexity scales exponentially in the search space dimension; optimizing over two 3D spaces ( and ) separately is significantly less costly than over the joint 6D space.
BB requires three major components: (1) a tessellation method for covering the optimization domain with subsets (see Sec. 4.1 and 5.1); (2) a branch/refinement procedure for subdividing any subset into smaller subsets (see Sec. 4.1 and 5.1); and (3) upper and lower bounds of the maximum objective on each subset to be used for pruning (see Sec. 4.2 and 5.2). BB proceeds by bounding the optimal objective in each subset, pruning those which cannot contain the maximum, subdividing the best subset to refine the bounds, and iterating. Note that in this work we select the node with the highest upper bound for subdivision. More nuanced strategies have been developed and could also be utilized [26, 31].
4 vMF Mixture Rotational Alignment
We model the distributions of surface normals as vonMisesFisher [19] mixture models (vMFMM) with means , concentrations , and positive weights , , for , with density
(2) 
While there are many techniques for inferring vMFMMs [3, 15, 45], we use a nonparametric method [45] that infers an appropriate automatically. The rotational alignment problem from Eq. (1) with this model becomes
(3) 
We obtain the following objective function by noting that the integral is the normalization constant of a vMF density with concentration :
(4) 
4.1 Cover and Refinement of the Rotation Space
In this section, we develop a novel tessellation scheme for the space of rotations, and show how to refine it in a way that guarantees convergence of BB for rotational alignment. We follow a similar approach to the geodesic grid tessellation of a sphere in 3D (i.e. ): as depicted in Fig. fig:tessellationS2, starting from an icosahedron, each of the triangular faces is subdivided into four triangles of equal size. Then the newly created triangle corners are normalized to unit length, projecting them onto the unit sphere.
In four dimensions we instead start with the analogue of the icosahedron, the 600cell [12] (shown in Fig. LABEL:fig:600cell), an object composed of 600 4D tetrahedra. We first generate its 120 vertices with the following algorithm [12, pp. 402–403]. Let . Then the (unnormalized) 120 vertices of the 600cell in 4D are even permutations of (96 vertices), all permutations of (8 vertices), and all permutations of (16 vertices). We then scale the 120 vertices to each have unit norm, representing a 3D quaternion rotation. Next, noting that the angle between any two connected tetrahedra vertices is , we iterate over all possible choices of 4 vertices, and only select those tetrahedra for which all pairwise angles are . This collection of tetrahedra, which are “flat” in 4D analogous to triangles in 3D, comprises a 4D object which approximates the 4D sphere, . Then, since the set of all quaternion rotations may be represented by any hemisphere of ( and
describe the same rotation), we define the “north” vector to be
, and only keep those tetrahedra for which at least one vertex has angle to the north vector. This results in 330 tetrahedra that approximate the 4D upper hemisphere in , i.e. the space of quaternion rotations. Note that this construction procedure is the same for any optimization on , so it can be performed once and the result may be stored for efficiency.One major advantage of the proposed tessellation is that it is exactly uniform at the 0th level and approximately uniform for deeper subdivision levels (Fig. fig:tessellationS2 shows the analogous nearuniformity for ). This generally tightens bounds employed by BB, leading to more efficient optimization. Another advantage is that this tessellation is a nearexact covering of the upper hemisphere of . Only 7% of rotation space is covered twice, meaning that BB wastes little time with duplicate searching. The widely employed AAtessellation scheme [32, 21, 39, 52], in contrast, uniformly tessellates a cube enclosing the axisangle space, a 3D sphere with radius , and maps that tessellation onto the rotation space. There are two major issues with the AA approach. First, it covers 46% of rotation space twice [32, 21] (see Fig. fig:tessellationAA). Second, it does not lead to uniform tessellation in rotation space. The reason for this is that the Euclidean metric in AA space is a poor approximation of the distance on the rotation manifold [32]. Fig. fig:tessellationAA shows the AA tessellation analog for , highlighting its significant nonuniformity. We empirically find that the tessellation leads to more efficient BB optimization than the AA tessellation (see results in Figs. LABEL:fig:bunnyFull and LABEL:fig:bunny_045).
We now discuss two properties of the proposed tessellation required by BB: 1) that it is a cover for the upper hemisphere of , guaranteeing that BB will search the whole space of rotations; and 2) that it is refinable, so BB can search promising subsets in increasingly more detail.
Cover Let the four vertices of a single tetrahedron from our approximation of be denoted , . Then, stacking them horizontally into a matrix , the projection of the tetrahedron onto is:
(5) 
In other words, is the set of unit quaternions found by extending the (flat in 4D) tetrahedron to the unit sphere using rays from the origin. For , this is displayed in the second row of Fig. fig:tessellationS2. The proposed set of 330 projected tetrahedra forms a cover of the upper hemisphere of .
Refinement Next, we require a method of subdividing any in the cover. Similar to the triangle subdivision method for refining the tessellation of , each 4D tetrahedron can be subdivided into eight smaller tetrahedra [33] as depicted in Fig. fig:tetrahedronSubdivision. The resulting six new vertices for the subdivided tetrahedra are scaled to unit length. As we have the freedom to choose one of three internal edges for subdivision, we choose the internal edge with the minimum angle between its unitnorm vertices. In other words, denoting for to be the three internal dot products,
(6) 
This process forms the eight new subdivided cover elements . For example, if , are the vertices of , then one of the subdivisions (corresponding to one of the “corner” subtetrahedra in Fig. fig:tetrahedronSubdivision) of would have vertices
(7) 
Selecting the internal edge via Eq. (6) is critical to our BB convergence guarantee in Sec. 4.4. If Eq. (6) is not used, the individual subsets
can become highly skewed due to repeated distortion from the unitnorm projection of the vertices, and refining
does not necessarily correspond to shrinking the angular range of rotations it captures. Since we use Eq. (6), however, Lemma 1 guarantees that subdividing shrinks its set of rotations appropriately:Lemma 1.
Let be the min dot product between vertices of any one at refinement level . Then
(8) 
This result (proof in the supplement) shows that the tetrahedra shrink and allow BB to improve its bounds during subdivision. Figure fig:minAngleBounds demonstrates the tightness of this bound, showing that converges to 0 as . We conjecture that the max dot product satisfies a similar recursion, , although this is not required for our convergence analysis. Fig. fig:minAngleBounds shows empirically that this matches the true max dot product, but we leave the proof as an open problem.
4.2 vMF Mixture Model Bounds
BB requires both upper and lower bounds on the maximum of the objective function within each projected tetrahedron , i.e. we need and such that
(9) 
For the lower bound , one can evaluate the objective at any point in (e.g. its center). For the upper bound , we use a quadratic upper bound on (see Fig. fig:fbound and the supplement for details), noting that for all , where
(10) 
whose computation is discussed in Sec. 4.3. This results in the upper bound where
(11) 
and is defined as the matrix for which for any quaternion (see the supplement for details). Writing as a linear combination of vertices of as in Eq. (5),
(12) 
Since , and we have the constraint , we can search over all possible combinations of components of being zero or nonzero. Thus we solve the optimization for given each possible subset of nonzero components of , and set
(13) 
For , we use a Lagrange multiplier for the equality constraint in Eq. (12
) and set the derivative to 0, yielding a small generalized eigenvalue problem of dimension
,(14) 
where is a dimensional vector, and subscript denotes the submatrix with rows and columns selected from . The condition that all elements of are nonnegative in Eq. (14) enforces that and thus corresponds to a solution that lies in . Note that if
is an eigenvector, so is
. If no satisfies , then we define .4.3 Computing and
To find the upper bound in Eq. (12), we require the constants and for each pair of mixture components . Given their definitions in Eq. (10), we have
(15) 
Since the inner optimization objective only depends on the rotation of by , we can reformulate the optimization as being over the set of 3D vectors such that for some . Thus, finding and is equivalent to finding the closest and furthest unit vectors in 3D to over the set of such vectors , shown in Fig. fig:closestpt. To solve this problem, let the vertices of be , , and define the matrix where . The inner optimization in Eq. (15) can be written as (for set ; for set )
(16) 
Showing that Eq. (16) is equivalent to solving the inner optimizations of Eq. (15) is quite technical and is deferred to the supplement. Again we search over all possible combinations of components of being zero or nonzero (we do not check the case since in this case the matrix below is rankdeficient). We thus solve the optimization for given each subset , of nonzero components, and set
(17) 
To solve for , we use a Lagrange multiplier for the equality constraint, and set derivatives to 0 to find that
(18)  
where  
(19) 
and is the matrix constructed from the set of columns in corresponding to . Note that is also defined to be if is not invertible. After solving for the value of via Eq. (17), we substitute it back into Eq. (15) to obtain or as desired.
4.4 Convergence Properties
We have now developed all the components necessary to optimize Eq. (4) via BB on . Theorem 1 (proof in the supplement) provides a bound on the worstcase search tree depth to guarantee BB terminates with rotational precision of degrees, along with the overall computational complexity. Note that the complexity of BB is exponential in , but since is logarithmic in (by Theorem 1, Eq. (20) and for ), the complexity of BB is polynomial in . Recall from Sec. 4.1 that for the 600cell is .
Theorem 1.
Suppose is the initial maximum angle between vertices in the tetrahedra tessellation of , and let
(20) 
Then at most refinements are required to achieve an angular tolerance of on , and BB has complexity .
5 Gaussian Mixture Translational Alignment
In this section, we reuse notation for simplicity and to highlight parallels between the translational and rotational alignment problems. We model the density of points in the two point clouds as Gaussian mixture models (GMMs) with means
, covariances , and weights , , for , with density(21) 
GMMs can be inferred in a variety of ways [29, 10]. Let be the optimal rotation corresponding to recovered using BB over . Then defining
(22) 
the translational optimization in Eq. (1) becomes:
(23) 
This is again a nonconcave maximization, motivating the use of a global approach. Thus, we develop a second BB procedure on to find the optimal translation.
5.1 Cover and Refinement of
We tessellate the space of translations, with rectangular cells. The initial tessellation is obtained by enclosing both point clouds with a single rectangular bounding box with diagonal length . For the refinement step, we choose to subdivide the cell into eight equalsized rectangular cells. Thus, the minimum diagonal of the rectangular cells at refinement level possesses a straightforward shrinkage property similar to Eq. (8),
(24) 
5.2 Gaussian Mixture Model Bounds
As in the rotational problem, the translational BB algorithm requires lower and upper bounds on the objective function in Eq. 23:
(25) 
For the lower bound , one can evaluate the objective at any (e.g. its center). For the upper bound , we use a linear upper bound on (see Fig. LABEL:fig:tfbound and the supplement for details), noting that for all , where
(26) 
whose computation is discussed in Section 5.3. This results in the upper bound , where
(27) 
This is a concave quadratic maximization over a rectangular cell . Thus, we obtain as the maximum over all local optima in the interior, faces, edges, and vertices of .
5.3 Computing and
Using the form of in Eq. (22), we have that
(28) 
Because of the concavity of the objective, can be obtained with the exact same algorithm as used to solve Eq. 27. can be obtained by checking the vertices of , as the minimum of a concave function over a rectangular cell must occur at one of its vertices.
5.4 Convergence Properties
We now have all the components necessary to optimize Eq. (23) via BB on . As in the rotational alignment case, we provide a characterization (Theorem 2, proof in the supplement) of the maximum refinement depth required for a desired translational precision , along with the complexity of the algorithm. Note that while the complexity of BB is exponential in , is logarithmic in (Theorem 2), so BB has polynomial complexity in .
Theorem 2.
Suppose is the initial diagonal length of the translation cell in , and let
(29) 
Then at most refinements are required to achieve a translational tolerance of , and BB has complexity .
6 Results and Evaluation
We evaluate BB (both with and without final local refinement [11]) on four datasets [13, 49, 40] compared to three global methods: an FTbased method [36], GoICP [52] ( trimming), and GOGMA [8]. To generate the vMFMMs and GMMs for BB, we cluster the data with DPvMFmeans [45] and DPmeans [29]
, and fit maximum likelihood MMs to the clustered data. To account for nonuniform point densities due to the sensing process, we weight each point’s contribution to the MMs by its surface area, estimated by the disc of radius equal to the fifth nearest neighbor distance. We use kNN+PCA
[54, 55] to extract surface normals. To improve the robustness of BB, it is run three times on each problem with scale values in DPvMFmeans (included in the timing results). The scale for DPmeans is manually selected to yield around mixture components. Using Theorems 1 and 2, we terminate rotational BB at and translational BB at for a rotational accuracy of and a translational accuracy of , where is defined in Eq. (24). All timing results include algorithmspecific preprocessing of the data. We used a 3GHz core i7 CPU and a GeForce GTX 780 GPU. While clustering via DPmeans and DPvMFmeans uses the GPU, we only use parallel CPU threads for the eight BB bound evaluations after each branch step.Stanford Bunny [49] Independent of the tessellation strategy, BB perfectly aligns the Stanford Bunny with a randomly transformed version of itself, as shown in Fig. LABEL:fig:bunnyFull. The results of aligning two partial scans of the Stanford Bunny with relative viewpoint difference are shown in Fig. LABEL:fig:bunny_045. BB’s initial alignment is close enough to allow ICP to converge to a perfect alignment. The proposed approach leads to a faster reduction in the bound gap, faster exploration, and a smaller number of active nodes, while reducing the computation time per iteration by an order of magnitude vs the AA tessellation. This shows conclusively that the proposed tessellation leads to more efficient BB optimization. Note that the AA tessellation starts at unexplored space because it covers the rotation space more than once as discussed in Sec. 4.1. In both cases BB finds the optimal translation within iterations.
Happy Buddha [13] This dataset consists of 15 scans taken at rotational increments about the vertical axis of a statue. This dataset is challenging, as the scans contain few overlapping points, and the surface normal distributions are anisotropic. We perform pairwise alignment of consecutive scans, and render the aligned scans together in one coordinate system (Fig. LABEL:fig:happyB). The only successful alignment is produced by BB+ICP. This shows the advantage of using surface normals for rotational alignment. Other methods using points (GoICP) or GMMs (GOGMA) have difficulty dealing with ambiguities due to the “flatness” of the scans.
Office Scan Figure LABEL:fig:rgbd demonstrates that BB+ICP finds accurate registrations on noisy, incomplete, cluttered and irregular point clouds as long as good surface normal estimates are available. This demonstrates the potential use of BB+ICP for loop closure detection.
Apartment Dataset [40] This dataset consists of LiDAR scans with an average overlap of . Figure LABEL:fig:apartment shows the BB+ICP aligned scans of the dataset. Table LABEL:tab:apartment compares the accuracy and inlier percentages defined by (C)oarse (m; ), (M)edium (m; ) and (F)ine (m; ) thresholds for all algorithms. For GoICP, we used scan points and an accuracy threshold of . We used the scale parameter of m for GMM computations in both GOGMA and BB.
Manmade environments such as this dataset exhibit “Manhattan World” (MW) symmetry in their surface normal distributions [44, 46]. We thus transform the rotation obtained via rotational BB by all MW rotations, and search over all using translational BB. Note that doing this is straightforward in the proposed decoupled BB approach, as opposed to a joint approach, e.g. GoICP and GOGMA.
Table LABEL:tab:apartment and Fig. LABEL:fig:CDF show that BB with searching over both scale and MW rotations leads to the best accuracy among all algorithms, with a 3x speedup over the best method, GOGMA (which uses a GPU). From the inlier percentages it is clear that FT and GoICP do not perform well. The CDFs in Fig. LABEL:fig:CDF show that accounting for MW symmetry (red, green) is important; ignoring it (blue) causes scans to be flipped by /, affecting the mean error strongly.
7 Conclusion
We introduced a BB approach to global point cloud alignment with convergence guarantees, based on a Bayesian nonparametric point cloud representation and a novel tessellation of rotation space. The method decouples translation and rotation via the use of surface normals, making it more efficient than previous joint approaches. Experiments demonstrate the robustness of the method to noisy real world data, partial overlap, and angular viewpoint differences. We expect that the proposed tessellation of will be useful in other rotational BB algorithms. All code is available at http://people.csail.mit.edu/jstraub/.
References
 [1] D. Aiger, N. J. Mitra, and D. CohenOr. 4points congruent sets for robust pairwise surface registration. In ACM TOG, volume 27, page 85, 2008.
 [2] C. Antoniak. Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. The Annals of Statistics, 1152–1174, 1974.
 [3] A. Banerjee, I. S. Dhillon, J. Ghosh, S. Sra, and G. Ridgeway. Clustering on the unit hypersphere using von MisesFisher distributions. JMLR, 6(9), 2005.
 [4] M. Bangert, P. Hennig, and U. Oelfke. Using an infinite von MisesFisher mixture model to cluster treatment beam directions in external radiation therapy. In ICMLA, 2010.
 [5] P. J. Besl and N. D. McKay. A method for registration of 3D shapes. TPAMI, 14(2):239–256, 1992.
 [6] P. Biber and W. Straßer. The normal distributions transform: A new approach to laser scan matching. In IROS, 2003.
 [7] M. Bosse and R. Zlot. Map matching and data association for largescale twodimensional laser scanbased SLAM. IJRR, 27(6):667–691, 2008.
 [8] D. Campbell and L. Petersson. Gogma: Globallyoptimal gaussian mixture alignment. In CVPR, June 2016.
 [9] R. J. Campbell and P. J. Flynn. A survey of freeform object representation and recognition techniques. Computer Vision and Image Understanding, 81(2):166–210, 2001.
 [10] J. Chang and J. W. Fisher III. Parallel sampling of DP mixture models using subclusters splits. In NIPS, 2013.
 [11] Y. Chen and G. Medioni. Object modeling by registration of multiple range images. In ICRA, 1991.
 [12] H. S. M. Coxeter. Regular polytopes. Courier Corporation, 1973.
 [13] B. Curless and M. Levoy. A volumetric method for building complex models from range images. In SIGGRAPH, 1996.
 [14] L. Devroye. A Course in Density Estimation. Birkhauser Boston Inc., 1987.
 [15] I. S. Dhillon and D. S. Modha. Concept decompositions for large sparse text data using clustering. Machine Learning, 42(12):143–175, 2001.
 [16] J. R. Driscoll and D. M. Healy. Computing Fourier transforms and convolutions on the 2sphere. Advances in Applied Mathematics, 15(2):202–250, 1994.
 [17] T. Ferguson. A Bayesian analysis of some nonparametric problems. The Annals of Statistics, 209–230, 1973.
 [18] M. Fischler and R. Bolles. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM, 24(6):381–395, 1981.
 [19] N. I. Fisher. Statistical Analysis of Circular Data. Cambridge University Press, 1995.
 [20] N. Gelfand, N. J. Mitra, L. J. Guibas, and H. Pottmann. Robust global registration. In Symposium on Geometry Processing, volume 2, page 5, 2005.
 [21] R. I. Hartley and F. Kahl. Global optimization through rotation space search. IJCV, 82(1):64–79, 2009.
 [22] P. Henry, M. Krainin, E. Herbst, X. Ren, and D. Fox. RGBD mapping: Using Kinectstyle depth cameras for dense 3D modeling of indoor environments. IJRR, 31(5):647–663, 2012.
 [23] B. K. Horn. Some notes on unit quaternions and rotation. 2001.
 [24] B. K. P. Horn. Extended Gaussian images. Proceedings of the IEEE, 72(12):1671–1686, 1984.
 [25] P. J. Huber. Robust statistics. Springer, 1981.
 [26] T. Ibaraki. Theoretical comparisons of search strategies in branchandbound algorithms. IJCIS, 5(4):315–344, 1976.
 [27] B. Jian and B. C. Vemuri. Robust point set registration using gaussian mixture models. PAMI, 33(8):1633–1645, 2011.
 [28] A. E. Johnson and M. Hebert. Surface matching for object recognition in complex threedimensional scenes. Image and Vision Computing, 16(9):635–651, 1998.

[29]
B. Kulis and M. I. Jordan.
Revisiting kmeans: New algorithms via Bayesian nonparametrics.
In ICML, 2012.  [30] A. H. Land and A. G. Doig. An automatic method of solving discrete programming problems. Econometrica: Journal of the Econometric Society, 497–520, 1960.
 [31] E. L. Lawler and D. E. Wood. Branchandbound methods: A survey. Operations research, 14(4):699–719, 1966.
 [32] H. Li and R. Hartley. The 3D3D registration problem revisited. In ICCV, 2007.
 [33] A. Liu and B. Joe. Quality local refinement of tetrahedral meshes based on 8subtetrahedron subdivision. AMS Math. Comp., 65(215):1183–1200, 1996.
 [34] M. Magnusson, A. Lilienthal, and T. Duckett. Scan registration for autonomous mining vehicles using 3DNDT. Journal of Field Robotics, 24(10):803–827, 2007.
 [35] M. Magnusson, A. Nüchter, C. Lörken, A. J. Lilienthal, and J. Hertzberg. Evaluation of 3D registration reliability and speeda comparison of ICP and NDT. In ICRA, 2009.
 [36] A. Makadia, A. Patterson, and K. Daniilidis. Fully automatic registration of 3D point clouds. In CVPR, 2006.
 [37] N. J. Mitra, A. Nguyen, and L. Guibas. Estimating surface normals in noisy point cloud data. IJCGA, 14:261–276, 2004.
 [38] R. A. Newcombe, A. J. Davison, S. Izadi, P. Kohli, O. Hilliges, J. Shotton, D. Molyneaux, S. Hodges, D. Kim, and A. Fitzgibbon. Kinectfusion: Realtime dense surface mapping and tracking. In ISMAR, 2011.
 [39] A. J. Parra Bustos, T.J. Chin, and D. Suter. Fast rotation search with stereographic projections for 3D registration. In CVPR, 2014.
 [40] F. Pomerleau, M. Liu, F. Colas, and R. Siegwart. Challenging data sets for point cloud registration algorithms. IJRR, 31(14):1705–1711, 2012.
 [41] S. Rusinkiewicz and M. Levoy. Efficient variants of the ICP algorithm. In 3D Digital Imaging and Modeling, 2001.
 [42] R. B. Rusu, N. Blodow, and M. Beetz. Fast point feature histograms (FPFH) for 3D registration. In ICRA, 2009.
 [43] J. Salvi, C. Matabosch, D. Fofi, and J. Forest. A review of recent range image registration methods with accuracy evaluation. Image and Vision Computing, 25(5):578–596, 2007.
 [44] J. Straub, N. Bhandari, J. J. Leonard, and J. W. Fisher III. Realtime Manhattan world rotation estimation in 3D. In IROS, 2015.

[45]
J. Straub, T. Campbell, J. P. How, and J. W. Fisher III.
Smallvariance nonparametric clustering on the hypersphere.
In CVPR, 2015.  [46] J. Straub, G. Rosman, O. Freifeld, J. J. Leonard, and J. W. Fisher III. A Mixture of Manhattan Frames: Beyond the Manhattan World. In CVPR, 2014.
 [47] Y. W. Teh. Dirichlet processes. In Encyclopedia of Machine Learning. Springer, New York, 2010.
 [48] Y. Tsin and T. Kanade. A correlationbased approach to robust point set registration. In ECCV, 2004.
 [49] G. Turk and M. Levoy. Zippered polygon meshes from range images. In SIGGRAPH, 1994.
 [50] G. Weiss, C. Wetzler, and E. Von Puttkamer. Keeping track of position and orientation of moving indoor systems by correlation of rangefinder scans. In IROS, 1994.
 [51] T. Whelan, M. Kaess, H. Johannsson, M. Fallon, J. Leonard, and J. McDonald. Realtime large scale dense RGBD SLAM with volumetric fusion. IJRR, 2014.
 [52] J. Yang, H. Li, and Y. Jia. GoICP: Solving 3D registration efficiently and globally optimally. In ICCV, 2013.
 [53] R. Webb. Stella software. http://www.software3d.com/Stella.php and https://en.wikipedia.org/wiki/600cell.
 [54] Meshlab. http://meshlab.sourceforge.net/. Accessed: 20161115.
 [55] Point cloud library. http://pointclouds.org/. Accessed: 20161115.