Fast and robust motion averaging via angle constraints of multi-view range scans

The motion averaging problem arises when recovering the motions for a set of scans in three-dimensional (3D) reconstruction. Motion averaging is taken as an optimisation problem, whose solution can be obtained via matrix decomposition in principal component analysis. The latest motion averaging methods utilize the Cauchy kernel function or overlap rate to construct the weighted matrix in the cost function, yielding weighted low-rank sparse matrix decomposition solutions. Nevertheless, the over-reliance on the initial value and high computational complexity remain two major problems. Here, we aim to improve the motion averaging performance via a fast and robust method based on the angle constraints of multiple scans. The additional term representing the angle constraint is added into the original cost function. By doing this, two beneﬁts are facilitated: 1) the iteration times for the optimisation process are reduced; and 2) it effectively mitigates the over-reliance on the initial value and thus caters for outliers and random noise of initial motions. Compared with state-of-the-art motion averaging methods, it is indicated that the proposed method is more efﬁcient and robust for numerical simulations of using both simulated motions and real 3D scanning data, given a large range of angle constraint values.


INTRODUCTION
Three-dimensional (3D) reconstruction is a hot topic in the field of computer vision, which is widely adopted in driverless technology, virtual reality technology, object recognition and tracking etc. The mission of 3D reconstruction is to move multiple views into one model through translation and rotation. The standard estimation of motion is based on iterative closest point (ICP) [1] or perspective geometry, which can obtain the relative motion of two different views. However, due to the influence of noise and outliers of initial motions, there are some errors in the recovery of these relative motions to the global coordinate system. Therefore, the major topic is to recover individual motions for a set of scans when a number of relative motion estimates between scan pairs are given. This problem is called motion averaging, also known as motion registration or synchronisation. A large number of alternative approaches have been developed to solve the motion averaging problem, which can be divided into the following categories: the iterative solution method based on Lie group structure and unit quaternion, the Govindu [2], Fusiello et al. [3] and Fantoni et al. [4] described rotation based on quaternions, which can avoid singularity and transform the solution of global motion into the minimisation objective function. Govindu [2] first proposed a motion averaging approach by singular value decomposition (SVD) and least-squares algorithm. Unfortunately, this method is nonrobust and only applicable to the ideal situation without noise. The iterative solution method based on Special Orthogonal group (SO(3)) and Special Euclidean group (SE(3)) represents the rigid motion in 3D space with a linear expression of six degrees of freedom. Hartley et al. [5] proposed the Weiszfeld method to decrease the influence of outliers. However, updating each vertex in the viewgraph resulted in slow convergence in the iteration. To improve the robustness, Govindu and Chatterjee [6][7][8] focused on large-scale 3D rotation of the Lie group structure. Chatterjee and Govindu [7] put forward the ℓ1-norm solution as the initial value of iteration before using the iteratively reweighted least-squares (IRLS) method; it was proved as a Quasi-Newton method in [8]. In addition, Guo et al. [9] estimated the overlapping percentage of each scan pair and proposed the weighted averaging method based on the Riemannian metric Lie group to enhance robustness. Wang and Singer [10] used the sum of unsquared deviations as a penalty function to enhance the accuracy and stability of motion recovery. However, this method focused on accuracy rather than efficiency, and they emphasized complete measurements and a specific noise model. Li et al. [11] formulated the cost function as subtraction operation of the rotation vector, which was solved by the weighted total least squares and the IRLS method. However, all of these least-squares methods based on the Lie group require exact initial parameters to prevent local minimum and have high computational complexity in iteration. Recently, Arrigoni et al. [12,13] have introduced the low-rank sparse (LRS) matrix decomposition into the motion averaging problem. According to the analysis of Martinec and Pajdla [14], the matrix-based approach is generally more effective than the minimisation method of quaternions. In the LRS matrix decomposition framework, the noise matrix is eliminated by decomposing the reconstruction matrix composed of the relative motion, and the global motion is restored by considering the missing value and abnormal data. In terms of robustness, Arrigoni et al. [13] proposed a robust LRS decomposition and matrix complement algorithm R-GoDec, which can eliminate noise, outliers and missing data in SE(3). However, the LRS matrix decomposition methods are sensitive to the sparsity of the reconstructed matrix and require an accurate initial value of relative motion. With the developing of computer vision, some up-to-date motion averaging methods were proposed. Huang et al. [15] proposed a neural network to learn the weights associated with each relative motion and reduced the heavily reliability on the quality of the input relative motions. Purkait et al. [16] built a graph-based neural network that learns the noise patterns from the data and replaces robust optimisation methods. However, the network was trained on a large number of synthetic graphs.
As for 3D reconstruction, the motion averaging method is wildly used in global structure from motion (SfM) and point clouds registration. For instance, Arrigoni et al. [12] proposed the EIG-SE(3)+IRLS approach, which utilized the Cauchy weight function and spectral decomposition, and the performance on real data of SfM and point clouds registration was outstanding. In global SfM, the translation vector is calculated by global rotation; thus, the quality of motion recovery is dependent on the rotation averaging algorithm. The mainstream in SfM is based on the Lie group method, and some details can be found in [2, 5-8, 10, 11]. In application of point clouds registration, the initial relative motions in the motion averaging problem are mostly estimated via the ICP algorithm. And some researchers ameliorated the cost function by adding weight of scan pairs or using different distances between corresponding points, such as state-of-the-art works in [3,9,12,13,17,20,26]. Jin et al. [17] proposed the weighted LRS method of multiview registration and reduce the sparsity of the reconstructed matrix by using the anti-symmetry of the relative motion and overlap rate of scan pairs. However, the weighted matrix of overlap rate is based on the trimmed iterative closest point (TrICP) algorithm; thus, generating the initial registration value in each iteration resulted in a large amount of computation. Besides, some methods based on the probabilistic model were emerged.
Mateo et al. [18] and Zach et al. [19] introduced the Bayesian framework into registration problems. Ma et al. [20] defined a Gaussian variational mixture model and used anisotropic components to weaken the effect of outliers in registration.
Consequently, motion averaging is a crucial part in multi-view 3D reconstruction. And the motion averaging method based on the Lie group, LRS matrix decomposition or others has heavy computing burden when solving the cost function iteratively. To sum up, all the above motion averaging methods have shown limitations in terms of efficiency, robustness or both. Therefore, the aim here is to effectively eliminate the over-reliance on the initial value and computational complexity so as to decrease the influence of outliers and random noise in initial motions and speed up runtime. This paper proposes a fast and robust motion averaging algorithm based on the angle constraints of multi-view range scans. Through the definition of angle constraint, we established the weighted LRS matrix decomposition model, which was added the additional term into the original cost function. Then, we solved the optimisation problem by the augmented Lagrange multiplier (ALM) method. The key innovations of the proposed method include the following.

High efficiency:
We put forward a novel idea, which considers the constraint of motion itself, thus reducing the computational complexity significantly. 2. Robustness: With the combination of angle constraint and weighted LRS matrix decomposition framework, the proposed method is robust to random noise and outliers when solving the cost function in motion averaging.
The rest of this paper is organized as follows. Section 2 introduces the motion averaging model of multiple views based on the LRS matrix decomposition. Then, the new motion averaging method is described in Section 3, and we add the angle constraint as an additional term of the motion into the cost function, which is based on the weighted LRS matrix. In Section 4, the proposed method is tested and developed in both motion averaging experiment by simulated data and 3D model reconstruction experiment by real data. Finally, the conclusions are briefly presented in Section 5.

LRS MATRIX DECOMPOSITION
This section describes the mathematical definition of motion registration based on the LRS matrix decomposition. In general, the motion registration problem can be solved by the robust principal component analysis (RPCA).

Problem statement
In this part, we adopt the recently introduced LRS framework [12,13] as our baseline. Suppose that there are N scans, which are acquired from multiple views. Let M i ∈ SE(3) represent the global rigid transformation, namely, the absolute motion of the ith scan where R i ∈ SO(3) denotes the rotation matrix and T i ∈ ℝ 3 denotes the translation vector; thus, the rank of M i is 4. The main objective is to optimize all of the initial global motion and improve the accuracy of motion recovery. Suppose that the first view is a global coordinate; the N − 1 remaining scans will be transformed to the global coordinate. Let M i j ∈ SE(3) denote the relative motion, which represents the rigid transform between the view of i and j , respectively where M i −1 represents the inverse of absolute motion M i . Obviously, it satisfies the constraint rank(M i j ) = 4.
In practical application, (2) is not satisfied due to the random noise and outliers in relative motions. In this case, let M i j ∈ SE(3) denote the theoretically estimate of M i j . Hence, the motion averaging problem can be defined as the following minimisation problem: where ‖ ⋅ ‖ F denotes the Frobenius norm of the matrix. In order to solve (3), three block matrices U, V and X are introduced as follows: where I denotes the 4 × 4 identity matrix, X is the reconstructed LRS relative motion matrix and the rank of X is also 4. Apparently, the constraints can be represented as X = UV and U T U = I 4 ; thus, X is a positive-semidefinite and symmetric matrix. Accordingly, we address the following problem in the noisy case: where E denotes the noisy sparse matrix,X denotes the estimation of relative motion matrix X ,  Ω denotes the projection of error matrix E onto the indicator matrix Ω and Ω implies whether the elements of reconstructed matrixX can be used.

Robust principle component analysis
The optimisation problem of motion registration mentioned in Section 2.1 can be solved by RPCA [21]. Therefore, we obtain an equivalent expression of (6) as follows: where rank(X ) means the rank of X and ‖ ⋅ ‖ 0 denotes the l 0 norm of matrix. Moreover, the problem mentioned above is not convex because of low-rank matrix X and sparse matrix E. The nuclear norm is the convex relaxation of rank and using l 1 norm to obtain the optimal convex approximation of l 0 norm. Thus, (7) can be reformulated as where U is satisfied U T U = I 4 . Theoretically, the arbitrary matrix recover algorithm can solve the problem in (8), such as the ALM method proposed by Lin et al. [22], the L1-ALM method proposed by Zheng et al. [23] and the Grasta algorithm proposed by He et al. [24].

THE NEW METHOD
In this section, we reformulate the cost function in terms of weighted LRS matrix and angle constraints. And the optimisation problem can be accurately and efficiently solved by the ALM method.

Angle constraints of relative motion
Suppose that the range scans from 1 to N are sequentially obtained by one object. Moreover, the transformation between i and j is known; let denote the rotation angle of relative motion M i j and let ⃗ u denote the rotation axis of M i j . Then, a coordinate system (X, Y, Z ) based on the ith scan is established, and the rigid rotation matrix R( , ⃗ u) ∈ SO(3) is represented as in (9), where u x , u y and u z denote the components of rotation axis in coordinate (X, Y, Z ).
According to Section 2.1, the relative motion M i j of any two views is satisfied, M i j = M i −1 M j , which can be expressed as (10) If we only consider the rotation and let R i j ∈ SO(3) represent the rotation matrix between i and j scan, we have The angle constraint of rigid motion is rewritten as As for the multi-view motion averaging problem, we introduce the angle constraint into the cost function, which is formulated by the LRS matrix. On the basis of (3) and (12), the cost function of the optimisation problem can be denoted as

Definition of weighted matrix
Dealing with outliers and random noise, the weighted adjacency matrix W is introduced into our framework. In different task, the weight updating strategy is different. The proposed approach has two applications. In the task of motion recovery, the weighted strategy based on the Cauchy kernel function can effectively reduce the influence of outliers in motion. Alternatively, in the multi-view point clouds registration, restoring all motions is the preparation for accurate registration. Therefore, the weighted strategy based on the overlap rate of two point sets can effectively improve the robustness of the algorithm.

Weighted matrix of the Cauchy kernel function
In order to recover the motions and improve the robustness, Arrigoni et al. [12] used the Cauchy kernel function to weight relative motion of two views. Compared with the ordinary leastsquares objective function, the Cauchy kernel function can suppress outliers when the residual error is too large.
Suppose that the relations of all views are known; let  denote the undirected weighted graph, which is stored in the symmetric adjacency matrix A = [w i j ], where w i j is non-negative and represents the measurement reliability of relative motion M i j . We use the Cauchy kernel function to update the weight matrix, namely where r i j = ‖R i j − R i −1 R j ‖ F denotes the residual of the rotational component, and c = 1.482med(|r − med(r )|) ⋅ 1 is a parameter in general, where med(⋅) is the median operator, r is the vectorisation of the residuals r i j and 1 = 2 is a tuning constant. Then, the weighted matrix can be expressed according to the adjacency matrix A, as follows: where ⊗ denotes the Kronecker product, and ]] 4×4 denotes a 4 × 4 matrix, whose all elements are 1.

Weighted matrix of overlap rate
In point clouds registration, Chetverikov et al. [25] proposed the TrICP algorithm to solve the registration problem of nonoverlapping region. The algorithm introduced an overlap percentage in the original ICP algorithm and predefined a threshold value thr . If the overlap ratio is higher than the threshold thr , then it is considered that the estimation of the relative motion is reliable. Based on the approach by Jin et al. [17], we estimate the relative motion of two point sets by TrICP and induct the overlap rate to the LRS decomposition framework. Thus, the weighted matrix is demonstrated as follows: , and A i j is defined in Section 3.2.1.

Motion averaging algorithm
In this work, we combine the angle constraint with the weighted LRS matrix framework, and the definition of angle constraint and weight matrix are mentioned in Section 3.2. With the rank relaxation method in RPCA, the following optimisation problem is obtained from (8) and (13): where ⊙ denotes the Hadamard product and W denotes the weight matrix. The initial value of weight matrix W is an N × N matrix consisting of 0 and 1, which is updated according to the update strategy described in Section 3.2 during iteration. From the above conditions X = UV and U T U = I 4 , it is not difficult to get an inference that the nuclear norm of X is equal to the nuclear norm of V . Therefore, we consider the following optimisation problem: Note that the optimisation problem, namely (18), can be solved based on the ALM method. Then, we demonstrate the augmented Lagrange function (U, V, E, Y, ) as where Y denotes the Lagrange multiplier, denotes the penalty term coefficient and R( , ⃗ u) represents the known rotation matrix of angle constraint. The Gauss-Seidel iterative algorithm is effective to solve the optimisation problem [23]. In each iteration, six steps are included as follows and U, V and E are alternately updated as the first to third step: (20) Solving (21) by the SVD method of (X − E k+1 2. Update V as: [see (22)] Generally, the rank minimisation problem is as follows: arg min and the solution of (23) is . The definition of soft threshold-ing operator S (Σ) is Consequently, the solution of (22)

Update E as
4. Update the Lagrange multiplier Y as 5. Update the as: k+1 = min( max , k ), where is the update rate of . 6. Calculate the objective function as where U V (i,i+1) is a 3 × 3 relative rotation matrix between view of i and (i + 1), and it is derived by matrix U k+1 multiply V k+1 . Iteration terminates while three iterative conditions is fulfilled as follows: Iteration stop condition 1: Iteration stop condition 2: Iteration stop condition 3: Then, we obtain the recovery relative motion as M = UV . In order to reduce the complexity of computation, in addition to determining the angle constraints in ALM, it is also essential to determine whether the angle constraint iteration needs to be terminated in the external loop. Therefore, the iteration stop condition of angle constraints can be described as If the iteration stop condition is not satisfied, we assign U * V to X , update the weight matrix W and continue the iteration with the ALM method. While the iteration is stop, the relative motion matrix M is obtained. The first to third rows of matrix M denote the global motion of every scans. Thus, the global motion of the ith view M i is the first to third rows and (4i − 3)th to 4ith columns of M .
In terms of computational complexity, due to the fact that the added constraints of view change are supposed to be a priori information, the iteration processes in (27) and (31) have played an important role in limiting the number of iterations. Compared with the ALM method without angle constraints, the speed of algorithm is improved without affecting the accuracy of results.

EXPERIMENTS
We evaluated and demonstrated the proposed method on both simulated and real data in accordance with accuracy and robustness. Moreover, the proposed approach was compared with two methods. In the simulated experiment, the proposed method was compared with the weighted LRS method, which was solved by L1-ALM [23] without angle constraint, and we called this approach WLRS_without_AC. In the real datasets experiment, the proposed method was compared with the TrICP_WLRS method proposed by Jin et al. [17]. Parameters in experiments are set as follows: = 1.2 and max = 10 20 . All of the experiments were implemented in MATLAB and performed on an Intel i5 core @ 2.20 GHz.

Simulated data
In this simulated experiment, we supposed the camera rotate 360 • around an object and obtained N motions in average. Using the random command such as rand (⋅) and randn(⋅) in MATLAB, we acquired a random set of translation vectors, random rotation noise matrices and translation noise vectors. The noise in rotation matrices follows a Gaussian distribution with a zero mean and a standard deviation R ∈ [0 • , 10 • ], and the noisy translation components are sampled from a Gaussian distribution with a zero mean and a standard deviation T ∈ [0, 0.1 m]. The weight matrix is updated with the Cauchy kernel function. We considered N = 100, and let p denote the ratio of outliers. In this way, we compared the simulation results with the relative motion and calculated the rotation error e R and translation error e T . Figure 1 illustrates the variation trend of rotation error e R and translation error e T with the added random noise R ∈ [0 • , 10 • ]  Figure 1(a), with the increase of R , the rotation error e R is gradually increased, while Figure 1(b) indicates a similar trend of the translation error. Therefore, the proposed method is effectively adapted to the random noise and outliers especially for the rotation error; the rotation error is decreased below three degrees when rotation noise R = 10 • and motion averaging is accurately obtained when 30% of initial motions are outliers.
In order to analyse the effectiveness of angle constraint, here, we only considered the error in rotational state not in translation direction. As shown in Figure 2 and Table 1, in the case of 5% outliers, the difference in rotation error with different methods is very small and under 0.1 • . However, with the increase of outliers, the accuracy of the proposed method is obviously better. In comparison, we found a significant difference in runtime and iteration times. The number of iterations is decreased more than one half, and the computational speed is approximately increased by twice in both p = 0.05 and p = 0.3 conditions. To sum up, the angle constraint in the proposed method can limit the number of iterations and reduce the complexity of the algorithm. Meanwhile, we can get almost the same precision under low noise and significantly higher precision when the ratio of outliers and random noise are large. Note that a known parameter in the angle constraint was introduced; thus, an impact analysis of angle constraint value is also crucial. Consequently, we discussed the variation of rotation error e R as the value of angle parameter was 0.4 , , 1.4 and 2 in the case of different outliers p. According to the results in Figure 3, the trend is similar in different angle constraint values . In contrast, their variation trends and transition nodes are highly correlated, namely, the variation of angle parameter have little effect on the accuracy of global motion recovery. Theoretically, different angle constraint values ∈ (0, 2 ] have the same limitation on the iteration; therefore, an accurate global motion can be obtained by given a large range of angle constraint values. However, at the late stage of Figure 3(a), referring to the critical parts of the blue curve located after seven, the line of = 0.4 was below others. This special case is most attributed to the little relative motion cumulative error in (19) when the random noise and outliers are relatively small. It denotes that when we calculate the objective function j is so small that the iteration stop condition is satisfied immediately. The proposed motion averaging method can significantly accommodate to a large proportion of random noise and outliers and is given a wildly range of angle constraint values. Therefore, the advantage of robustness is more obvious. In addition, it is potentially to be applied in SfM, particularly on large-scale datasets.

Real data
We applied our approach to multi-view point clouds registration and tested on the Stanford 3D Scanning Datasets [26]. The Dragon dataset includes 15 point sets, which scanned around the model. We calculate the objective function value V obj to evaluate the accuracy of registration as follows: where ( i , M i ) implies the distance of the nearest neighbour point searched by the method of K-D tree [17], and more details of V obj can be found in [27]. In practical applications, smaller V obj means more precise 3D point clouds reconstruction. In this experiment, we established the LRS framework with angle constraint and adopted the same registration and update the weighted matrix method in [17], namely, the TrICP algorithm and weighted matrix based on the overlap rate. In this part, a set of random noise was added in each initial relative motion. As displayed in Figures 4 and 5, the value of V obj in the proposed method is almost the same as the method   proposed by Jin et al. [17], but the iteration times and runtime in different methods exist a specific link. And the trend of iteration and runtime (seconds) is similar with the increase of rotation noise, but our approach is significantly more effective than the TrICP_WLRS method. The random noise we added was indeterminate; thus, we repeated experiments to verify the robustness of the proposed method. As for the TrICP_WLRS method, the experiment was repeated 100 times for each value of R ∈ [2 • , 16 • ] and took the average of V obj . Alternatively, in the proposed method, we generated 10 groups of random noise, which included 14 noise rotation matrixes in each group, then repeated 10 times for each value of R and took the average of V obj . The comparison results are shown in Table 2. As can be seen from Table 2, the proposed method can reduce the sensitivity to noise and improve the success rate of registration while ensuring the registration result V obj . Especially, when rotation noise R is above 10 • , the accuracy rate and registration success rate in the proposed method are better than those in the TrICP_WLRS method.

CONCLUSION
In order to solve the over-reliance on the initial value and high computational complexity problems of the weighted LRS matrix decomposition motion averaging methods, which utilize the Cauchy function or overlap rate to construct the weighted matrix in the cost function, we proposed a novel motion averaging method and thus effectively catered for noise and lessen the execution time. This paper first derived the angle constraint theoretically where at least one exact rotation is known. In the optimisation process, we added the angle constraint into the cost function as an additional term. Finally, we derived a fast and robust iterative solution for the cost function based on the ALM method. Compared with WLRS_without_AC and TrICP_WLRS algorithms, the experiments on both simulated motions and real 3D scanning datasets demonstrate that the proposed method can improve the efficiency and robustness of motion averaging. The proposed method can not only be applied to the 3D reconstruction methods by point clouds registration and SfM, but also provide theoretical basis and efficient solutions for potential engineering applications such as simultaneous localisation and mapping, satellite navigation and tracking and pose estimation.