PERFORMING 3D SIMILARITY TRANSFORMATION WITH LARGE ROTATION ANGLES USING CONSTRAINED MULTIVARIATE TOTAL LEAST SQUARES Realizando transformação de similaridade 3D com grandes ângulos de rotação usando quadrados mínimos multivariados restritos

3D similarity transformation is frequently encountered operation in the field of geodetic data processing, and there are many applications that involve large rotation angles. In previous studies, the errors of the coefficient matrix were usually neglected and a least squares algorithm was applied to calculate the transformation parameters. However, the coefficient matrix is composed of the point coordinates in source coordinate system, i.e., the coefficient matrix is also contaminated by errors. Therefore, a total least squares algorithm should be applied. In this paper, a new method is proposed to address the 3D similarity transformation problem with large rotation angles. Firstly, the scale factor and rotation matrix are put together as the parameter matrix to avoid the rank-defect problem. Then, the translation vector is removed and the multivariate model is constructed. Finally, the constraints are introduced according to the properties of the parameter matrix and the constrained multivariate total least squares algorithm is derived to obtain the transformation parameters. The experimental results show that the proposed method has a high computational efficiency.


Introduction
The 3D coordinate transformation is a fundamental problem in the field of geodesy, photogrammetry and map projection etc. The 3D coordinate transformation model with seven parameters (i.e., scale factor, three translation parameters and three rotation angles), also known as similarity transformation, is the most frequently employed. When the rotation angles are small, the Bursa-Wolf model, Molodensky-Badekas model, Veis model and Vaniček-Wells model can be used in the 3D similarity transformation problem (Rapp 1993). Yang (1999) employed the robust least squares in geodetic datum transformation so as to resist the influence of gross errors. You and Hwang (2006) utilized the least squares collocation in geodetic datum transformation so as to compensate the systematic distortion. However, the 3D similarity transformation with large rotation angles is among the frequently encountered operation in a practical project, e.g., the transformation between physical coordinate system and image coordinate system in photogrammetry. Therefore, it is necessary to develop a solution for the 3D similarity transformation with large rotation angles. In Chen et al. (2004), the nine elements of the rotation matrix were used to replace the three rotation angles. The coordinate transformation model with thirteen parameters was established, which can be used to deal with the 3D similarity transformation with large rotation angles, and the constraints were introduced because the rotation matrix is orthogonal. Then, the constrained least squares (CLS) algorithm was applied to estimate the transformation parameters. When both the coefficient matrix and the observation vector are affected by random errors, the total least squares (TLS) algorithm with errors-in-variables (EIV) model should be employed, which was termed by Golub and van Loan (1980). Numerous scholars have paid attention to TLS algorithm. De Moor (1993) considered the structural characteristics of the coefficient matrix, the structured TLS (STLS) algorithm was studied. Yan and Fan (2001) proposed the mixed LS-TLS method, in which the constant columns were separated from the coefficient matrix. However, the STLS algorithm and mixed LS-TLS algorithm are not satisfactory in measurement data processing. For this reason, Xu et al. (2012) extracted the random elements from the coefficient matrix, the weighted total least squares (WTLS) algorithm with partial EIV (PEIV) model was developed. The quantities that need to be estimated were reduced so that the computational efficiency was improved. Subsequently, Shi et al. (2015) improved the WTLS algorithm with PEIV model, the computational efficiency was further improved. In Schaffrin and Felus (2008), the observation vector and parameter vector were written in matrix form and developed the multivariate total least squares (MTLS) algorithm, but for the case of identity weight matrix. Fang (2011) formulated the weighted MTLS (WMTLS) algorithm, and the weight matrix was universal. Zeng et al. (2015) derived the WTLS algorithm with inequality constraints. Zhao and Gui (2017) studied the outlier detection for PEIV model. Because different observation types usually lead to a model with heterogenous observations, variance component estimation for the EIV model was also derived (Xu et al. 2014;Wang and Xu 2016).
The TLS algorithm has been widely applied in the 3D coordinate transformation as well (Fang 2014;Fang 2015). Lu et al. (2014) took the gross errors in observed coordinates into consideration, and constructed the robust weighted total least squares (RWTLS) algorithm. In Laari et al. (2016) and Okwuashi and Eyoh (2012), the TLS solution based on singular value decomposition was applied in the Bursa-Wolf and Molodensky-Badekas model. Because the aim of the 3D coordinate transformation is to obtain the coordinates of non-common points in the target coordinate system, Li et al. (2012) formulated the total least squares collocation to solve the 3D coordinate transformation problem. Thus, the errors of common points in both the source and target coordinate system could be taken into consideration. Furthermore, both the errors of non-common points in the source coordinate system and the correlations between the coordinates of common and non-common points were also taken into consideration. Both the computation of the transformation parameters and the transformation of the coordinates of non-common points were integratively implemented. In Lu et al. (2015), the weighted total least squares (WTLS) algorithm based on the nonlinear Gauss-Helmert model (Neitzel 2010) was applied in the 3D similarity transformation with large rotation angles. However, the method of Lu et al. (2015) has potential numerical instability because of the existing of the rank-defect problem (Yang et al. 2017). For this reason, Chang (2016) derived the closed-form least-squares solution, but the method would only be applied in the case with rotational invariant covariance structure. Mercan et al. (2018) provided an efficient solution to 3D similarity transformation, which is based on quaternions. In Amiri-Simkooei (2018), the least-squares variance component estimation (LS-VCE) was introduced into EIV model for handling the 3D similarity transformation with hard constraints.
In this paper, the nine elements in the rotation matrix are treated as the unknown parameters. Thus, the complex formula derivation and tedious computation of trigonometric functions are avoided. The scale factor and rotation matrix are put together as the parameter matrix to avoid the rank-defect problem. Then, the translation parameters are eliminated and the multivariate model is constructed. Although the multivariate model had been constructed for the 3D coordinate transformation by Felus and Burtch (2009), the orthogonal Procrustes algorithm was applied so that their method is only suitable for the condition where the errors of observed coordinates are symmetrical and independent. Here, we develop the constraints according to the properties of the parameter matrix, and the constrained MTLS (CMTLS) algorithm is derived to determine the parameter matrix. Because we eliminate the translation parameters, the constants that result in a large number of computation will not appear in the coefficient matrix. Furthermore, we do not need to obtain the Jacobian matrix whose calculation is time-consuming for calculating the cofactor matrix of the coefficient matrix. Therefore, the proposed algorithm is computationally cheap. The main contributions of our work are summarized as follows: The random errors of both the coefficient matrix and observation matrix are taken into consideration. Thus, the more reliable posteriori standard deviation can be obtained.
The translation vector is eliminated and calculated separately. The multivariate model is constructed to estimate the parameters. As a consequence, we avoid calculating the Jacobian matrix, and the coefficient matrix will not contain a large amount of constants. This reduce the computation time of the algorithm.

3D Similarity Transformation Model with Large Rotation Angles
In the 3D similarity transformation model with large rotation angles, the nine elements in the rotation matrix are treated as the unknown parameters in addition to the three translation and one scale parameters. The nine elements need to be estimated, rather than the three rotation angles (Lu et al. 2015). Let If one wants to estimate the three rotation angles, the rotation matrix needs to be linearized, i.e., seeking the partial derivative of the rotation matrix with respect to the three rotation angles. The formula is rather complicated, and the tedious computation of the trigonometric functions will also be introduced. Therefore, Chen et al. (2004) directly estimated the nine elements in the rotation matrix. Then, to preserve the orthogonality of the rotation matrix, the following constraints should be applied For n number of common points, we have where n 1 denotes the 1 × n vector of ones and ⊗ denotes the Kronecker product operator.

Designing a matrix G so that
We need to make The rotation matrix and scale factor are put together as the parameter matrix. That is The aim is to avoid the rank-defect problem. The reason for the rank-defect problem was present in Yang et al. (2017). Let the coefficient matrix A = GS and the observation matrix Y = GT. Then, equation (9) can be rewritten in the following form According to the Kronecker product and the property of the vectorized operator, we obtain where ) (⋅ vec is the mathematical operator to convert a matrix to a column vector by stacking one column of this matrix underneath the previous one, and 3 I denotes a 3 3× identity matrix. According to the cofactor propagation law, the cofactor matrices of the coefficient matrix and observation matrix are where T Q is the cofactor matrix of the coordinates in the target system, i.e., the cofactor matrix of T . S Q is the cofactor matrix of the coordinates in the source system, i.e., the cofactor matrix of S .
Because the rotation matrix is an orthogonal matrix, we can derive Then, according to equation (15) As we can see, the diagonal elements of In order to adjust the nonlinear model, we need to linearize equations (10) and (18). Whilst the random errors in the coefficient matrix and observation matrix are taken into consideration are the elements of 0 ξ , which is the initial value of the parameter matrix. A E is the random error matrix of the coefficient matrix, and Y E is the random error matrix of the observation matrix.
Although ξ E δ A is the second-order smaller quantity, the iterative computation could diverge or converge to a wrong stationary point if it is neglected (Pope 1972). However, this pitfall can be avoided by using the estimate (19), as suggested by Li et al. (2012). Noting Next, the CMTLS algorithm will be derived. We start with the objective function where λ denotes a matrix of Lagrange multipliers, u denotes a vector of Lagrange multipliers.
The solution of this target function is derived via the Lagrange necessary condition, namely where hats indicate "estimated" quantities.
From equations (26) and (27), we can obtain Inserting equations (31) and (32) into equation (29), we derive Inserting equation (33) into equation (28), we derive Combining equation (34) and equation (30), we can obtain Thus, we get the CMTLS solution Then, the updated parameter matrix is calculated as After the parameter matrix is calculated, we can take the value of any diagonal elements of the matrix T ξξ to calculate the scale factor. For example where T ξξ (1,1) denotes the first diagonal element of T ξξ .
Thus, the rotation matrix is given by Because the observation matrix Y is an 3 ) 1 ( × − n matrix and the parameter matrix ξ is a 3 x 3 matrix, according to Schaffrin and Felus (2008) and Zeng et al. (2015), the standard deviation is estimated as In addition, we can calculate the estimated random error vectors of S and T with where S Ê is the estimated random error matrix of S and T Ê is the estimated random error matrix of T .
After the scale factor is computed by using equation (38) and the rotation matrix is computed by using equation (39), the adjusted coordinates are used to calculate the translation vector of each point as follows Here, we observe that the calculated translation vectors of all the points are the same, so any row ∆ i of ∆ can be considered as the translation vector.
The calculation step is summarized as follows: Step 1: The prepared data include the coordinates in the source system, the coordinates in the target system, the cofactor matrix S Q of the coordinates in the source system, and the cofactor matrix T Q of the coordinates in the target system.
Step 2:Compute the initial value of the parameter matrix 0 ξ and set the initial random error matrix 0 A E of the coefficient matrix as zero matrix.
Step 5: Compute the scale factor and rotation matrix by equations (38) and (39), and compute the translation vector using equation (43).

Three-dimensional coordinate transformation
In this section, the simulation experiment is performed to illustrate the advantages of the derived CMTLS algorithm. Matlab is used to perform the experiment. We randomly choose ten points as the points in the source system. The true values of the scale factor, translation parameters and rotation angles are set as Thus, we can calculate the coordinates of the ten points in the target system. The true values of the coordinates in the two coordinate systems are shown in Figure 1. The blue points are the points in the source coordinate system and the red points are the points in the target coordinate system. The random errors with different variances are added to the true values of the coordinates. The coordinates of point 1 to point 5 in the target system contain the random errors with zero mean and standard deviation of 0.03m. The coordinates of point 6 to point 10 in the target system contain the random errors with zero mean and standard deviation of 0.06m. The coordinates of point 1 to point 5 in the source system contain the random errors with zero mean and standard deviation of 0.09m. The coordinates of point 6 to point 10 in the source system contain the random errors with zero mean and standard deviation of 0.12m.
At first, the WMTLS (Fang 2011) and CMTLS are implemented to obtain the parameter matrix. The priori standard deviation is set as 0.03 to determine the cofactor matrices S Q and T Q . The matrix T ξξ is calculated according to the calculated parameter matrix, as shown in Table 1.  Because the rotation matrix is an orthogonal matrix, the diagonal elements of the matrix T ξξ should be equal and the off-diagonal elements should be equal to 0. From the calculated results of Table 1, we note that the results of the CMTLS are consistent with the constraints. However, in the results calculated by the WMTLS, the diagonal elements are not equal, and non-zero values exist in the off-diagonal elements. Therefore, the results of the CMTLS are better and the additional constraints are necessary in this context.
Next, the CWLS, CWTLS and CMTLS are implemented. The calculated results of the three algorithms are compared to perform analysis. In order to have an insight to the three algorithms, the CWLS and CWTLS for the 3D similarity transformation are presented in Appendix. We carried out 1000 simulation experiments and randomly extract the calculated results of one simulation experiment. The results of the CWLS are as follows:

R
As we can see, the results of the CWTLS are the same as those of the CMTLS. It indicates that the CMTLS is equivalent to the CWTLS in terms of the transformation parameter estimates. The results of the CWLS are slightly different from the results of the CWTLS and CMTLS. This is because the CWLS ignores the errors of the coefficient matrix. However, the estimated standard deviations are different. As shown in Figure 2, the standard deviations calculated by the CWLS are obviously higher than those calculated by the CWTLS and CMTLS. The mean value of the standard deviations calculated by the CWLS for 1000 simulations is 0.0787, which is higher than the priori standard deviation. The mean values of the standard deviations calculated by the CWTLS and CMTLS for 1000 simulations are all 0.0296, which is close to the priori standard deviation. Therefore, the standard deviations obtained by the CWTLS and CMTLS are more reasonable. This advantage is significant for the statistical analysis of the results. The standard deviations calculated by the CWTLS and CMTLS are useful for gross error detection and hypothesis testing etc. Although the calculated results of the CWTLS are the same as those of the CMTLS, the computational efficiency of two methods is different. The cofactor matrix of the coefficient matrix cannot be directly obtained when implementing the CWTLS. We need to compute the jacobian matrix of the coefficient matrix with respect to the coordinates in the source coordinate system. Then, we can obtain the cofactor matrix of the coefficient matrix according to the cofactor propagation law. In contrast, the cofactor matrices of the coefficient matrix and observation matrix can be obtained without extra computation within the CMTLS, because the matrix G is given and easy to compile code, as shown in equations (13) and (14). Therefore, compared with the CWTLS, the CMTLS is easier to implement. In addition, for the CWTLS, there are numerous constants in the coefficient matrix, which result in a large computational burden.
In order to compare the computational efficiency of the CWTLS and CMTLS, we simulate five common points, ten common points, fifteen common points and twenty common points, respectively. The computation time is plotted as a function of the number of the common points. As shown in Figure 3, the computation time is obviously less when implementing the CMTLS. With the increasing number of the common points, the growth of the computation time is also slow. However, when implementing the CWTLS, the growth of the computation time is significant. Hence the CWTLS algorithm is relatively time-consuming. For big data processing, the computational efficiency is very important. In some practical project, the point coordinates are obtained by control network adjustment, so the cofactor matrix of the point coordinates should be fully-populated, i.e., the cofactor matrix is a non-diagonal matrix. For this reason, a case with non-diagonal cofactor matrices S Q and T Q is presented. The data are the twenty common points used before. The random errors with standard deviation of 0.01m are added to the true values of the coordinates. The priori standard deviation is set as 0.01. Thus, the cofactor matrix of each point is the identity matrix. We assume that the point coordinates are correlated. The empirical covariance function in Li et al. (2013) is employed to construct the cofactor matrices. The correlation coefficient of two arbitrary points is where ij d is the distance between points i and j , and 0 d is a constant to govern the correlation strength, which is set as 1000 in this experiment. The maximum, minimum and mean correlation coefficients are respectively 0.9412, 0.0513 and 0.2813. Note that the correlation coefficient is only computed between two points in the same coordinate system.

R
It can be seen that the results of the CMTLS and CWTLS are the same, and the results of the CWLS are slightly different from those of the CMTLS and CWTLS. The estimated standard deviation and the computation time are listed in Table 2. As we can see, the standard deviations calculated by the CWTLS and CMTLS are closer to the priori standard deviation, while the one calculated by the CWLS is much far away from the priori standard deviation. We can also see that the computation time of the CMTLS is significantly smaller than that of the CWTLS. For the case with non-diagonal cofactor matrices, the CMTLS still yields the more reliable standard deviation in comparison with the CWLS and has a significantly better computational efficiency in comparison with the CWTLS. In our algorithm, the scale factor and rotation matrix are put together as the parameter matrix, and the translation vector is calculated by equation (43), so our algorithm can only get the variances of the parameters (the elements in the parameter matrix) rather than the variances of the transformation parameters. This is a drawback of the proposed algorithm.

Point cloud registration
In order to indicate the importance of computational efficiency, a case involved point cloud registration is provided. The iterative closest point (ICP) algorithm (Besl and McKay, 1992) is the most popular registration method. It alternates between two steps until the convergence condition is reached. The first step is to establish correspondences, and the second step is to estimate the transformation. Here, the CWTLS and CMTLS algorithms are used for the second step. The TLS adjustment has already been used for the ICP-like algorithm by Tao et al. (2018) and Ge and Wunderlich (2016).
The data are the point clouds of the "Man head" model, as shown in Figure 4(a). The two point clouds respectively contain 11283 and 10862 points. We down-sample them to 25%, 50% and 75% of their original points. The computation time required by the CWTLS and CMTLS algorithms is recorded under different sampling rates, which is shown in Figure 4(b). The registration result without down-sampling is also illustrated in Figure 4(a). As we can see, the computation time of the CWTLS algorithm is significantly more than that of the CMTLS algorithm. For the point cloud registration, it often involves a large number of points, so the computation efficiency is crucial.

Conclusion
The 3D similarity transformation with large rotation angles is studied in this paper. A new method is developed to obtain the transformation parameters with multivariate model. Considering the property of the parameter matrix, the CMTLS algorithm is derived in this paper. Through the experiments, we can conclude: Because the rotation matrix is an orthogonal matrix, the parameter matrix, which is the product of the scale factor and transposed rotation matrix, should satisfy the proper constraints. By comparing the results of the WMTLS and CMTLS, the matrix T ξξ calculated by the CMTLS is more consistent with the constraints. Hence the constraints are necessary.
The CMTLS algorithm and CWTLS algorithm take the errors of the coordinates in the two coordinate systems into consideration. The calculated standard deviations are more reliable, while the standard deviation calculated by the CWLS are greatly bigger than the priori standard deviation.
The CWTLS algorithm is more time-consuming due to a large number of constants existing in the coefficient matrix and the calculation of jacobian matrix. Hence, the computational efficiency of the CMTLS is higher compared to the CWTLS, so it is useful to process big data.

APPENDIX
The CWLS algorithm and CWTLS algorithm for the 3D similarity transformation are presented in this section. where 3 5× 0 is a 3 5× matrix whose elements are all zero, the meanings of other letters are the same as those mentioned above.
If the errors of the coefficient matrix B are neglected, the CWLS algorithm can be used to solve equation (A3). The cofactor matrix of L is T Q when implementing the CWLS algorithm. If the errors of the coefficient matrix are considered, the CWTLS algorithm can be used to solve equation (A3). The cofactor matrix of the coefficient , and the cofactor matrix of L is still T Q . M is the jacobian matrix of the coefficient matrix with respect to the coordinates in the source coordinate system. The CWLS algorithm and CWTLS algorithm have been studied extensively. Therefore, no longer described here. It is worth noting that the translation vector is not eliminated and the parameter matrix î is still the product of the scale factor and transposed rotation matrix.