Total Least-Squares Iterative Closest Point Algorithm Based on Lie Algebra

: In geodetic surveying, input data from two coordinates are needed to compute rigid transformations. A common solution is a least-squares algorithm based on a Gauss–Markov model, called iterative closest point (ICP). However, the error in the ICP algorithm only exists in target coordinates, and the algorithm does not consider the source model’s error. A total least-squares (TLS) algorithm based on an errors-in-variables (EIV) model is proposed to solve this problem. Previous total least-squares ICP algorithms used a Euler angle parameterization method, which is easily a ﬀ ected by a gimbal lock problem. Lie algebra is more suitable than the Euler angle for interpolation during an iterative optimization process. In this paper, Lie algebra is used to parameterize the rotation matrix, and we re-derive the TLS algorithm based on a GHM (Gauss–Helmert model) using Lie algebra. We present two TLS-ICP models based on Lie algebra. Our method is more robust than previous TLS algorithms, and it suits all kinds of transformation matrices.


Introduction
An iterative closest point (ICP) algorithm [1] estimates a pose transformation matrix using coordinates from a source model and target model. This technique is frequently used in ego-motion estimation [2], geodesy, computer vision, and system identification [3]. In this paper, the ICP algorithm is only concerned with the estimation of the pose transformation matrix. This means we have known correspondences between two point clouds, and there are no outliers in the source model and target model. ICP forms a cost function using two point clouds' coordinates and searches for a rigid transformation to minimize its residual error. In [4], Horn proposed an analytical solution to the ICP problem using unit-quaternion. This approach only considers target errors and assumes the error is isotropic, identical, and independent Gaussian noise. If errors from the target model are not isotropic and identical, then we can use a nonlinear optimization method [5] to obtain a rigid transformation matrix. However, errors do not only exist in the target model in practice. A more accurate model assumes that errors exist in both the source model and target model; this is called an errors-in-variables (EIV) model. Golub and Van Loan [6] proposed a total least-squares (TLS) method, especially for solving EIV problems. If the cost function of EIV is linear in the estimated parameter, we can use a Newton-Gauss TLS algorithm [7] or Lagrange TLS algorithm [8] to solve an EIV model. However, the cost function of an ICP algorithm is nonlinear because of the rotation matrix. New TLS algorithms were designed to solve nonlinear ICP problems [9,10]. The Felus [9] uses a Procrustes approach to obtain a rotation matrix to avoid nonlinearity in the TLS cost function. Then the scale s and translation t are obtained using a linear EIV model. Felus assumes that errors in the target model and source model are the same, so a rotation matrix from the Procrustes approach satisfies the EIV model. In [10], Brun fixes the translation and estimates the rotation matrix using [10]. After obtaining the rotation matrix, the translation matrix is estimated using a linear EIV model. [10] cannot guarantee the optimality and effectiveness of the algorithm because of its alternating optimization strategy. [9,10] use approximation methods and Euler angle parameterization to solve the nonlinear TLS problem. In [11], Ohta constructs a cost function using elements of the rotation matrix directly. This method is convenient for obtaining results using existing methods such as [7,12]. However, it minimizes the cost function in Euclidean space, which does not conform to reality because the rotation matrix belongs to a special orthogonal group. We cannot promise that rotation matrix constructed by elements is identity orthogonal matrix. We must minimize the Frobenius norm to project the rotation matrix of the TLS algorithm to a special orthogonal group. This does not guarantee that the final rotation matrix is the global optimal solution. Previous papers focus on how to construct a robust cost function and how to deal with the weight matrix. The Euler angle is a commonly used variable to parameterize the rotation matrix [13][14][15]. However, the Euler angle loses its freedom in the case of a gimbal lock and is not suitable for interpolating in nonlinear optimization. In practical situations, measurement data cannot promise that a gimbal lock does not occur. To solve this problem, we introduce Lie algebra to the TLS problem. We use a variant of the Gauss-Helmert [16] model, which is a general method to solve nonlinear TLS problems without assumptions, and Lie algebra is used to parameterize the rotation matrix without the gimbal lock problem.
The rest of this paper is organized as follows. In Section 2, we introduce definitions and properties of Lie algebra and Lie groups. Then we introduce a variant of the Gauss-Helmert model that is used in the ICP model. We specify the Jacobian matrix used in the Gauss-Helmert model. In Section 3, we use a numerical example to verify our TLS algorithm, and we compare it to least-squares results. We summarize our conclusions in Section 4.

ICP Algorithm Based on Lie Algebra
In Section 2.1, we describe the fundamental ICP problem and introduce the TLS model to solve it. In Section 2.2, the GHM model is used to solve the nonlinear TLS problem. In Section 2.3, we introduce Lie groups and Lie algebra to parameterize the rotation matrix. In Section 2.4, we derive the Jacobian matrix for use in the GHM model and parameter updating method.

Iterative Closest Points
We use ICP to compute a rigid transformation matrix between two frames. We measure points' coordinates from two frames, and their correspondences are known. Measurements are assumed without data, and we can obtain the equation: where X i is a 3D point from the target model without noise, and X i is a 3D point of the source model without noise. X i = X i + ∆X i , where X i is the measurement of the 3D point from the target model, and ∆X i is the measurement error in the target model. X i = X i + ∆X i , where X i is the measurement of the 3D point from the source model, and ∆X i is the measurement error in the source model. The aim of ICP is to estimate the rotation matrix R and t using measurements {X i } and {X i }. The least-squares (LS) method thinks that ∆X i = 0, so the optimization equation is: Appl. Sci. 2019, 9, 5352 3 of 12 where n is the total number of measurements from the target model and source model. We can see from Equation (2) that the LS approach only minimizes residual error in the target model. TLS considers the error from both the target model and source model, and the optimization equation is: We do not only minimize Equation (3); we must guarantee that X i + ∆X i = R(X i + ∆X i ) + t. We add a Lagrange multiplier to Equation (3) and the optimization equation becomes: where λ i is a real number. Equation (4) is a constrained nonlinear optimization problem; we use the GHM model to solve this problem.

Variant of the Gauss-Helmert Model
First, we transform Equation (4) to: where e = vec(V A ) vec(V y ) 6n * 1 , and the operator "vec" is to stack one column of a matrix under the other from left to right. p is the Fisher information matrix and is equal to the identity matrix in this paper: x 1 + ∆x 1 y 1 + ∆y 1 z 1 + ∆z 1 x 2 + ∆x 2 y 2 + ∆y 2 z 2 + ∆z 2 . . . . . . . . . x n + ∆x n y n + ∆y n z n + ∆z n x 1 + ∆x 1 y 1 + ∆y 1 z 1 + ∆z 1 x 2 + ∆x 2 y 2 + ∆y 2 z 2 + ∆z 2 . . . . . . . . . x n + ∆x n y n + ∆y n z n + ∆z n l n is an n × 1 vector of ones. The operator ⊗ denotes the Kronecker-Zehfuss product, and related to the estimated parameter ξ and random error e. We use a first-order Taylor expansion to Ψ and obtain a linear equation with respect to ξ and e: where J e = ∂Ψ ∂e T , J ξ = ∂Ψ ∂ξ T , δe = e − e 0 , and δξ is the update step of the estimated parameter. We substitute Equation (8) into Equation (5) to obtain: We take the partial derivative of Equation (9) with respect to e and the estimated parameter ξ to obtain: Appl. Sci. 2019, 9, 5352 4 of 12 where B 0 = [R ⊗ I n −I 3 ⊗ I n ] and A 0 = J ξ . The derivations of these two parameters are given in Section 2.4. We combine Equations (8)- (11) and obtain iterative results: where ω 0 = −J e e + Ψ(e 0 , ξ 0 ). The initial value of the random error is e 0 = 0, and ξ 0 = 0 0 0 0 0 0 T . We can use Equations (12)- (14) to update ξ.

Lie Algebra and Lie Group
Euler angle constructs rotation matrix as following equation: where θ z , θ y , and θ x are, respectively, the rotation angles around the z-, y-, and x-axes. If θ y = π 2 , then the Euler angle is unsuitable for interpolation because of the gimbal lock problem. When the second rotation angle is 90 degrees, the rotation matrix constructed by the Euler angle loses z-and x-axis freedom. In this paper, we introduce Lie algebra and a Lie group to replace the Euler angle.
A group is a set of elements together with an operation that combines any two of its elements to form a third element also in the set. A Lie group is also a differential manifold, with the property that the group operations are smooth [17]. A Lie group must satisfy the four properties listed in Table 1. Table 1. Properties for Lie group.

SO (3) SE (3)
Closure A rotation matrix is a special orthogonal group, denoted by SO (3), which also belongs to a Lie group. A rigid transformation is in a special Euclidean group, denoted by SE (3), which also belongs to a Lie group. Every Lie group is associated with a Lie algebra. For example, a rotation matrix is associated with a three-dimensional vector φ, which is the angle-axis representation of the rotation matrix R. The rigid transformation matrix R t 0 1 is associated with the six-dimensional vector ξ. SO where φ ∧ is the antisymmetric matrix of vector φ. The Lie bracket of ξ is: Appl. Sci. 2019, 9, 5352 5 of 12 The operator "∧" has different definitions for the Lie algebra of SO (3) and We can see from Equation (8) that a nonlinearity factor only exists in a rotation matrix. According to SO (3) and SE (3), we have two strategies to parameterize a rigid transformation matrix. The first is to use SO (3) and a translation vector to parameterize the transformation matrix. The second uses SE (3) directly to parameterize the transformation matrix. We will compare the performance of these two methods and come to a conclusion which method is more suitable for the TLS-ICP algorithm. Here, we provide two parameterization methods to calculate J ξ .

Method 2:
The estimated parameter is ξ = [φ, ρ]: where θ is the norm of φ, and a is the unit vector of φ. After obtaining δξ, we should update the estimated parameter ξ. If the parameter is in Euclidean space, such as a Euler angle, then the parameter can be updated using ξ k+1 = ξ k + δξ. Here, we notice that Lie algebra does not belong to Euclidean space, and the operators "+" and "-" are undefined. We define a new plus operator, "⊕," which conforms to the Lie algebra update rule. If we use Lie algebra to parameterize a rotation matrix or transformation matrix, the parameter is updated as ξ k+1 = ξ k ⊕ δξ. Section 2.4 provides details of the update formula.
where e is linear in Equation (5), and J e is easily obtained as: where I n is an identity matrix of dimension n. Then we need to calculate J ξ = ∂ψ ∂ξ T . Method 1: If we use parameterization method 1, then the Jacobian matrix is: ∂ψ ∂ t 1 t 2 t 3 is easily obtained from Equation (20) as: ∂ψ ∂ φ 1 φ 2 φ 3 can be simplified as follows: According to deductions from Appendix A, we know that: We let R(X i + ∆X i ) = f i g i h i T , so we can obtain the following equation: The elements in Equation (26) can be obtained from Equation (25). Thus, ∂ψ ∂ φ 1 φ 2 φ 3 is equal to: The estimated parameters in method 1 are the Lie algebra of SO (3) and the translation vector.
Thus, the update parameter is equal to [ ∆φ ∆t ] T . The translation vector is in Euclidean space, so it is easy to update the new translation variable: t k+1 = t k + ∆t. ∆φ belongs to the Lie algebra of SO (3), and the updated formula is: R k is the kth step rotation matrix. The operator "∧" is defined in Section 2.4. The "ln() ∨ " operator transforms the Lie group of SO (3) to the Lie algebra of SO (3). We do not provide specific expressions for "ln() ∨ " here; these can be found in [17]. In Appendix A, we deduce the Jacobian of Equation (25) using left multiplication so that update Formula (30) also uses left multiplication. Thus, the method 1 parameter update formula is: Method 2: If we use parameterization method 2, the Jacobian matrix is: From Appendix B, we know that: The left side of Equation (33) uses homogeneous coordinates, and the right side uses inhomogeneous coordinates. We transform homogeneous coordinate to inhomogeneous coordinates by default. We let T(X i + ∆X i ) = f i g i h i T , so we can obtain the following equation: The elements in Equation (34) can be obtained by Equation (35). Thus, we have: The estimated parameters in method 2 are the Lie algebra of SE (3). The update parameter is equal to ∆ξ, which belongs to the Lie algebra of SE (3), and the updated formula is: where T k is the kth step rigid transformation matrix. The Operator "∧" is defined in Section 2.4. The operator "ln() ∨ " transforms the Lie group of SE (3) to the Lie algebra of SE (3). We provide no specific expressions for "ln() ∨ ." These can be found in [17]. In Appendix B, we deduce the Jacobian of Equation (33) using left multiplication, so that update Equation (36) also uses left multiplication. Thus, the method 2 parameter updated formula is:

Experiments
In this section, we use the dataset from [9]. This numerical example is from an actual experiment of fitting two surfaces that are surveyed in two different reference systems. Four control points are chosen in two different coordinate systems. We assume that we have a known scale equal to 2.1, and we multiply the source coordinate by 2.1. Data from [9] are shown in Table 2. We compare the LSand TLS-ICP algorithms using simulated data and demonstrate that the TLS-ICP algorithm is more accurate. Then we use real measurement data from Table 2 and demonstrate that the TLS algorithm has a smaller residual error than the LS method.  1  63  84  21  290  150  15  2  210  84  21  420  80  2  3  210  273  21  540  200  20  4  63  273  21  390  300  5 First, we compare results from TLS and LS using simulated data. A true translation vector is set to (190, 110, −15), and a rotation matrix is constructed by the Euler angle. The true rotation angles around the z-, y-, and x-axes are set to 45 degrees, 90 degrees, and 60 degrees, respectively. Then, we generate controlled points in target coordinates according to the true rotation matrix, true translation, and controlled points in the source coordinate. We add noise in the source and target coordinates according to two distributions: normal and uniform. We classify noise into eight classes, as shown in Table 3. The normal distribution is a Gaussian distribution with zero mean and a covariance. We use the uniform distribution to generate samples in some range, and every sample is generated by equal chance. In the first row of Table 3, 0.1 is the variance of the Gaussian distribution, and a noise number of 5 means that the uniform distribution generates samples ranging from −5 to 5. The rotation error is: Error R = acos(max(min(0.5 * (∆R(0, 0) + ∆R(1, 1) + ∆R(2, 2)), 1), −1), where ∆R = R T true R TLS , R TLS is the rotation matrix from the TLS algorithm, and R true is the true rotation matrix. The error measure in Equation (38) describes how close ∆R is to the identity matrix. If ∆R is closer to the identity matrix, then our algorithm is more accurate. The translation error is evaluated as Error t = t true − t TLS , where t TLS is the translation from the TLS algorithm, and t true is the true translation.
If noise is generated by a normal distribution, then the information matrix p is set to the inverse of the variance. If noise is generated by a uniform distribution, then the information matrix p is set to the identity matrix. Figures 1 and 2 show the rotation error and translation error, respectively, of least-squares, method 1, and method 2. We see from these figures that larger variances between the x-, y-, and z-axes lead to a larger error of the least-squares method. This is because the rotation matrix is obtained by the Procrustes approach [9], which assumes that noise in the x-, y-, and z-axes is equal. In this situation, the LS algorithm is nearly as accurate as the TLS algorithm, and we can see this result in the first column of Figures 1 and 2. However, the TLS algorithm performs much better than the LS algorithm when variances between the x-, y-, and z-axes become large, as shown in the fourth column and the eighth column in Figures 1 and 2. Method 1 is more accurate than method 2 in estimating the rotation matrix, but method 2 is more accurate than method 1 in estimating the translation vector. This phenomenon is caused by different parameterization approaches in these two methods. The parameterization of method 2 uses the Lie algebra of SE (3), which is a six-dimensional vector. The first three elements in this vector describe the translation and rotation matrices. The last three elements only describe the rotation matrix. The first three elements are optimized during the TLS iterative process, and then we recover the translation vector using the last three elements in the Lie algebra vector. The benefit of this is that the Jacobian matrix conditional number is not so large, which means that the estimated parameters equally decide the direction of optimization. In Method 1, the Jacobian matrix of φ is −(Rp) ∧ , and the Jacobian matrix of t is always the identity matrix. If p is very large, then the Jacobian matrix of φ mainly decides the direction of optimization, so that method 1 has a smaller rotation error than method 2.
In this situation, the LS algorithm is nearly as accurate as the TLS algorithm, and we can see this result in the first column of Figures 1 and 2. However, the TLS algorithm performs much better than the LS algorithm when variances between the x-, y-, and z-axes become large, as shown in the fourth column and the eighth column in Figures 1 and 2. Method 1 is more accurate than method 2 in estimating the rotation matrix, but method 2 is more accurate than method 1 in estimating the translation vector. This phenomenon is caused by different parameterization approaches in these two methods. The parameterization of method 2 uses the Lie algebra of SE (3), which is a six-dimensional vector. The first three elements in this vector describe the translation and rotation matrices. The last three elements only describe the rotation matrix. The first three elements are optimized during the TLS iterative process, and then we recover the translation vector using the last three elements in the Lie algebra vector. The benefit of this is that the Jacobian matrix conditional number is not so large, which means that the estimated parameters equally decide the direction of optimization. In Method 1, the Jacobian matrix of is −( ) ∧ , and the Jacobian matrix of t is always the identity matrix. If p is very large, then the Jacobian matrix of mainly decides the direction of optimization, so that method 1 has a smaller rotation error than method 2.  Table 3.  Table 3.
result in the first column of Figures 1 and 2. However, the TLS algorithm performs much better than the LS algorithm when variances between the x-, y-, and z-axes become large, as shown in the fourth column and the eighth column in Figures 1 and 2. Method 1 is more accurate than method 2 in estimating the rotation matrix, but method 2 is more accurate than method 1 in estimating the translation vector. This phenomenon is caused by different parameterization approaches in these two methods. The parameterization of method 2 uses the Lie algebra of SE (3), which is a six-dimensional vector. The first three elements in this vector describe the translation and rotation matrices. The last three elements only describe the rotation matrix. The first three elements are optimized during the TLS iterative process, and then we recover the translation vector using the last three elements in the Lie algebra vector. The benefit of this is that the Jacobian matrix conditional number is not so large, which means that the estimated parameters equally decide the direction of optimization. In Method 1, the Jacobian matrix of is −( ) ∧ , and the Jacobian matrix of t is always the identity matrix. If p is very large, then the Jacobian matrix of mainly decides the direction of optimization, so that method 1 has a smaller rotation error than method 2. Figure 1. The rotation error is compared between least-squares and our proposed methods. The blue, green, and yellow bars are the rotation errors of least-squares, method 1, and method 2, respectively. The x-axis of this figure is a different noise class shown in Table 3.  Figure 2. The translation error is compared between least-squares method and our proposed methods. The blue, green, and yellow bars are the translation errors of least-squares, method 1, and method 2, respectively. The x-axis of this figure is a different noise class shown in Table 3.
Here, we directly use data from [12], which is an actual experiment of fitting two surfaces surveyed in two reference systems. Data are shown in Table 2, and the result is shown in Table 4. We give results of LS, TLS based on Euler angle, method 1, and method 2. Our two methods have much smaller sums of residual errors than the LS method. This is mainly because the TLS algorithm considers noise from the source and target coordinate systems. Further discussion of errors and the model are beyond the scope of this paper; refer to Grafarend [18]. Compared to the Euler angle method, our method has smaller SSE after optimization, which proves that Lie algebra is more suitable for the TLS-ICP algorithm. We should note that our method also can be used in similar transformation estimation. The difference is that the Jacobian of scale s is needed for additional calculations. TLS-ICP algorithm often uses 4-8 points to calculate transformation between two coordinates in geological mapping. Thus, in previous experiments, datasets from [9] are used to evaluate accuracy. However, in the field of three-dimensional reconstruction, ICP algorithm also can be used to calculate robot pose. Thus, the cost time of ICP is very important. In our experiment, all of codes run on a laptop with 2.7 GHz quad cores in Ubuntu. We use 100 points, 500 points, and 1000 points, which are copied from [9], to evaluate the cost time of every step in our algorithm. The cost time of every iterative step is 0.2 ms using four points. The cost time of every iterative step is 380 ms using 100 points. The cost time of every iterative step is 46 s using 500 points. The cost time of every iterative step is 369 s using 1000 points. In our experiment, the process of our algorithm converges to final result after three or four steps.

Conclusions
In this paper, we present two total least-squares ICP algorithms based on Lie algebra using a GHM framework. These methods can estimate a rigid transformation matrix from adjusted source coordinates and adjusted target coordinates. Our methods perform better than the common least-squares ICP algorithm because they consider random noise from the source coordinate, which results in a smaller sum of squared errors. The most important contribution of this paper is that our methods are more robust than [12,19], which utilize the Euler angle to describe the rotation matrix. The nonlinearity in the cost function due to the rotation matrix results in there being no analytical solution to this problem. Thus, an iterative optimization algorithm is employed to solve the GHM model. The Euler angle is unsuitable for interpolation because of the gimbal lock problem. We replace the Euler angle with Lie algebra to execute the TLS optimization algorithm and give deduction of our method in detail.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A
Target: We want to calculate Jacobian matrix of Rp with respect to Lie algebra of SO (3). Here we use left multiplication to calculate Jacobian matrix: We use first order Taylor series expansion to exp(δφ ∧ ) and we obtain:

Appendix B
Target: We want to calculate Jacobian matrix of Tp with respect to Lie algebra of SE (3). Here we use left multiplication to calculate Jacobian matrix: We use first order Taylor series expansion to exp(δξ ∧ ) and we obtain: