Globally optimal point set registration by joint symmetry plane fitting

The present work proposes a solution to the challenging problem of registering two partial point sets of the same object with very limited overlap. We leverage the fact that most objects found in man-made environments contain a plane of symmetry. By reflecting the points of each set with respect to the plane of symmetry, we can largely increase the overlap between the sets and therefore boost the registration process. However, prior knowledge about the plane of symmetry is generally unavailable or at least very hard to find, especially with limited partial views, and finding this plane could strongly benefit from a prior alignment of the partial point sets. We solve this chicken-and-egg problem by jointly optimizing the relative pose and symmetry plane parameters, and notably do so under global optimality by employing the branch-and-bound (BnB) paradigm. Our results demonstrate a great improvement over the current state-of-the-art in globally optimal point set registration for common objects. We furthermore show an interesting application of our method to dense 3D reconstruction of scenes with repetitive objects.


Introduction
The alignment of two point sets is a fundamental geometric problem that occurs in many computer vision and robotics applications. In computer vision, the technique is used to stitch together partial 3D reconstructions in order to form a more complete model of an object or environment [23]. In robotics, point set registration is an essential ingredient to simultaneous localization and mapping with affordable consumer depth cameras [20] or powerful 3D laser range scanners [5]. The general approach for aligning two point sets does not require initial correspondences. It is given by the Iterative Closest Point (ICP) algorithm [28], a local search strategy that alternates between geometric correspondence establishment (i.e. by simple nearest neighbour search) and Procrustes alignment. The iterative procedure depends on a sufficiently accurate initial guess about the relative transformation (e.g. an identity transformation in the case of incremental ego-motion estimation).
The present work is motivated by a common problem that occurs when performing a dense reconstruction of an environment which contains multiple instances of the same object. An example of the latter is given by a room in which the same type of chair occurs more than once. Let us assume that the front-end of our reconstruction framework encompasses semantic recognition capabilities which are used to segment out partial point sets of objects of the same class and type [14]. There is a general interest in aligning those partial object point sets towards exploiting their mutual information and completing or even improving the reconstruction of each instance. The difficulty of this partial point set registration problem arises from two factors: • The relative pose between the different objects is arbitrary and unknown.
• Since the objects are observed under an arbitrary pose and with potential occlusions, the measured partial point sets potentially have very little overlap.
Our contribution focuses on the registration of only two partial point sets, which appears as a worthwhile starting point given that the registration of more than two point sets may be broken down into many pair-wise alignments. The plain ICP algorithm is only a local search strategy that depends on a sufficiently accurate initial guess, which is why it may not serve as a valid solution to our problem. A potential remedy is given by the globally optimal ICP algorithm presented by Yang et al. [26]. However, the algorithm still depends on sufficient overlap in the partial point sets, which is not necessarily a given (50% is reported as a requirement for high success rate).
The core idea of the present work consists of exploiting the fact that the majority of commonly observed objects contain a plane of symmetry. By reflecting the points of each partial point set with respect to the plane of symmetry, we may effectively increase the overlap between the two sets and vastly improve the success rate of the registration process. However, given that each point set only observes part of the object, the identification of the plane of symmetry in each individual point set appears to be an equally diffi-cult and ill-posed problem than the partial point set registration problem itself. It is only after the aligning transformation is found that symmetry plane detection would become a less challenging problem. In conclusion, the solution of each problem strongly depends on a prior solution to the other. We therefore present the following contributions: • We solve this chicken-and-egg problem by a joint solution of the aligning transformation as well as the symmetry plane parameters. To the best of our knowledge, our work is the first to address those two problems jointly, thus leading to a vast improvement over the existing state-of-the-art in globally optimal point set registration, especially in situations in which two point sets contain very limited overlap.
• We immediately provide a globally optimal solution to this problem by employing the branch-and-bound optimisation paradigm. Our work implicitly provides the first solution to globally optimal symmetry plane estimation in a single point set, or-more generallysymmetry plane detection across two point sets.
• We show a meaningful application of our algorithm in a dense 3D reconstruction scenario in which multiple instances of the same object occur.

Related Work
Despite strong mutual dependency, our method is the first to perform joint point set alignment and symmetry plane estimation. Our literature review therefore cites prior art on those two topics individually.

Point-set Registration
The Iterative Closest Point (ICP) [6,3,28] algorithm is a popular method for aligning two point sets. It does not depend on a prior derivation of point-to-point correspondences, and simply aligns the two sets by iteratively alternating between the two steps of finding nearest neighbours (e.g. by minimising point-to-point distances), and computing the alignment (e.g. using Arun's method [2]). To improve robustness of the algorithm against occlusions and reduced overlap, the method has been extended by outlier rejection [28,12] or data trimming [7] techniques. However, the classical ICP algorithm remains a local search algorithm for which the convergence depends on initial guess and sufficient overlap between both point sets.
An entire family of alternative approaches relies on the idea of expressing both point sets by a Gaussian Mixture Model (GMM) and aligning the latter using Gaussian Mixture Alignment (GMA). Notable GMA-based techniques for rigid and non-rigid registration are given by the robust point matching algorithm by Chui and Rangarajan [8], the coherent point drift strategy by Myronenko and Song [19], kernel correlation by Tsin and Kanade [24], and the GMM-Reg algorithm by Jian and Vemuri [15]. While GMA is advertised by improved robustness against poor initialisations, noise, and outliers, another key advantage with respect to point-based methods is given by a closed-form expression to evaluate the quality of the alignment (i.e. the method does not depend on an alternating search for nearest neighbours). However, the listed GMA algorithms remain local search algorithms, which makes them inapplicable in scenarios in which no prior about the relative transformation is given upfront.
In contrast, globally optimal algorithms avoid local minima by searching the entire space of relative transformations, often using the branch-and-bound paradigm [17,21,22]. Yang et al. [27,26] propose the Go-ICP algorithm, which applies branch-and-bound to the ICP problem to find the globally optimal minimum of the sum of L2-distances between nearest neighbours from two aligned point sets. The method is accelerated by using local ICP in the loop. However, missing robustness of the cost-function causes the method to remain sensitive with respect to occlusions and partial overlaps. Campbell et al. [4] finally devise GOGMA, a branch-and-bound variant in which the objective of minimising point-to-point errors is again replaced by the convolution of GMMs.
Inspired by those recent advances, we also employ the nested branch-and-bound strategy integrated with local ICP to find a globally optimal alignment. However, in contrast to all prior art, we are the first to jointly fit a plane of symmetry, which leads to a large improvement in partial scan alignment for common symmetric objects.

Symmetry Plane Estimation
For a comprehensive review of symmetry detection, we kindly refer the reader to Liu et al.'s review [18]. Here we only focus on the problem of symmetry plane fitting with missing data. The most straightforward solution is given by employing the RANdom SAmple Consensus (RANSAC) algorithm proposed by Fischler and Bolles [11], a wellknown algorithm for robust model fitting for outlier affected data. In the context of shape matching, the basic idea is to extract sparse characteristic points and match them between both sets. We then choose a random subset of correspondences and derive a hypothesis for the global transformation induced by these samples. The alignment quality is finally evaluated by the matching error between the two shapes. The method can be easily applied for detecting a plane of symmetry in a single point set by hypothesising the latter to be orthogonal to the axis connecting a correspondence. For example, Cohen et al. [10] detect symmetries in sparse point clouds by using appearance-based 3D-3D point correspondences in a RANSAC scheme. The detected symmetries are subsequently explored to eliminate noise from the point-clouds. Xu et al. [25] present a voting algorithm to detect the intrinsic reflectional symmetry axis. Using the axis as a hint, a completion algorithm for missing geometry is again shown. Jiang et al. [16] on the other hand propose an algorithm to find intrinsic symmetries in point clouds by using a curve skeleton. A set of filters then produces a candidate set of symmetric correspondences which are finally verified via spectral analysis. Although these works show results on partial data, the amount of missing data is typically small. Inspired by the work of [9] which detects symmetry by registration, we propose to detect the symmetry plane alongside partial point set registration, thus leading to improved performance in situations with limited overlap.

Preliminaries
We start by introducing the notation used throughout the paper and review the basic formulation of the ICP problem as well as planar reflections.

Notations and Assumptions
Let us denote the two partial object point sets by (sometimes called the model and data point sets, respectively). The goal pursued in this paper is the identification of a Euclidean transformation given by the rotation R and the translation t that transforms the points of Y such that they align with the points of X . If X and Y contain points in the 2D plane, R and t form an element of the group SE(2). If X and Y contain 3D points, R and t will be an element of SE(3). Note that alignment denotes a more general idea rather than just the minimisation of the sum of distances between each point of Y and its closest point within X . The point sets have different cardinality and potentially observe very different parts of the object with only very little overlap. This motivates our approach that takes object symmetries into account.

Registration of Two Point Sets
The standard solution to the point set registration problem is given by the ICP algorithm [28], which minimizes the alignment error given by where e r i (R, t|y i ) is the per-point residual error for y i , and x j is the closest point to y i within X , i.e.
Given an initial transformation R and t, the ICP algorithm iteratively solves the above minimization problem by alternating between updating the aligning transformation with fixed x j (i.e. using (1)), and updating the closest-point matches x j themselves using (2). It is intuitively clear that the ICP algorithm only convergence to a local minimum.

Modelling and Identifying Symmetry
Symmetry is modelled by a reflection by the symmetry plane. Let us define the symmetry plane by the normal n and depth-of-plane d such that any point x on the plane satisfies the relation x T n + d = 0. Let x now be a (2D or 3D) point from X . The reflection plane reflects x to a single reflected pointx given bŷ The term in parentheses is the signed distance between x and the reflection plane. The subtraction of 2n times this distance reflects the point to the other side of the plane. The problem of symmetry identification may now be formulated as a minimisation of the symmetry distance defined by where e s i (n, d|x i ) is the per-point residual error for x i , and x j is the nearest point tox i in X , i.e.
It is intuitively clear that the symmetry plane fitting problem may also be solved via ICP, the only difference being the parameters over which the problem is solved (i.e. n and d).

Transformed Symmetry Plane Parameters
Let us still assume that n and d are the symmetry plane parameters of a point set X , and R and t are the parameters that align a point set Y with X . Each transformed point Ry + t and its transformed, symmetric equivalent Rŷ + t must still fulfill the original reflection equation (3): Cancelling t and multiplying by R T on either side, we easily obtainŷ By comparing to (3), it is obvious thatn = R T n andd = t T n + d must represent the symmetry plane parameters for the original, untransformed set Y.

Alignment and Symmetry as a Joint Optimization Problem
We now introduce our novel optimization objective which jointly optimizes an aligning point set transformation as well as the plane of symmetry. The objective is then solved in a branch-and-bound optimization paradigm, for which we introduce both the domain parameterization as well as the derivation of upper and lower bounds.

Objective Function
We still assume that our two partial object point sets are , and that the symmetry plane is represented by the normal n and depth d.
, N } to be the aligned sets in either direction. The symmetry fitting objective of X employs where the difference to (4) is given by the fact that x j is now the nearest neighbour from the set X Y r . Similarly, using equation (7), the symmetry objective function for Y employs where y j is now chosen as the nearest neighbour from the set X r Y. The final registration error itself employs where x j is chosen as the nearest neighbour from the set X X s and weight w i is used to take the set X as aligned points first. The overall objective function becomes Direct optimization over R, t, n, and d using traditional ICP would easily get trapped in the nearest local minimum. We therefore propose to minimize the energy objective using the globally optimal branch-and-bound paradigm, an exhaustive search strategy that branches over the entire parameter space. In order to speed up the execution, the method derives upper and lower bounds for the minimal energy on each branch (i.e. sub-volume of the optimization space), and discards branches for which the lower bound remains higher than the upper bound in another branch. In the remainder of this section, we will discuss the two important questions of (i) how to parametrize and branch the parameter space, and (ii) how to find concrete values for the upper and lower bounds.

Domain Parameterization
2D problem: For a 2D point set, the domain of all feasible alignment parameters is represented by an angle r to parameterize the planar rotation R r = cos(r) − sin(r) sin(r) cos(r) , and a 2D translation vector t. The space of all 2D rotations can therefore be compactly represented by the interval [−π, +π]. For the translation, we assume that the optimal translation must lie within a cube defined by the interval [− , + ] 2 . The symmetry plane in the 2D case becomes a line which is parameterized by an angle α defining the normal vector n = [cos(α) sin(α)] T , and a scale d. Given the dual representation of line normal vectors, α is constrained to the interval [− π 2 , + π 2 ], and d lies in the interval [−ε, +ε]. 3D problem: The disadvantage of BnB is that its complexity grows exponentially in the dimensionality of the problem. We therefore take prior information about the 3D point sets into account that helps to decrease the dimensionality. We make the assumption that most objects are standing upright to the ground plane. The rotation between the different partial point sets is therefore still constrained to be a 1D rotation about the vertical axis, and the normal vector n remains a 1D variable the horizontal plane. The translation however becomes a 3D vector over the interval In conclusion, the 2D alignment problem has 5 degrees of freedom, whereas the 3D problem turns into a 6dimensional estimation problem.

Derivation of the Upper and Lower Bounds
The basic idea of BnB is to partition feasible set into convex sets which means it crucially depends on upper and lower bounds for the objective L2-energy on a given interval. Given that our objective (12) consists of a sum of many squared and positive sub-energies, lower and upper bounds on the overall objective energy can be derived by calculating individual upper and lower bounds on the alignment and symmetry residuals. Using the above introduced parametrizations, the upper bound E and the lower bound E of the optimal, joint L 2 registration and symmetry cost E * on a given interval of variables are therefore given as Upper bounds on an interval are easily given by the energy of an arbitrary point within the interval. Given an interval centered at {r 0 , t 0 , α 0 , d 0 }, upper bounds are therefore easily defined as e s i (α, d|xi) = e s i (α0, d0|xi) e r i (r, t|yi) = e r i (r0, t0|yi) e s i (r, t, α, d|yi) = e s i (r0, t0, α0, d0|yi).
The remainder of this section discusses the derivation of the lower bounds.
We can similarly derive a translation uncertainty radius γ t . For a translation volume with half side-length σ t centered at t 0 , we have Note that in the 2D case, we have γ t . = √ 2σ t . The lower bound of the registration term in equation (10) For more details, we kindly refer the reader to [13,26].
Lower bound of symmetry term e s i (α, d|x i ): Assuming that the normal is defined by an α-interval of half-length σ α and with center α 0 , we have For Now let x j ∈ X ∪Y r be the closest point to (x i −2n(x T i n+ d)) , and let x j0 ∈ X ∪ Y r be the closest point to (x i − 2n 0 (x T i n 0 + d 0 )). The lower bound is derived as follows: We furthermore have (22) and Substituting (22) and (23) in (21), we finally obtain Lower Bound of symmetry term e s i (r, t, α, d|y i ): By usingn = R T n andd = t T n + d, we analogously derive By substitutingn = R T n, similar to 22, the first term gives By also substitutingd = t T n + d, the second term gives where we have used t T 0 n 0 − t T n = t T 0 n 0 − t T 0 n + t T 0 n − t T n ≤ t 0 n 0 − n + n t 0 − t . Substituting (26) and (27) in (25), we finally obtain e s i (r, t, α, d|y i )

Implementation
Similar to prior art [26], we improve the algorithm's ability to handle the dimensionality of the problem by installing a nested BnB paradigm.

Nested BnB
We install a nested BnB scheme in which the outer layer searches through the space C rα of all angular parameters, while the inner layer optimizes over the space C td of translation and depth. While finding the bounds in a sub-volume Algorithm 1 Main Algorithm: BnB search for optimal registration and symmetry parameters Input: Data and model point set; threshold τ ; initial intervals C rα and C td . Output: Globally minimal error E * and the optimal r * , t * , α * , d * Put C rα into priority queue Q rα . Set E * = +∞. loop Read out interval with lowest E rα from Q rα . Quit the loop if E * − E rα < τ . Divide interval into 4 sub-intervals. for each sub-interval C rα do Compute the corresponding optimal t and d by calling Algorithm 2 with R r0 and n 0 (zero rotation and normal uncertainty). Compute E rα and E rα for C rα with the optimal t, d. if E rα < E * then Run ICP with initialization of (r 0 , t 0 , α 0 , d 0 ) Update E * , r * , α * , t * , and d * with the results of ICP. end if Discard C rα if E rα ≥ E * ; otherwise put it into Q rα end for end loop of the angle space, the algorithm calls the inner BnB algorithm to identify the optimal translation and scale. One important approximation that accelerates the execution is that-whenever the bounds in a sub-volume are derivedthe uncertainty of the non-optimized parameters of that particular layer are set to zero. Detailed descriptions are given in Algorithm 1 (the Main Algorithm) and Algorithm 2 (the Inner BnB).

Integration with Local ICP
Within the outer layer, whenever BnB identifies an interval C rα with an improved upper bound, we will execute a conventional local ICP algorithm starting from the center of the C rα and taking t * and d * as an initial value. Once ICP converges to a local minimum with a lower function value, the new value is used to further reduce the upper bound. The technique is inspired by Yang et al. [27].

Trimming for Outlier Handling
A general problem with partially overlapping point sets is that-even at the global optimum-some points may simply not have a correspondence, and should hence be treated as outliers. Although the addition of symmetry and point reflections already greatly alleviates this problem, we still add the strategy proposed in Trimmed ICP [7] for robust pointset registration. More specifically, in each iteration, only Algorithm 2 BnB search for optimal translation and depth given rotation and normal Input: Data and model point set; threshold τ ; initial intervals C td ; Currently lowest error E * Output: Minimal error E * and the optimal t * , d * Put C td into priority queue Q td . loop Read out interval with lowest E td from Q td . Quit the loop if E * − E td < τ . Divide interval into 8(2D)/16(3D) sub-intervals. for each sub-interval C td do Compute E td and E td for C td with the r 0 , n 0 . if E td < E * then Update E * and t * , d * . end if Discard C td if E td ≥ E * ; otherwise put it into Q td end for end loop a subset of the matched data points with smallest point-topoint distances are used for motion computation. In this work, we choose a 70%-subset for both symmetry and registration residuals.

Experiments
We now report our experimental results on both synthetic and real data. In all our experiments, we pre-normalized the pointsets such that all points are within the domain of [−1, 1]. The parameter ε is set to 0.5 (cf. Section 4.2). We run experiments on both 2D and 3D data, and compare our results against the open-source implementation of Go-ICP complemented by the ransac-based symmetry detection method presented in [10] and applied to a fusion of the aligned point sets. For Go-ICP, the stopping criterion τ is set to 0.001 · 0.7 · N , and for our method it is set to 0.001 · 0.7 · (M + 2N ) where 0.7 is the trimming ratio, M, N are the number of points in X and Y, respectively. For [10], the number of iterations is limited to 5000.

Performance on 2D Synthetic Data
Each experiment is generated by taking an image that contains a symmetrical object, and using the Sobel edge detector to extract the object contour points. To evaluate the performance, we randomly divide the contour points into two subsets with a defined and structured overlap. To conclude, Y is transformed by a random rotation and translation drawn from the intervals ±180 degrees and ±0.5. the centre is the result of the 2D version of Go-ICP, and the right one shows the ground truth alignment. Figure 2 shows all mean and median errors of all optimisation parameters over all evaluation results. As can be observed, our method has a substantially better ability to deal with partial overlaps compared to Go-ICP, thus leading to vastly improved symmetry plane parameters as well.
Outlier handling: To test resilience against outliers, we repeat the same experiment but add up to 30% outliers to both X and Y. Figure 3 indicates an example result, and Figure 4 again illustrates the mean and median errors over all experiments. While the registration error starts to increase earlier and the average angular errors tend to be higher, it can still be concluded that our method significantly outperforms Go-ICP followed by symmetry plane detection using [10].

Performance on 3D Synthetic Data
We choose 24 symmetrical CAD models from ShapeNet [1], 3 different types from 8 classes. There 8 objects contain more than one symmetry plane. For each object, we generate 16 depth images (with occlusions) from random views around the object, and producing 109 pairs of point sets for each object instance with varying overlap, and finally a total of about 2600 point-set registration    experiments. Our result is shown in Figure 5, which again illustrates mean and median errors of all estimated quantities as a function of the overlap between the sets. In 3D, it is natural that the camera captures only a small part of the object, and that in turn even the fusion of both aligned point scans may not enable a stable estimation of the symmetry plane, hence the slightly increased mean error. However, especially the median error is still lower compared to Go-ICP, thus confirming the mutual benefits of our joint estimation paradigm. Figure 6 shows a few concrete examples for which the 3D objects only contain a single symmetry plane and for which our method is largely outperformed. Figure 7 shows failure examples where the partial point sets lead to an ambiguity in the symmetry planes. In particular, in the first example, Go-ICP still works as the overlap ratio is sufficient for the registration.

Experiments on Real Data
Our last experiment is an exciting application to real data that goes back to the initial motivation in the introduction. Figure 8 shows two different depth images captured by a Kinect camera, each one containing three instances of the same object under different orientations. By pairwise alignment of partial object scans, the mutual information is transferred thus leading to more complete perception of each individual object. Note that we use simple ground plane fit- mean translation error of our method median translation error of our method mean translation error of 3D Go-ICP median translation error of 3D Go-ICP mean depth error of our method median depth error of our method mean depth error of ransac-based method median depth error of ransac-based method Figure 5. Mean and median errors of 3D registration compared against 3D Go-ICP followed by ransac-based symmetry detection [10]. Note that no outliers are added, but-as for the case of the results in 2D-a similar trend has been confirmed for up to 30% outliers.
ting and depth discontinuity-aware point clustering within object bounding boxes to isolate the partial object scans. With known position of the ground plane, we then transform the whole scene to be orthogonal to the ground plane and meet the assumption that all objects are placed upright and-in terms of relative rotation-differ only by an angle about the vertical axis. In Figure 8, the first column shows the original scan in different orientations, the second one the partial object measurements, the third one the completion obtained by using Go-ICP as an alignment algorithm and ransac-based symmetry detection, and the last one the result obtained by using our algorithm. As can be observed, our joint alignment strategy outperforms the comparison method, and achieves meaningful shape completion.

Discussion
Symmetry detection and point set alignment over sets with small overlap are challenging problems if handled separately. Our work demonstrates a substantial improvement in both accuracy and success rate of the alignment by solving those two problems jointly. The information gained from estimating symmetry and reflecting points notably makes up for otherwise missing correspondences. However, our current implementation is not competitive in terms of running time, hence we are working on a parallel implementation. We furthermore plan to extend the algorithm to multi-point set registration, and improve its ability to deal with the situation of multiple symmetry planes.