Least-squares registration of point sets over SE (d) using closed-form projections

TL;DR

This paper introduces a novel least-squares approach for registering multiple point sets in d-dimensional space by reducing the problem to a positive semidefinite cone optimization, solved efficiently with ADMM and closed-form projections.

Contribution

It presents a new nonconvex-to-convex reduction for point set registration and a solver with closed-form subproblems, improving efficiency and stability over existing methods.

Findings

The proposed method converges reliably with appropriate parameters.

It achieves faster registration with comparable or better accuracy.

Demonstrated effectiveness in 2D shape matching and 3D multiview registration.

Abstract

Consider the problem of registering multiple point sets in some $d$-dimensional space using rotations and translations. Assume that there are sets with common points, and moreover the pairwise correspondences are known for such sets. We consider a least-squares formulation of this problem, where the variables are the transforms associated with the point sets. The present novelty is that we reduce this nonconvex problem to an optimization over the positive semidefinite cone, where the objective is linear but the constraints are nevertheless nonconvex. We propose to solve this using variable splitting and the alternating directions method of multipliers (ADMM). Due to the linearity of the objective and the structure of constraints, the ADMM subproblems are given by projections with closed-form solutions. In particular, for $m$ point sets, the dominant cost per iteration is the partial eigendecomposition of an $md \times md$ matrix, and $m-1$ singular value decompositions of $d \times d$ matrices. We empirically show that for appropriate parameter settings, the proposed solver has a large convergence basin and is stable under perturbations. As applications, we use our method for $2$D shape matching and $3$D multiview registration. In either application, we model the shapes/scans as point sets and determine the pairwise correspondences using ICP. In particular, our algorithm compares favorably with existing methods for multiview reconstruction in terms of timing and accuracy.

Tables

Table 1. TABLE I: Breakup of the complexity of Algorithm 1 ( k 𝑘 k in step 4 4 4 is the number of Arnoldi iterations).

Operation	Complexity
$1.L^{†}$	$O\left(m^3\right)$
$2.B{L}^{†}{B}^{\top }$	$O\left(m^{2}d\right)$
$3.C=D-B{L}^{†}{B}^{\top }$	$O\left(m^2{d}^2\right)$
$4.G\leftarrow {\Pi }_{\Omega }\left(H-{\rho }^{-1}\left(C+\Lambda \right)\right)$	$O\left(k{m}^2{d}^2\right)$	\rdelim}20.1mm[recurring cost]
$5.H\leftarrow {\Pi }_{\Theta }\left(G+{\rho }^{-1}\Lambda \right)$	$\left(m-1\right)O\left(d^3\right)$
$6.$final rounding	$O\left(d^3\right)$

Table 2. TABLE II: Comparison of rotational error of proposed algorithms with that in Ahmed and Chaudhury ( 2017 ) for 10 10 10 point sets. Rows I, II and III correspond to the ground truth, results from Ahmed and Chaudhury ( 2017 ) , and results from our method (see text for details).

	Determinant of optimal transforms										Error
I	+1	+1	+1	+1	+1	+1	+1	+1	+1	+1	$0°$
II	+1	+1	+1			+1	+1	+1	+1	+1	$59.4°$
III	+1	+1	+1	+1	+1	+1	+1	+1	+1	+1	$5.23°$

Table 3. TABLE III: Registration results from known correspondences. Left to right: ground truth, before registration, and after registration. Six point sets were used in each case (which are color coded). We used ρ = 10 𝜌 10 \rho=10 for all the experiments.

Table 4. TABLE IV: Visual comparison of reconstructions based on their cross-sections.

Model	Scans	Proposed	MAICP

Equations

$$\mathbf{R}\in \mathbb{O}(d)\text{ and }\det (\mathbf{R})=1.$$

$$\{{\mathbf{x}}_{ij}^k:1\le k\le n_{ij}\}\text{and}\{{\mathbf{x}}_{ji}^k:1\le k\le n_{ij}\}.$$

$${\mathbf{R}}_i{\mathbf{x}}_{ij}^k+{\mathbf{t}}_i={\mathbf{R}}_j{\mathbf{x}}_{ji}^k+{\mathbf{t}}_j(i\sim j).$$

$$\min \sum _{i\sim j}\sum _{k=1}^{n_{ij}}\Vert {\mathbf{R}}_i{\mathbf{x}}_{ij}^k+{\mathbf{t}}_i-{\mathbf{R}}_j{\mathbf{x}}_{ji}^k-{\mathbf{t}}_j{\Vert }_2^2.$$

$$\mathbf{R}=\left[{\mathbf{R}}_1\cdots {\mathbf{R}}_m\right]\text{and}\mathbf{T}=\left[{\mathbf{t}}_1\cdots {\mathbf{t}}_m\right],$$

$$\min \sum _{i\sim j}\sum _{k=1}^{n_{ij}}\Vert \mathbf{R}{\mathbf{d}}_{ij}^k+\mathbf{T}{\mathbf{e}}_{ij}{\Vert }_2^2,$$

$${\mathbf{d}}_{ij}^k=({\mathbf{e}}_i\otimes \mathbf{I}){\mathbf{x}}_{ij}^k-({\mathbf{e}}_j\otimes \mathbf{I}){\mathbf{x}}_{ji}^k,$$

$$\mathbf{L}=\sum _{i\sim j}{n}_{ij}{\mathbf{e}}_{ij}{\mathbf{e}}_{ij}^{\top }\text{and}\mathbf{B}=\sum _{i\sim j}\sum _{k=1}^{n_{ij}}{\mathbf{d}}_{ij}^k{\mathbf{e}}_{ij}^{\top },$$

$$⟨\mathbf{C},{\mathbf{R}}^{\top }\mathbf{R}⟩=\sum _{i,j=1}^m\text{trace}({\mathbf{C}}_{ij}{\mathbf{R}}_i^{\top }{\mathbf{R}}_j),$$

$$\mathbf{D}=\sum _{i\sim j}\sum _{k=1}^{n_{ij}}{\mathbf{d}}_{ij}{\mathbf{d}}_{ij}^{\top }\text{and}\mathbf{C}=\mathbf{D}-\mathbf{B}{\mathbf{L}}^{†}{\mathbf{B}}^{\top }.$$

$${\mathbf{C}}_{ij}(p,q)=\mathbf{C}\left((i-1)d+p,(j-1)d+q\right),$$

$$\underset{\mathbf{R}}{\min }⟨\mathbf{C},{\mathbf{R}}^{\top }\mathbf{R}⟩\text{s.t.}{\mathbf{R}}_1,\dots ,{\mathbf{R}}_m\in \mathbb{SO}(d).$$

$${\mathbf{G}}_{ij}={\mathbf{R}}_i^{\top }{\mathbf{R}}_j(i,j\in [m]).$$

$$\mathbf{G}=\left[\begin{array}{cccc}{\mathbf{G}}_{11}& {\mathbf{G}}_{12}& \dots & {\mathbf{G}}_{1m}\\ {\mathbf{G}}_{21}& {\mathbf{G}}_{22}& \dots & {\mathbf{G}}_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{G}}_{m1}& {\mathbf{G}}_{m2}& \dots & {\mathbf{G}}_{mm}\end{array}\right].$$

$$\underset{\mathbf{G}}{\min }$$

$$\mathbf{G}\in {\mathbb{S}}_{+}^{dm},\mathrm{rank}(\mathbf{G})\le d,$$

$${\mathbf{G}}_{ii}=\mathbf{I},i\in [m],$$

$${\mathbf{G}}_{i,i+1}\in \mathbb{SO}(d),i\in [m-1].$$

$$\mathrm{rank}(\mathbf{G})\le d\text{and}{\mathbf{G}}_{i,i+1}\in \mathbb{SO}(d)$$

$$\Omega=\Big{\{}\mathbf{X}\in\mathbb{S}^{dm}_{+};\ \mathrm{rank}(\mathbf{X})\leq d\Big{\}},$$

$$\Theta=\Big{\{}\mathbf{X}_{ii}=\mathbf{I},i\in[m];\ \mathbf{X}_{j,j+1}\in\mathbb{SO}(d),j\in[m-1]\Big{\}}.$$

$$\underset{\mathbf{G},\mathbf{H}}{\min }$$

$$\mathbf{G}\in \Omega ,$$

$$\mathbf{H}\in \Theta ,$$

$$\mathbf{G}-\mathbf{H}=0.$$

$${\mathcal{L}}_{\rho }(\mathbf{G},\mathbf{H},\mathbf{\Lambda })=⟨\mathbf{C},\mathbf{G}⟩+⟨\mathbf{\Lambda },\mathbf{G}-\mathbf{H}⟩+\frac{\rho }{2}\Vert \mathbf{G}-\mathbf{H}{\Vert }_{\mathrm{F}}^2,$$

$${\mathbf{G}}^{k+1}=\underset{\mathbf{G}\in \Omega }{\mathrm{argmin}}{\mathcal{L}}_{\rho }(\mathbf{G},{\mathbf{H}}^k,{\mathbf{\Lambda }}^k),$$

$${\mathbf{H}}^{k+1}=\underset{\mathbf{H}\in \Theta }{\mathrm{argmin}}{\mathcal{L}}_{\rho }({\mathbf{G}}^{k+1},\mathbf{H},{\mathbf{\Lambda }}^k),$$

$${\mathbf{\Lambda }}^{k+1}={\mathbf{\Lambda }}^k+\rho ({\mathbf{G}}^{k+1}-{\mathbf{H}}^{k+1}).$$

$$\frac{\rho }{2}\Vert \mathbf{G}-({\mathbf{H}}^k-{\rho }^{-1}(\mathbf{C}+{\mathbf{\Lambda }}^k)){\Vert }_{\mathrm{F}}^2+\text{ constant },$$

$$\mathbf{G}^{k+1}=\Pi_{\Omega}\big{(}\mathbf{H}^{k}-\rho^{-1}(\mathbf{C}+\bm{\Lambda}^{k})\big{)},$$

$${\Pi }_{\Omega }(\mathbf{A})=\underset{\mathbf{X}\in \Omega }{\mathrm{argmin}}\Vert \mathbf{X}-\mathbf{A}{\Vert }_{\mathrm{F}}^2.$$

$${\mathbf{H}}^{k+1}={\Pi }_{\Theta }({\mathbf{G}}^{k+1}+{\rho }^{-1}{\mathbf{\Lambda }}^k).$$

Full text

I Introduction

Surface reconstruction from multiview scans has applications in computer-aided design, computer graphics, computer vision, medical imaging, and virtual reality. A range scanner is typically used to scan the object from different views. The scans are extracted by fixing the object and moving the scanner around, or by fixing the scanner and rotating the object on a turntable. Each scan is represented using a mesh, composed of vertices, faces, and normals (Levoy et al. (2005)). The vertices are simply the sampled points, the faces encode the connectivity between vertices, and the normals can be used to estimate the curvature (Guennebaud et al. (2008)). As with most existing reconstruction methods, we will use just the vertex information, i.e., each scan is represented as a point set. The computational problem in this regard is to register the scans using translations and rotations. The problem has two components, namely, finding the point-to-point correspondences between scans and determining the alignment of each scan using the found correspondences. The former is harder to deal with due to its combinatorial nature (Huang and Guibas (2013)).

On the other hand, if the point-to-point correspondence is known for two scans, then they can be registered simply using singular value decomposition (Umeyama (1991)). However, the correspondences are not known in practice and have to be inferred from the scan data. A natural approach is to solve the correspondence and registration problems in an alternating manner. This approach is used in classical ICP (Iterative Closest Point) and its variants (Besl and McKay (1992); Rusinkiewicz and Levoy (2001); Zinßer et al. (2003)), where the correspondences are obtained by aligning the scans and matching nearest points. The two steps are repeated until convergence. However, ICP requires a good estimate for the initial alignment. Deformations (scaling/stretching) in the scans further deteriorates the performance of ICP. To address this, Ying et al. (2009a, b) proposed an unified mathematical model for registration, by combining ICP with optimization over Lie groups. Classical ICP is also sensitive to outliers and partial overlaps between scans. Robust variants of ICP have been proposed that can reject outliers (Chetverikov et al. (2002); Zinßer et al. (2003); Phillips et al. (2007); Du et al. (2011); Zhang et al. (2011); Dong et al. (2014); Yang et al. (2016)). However, this involves pruning or reweighting the correspondences, which can be computation intensive. A more efficient alternative based on sparsity inducing norms was explored in Bouaziz et al. (2013).

The situation is more difficult with multiple scans. A straightforward approach is sequential registration in which multiple scans are aligned one at a time using ICP. However, sequential registration is prone to error propagation, and the situation gets worse in the presence of outliers. A more robust alternative is to take into account the pairwise correspondences between scans (obtained using ICP), integrate them in a single objective (cost function), and jointly optimize the transforms with respect to the objective. This way we can prevent error propagation and distribute the registration error across the scans. Our approach is based on this so-called global registration.

I-A Global Registration

There exist several methods for global registration in the literature. These can be divided into two broad classes depending on whether the registration is performed in frame or point space. Frame-space methods use the relative transform between scans (Sharp et al. (2002); Fusiello et al. (2002); Torsello et al. (2011); Govindu (2004); Govindu and Pooja (2014)), whereas point-space methods work directly with the local coordinates (Pennec (1996); Bergevin et al. (1996); Benjemaa and Schmitt (1999); Williams and Bennamoun (2001)). In one of the earliest frame-space methods (Sharp et al. (2002)), the view graph is decomposed into cycles such that the optimal transforms for each cycle can be computed in closed-form. In Fusiello et al. (2002), the rotations are parameterized using quaternions and the registration error is distributed evenly across all scans. In Torsello et al. (2011), dual quaternions are used to represent the transforms and geodesic averaging is used to denoise the relative transforms. In Govindu (2004) and Govindu and Pooja (2014), the authors used Lie-algebraic methods for averaging relative transforms. In Arrigoni et al. (2016b) and Bernard et al. (2015), it is shown that the null space of an appropriate matrix (constructed from the relative transforms) can be used to extract cycle-consistent transforms. In Arrigoni et al. (2016a), the authors proposed to use low-rank matrix completion for global registration. The pioneering point-space methods include Pennec (1996); Bergevin et al. (1996); Benjemaa and Schmitt (1999); Williams and Bennamoun (2001). The authors in Pennec (1996) alternately computed an average shape from the scans and aligned the scans against this shape. In Bergevin et al. (1996), the scans are repeatedly selected and registered against other scans in a global reference frame. The method in Bergevin et al. (1996) is enhanced in Benjemaa and Schmitt (1999) by accelerating the correspondence search using multi-z buffer technique, and by updating the surface positions immediately after the transform computation in each iteration. In Raghuramu (2015), the registration error is measured using the $\ell_{1}$ norm instead of the $\ell_{2}$ norm. A generalization of two-view ICP for multiple views was proposed in Fantoni et al. (2012). In Toldo et al. (2010), the authors combined ICP with general procrustes analysis for multiset registration. In Evangelidis et al. (2014), the point sets are modeled using the Gaussian mixture model and registration is performed using the EM algorithm.

I-B Contribution

In this paper, we consider the abstract problem of registering $m$ point sets over the group of rotations and translations, the special Euclidean group $\mathbb{SE}(d)$, given the local coordinates of points in each point set. We assume that there exist sets with common points, and that the point-to-point correspondences are known for such sets. We consider a least-squares formulation for this problem, where the variables are the $\mathbb{SE}(d)$-transforms associated with the point sets. Remarkably, while this problem is inherently nonconvex, its solution (global optimum) can be computed using SVD when $m=2$ (Umeyama (1991)). Unfortunately, such a solution is not available when $m\geq 3$, and one must resort to some iterative method. In this regard, the precise contributions are as follows:

• •

We observe that the constrained least-squares optimization can be reduced to an SDP (semidefinite program), where the objective is linear but the constraints are nonconvex. The solution of the original least-squares problem can be derived from the SDP solution via a linear map.

• •

We propose to solve the SDP using variable splitting and ADMM (alternating directions method of multipliers). In particular, the variable splitting is such that the resulting subproblems are given by matrix projections with closed-form solutions. For $m$ point sets, the per-iteration cost is essentially the partial eigendecomposition of an $md\times md$ matrix, and $m-1$ SVD-based projections onto the rotation group $\mathbb{SO}(d)$.

• •

We apply the proposed algorithm to the motivating problem, namely, the registration of multiview scans extracted from a three-dimensional surface. We model each scan as a point set and determine the pairwise correspondences using Picky-ICP Zinßer et al. (2003). Since the dominating computation per iteration is a partial eigendecomposition, we can scale our iterative solver to practical problems involving large number of scans. To demonstrate the applicability of the algorithm beyond multiview registration, we also use it for 2D shape matching.

Based on toy examples and simulated multiview data, we empirically demonstrate that the proposed solver has a large basin for global convergence and exhibits fast convergence for appropriate parameter settings. We present simulation results on scans from the Stanford repository and compare them with existing multiview methods. We observe that our algorithm is generally more robust to noise (both in the coordinates and the correspondences) and is quite fast.

I-C Related Work

The present contribution is an extension of Ahmed and Chaudhury (2017) where the registration is performed over the Euclidean group $\mathbb{E}(d)$, which includes reflections along with translations and rotations. Since the optimization is performed over a subgroup of $\mathbb{E}(d)$ in this paper, namely $\mathbb{SE}(d)$, we are required to introduce additional constraints. The novelty in this regard is the proposed variable splitting, whereby we are able to add more constraints and yet retain the efficiency of the ADMM solver in Ahmed and Chaudhury (2017). Similar to Ahmed and Chaudhury (2017), the subproblems of the proposed ADMM solver are given by closed-form projections. In the context of multiview registration, the ADMM algorithm in Ahmed and Chaudhury (2017) works well for small noise. However, when the noise is large, the optimal transforms returned by the algorithm often include reflections. This is possible since the domain is $\mathbb{E}(d)$, whereby both rotations and reflections are allowed. The present proposal fixes this problem by forcing the optimal transforms to be in $\mathbb{SE}(d)$.

We note that least-squares formulations of global registration have been considered in Williams and Bennamoun (2001); Krishnan et al. (2005); Chaudhury et al. (2015). The optimization is performed over $\mathbb{SE}(d)$ in Williams and Bennamoun (2001); Krishnan et al. (2005), and over $\mathbb{E}(d)$ in Chaudhury et al. (2015). The optimal rotations are iteratively computed using SVD in Williams and Bennamoun (2001), while Gauss-Newton iterations are performed on the $\mathbb{SO}(3)$ manifold to compute the optimal rotations in Krishnan et al. (2005). The solver in Krishnan et al. (2005) was computationally enhanced in Bonarrigo and Signoroni (2011), and later it was shown in Mateo et al. (2014) that the problem can be posed within a Bayesian framework. In a different direction, it was observed in Chaudhury et al. (2015) that least-squares registration over $\mathbb{E}(d)$ can be reduced to a rank-constrained SDP; this can further be relaxed to a convex SDP with provable tightness and stability guarantees.

We note that though ADMM has found wide applications in convex programming Boyd et al. (2011), the fact that it works well with nonconvex problems was reported more recently; e.g., see Chartrand and Wohlberg (2013); Miksik et al. (2014); Wang et al. (2018); Diamond et al. (2017). Moreover, while ADMM comes with with strong convergence guarantees for convex problems, analysis of nonconvex ADMM is still in its infancy (Wang et al. (2018, 2016)).

I-D Organization

The paper is organized as follows. In Section II, we give the mathematical description of the registration problem and the algorithmic solution. We empirically analyze the algorithm in Section III using toy examples and simulated data. In Section V, we present results on multiview registration using scans from the Stanford repository and compare our results with existing methods. We conclude the paper in Section VI.

I-E Notation

We use $[m]$ to denote $\{1,\ldots,m\}$. The Euclidean norm is denoted by $\lVert\bm{x}\rVert_{2}$. Symmetric matrices in $\mathbb{R}^{k\times k}$ are represented by $\mathbb{S}^{k}$, and the subset of positive semidefinite matrices by $\mathbb{S}^{k}_{+}$. We use the standard inner-product on $\mathbb{R}^{k\times k}$ given by $\langle\mathbf{A},\mathbf{B}\rangle=\text{trace}(\mathbf{A}^{\top}\!\mathbf{B})$; the norm induced by this inner-product is $\|\mathbf{A}\|_{\mathrm{F}}=\langle\mathbf{A},\mathbf{A}\rangle^{1/2}$. The orthogonal group $\mathbb{O}(d)$ consists of matrices $\mathbf{U}\in\mathbb{R}^{d\times d}$ such that $\mathbf{U}^{\top}\!\mathbf{U}=\mathbf{I}$, where $\mathbf{I}\in\mathbb{R}^{d\times d}$ is the identity matrix. The special orthogonal (or rotation) group $\mathbb{SO}(d)$ consists of matrices $\mathbf{R}\in\mathbb{R}^{d\times d}$ such that

[TABLE]

The special Euclidean group $\mathbb{SE}(d)$ consists of orientation-preserving rigid motions in $d$-dimensions, which are simply rotations and translations. An element of $\mathbb{SE}(d)$ is represented using a pair $(\mathbf{R},\bm{t})$, where $\bm{t}\in\mathbb{R}^{d}$ and $\mathbf{R}\in\mathbb{SO}(d)$.

II Registration over $\mathbb{SE}(d)$

We generally consider the global registration of point sets over $\mathbb{SE}(d)$, though we are mainly interested in $d=2$ and $d=3$. Suppose we have $m$ point sets $\mathcal{P}_{1},\ldots,\mathcal{P}_{m}\subset\mathbb{R}^{d}$. Two point sets are said to overlap if they have at least one common point. In general, not every pair of sets overlap. For $i,j\in[m]$, we use $i\sim j$ to mean that $\mathcal{P}_{i}\cap\mathcal{P}_{j}$ is not empty, and we denote the number of common points in this case by $n_{ij}$. Moreover, we denote the local coordinates of the common points in $\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ as

[TABLE]

In the ideal case, the hypothesis is that these points are related to each other via a transform from $\mathbb{SE}(d)$. In particular, if we associate the transform $(\mathbf{R}_{i},\bm{t}_{i})\in\mathbb{SE}(d)$ with $\mathcal{P}_{i}$, then

[TABLE]

Needless to say, the transforms are specified relative to some global reference frame.

II-A Least-Squares Formulation

In practice, the consistency relations (2) only hold approximately due to various imperfections. The task of computing the transforms can be posed as an optimization problem in this case. In particular, we consider the least-squares formulation

[TABLE]

By introducing the matrix variables

[TABLE]

we can express (3) as

[TABLE]

where $\bm{e}_{ij}=\bm{e}_{i}-\bm{e}_{j}$, $\bm{e}_{i}\in\mathbb{R}^{m}$ is the all-zeros vector with one at the $i$-th position, and

[TABLE]

where $\otimes$ is the Kronecker product. For fixed $\mathbf{R}$, the minimum of (5) is attained when $\mathbf{T}=-\mathbf{R}\mathbf{B}\mathbf{L}^{\dagger}$, where

[TABLE]

and $\mathbf{L}^{\dagger}$ is the Moore-Penrose pseudo-inverse of $\mathbf{L}$. Substituting $\mathbf{T}=-\mathbf{R}\mathbf{B}\mathbf{L}^{\dagger}$, the objective in (5) becomes

[TABLE]

where

[TABLE]

In (V-B), we use $\mathbf{C}_{ij}\in\mathbb{R}^{d\times d}$ to denote the $(i,j)$-th block (sub-matrix) of $\mathbf{C}$:

[TABLE]

where $i,j\in[m]$ and $p,q\in[d]$. In summary, we have reduced (5) to

[TABLE]

The reduction is in the following sense: if the minimizer of (7) is $\mathbf{R}^{*}$, then the minimizer of (5) is $(\mathbf{R}^{*},\mathbf{T}^{*}),\ \mathbf{T}^{*}=-\mathbf{R}^{*}\mathbf{B}\mathbf{L}^{\dagger}$. Note that the domain of (7) is $\mathbb{SO}(d)\times\cdots\times\mathbb{SO}(d)$, which is not convex. However, $\mathbb{SO}(d)$ has the structure of a smooth Riemannian manifold (Absil et al. (2009)), and this is be used to design efficient numerical solvers in Krishnan et al. (2005).

II-B Nonconvex SDP

By a change-of-variables, we can express (7) as an SDP, where the objective is linear but the constraints are nonconvex. In particular, this is done using the Gram matrix $\mathbf{G}\in\mathbb{R}^{dm\times dm}$ of the variables $\mathbf{R}_{1},\ldots,\mathbf{R}_{m}\in\mathbb{SO}(d)$, given by

[TABLE]

We reiterate that $\mathbf{G}_{ij}$ denotes the $(i,j)$-th block of $\mathbf{G}$, i.e.,

[TABLE]

Proposition II.1**.**

The Gram matrix has the following properties:

$(\mathrm{P1})$* $\mathbf{G}\in\mathbb{S}^{dm}_{+}$,*

$(\mathrm{P2)}$* $\mathrm{rank}(\mathbf{G})\leq d$,*

$(\mathrm{P3})$* $\mathbf{G}_{ii}=\mathbf{I},\ i\in[m]$, and*

$(\mathrm{P4})$* $\mathbf{G}_{i,i+1}\in\mathbb{SO}(d),\ i\in[m-1]$.*

Proof.

Since we can write $\mathbf{G}=\mathbf{R}^{\top}\!\mathbf{R}$, $(\mathrm{P1})$ is clear. Moreover, since $\mathbf{R}$ has full rank, and $\mathrm{rank}(\mathbf{G})=\mathrm{rank}(\mathbf{R})$, we in fact have equality in $(\mathrm{P2})$. Finally, note that $\mathbf{G}_{i,i+1}=\mathbf{R}_{i}^{\top}\mathbf{R}_{i+1}$, i.e., $\mathbf{G}_{i,i+1}$ is the product of rotations. In particular, $\mathbf{G}_{i,i}=\mathbf{R}_{i}^{\top}\mathbf{R}_{i}=\mathbf{I}$. This establishes $(\mathrm{P3})$ and $(\mathrm{P4})$. ∎

Apart from ($\mathrm{P1}$)-($\mathrm{P4}$), we can of course list other properties of $\mathbf{G}$. But the key observation is that ($\mathrm{P1}$)-($\mathrm{P4}$) are its essential properties, namely, they are sufficient to characterize $\mathbf{G}$ as a Gram matrix of rotations (see Appendix for the proof).

Theorem II.2**.**

If $\mathbf{G}$ satisfies $\mathrm{(P1)}\mbox{-}\mathrm{(P4)}$, then $\mathbf{G}$ is given by (8) for some $\mathbf{R}_{1},\ldots,\mathbf{R}_{m}\in\mathbb{SO}(d)$.

At this point, we note that Theorem II.2 remains valid if we use $\mathrm{rank}(\mathbf{G})=d$ in ($\mathrm{P2}$). The reason why we use the present formulation is that the set $\{\mathbf{G}\in\mathbb{S}^{dm}:\mathrm{rank}(\mathbf{G})\leq d\}$ is closed in $\mathbb{S}^{dm}$, i.e., it contains all its limit points. In contrast, the set $\{\mathbf{G}\in\mathbb{S}^{dm}:\mathrm{rank}(\mathbf{G})=d\}$ is not closed. For example, the sequence of matrices $\mathbf{I},(1/2)\mathbf{I},(1/3)\mathbf{I},\ldots$ are of rank $d$, but their limit is the zero matrix which has zero rank. The algorithmic implication of this technical point will be evident in the next section, where we will be required to compute the projection on a set. If the set is not closed, then the projection might not be defined.

Based on Theorem II.2, we substitute $\mathbf{G}=\mathbf{R}^{\top}\!\mathbf{R}$ in (7) and consider the following problem:

[TABLE]

Note that the objectives in (7) and (9) are identical. Moreover, following Theorem II.2, there exists a one-to-one correspondence between the domains of (7) and (9). Therefore, we have the following important observation.

Theorem II.3**.**

Problems (7) and (9) are equivalent, that is, $\mathbf{G}^{\star}$ is a minimizer of (7) if and only if $\mathbf{G}^{\star}$ a minimizer of (9).

The optimization in (9) is a nonconvex SDP (Diamond et al. (2017)), where the variable is positive semidefinite and the objective is linear. The constraint $\mathbf{G}_{ii}=\mathbf{I}$ is affine. However, the constraints

[TABLE]

are nonconvex. We develop an ADMM-based solver for (9). This is based on two crucial observations—the objective in (9) is linear, and the projections onto the feasible sets in (10) can be computed in closed-form.

II-C Variable Splitting and ADMM

Following the success of ADMM for convex programming (Boyd et al. (2011)), the ADMM framework has been extended to several nonconvex problems (Chartrand and Wohlberg (2013); Miksik et al. (2014); Diamond et al. (2017)). Preliminary theoretical results concerning the validity of such formal extensions have also been reported (Wang et al. (2018, 2016)). We propose an ADMM solver for (9) using variable splitting. In particular, we define the following subsets of $\mathbb{S}^{dm}$:

[TABLE]

and

[TABLE]

Simply stated, $\Theta$ consists of symmetric matrices whose diagonal blocks are identity and super-diagonal blocks are rotations. Note that we can equivalently write (9) as

[TABLE]

The purpose of this splitting is to group the constraints into two distinct classes, albeit at the expense of a linear constraint. For fixed $\rho>0$, the augmented Lagrangian for (11) is

[TABLE]

where $\bm{\Lambda}\in\mathbb{S}^{dm}$ is the dual variable for the constraint $\mathbf{G}-\mathbf{H}=\mathbf{0}$ (Boyd et al. (2011)). Starting with initializations $\mathbf{H}^{0}$ and $\bm{\Lambda}^{0}$, the ADMM solver uses the following sequence of updates for $k\geq 0$:

[TABLE]

[TABLE]

[TABLE]

Note that we can write the objective in (12) as

[TABLE]

where the constant does not depend on $\mathbf{G}$. Therefore,

[TABLE]

where $\Pi_{\Omega}(\mathbf{A})$ is the projection of $\mathbf{A}\in\mathbb{S}^{dm}$ onto $\Omega$:

[TABLE]

Similarly, it can be verified that

[TABLE]

In summary, we were able to express the primal subproblems as orthogonal projections. This is where the linearity of the objective in (11) plays an important role. At this point, note that both $\Omega$ and $\Theta$ are closed sets. Hence, $\Pi_{\Omega}$ and $\Pi_{\Theta}$ are well-defined. However, if we had defined $\Omega$ to be the set of rank-$d$ matrices, then $\Pi_{\Omega}(\mathbf{A})$ would not be defined if the rank of $\mathbf{A}$ were strictly less than $d$.

The complete algorithm is summarized in Algorithm 1. We now explain how steps 1 and 1 can be computed in closed-form.

II-D Matrix Projections

Recall that $\Omega$ and $\Theta$ are closed in $\mathbb{S}^{dm}$. As a result, projections (15) and (17) are well-defined. In particular, by adapting the Eckart-Young theorem Eckart and Young (1936) for positive semidefinite matrices, we can deduce the following.

Theorem II.4**.**

Let the eigendecomposition of $\mathbf{A}\in\mathbb{S}^{dm}$ be

[TABLE]

where $\lambda_{1}\geq\cdots\geq\lambda_{dm}$ are its eigenvalues, and $\bm{u}_{1},\ldots,\bm{u}_{dm}$ the corresponding eigenvectors. Then

[TABLE]

For completeness, we sketch the proof of Theorem II.4 in the Appendix. In practice, we can efficiently compute $\Pi_{\Omega}(\mathbf{A})$ using the power method (Golub and Van Loan (2012)), since we only require the top-$d$ eigenvalues/eigenvectors of $\mathbf{A}$. The relevance of property ($\mathrm{P2}$) is now evident. Namely, $\Pi_{\Omega}(\mathbf{A})$ might not be defined if instead of requiring the matrices in $\Omega$ to have rank at most $d$, we insist that the rank be exactly $d$. In particular, if $\mathbf{A}$ has less than $d$ positive eigenvalues, then $\Pi_{\Omega}(\mathbf{A})$ does not exist as the minimum in (16) is not attained in this case.

To compute $\Pi_{\Theta}(\mathbf{A})$, note that if $\mathbf{X}\in\Theta$, then we can write

[TABLE]

Following the definition of $\Theta$, it is then evident that

[TABLE]

Simply stated, we can determine $\Pi_{\Theta}(\mathbf{A})$ by setting the diagonal blocks of $\mathbf{A}$ to $\mathbf{I}$, projecting the super- and sub-diagonal blocks of $\mathbf{A}$ onto $\mathbb{SO}(d)$, and keeping the remaining blocks of $\mathbf{A}$ unchanged. Of course, since the resulting matrix is required to be symmetric, we simply project the $m-1$ super-diagonal blocks of $\mathbf{A}$ and set the sub-diagonal blocks using symmetry. In this regard, we record the following result.

Theorem II.5**.**

Let the SVD of $\mathbf{H}\in\mathbb{R}^{d\times d}$ be

[TABLE]

where $\sigma_{1}\geq\cdots\geq\sigma_{d}\geq 0$ and $\mathbf{U},\mathbf{V}\in\mathbb{O}(d)$. Then

[TABLE]

A proof of this result using the theory of Lagrange multipliers can be found in Umeyama (1991); Kabsch (1976). We provide a somewhat different proof in the Appendix.

We note that the optimum $\mathbf{H}^{\ast}$ obtained using Algorithm 1 has rank $d$ or more (this follows from the definition of $\Theta$). If the rank of $\mathbf{H}^{\ast}$ is exactly $d$, then it follows from Theorem II.2 that we can factor it as $\mathbf{H}^{\ast}=\mathbf{R}^{\top}\!\mathbf{R}$, where $\mathbf{R}\in\mathbb{R}^{d\times dm}$. This can be done using the eigendecomposition of $\mathbf{H}^{\ast}$; see Chaudhury et al. (2015). In particular, if we partition $\mathbf{R}$ into $m$ blocks of size $d\times d$, then each block would be a rotation matrix. However, if the rank of $\mathbf{H}^{\ast}$ is greater than $d$, we have to use some form of “rounding” to extract the rotations (which are no longer uniquely defined). In the present case, we have use spectral rounding (see Section $2.3$ in Chaudhury et al. (2015)).

II-E Complexity Analysis

The breakup of the computation complexity for the proposed method is given in Table I. The one-time cost of building the matrix $C$ is $O(m^{2}(m+d+d^{2}))=O(m^{3})$, since $m\gg d$ in practice. The cost per iteration is dominated by the projections steps $4$ and $5$ in Table I. For step $5$, we need to find the top $d$ eigenvectors of an $md\times md$ matrix. This can be done efficiently using Arnoldi iterations (Golub and Van Loan (1996)), where the per-iteration cost is $O(m^{2}d^{2})$. On the other hand, we need to compute the SVD of $m-1$ matrices in step $5$, each of size $d\times d$. This can be done efficiently using QR-type algorithms (Golub and Van Loan (1996)). Finally, a single SVD is required in step $6$ (one time cost). In summary, the per-iteration cost of Algorithm 1 is essentially the partial eigendecomposition of a $md\times md$ matrix, and the SVD of $m-1$ matrices of size $d\times d$. As a result, we can scale up our algorithm to problems involving large number of point sets.

III Numerical Results

In this section, we present numerical results using simulated point sets for which the point-to-point correspondences are known. This allows us to factor out the task of finding the correspondences, which will be addressed in Section V. Moreover, this allows us to manipulate the true correspondences and access their impact on the registration. In particular, we empirically investigate the following questions:

• •

Are the iterates in Algorithm 1 convergent?

• •

If so, do they converge to the global minimum of (9), or simply to some saddle point?

• •

How does the initialization and $\rho$ affect convergence?

• •

Is the algorithm stable to perturbations of coordinates and correspondences?

The first three questions can be theoretically resolved for convex problems (Boyd et al. (2011)). However, the situation is more complicated with nonconvex problems, and only some preliminary results have been established (Wang et al. (2018, 2016)). Similar to (Chartrand and Wohlberg (2013); Miksik et al. (2014); Wang et al. (2018); Diamond et al. (2017)), we will demonstrate the practical utility of the ADMM solver. In fact, as with the above papers, we will show that the solver works well in practice.

To initialize Algorithm 1, we set $\bm{\Lambda}$ as the zero matrix. The solution of the spectral relaxation of (9) is used for $\mathbf{H}$ (Chaudhury et al. (2015)). The latter requires us to compute the smallest $d$ eigenvectors of $\mathbf{C}$. We will sometimes initialize $\mathbf{H}$ with the all-identity matrix (a matrix of size $md\times md$ whose $m^{2}$ blocks are the $d\times d$ identity matrix) for some experiments, which is weaker than the spectral initialization.

To address the second question, we somehow need to compute the global minimum of (9). An useful observation is that this can be done for two point sets. In this case, the minimizer is simply given by the Gram matrix of the optimal rotations obtained using Umeyama’s algorithm Umeyama (1991). This will be particularly useful for analyzing the stability of our algorithm, since it is non-trivial to determine the minimum of (9) in the presence of noise. On the other hand, note that the minimum of (9) is simply zero with simulated data (without noise). In this case, the minimizer is the Gram matrix of the ground-truth rotations.

Since Umeyama’s original formulation (Umeyama (1991)) includes scaling (along with rotation and translation), we describe the solution over $\mathbb{SE}(d)$ for completeness. For $m=2$, we can write (3) as

[TABLE]

where $n$ is the number of common points, and $(\bm{x}^{k})$ and $(\bm{y}^{k})$ are the respective local coordinates.

The problem can be simplified by fixing one point set and computing the relative rotation and translation of the other. In particular, we can take $\mathbf{R}_{1}=\mathbf{I}$ , $\bm{t}_{1}=\bm{0}$ and $\mathbf{R}_{2}=\mathbf{R}$ , $\bm{t}_{2}=t$. That is, we consider the following simplification of (18):

[TABLE]

As shown in Umeyama (1991), the minimizers of (19) are

[TABLE]

where

[TABLE]

and

[TABLE]

For our first experiment, we consider the example in Umeyama (1991) involving two point sets. Each set has three points whose coordinates are $\{(0,0),(1,0),(0,2)\}$ and $\{(0,0),(-1,0),(0,2)\}$; one point set is simply a reflection of the other (see Figure 1). Clearly, the minimum of (19) cannot be zero in this case, since the points cannot be perfectly aligned using just translations and rotations. In fact, the optimal value corresponding to (20) is $3.7185$ (up to four decimal places).

We next solve the above problem using Algorithm 1, where we work with formulation (18). Note that the algorithm does not make use of Umeyama’s solution; it iteratively computes the solution starting from some initialization. Therefore, this simple problem is nonetheless non-trivial for Algorithm 1. For different values of $\rho$, the objective at each iteration is shown in Figure 2, where we initialize using $\mathbf{H}=\mathbf{I}$. Note that, while the convergence speed changes with $\rho$ (as in convex ADMM), the objective asymptotically converges to the minimum (indicted in the figure) obtained using Umeyama’s formula. The solution obtained using our algorithm is depicted in Figure 1, which agrees perfectly with the result in Umeyama (1991).

We next experiment with three-dimensional models from the Stanford repository Levoy et al. (2005). To extract overlapping point sets from a given model (ground-truth), we follow the process in Evangelidis et al. (2014). We first center the model by subtracting its centroid, and then rotate it about the $x$-axis by $m$ different angles. For a fixed rotation, the points above the $x\mbox{-}y$ plane are formed into a point set. After creating $m$ such point sets, we randomly rotate and translate them. Obviously, we know the exact correspondences in this case.

We test the robustness of our algorithm by (i) corrupting the coordinates with additive Gaussian noise, and (ii) introducing false correspondences (outliers). The goal is to mimic real world scenarios in which the scanned measurements are invariably noisy and some of the correspondences are wrongly estimated. To objectively assess the reconstruction quality, we use the rotation error from Arrigoni et al. (2016a). If $(\mathbf{R}_{i})$ are the true rotations and $(\hat{\mathbf{R}}_{i})$ are the rotations obtained using our algorithm, we first remove the global rotation by fixing one of the rotations (say, the first one) and performing the corrections $\mathbf{R}_{i}^{\prime}=\mathbf{R}_{1}^{\top}\mathbf{R}_{i}$ and $\hat{\mathbf{R}}_{i}^{\prime}=\hat{\mathbf{R}}_{1}^{\top}\hat{\mathbf{R}}_{i}$. The rotation error is defined as

[TABLE]

where

[TABLE]

is the geodesic distance on $\mathbb{SO}(3)$ (Arrigoni et al. (2016a)).

To perturb the coordinates, we add isotropic Gaussian noise of variance $\sigma^{2}$ to the points in each set. For the correspondence noise, we first fix some $\eta\in[0,1]$ and randomly shuffle $\eta$-fraction of the known correspondences (these are the outliers), keeping other correspondences fixed. The results obtained on the models of Bunny and Buddha are shown in Figures 3 and 4. We have used $10$ scans per model for the experiments. Notice that the rotation error scales gracefully with increase in $\sigma$ and $\eta$. Moreover, thanks to the global registration framework, we can obtain accurate reconstructions even with $50\%$ outliers. This is further demonstrated using some visual results in Figure 5. Some more visual results for Bunny, Buddha, Armadillo and Dragon are shown in Table III. In each case, we have used six point sets. The point sets before and after registration are shown in the figure along with the ground truth.

As in Umeyama’s example, the ADMM objective was found to converge to the global minimum (zero) in the noiseless case. Since it is not possible to determine the global minimum in general for multiple point sets, we cannot decide on optimality in the presence of noise. However, the iterates were found to converge to a fixed point in such cases.

Going back to the questions posed at the start of this section, we have empirically noticed the following:

• •

The iterates of Algorithm 1 converge to a fixed point for any arbitrary initialization when $\rho>0$. The convergence is generally fast if $\rho\in[1,10]$ and spectral initialization is used.

• •

For special cases where we know the global minimum of (9), the objective remarkably converges to the global minimum if we use the spectral initialization. Generally, the algorithm always seemed to work in the noiseless setting.

• •

The algorithm behaves stably with perturbations in the coordinates and correspondences.

For completeness, we compare the proposed algorithm with Ahmed and Chaudhury (2017), where the optimization in (7) is performed over $\mathbb{O}(d)$. We extract ten point sets from Bunny and randomly shuffle $60\%$ of the correspondences. We feed this data into the algorithms and check the determinants of the computed transforms. The results of an experiment are presented in Table II. Notice that the transforms computed by our algorithm are indeed rotations, whereas the algorithm in Ahmed and Chaudhury (2017) returns a mix of rotations and reflections. Moreover, the reconstruction error is much higher for the latter (see Figure 6).

IV Shape Matching

We now apply our method for matching $2$D shapes (Belongie et al. (2002)). Without getting into a rigorous analysis, we simply present few results to demonstrate that our method can be used for shape matching. In particular, we wish to match a collection of shapes of a $2$D model without having to compute all pairwise similarities. We have used the models plane, car and bicego from the hmm-gdb database111Download from http://visionlab.uta.edu/shape_data.htm (accessed on $13$ Nov, $2018$).. The scans are first centered and arbitrarily labeled. Picky-ICP is then used to estimate the correspondences between successive pairs of scans. Finally, the scans are registered using our method. The results are shown in Figure 8. Notice that the scans do not match perfectly. This is expected since the scans are originally deformed and also because we use only rigid transforms. Nevertheless, the overall matching appears to be reasonably good.

V Multiview Registration

We next use the proposed algorithm for the registration of 3D scans (Sharp et al. (2002)). The important consideration here is that the point-to-point correspondences need to be estimated from the scan data. We propose to use two-scan registration for the same which is discussed next.

V-A Correspondence Estimation

As mentioned, the scans extracted from a model are represented using a mesh (Levoy et al. (2005)). We treat each scan as a point cloud, where the points are simply the mesh vertices. We first determine which pairs of scans overlap and the correspondences between them. This information is supplied to our registration algorithm. Note that, while we determine the correspondences in a pairwise manner, the registration algorithm takes all the pairwise correspondences into account. To find pairwise correspondences, we can use ICP or its fast variant (Besl and McKay (1992); Rusinkiewicz and Levoy (2001)). However, we noticed in our experiments that these are sensitive to outliers. After trying different methods, we found that the correspondences obtained using Picky-ICP Zinßer et al. (2003) give good reconstructions for our registration algorithm. In ICP, multiple points from scan $\mathcal{P}_{i}$ are often assigned to a single point from some target scan $\mathcal{P}_{j}$. However, the correspondences between $\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ should ideally be one-to-one. In Picky-ICP, the correspondences between $\mathcal{P}_{i}$ and $\mathcal{P}_{j}$ are first estimated as in ICP. Multiple assignments are then resolved by selecting the point in $\mathcal{P}_{i}$ (among several candidates) that is closest to the matching point in $\mathcal{P}_{j}$ (ties are randomly broken). Let $(d_{k})$ be the distances between corresponding points and $s$ be their standard deviation. Pairs for which $d_{k}$ is within a certain factor of $s$ are considered as overlapping points, and the remaining points are discarded. In our case, we set the factor as three. We simply used PickyICP as a black box and have not engineered anything on our own. Needless to say, if a better method is used for finding correspondences, the performance of our registration is expected to improve.

V-B Comparisons

We report results on four datasets from four datasets: Bunny (Turk and Levoy (1994)), Buddha, Dragon (Curless and Levoy (1996)) and Armadillo (Krishnamurthy and Levoy (1996)). We also compare with some recent methods for multiview registration: MAICP (Govindu and Pooja (2014)), LRS (Arrigoni et al. (2016a)), and JRMPC (Evangelidis et al. (2014)). These methods have already been demonstrated to perform better than Sharp et al. (2002); Torsello et al. (2011); Benjemaa and Schmitt (1999); Williams and Bennamoun (2001); Bernard et al. (2015). Codes for LRS and JRMPC are available online and that of MAICP was provided by the authors. All the competing methods were run using default parameters.

Comparison 1. We first compare the reconstructions using the scans from Stanford dataset. For a fair comparison, we have initialized all methods using PickyICP Zinßer et al. (2003). The reconstruction accuracy is assessed using the error metric in (21). The results are presented in Table V. It is evident that our method performs much better than LRS and JRMPC. The performance is generally comparable to MAICP. Note that the execution time for our method is significantly less. This aspect is especially important when registering several scans. For a visual comparison, the cross-sections of the reconstructions are compared in Table IV. Since the rotation errors for LRS and JRMPC are large, we have only shown the cross-sections for MAICP and our method in Table IV. The number of scans and angular differences are as follows: Bunny ($12$ scans, $30$ degree), Armadillo ($12$ scans, $30$ degree), Buddha ($15$ scans, $24$ degree), and Dragon ($15$ scans, $24$ degree). Notice that the cross-section for Buddha is much better for our reconstruction. Following JRMPC, LRS and MAICP, we have tried comparing the convergence rate (of the rotation error) using the following protocol:

read the full scans. 2. 2.

initialize all methods using PickyICP. 3. 3.

run LRS, JRMPC, MAICP and our method on the scans. 4. 4.

record the rotation error at each iteration.

The results are shown in Figure 9. Notice that the rotation error decreases quickly for our method and MAICP. The rotation error at convergence is comparable for MAICP and our method. However, the timing is significantly lower in our case. Note that $100$ iterations were used for all methods for comparison, though fewer iterations are required for our method and MAICP (cf. Figure 9). Thus only $10$ iterations were run MAICP and proposed method in Table V.

Comparison 2. To test robustness, we perturb the individual scans by Gaussian noise and register the scans with all the methods to test for robustness. Some typical results are reported in Figure 7. The performance of the proposed method is generally comparable to MAICP. But again, the timings are significantly less for our method.

Comparison 3. Finally, we carry out experiments with random noise (in the rotations and coordinates) using the following protocol:

read the model. 2. 2.

center the model by subtracting its centroid and rotate it about the $x$-axis by angles $\theta,2\theta,\ldots$ for some fixed $\theta$ (this mimics the collection of point sets using a turn table). For each rotation, the points above the $x$-$y$ plane are formed into a point set. 3. 3.

perturb the coordinates with additive Gaussian noise.

The trials are carried out $25$ times for each noise setting and the resulting errors are averaged. The angle of rotation $\theta$ determines the overlap between successive scans. For all experiments, we fixed $\theta$ to be $15$ degrees. We generate $24$ scans for Dragon and Buddha. Standard deviation for Gaussian noise is $0.01$, while the true rotations are randomly perturbed in the range $[0\mbox{-}1]$ degrees. Reconstruction results from MAICP and our method are shown in Figure 10 for one of the trials. PickyICP is used for initialization in both methods. The averaged rotation error is slightly better for Dragon in our case (our: $1.39$, MAICP: $1.45$). The visual quality also appears to be better for our method (see Figure 10). The timings are $1.49$ sec and $23.2$ min for our method and MAICP. The results are almost identical for Buddha (see Figure 10). Average rotation errors for our method and MAICP are $1.486$ and $1.485$. However, the timing is significantly lower for our method (our: $1.39$ sec, MAICP: $21.9$ min).

We next simulate scenarios where scans covering the entire model may not be available. For this purpose, we generate $20$ scans from Buddha and Armadillo. The scans for Armadillo are perturbed with Gaussian noise having standard deviation $0.03$. No Gaussian noise is added for Buddha in order to verify reconstruction with just rotation noise. Each scan is randomly rotated in the range $[0\mbox{-}2]$ degrees. As earlier, we use $25$ trials and average the rotation error. PickyICP is used to initialize both methods. Average rotation errors for Buddha are $1.9$ and $2.2$ for our method and MAICP. For Armadillo, the corresponding rotation errors are $1.59$ and $1.63$ respectively. The registration results for one of the trials is shown in Figure 11. Notice that our reconstruction result is much better than MAICP for both models. The possible reason for this is that MAICP uses correspondences between the first and and the last scan. Since full scans are not considered, there are very few overlapping points between them. Again, our method is significantly faster: $1.06$ sec and $0.43$ sec for Buddha and Armadillo compared to $18.4$ min and $17.6$ min for MAICP.

VI Conclusion

We proposed a novel optimization algorithm for registering multiple points sets in a globally consistent fashion using rotations and translations. We empirically analyzed the algorithm and showed that it works well on both simulated and real scan data. An intriguing property of the proposed ADMM solver is that it converges to the global minimum with a good initialization. We could certify this for cases where the global minimum can be computed by some other means, namely, for clean data (where the minimum is zero) and for two scans. We applied the proposed algorithm for 2D shape matching and 3D multiview registration. For the later we compared its performance with recent methods. The reconstruction accuracy of the proposed algorithm is comparable to MAICP but we are much faster. The results also suggest that the overall algorithm is quite robust to noise.

VII Appendix

VII-A Proof of Theorem II.2

It follows from ($\mathrm{P3}$) that $\mathrm{rank}(\mathbf{G})\geq d$. Along with ($\mathrm{P2}$), we conclude that $\mathrm{rank}(\mathbf{G})=d$. Therefore, using ($\mathrm{P1}$) and the spectral theorem for symmetric matrices, we can write

[TABLE]

where $\lambda_{i}>0$ and $\bm{u}_{1},\ldots,\bm{u}_{d}$ is a orthonormal basis of $\mathbb{R}^{dm}$. Define $\mathbf{R}\in\mathbb{R}^{d\times dm}$ to be

[TABLE]

and let $\mathbf{R}=[\mathbf{R}_{1}\cdots\mathbf{R}_{m}]$, where $\mathbf{R}_{i}\in\mathbb{R}^{d\times d}$. By construction, $\mathbf{G}=\mathbf{R}^{\top}\mathbf{R}$ and, in particular,

[TABLE]

Therefore, it follows from ($\mathrm{P3}$) that $\mathbf{R}_{i}^{\top}\mathbf{R}_{i}=\mathbf{I}$. Furthermore, we conclude from ($\mathrm{P4}$) that for $i=1,\ldots,m-1$,

[TABLE]

We deduce that for $i\in[m]$, $\mathrm{det}(\mathbf{R}_{i})=1$ or $\mathrm{det}(\mathbf{R}_{i})=-1$. In the latter case, we simply pick a global reflection $\mathbf{Q}$ with $\mathrm{det}(\mathbf{Q})=-1$, and reassign $\mathbf{R}_{i}\leftarrow\mathbf{Q}\mathbf{R}_{i}$. This gives us the desired $\mathbf{R}_{1},\ldots,\mathbf{R}_{m}\in\mathbb{SO}(d)$ such that (8) holds.

VII-B Proof of Theorem II.4

We can write $\mathbf{A}=\mathbf{U}\Lambda\mathbf{U}^{\top}$, where

[TABLE]

Let $\Psi$ denote matrices in $\mathbb{S}_{+}^{dm}$ with rank at most $d$. Any $\mathbf{X}\in\Psi$ can be represented as $\mathbf{X}=\mathbf{V}\Gamma\mathbf{V}^{\top}$,

[TABLE]

where $\mu_{1}\geq\cdots\geq\mu_{dm}\geq 0$ and at most $d$ of these are positive.

Note that $\lVert\mathbf{X}-\mathbf{A}\rVert_{\mathrm{F}}=\lVert\mathbf{K}\Gamma\mathbf{K}^{\top}-\Lambda\rVert_{\mathrm{F}}$, where $\mathbf{K}=\mathbf{U}^{\top}\mathbf{V}$. As a result,

[TABLE]

where $\Phi$ denotes non-negative diagonal matrices with rank at most $d$. For fixed $\Lambda\in\Phi$, it can be shown that the minimum over $\mathbf{K}\in\mathbb{O}(d)$ is attained when $\mathbf{K}=\mathbf{I}$, that is, when $\mathbf{V}=\mathbf{U}$. In particular, this reduces (22) to

[TABLE]

The minimizer of (23) is given by the projection $\Gamma=\Pi_{\Phi}(\Lambda)$, that is, $\mu_{i}=\max(0,\lambda_{i})\text{ for }1\leq i\leq d$, and $\mu_{i}=0\text{ for }d+1\leq i\leq dm$.

VII-C Proof of Theorem II.5

The projection problem in question is

[TABLE]

where we have used the fact that $\mathbf{X}\in\mathbb{SO}(d)$. Since $\mathbb{SO}(d)$ is compact, there exists $\mathbf{X}_{0}\in\mathbb{SO}(d)$ such that

[TABLE]

We claim that $\mathbf{P}=\mathbf{A}^{\top}\!\mathbf{X}_{0}\in\mathbb{S}^{d}$. Indeed, consider an arbitrary anti-symmetric matrix $\mathbf{M}\in\mathbb{R}^{d\times d}$ such that $\mathbf{M}^{\top}=-\mathbf{M}$, and define

[TABLE]

Since $e^{t\mathbf{M}}\in\mathbb{SO}(d)$, it follows that for all $t\in\mathbb{R}$,

[TABLE]

Hence,

[TABLE]

Since (24) holds for any anti-symmetric $\mathbf{M}$, it easily follows that $\mathbf{P}\in\mathbb{S}^{d}$.

Having shown that $\mathbf{P}$ is symmetric, let $\mathbf{P}=\mathbf{Q}\Lambda\mathbf{Q}^{\top}$, where $\Lambda$ is a diagonal matrix and $\mathbf{Q}\in\mathbb{O}(d)$. Then

[TABLE]

where $\mathbf{K}=\mathbf{X}_{0}^{\top}\!\mathbf{U}\in\mathbf{O}(d)$. Since $\mathbf{P}$ and $\mathbf{P}^{2}$ commute, they can be diagonalized in the same basis. In particular,

[TABLE]

In terms of this representation, we have

[TABLE]

Since $\sigma_{i}\geq 0$, (25) is maximum when each $s_{i}=1$. However, since $\mathbf{X}_{0}\in\mathbb{SO}(d)$, it follows that

[TABLE]

Therefore, if $\mathrm{det}(\mathbf{A})\neq 0$, we must have

[TABLE]

since $\mathbf{U}\mathbf{V}^{\top}\in\mathbb{O}(d)$. If $\mathrm{det}(\mathbf{U}\mathbf{V}^{\top})=1$, then (25) is maximum when $s_{i}=1$ for all $i$. However, if $\mathrm{det}(\mathbf{U}\mathbf{V}^{\top})=-1$, then it follows from (26) that the maximizer is $s_{i}=1\text{ for }1\leq i\leq d-1$, and $s_{d}=-1$ (this is where we use the fact that $\sigma_{i}\geq\sigma_{d}$ for all $i$). The case $\mathrm{det}(\mathbf{A})=0$ can be worked out similarly.