WTLS iterative algorithm of 3D similarity coordinate transformation based on Gibbs vectors

The 3D similarity coordinate transformation is fundamental and frequently encountered in many areas of work such as geodesy, engineering surveying, LIDAR, terrestrial laser scanning, photogrammetry, machine vision, etc. The algorithms of 3D similarity transformation are divided into two categories. One is a closed-form algorithm that is straightforward and fast. However, it cannot provide the accuracy information for the transformation parameters. The other category of algorithm is iterative, and this can offer the accuracy information for the transformation parameters. However, the latter usually needs a good initial value of the unknown. Considering the accuracy information for transformation parameters is essential or indispensable from the viewpoint of uncertainty, this contribution proposes a weighted total least squares (WTLS) iterative algorithm of the 3D similarity coordinate transformation based on Gibbs vectors. It is fast in terms of fewer iterations, reliable and does not need good initial values of transformation parameters. Two cases including the registration of LIDAR points with big rotation angles and a geodetic datum transformation with small rotation angles are demonstrated to validate the new algorithm.


Introduction
The 3D similarity coordinate transformation aims to align the coordinate of points in different coordinate systems into a common coordinate system. It usually implements the transformation based on 3D seven-parameter similarity transformation utilizing a set of control points with the coordinates both in the source coordinate system and in the target coordinate system. This work is very popular in many fields, such as geodesy, engineering surveying, LIDAR, terrestrial laser scanning, photogrammetry, machine vision, etc. (Besl and McKay 1992;Crosilla and Beinat 2002;Horn 1987;Jaw and Chuang 2008;Kashani 2006;Krarup 1985;Marx 2017;Paffenholz and Bae 2012;Walker et al. 1991;Wang et al. 2014;Závoti and Kalmár 2016;Zeng 2014;Zeng et al. 2018).
As far as the algorithm of 3D similarity coordinate transformation is concerned, much related literature has been published. The algorithms proposed in that literature can be divided into two categories. One is the closedform or analytical algorithm. This category of algorithms utilizes the technologies such as singular value decomposition (SVD) also known as the Procrustes algorithm, e.g., Umeyama (1991), Crosilla and Beinat (2002), Grafarend and Awange (2003), eigenvalue-eigenvector decomposition (a quaternion-based algorithm), e.g., Horn (1987), Shen et al. (2006), orthonormal matrix, e.g., Horn et al. (1988), Zeng (2015). Additionally, Leick and Van Gelder (1975), Zeng and Yi (2010) presented, respectively, a stepwise algorithm based on the physical meaning of the similarity transformation. Wang et al. (2014) proposed a closed-form pairwise registration algorithm of point clouds utilizing the dual quaternion. The closed-form Open Access *Correspondence: zenghuaien_2003@163.com 1 Hubei Key Laboratory of Construction and Management in Hydropower Engineering, China Three Gorges University, Yichang 443002, China Full list of author information is available at the end of the article algorithms are fast and efficient because the parameters are recovered directly by exact formulae without solution of nonlinear equations. However, this type of algorithm cannot deal with the general weight matrix of observations (i.e., the coordinates as pseudo-observations) as well as the accuracy estimation of transformation parameters. The other one is the iterative algorithm, which utilizes the numerical (iterative) computation technique to find the unknowns. Lots of iterative algorithms have been presented at present, e.g., Zeng and Tao (2003), Chen et al. (2004), El-Habiby et al. (2009), Zeng and Yi (2011), Zeng et al. (2016), Kurt (2018), Zeng et al. (2019). This type of algorithm usually requires a good initial value of unknown and iterative computation (e.g., Zeng and Tao 2003;Zeng and Yi 2011). But, in some situations, e.g., registration of LIDAR points due to the arbitrary size of rotation angles, it is difficult or even unlikely to obtain a good initial value of parameter. As a result, the algorithm needs many iterations or falls into a local minimum or diverges. Sometimes initial values are not required if global optimization algorithms are performed (e.g., Xu 2003); unfortunately it may lead to much more computation time and burden. The advantage of iterative algorithms over closed-form algorithms is the former can supply the accuracy information of unknowns, which is essential from the viewpoint of uncertainty.
In recent years, the errors-in-variables (EIV) model has evoked a lot of research interest. It considers the errors in all variables; for instance, in the 3D similarity transformation, the errors of target coordinates and source coordinates of control points are considered. Thus it is more reasonable than the traditional Gauss-Markov (GM) model which just considers the errors of target coordinates of control points. Golub and Van Loan (1980) presented the approach dealing with the EIV model and named it total least squares (TLS) technique. Afterward, a few types of TLS algorithms of 2D or 3D similarity coordinate transformation have been proposed. Teunissen (1988) derived a closed-form solution for 2D similarity transformation in the EIV model, i.e., symmetric 2D similarity transformation. Goryn and Hein (1995) presented a TLS solution for the 3D rigid transformation based on the Procrustes algorithm considering homogenous and uncorrelated errors. Schaffrin and Felus (2008) proposed TLS approaches to empirical coordinate transformation, which improved the accuracy in 2D affine transformations. Felus and Burtch (2009) proposed a closed-form weighted TLS (WTLS) solution to the 3D similarity transformation with pointwise weights and uncorrelated errors among points based on the Procrustes algorithm. Neitzel (2010) proposed a TLS solution to 2D similarity transformation within the nonlinear Gauss-Helmert model. Fang (2015) presented a WTLS algorithm with constraints and universal formula for geodetic transformation. Chang (2015) presented a rigorous solution without presupposed fixing of the scale parameter (which some authors assume to be 1). Mahboub (2016) proposed a WTLS solution to 3D symmetrical similarity transformation without linearization; however, more iterations may be needed than the linearized model method. Mercan et al. (2018) proposed a weighted similarity transformation based on quaternions. The algorithm has seven unknowns including the translation parameters and scaled quaternion, and the iterations may lead to divergence.
As mentioned above, the closed-form solution to 3D similarity coordinate transformation is straightforward and fast, but cannot provide the accuracy of recovered transformation parameters. Iterative solutions to 3D similarity coordinate transformation can deal with general weight and provide the accuracy of recovered transformation parameters; however, it usually needs a good initial value of unknowns. Therefore, this contribution intends to present a new WTLS iterative solution to 3D similarity coordinate transformation based on Gibbs vectors, which is fast and does not need a relatively good initial value of the unknowns; in other words, is not sensitive to the initial value of unknowns.
In the next section, firstly 3D similarity transformation is introduced and the model of similarity transformation based on Gibbs vectors is established. Secondly in order to simplify the transformation model for improvement of computation performance, the translation parameters are derived in the WTLS sense and a WTLS iterative algorithm of 3D similarity transformation based on Gibbs vectors is proposed with detailed derivation. Two cases including the registration of LIDAR points and geodetic datum transformation are studied; the results show the new algorithms are fast and reliable. Finally, conclusions are drawn in the last section.

Formulation and WTLS iterative algorithm of 3D similarity transformation
In this section, the basic EIV model of 3D similarity transformation is introduced, and in order to avoid transcendental function and functional constraints, the model of similarity transformation based on Gibbs vectors is established. Lastly, a WTLS iterative algorithm of 3D similarity transformation based on Gibbs vectors is put forward.

3D similarity transformation in EIV model without functional constraints
Assume that a set of control points are given with their coordinates in both source and target coordinate systems.
Then the 3D similarity transformation in EIV model can be expressed as subject to T are the 3D coordinate vectors of control point i ( i = 1, 2, . . . , n ) in the source coordinate system and the target coordinate system (denoted with superscript o and t ), respectively. e o i and e t i are the error vectors of p o i and p t i , respectively. And superscript T is transpose of matrix, det represents determinant computation of matrix, I 3 is a 3 × 3 identity matrix. denotes the scale factor, t = t x t y t z T is the vector of three translation parameters and R represents the rotation matrix of size 3 × 3, which is traditionally expressed by three rotation angles. Assume R is produced by counterclockwise coordinate-frame rotations θ x , θ y , θ z ( θ x applied first, θ z applied last). Then R can be represented by rotation angles as If R is known, and the 3 main diagonal elements are positive, θ x , θ y , θ z can be calculated quickly by (3) as where R ij is the element of R in the row i and column j.
For n control points, the 3D similarity transformation in EIV model is constructed easily from Eq. (1). where . . e tT n T , and . I n is the identity matrix of size n × n , the symbol ⊗ means the Kronecker product. The objective of 3D similarity transformation in EIV model is to recover the seven transformation parameters in the principle of total least squares, i.e., where W of size 3n × 3n is a weight matrix of observations, adopting the point (i.e., row or column) weight in (4) e tT We t + e oT We o = min, Felus and Burtch (2009) as follows, since the point weight is more reasonable than the identity matrix weight usually employed in practice.
where w i is the weight of point i ( i = 1, 2, . . . , n ). Note that the components of R in Eq. (3) are trigonometric functions of rotation angles, which causes a high computation burden in the solution of parameters. In some particular situations, for instance in geodesy, the small rotation angles are small usually at the level of seconds; the elements of rotation matrix can be reduced to rotation angles (in radians) or constants (one or zero) and the 3D similarity transformation model is simplified to a linear one. Evidently this reduced treatment causes model error and then is not applicable to big rotation angle cases. In this contribution, the nonlinear 3D similarity transformation model is adopted to maintain validity for any rotation angles.
Apart from using rotations, there are other representation of R , such as unit quaternion (e.g., Horn 1987;Shen et al. 2006;Zeng and Yi 2011;Závoti and Kalmár 2016), dual quaternion (e.g., Walker et al. 1991;Wang et al. 2014;Zeng et al. 2018Zeng et al. , 2019, direction cosine matrix (e.g., Chen et al. 2004;Wang et al. 2018), and the Rodrigues matrix or Gibbs vector (e.g., Zeng and Yi 2010;Závoti and Kalmár 2016;Zeng et al. 2016;Kurt 2018). The first three approaches to the representation of R introduce functional constraints, for instance the norm of quaternion is unity as the unit quaternion representation is concerned. Additional functional constraints increase the complexity of formulation and computation. However, the 4th approach is free of transcendental function and functional constraints, and then improves the speed of computation. For this sake, this paper adopts the Rodrigues matrix or Gibbs vector to represent the rotation matrix R . Due to the orthogonality of rotation matrix, it can be represented by a Rodrigues matrix as It is easy to see that the rotation matrix R is expressed in terms of polynomials without trigonometric functions and functional constraints.

Deduction of translation parameters and simplified 3D similarity transformation model
The similarity transformation parameters in preceding EIV model based on the Rodrigues matrix or Gibbs vector theoretically can be recovered by the TLS derivation; however, it is found in computation practice that the normal equation is singular because the numerical magnitude of translation parameter t and those of the Gibbs vector v as well as scale factor probably differ greatly from each other, which prevents a solution to the equation. For this reason, we seek to eliminate the translation parameters t from the original model.
Lagrangian extremum principle is a reliable way to solve the 3D similarity transformation problem in the sense of TLS. In order to obtain the explicit solution of t regarding and R as known parameters, the Lagrangian extremum problem which just considers t as unknown is established as where k is a vector of the Lagrangian multiplier. In order to obtain the solution, partial derivations with respect to all variables should be set to zero. From the partial derivations with respect to e o , one gets From the partial derivations with respect to e t , one gets Inserting Eqs. (12) and (13) into the partial derivations with respect to k which is identical to Eq. (5) and arranging the terms, one obtains substituting B = I n ⊗ R into Eq. (15), the following formula can be easily derived: (14) with k T C = 0 from the partial derivative of t, one obtains the solution of t as Inserting Eq. (16) into Eq. (17), one obtains the following formula and its proof is given in Additional file 1: Appendix S1: where and are the weighted coordinates of the barycenter in the target coordinate system and the source coordinate system, respectively.

thus Eq. (21) is rewritten as
It is easily proved that and the proof is given in Additional file 2: Appendix S2. Then Eq. (23) is re-expressed as where and are centralized coordinate vectors, i.e., the difference of original coordinates and the weighted coordinates of the barycenter, in the target coordinate system and source coordinate system, respectively. D is a centering matrix.
It is obvious to see that the translation parameter t is eliminated.

WTLS iterative algorithm of 3D similarity transformation in errors-in-variables model
Linearizing Eq. (25) by Taylor's formula, one obtains where superscript j denotes the approximation from iteration j of the corresponding variable, dx = d da db dc T is the correction vector of unknowns x = a b c T , and Next, we can utilize an iterative algorithm to compute the final solution of x . Set the initial value x 0 = 1 0 0 0 T , the initial value of the unknown e o as 0 0 · · · 0 T 3n and compute the dx , if the absolute values of all elements of dx are less than a small tolerance value τ (e.g., 1.0 × 10 −10 ), stop the iteration, otherwise update the x with previous iterative value of x plus dx and repeat the iterative computation until the stop criterion is satisfied. After iteration, the variance factor of unit weight is estimated as or if A j dx is a sufficiently small term, the approximate formula as follows is adopted: By Eq. (40) the estimated covariance matrix of transformation parameters x is derived as Similarly by Eq. (17) the estimated covariance matrix of transformation parameters t is derived as To sum up, the above presented WTLS iterative algorithm is listed in Table 1.

Numerical analysis and discussion
Two cases are studied. One is the registration of LIDAR points involving big rotation angles and the weight matrix is not considered, i.e., an identity matrix is adopted. The other one is 3D coordinate transformation in geodesy, which encounters small rotation angles and point-weighting is adopted.

Case 1 (registration of LIDAR points)
This data are adopted from Wang et al. (2014). The point features extracted from two neighboring LIDAR point clouds are listed in Table 2. In order to verify the new algorithm, the 18 points are divided into control points and check points. The first 10 points in Table 2 are arbitrarily chosen as control points and the last 8 points are check points. The distribution of all points is depicted in Fig. 1. In Fig. 1 plus signs denote control points and cross signs denote check points. The weight matrix adopted is an identity matrix in this case; in other words, the weights of all points are the same, and the correlation among the three coordinates of each point is not considered.
The transformation parameters in the WTLS sense are computed by the new algorithm and Algorithm 4 presented in Felus and Burtch (2009), henceforth described as the Felus-Burtch Algorithm. All results are listed in Table 3. It is seen from Table 3, that the new algorithm needs only 6 iterations to converge. And for two algorithms, the results are exactly the same. Thus the new algorithm is fast and reliable. Additionally, the Felus-Burtch Algorithm does not provide the accuracies of recovered transformation parameters; however, the new algorithm does. It gives the accuracy information which is shown in Table 3. For more details, refer to the covariance matrix of recovered transformation parameters offered by the new algorithm, which is listed in Table 4. The predicted errors of coordinates of control points by the new algorithm are listed in Table 5. It is worth mentioning that the formula to compute E Y in Algorithm 4 presented from Felus and Burtch (2009) is wrong. The first term namely the inverse of weight matrix should be removed from the formula. If the LS transformation parameters and predicted error of coordinates of control Table 1

WTLS iterative algorithm of 3D similarity transformation in EIV model
Initiation: Input 3D coordinates vector p o , p t of control points and pointwise weight matrix W , set the initial value of the unknown x as 1 0 0 0 T and the initial value of the unknown e o as 0 0 · · · 0 T 3n Iterative computation: Step 1. Calculate the centering matrix D by Eq. (22), and then calculate the centralized coordinate vectors p o , p t by Eqs. (27) and (26) Step 2. Compute R by Eq. (8) or Eq. (10), W 1 by Eq. (16), and A j by Eqs. (29) to (37) Step 3. Compute dx by Eq. (40). If absolute values of all the elements of dx are less than a threshold τ (1.0 × 10 −10 is set in the paper), turn to Step 4, otherwise turn to Step 2 with the undated value x = x + dx and e o by Eq. (41) and Eq. (43) Step 4. Compute R by Eq. (8) or Eq. (10). If rotation angles θ x , θ y , θ z are needed, compute them by Eq. (4). Compute t by Eqs. (18) to (20) Step 5. Estimate the variance factor of unit weight σ by Eq. (44) or Eq. (45), estimate the covariance matrix of transformation parameters x and t by Eqs. (46) and (47), respectively Output the scale factor (the first element of x ), the Gibbs vector v = a b c T (the latter three elements of x ), the rotation matrix R and the translation parameter t . If required, output rotation angles θ x , θ y , θ z . Output the variance factor of unit weight σ and the covariance matrixes of transformation parameters x and t points are computed with the algorithm in Wang et al. (2014), one finds that LS variance factor of unit weight is 0.0234 m, which is bigger than that of WTLS. Therefore, the WTLS solution is better than the LS solution from the viewpoint of variance factor of unit weight. Next, the new algorithm in this contribution and the Felus-Burtch Algorithm are used to calculate the target coordinate of check points once the transformation parameters are recovered by the 3D similarity transformation model, i.e., Eq. (1). The errors of computed coordinates of check points are obtained by subtracting the computed coordinates of check points from the known coordinates of check points. The errors are listed in Table 6.
Lastly, in order to test the dependence of the new algorithm on the initial values of parameters, six sets of initial values of parameters are adopted which are listed in Table 7. From Table 7, it is seen that the maximum biased angle from the correct solution is 2.5°, 18.9°, 29.4°, 44.5°, 59.4°, 74.9° for Set 1, Set 2, Set 3, Set 4, Set 5, Set 6, respectively. Set 3 gives the default initial values of parameters. The corresponding initial values of the Gibbs vector elements are computed from the initial values of rotation angles in Table 7. Then the new algorithm is performed with six sets of initial values of parameters, respectively. Results show that the algorithm successively finds the correct solutions with six sets of initial values of parameters. Number of iterations is 5, 5, 6, 6, 8, and 8 for Set 1, Set 2, Set 3, Set 4, Set 5, Set 6, respectively. It is evidently seen that with the increase of maximum biased angle from Set 1 to Set 6, the number of iterations increases slowly from 5 to 8. Thus the new algorithm is not sensitive to the initial values of parameters. Further, the iterative process is drawn in Fig. 2. From this figure, it is seen that the convergence rate is fast since that the logarithm (base 10) of objective function, i.e., log 10 e oT We o + e tT We t linearly decreases for all six sets of initial values. So the new algorithm is very fast in terms of fewer iterations. It is worth noting that for this big rotation angle case, no good initial values of parameters are required while usually good initial values of parameters are needed, as shown in Zeng and Yi (2011).

Case 2 (geodetic datum transformation)
The data are from Grafarend and Awange (2003). The 3D coordinates of 7 points in the source coordinate system and the target coordinate system, respectively, are given in Table 8. The 3D space distribution of 7 points in the source coordinate system is shown in Fig. 3. In    Fig. 3 based on the consideration that usually, in order to obtain good transformation parameters, control points should be distributed evenly and enclose as much of the surveying area as possible. The blue asterisk denotes the control points, and the red circle denotes the check points in Fig. 3. The selected four control points are shown with bold font in Table 8. The weight matrix adopted is a pointwise one, i.e., for each point, it has isotropic weight and its different coordinate components are not correlated, and for all points, they are independent of each other. The pointwise weights are listed in Table 9. The recovered transformation parameters in the WTLS sense by the new algorithm and the Felus-Burtch Algorithm are listed in Table 10. It is seen from Table 10 that the new algorithm needs only 2 iterations to complete computation. The results are identical for the two algorithms if the decimal rounding is ignored. Thus the new algorithm is fast and reliable for this case. As mentioned in the preceding case, the new algorithm can provide the accuracies of recovered transformation parameters while the Felus-Burtch Algorithm cannot. The accuracy information is shown in Table 10. For   Table 11. The predicted errors of coordinates of control points by the new algorithm are listed in Table 12. Table 12 shows as far as absolute values are concerned, the predicted errors of coordinates of the new algorithm control points in the source coordinate system are approximately identical to those in the target coordinate system. It is because the point weights are the same for the two systems, the scale factor is very close to 1 and the rotation angles are very small resulting in almost an identity rotation matrix. Next, the new algorithm and the Felus-Burtch Algorithm are employed to calculate the target coordinates of check points. Then errors of computed coordinates of the check points are computed. The errors are the same for the two algorithms and are listed in Table 13.

Conclusions
The paper presents a WTLS iterative algorithm of the 3D similarity coordinate transformation based on the Gibbs vector considering that transformation model is compact and without redundant constraints of parameter, such as the norm of the unit quaternion is 1. The