Analytical algorithm of weighted 3D datum transformation using the constraint of orthonormal matrix

Based on the Lagrangian extremum law with the constraint that rotation matrix is an orthonormal matrix, the paper presents a new analytical algorithm of weighted 3D datum transformation. It is a stepwise algorithm. Firstly, the rotation matrix is computed using eigenvalue-eigenvector decomposition. Then, the scale parameter is computed with computed rotation matrix. Lastly, the translation parameters are computed with computed rotation matrix and scale parameter. The paper investigates the stability of the presented algorithm in the cases that the common points are distributed in 3D, 2D, and 1D spaces including the approximate 2D and 1D spaces, and gives the corresponding modified formula of rotation matrix. The comparison of the presented algorithm and classic Procrustes algorithm is investigated, and an improved Procrustes algorithm is presented since that the classic Procrustes algorithm may yield a reflection rather than a rotation in the cases that the common points are distributed in 2D space. A simulative numerical case and a practical case are illustrated.


Background
Three-dimensional datum transformation is a frequently used work in geodesy, engineering surveying, photogrammetry, mapping, geographical information science (GIS), machine vision, etc., e.g., Aktuğ (2009), Akyilmaz (2007), El-Mowafy et al (2009), Ge et al (2013, Han and Van Gelder (2006), Horn (1986), Kashani (2006), Neitzel (2010), Paláncz et al (2013), Soler (1998), Soler and Snay (2004), Soycan and Soycan (2008), Zeng (2014). Usually, in order to compute the transformed coordinate, the transformation parameters in the transformation model (e.g., seven-parameter similarity transformation, see Aktuğ 2012, Leick 2004, Leick and van Gelder 1975 need to be solved with several control points in advance. So far, a large number of algorithms for recovering the parameters have been presented, which can be divided into two classes. One is the numerical iterative algorithm, and the other one is analytical algorithm. The former needs the initial parameter values, linearization, and iterative computation, e.g., Zeng and Tao (2003), Chen et al. (2004), Zeng and Huang (2008), El-Habiby et al. (2009), Zeng and Yi (2011), etc. In the case that the rotation angle is large, the initial values are difficult even never to be obtained in advance, and consequently, it leads to the failure of solution (see Zeng and Yi 2011). We should note that if global optimization algorithms are used, then no initial values are required (see, e.g., Xu 2003a, Xu 2003b, Xu 2002. In contrast, the latter does not involve the initial parameter values, linearization, and iterative computation, and can give the exact solution quickly. However, because of the complexity of mathematical derivation, only several analytical algorithms have been put forward. Grafarend and Awange (2003) presented the Procrustes algorithm which utilized the singular value decomposition technique. Shen et al. (2006) presented a quaternion-based algorithm which utilized the quaternion property and eigenvalue-eigenvector decomposition. Han (2010) presented a stepwise approach to individually calculate the transformation parameters by the physical interpretation of similarity transformation. Zeng and Yi (2010) presented a new analytical algorithm based on the good properties of Rodrigues matrix and Gibbs vector.
The present study is organized as follows. In the Methods section, a new analytical algorithm to weighted 3D datum transformation is derived in detail, based on the Lagrangian extremum law with the constraint that the rotation matrix is an orthonormal matrix. In the meanwhile, its stability is discussed when the distribution of 3D control points degenerates into 2D (planar) or even 1D (collinear). The presented algorithm and classic Procrustes algorithm are compared, and an improved Procrustes algorithm is presented since that the classic Procrustes algorithm may yield a reflection rather than a rotation in the cases that the common points are distributed in 2D space. In the Results and discussion section, a simulative numerical case and a practical case are given to demonstrate the presented algorithm, classic Procrustes algorithm, and improved Procrustes algorithm. Lastly, conclusions are made in the Conclusions section.

Presentation of the basic algorithm
The seven-parameter similarity transformation model can be expressed as subject to where a i ¼ X i Y i Z i ½ T and b i ¼ x i y i z i ½ T i = 1, 2, ⋯, n are the 3D coordinates of a common point in target and source coordinate systems of transformation, tagged as system A and system B, respectively. Superscript T stands for transpose, I 3 denotes a 3 × 3 identity matrix, and det means computing the determinant of matrix. λ denotes the scale parameter, t ¼ ΔX ΔY ΔZ ½ T denotes the three translation parameters, and R denotes the 3 × 3 rotation matrix, which contains the three rotation angles. Supposing R is formed by rotating angles α, β, and γ counterclockwise around the Cartesian X, Y, and Z axes, respectively, then R can be expressed by rotation angles as R ¼ cos γ cos β sin γ cos α þ cos γ sin β sinα sin γ sin α− cos γ sin β cos α − sin γ cos β cos γ cos α− sin γ sin β sin α cos γ sin α þ sin γ sin β cos α sin β − cos β sin α cos β cos α Using (3), the rotation angles α, β, and γ can be computed if R is recovered as where R ij is the element of R in the ith row and jth column. Introducing the following matrix form of the coordinates as then Eq.
(1) is rewritten as where 1 n ¼ 1 ⋯ 1 ½ |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} n , i.e., a row vector with n elements and all elements are 1. It is obvious that in order to determine the seven parameters, the number n of common points must be greater than or equal to 3. Considering the coordinates include errors, Eq. (7) is transformed as where E is the transformation error matrix. The criterion of the least squares can be constructed by the Lagrangian extremum law with the constraint of Eq.
(2), i.e., orthonormal matrix as follows. It is worthy of note that the constraint det(R) = + 1 is not imposed, since it can be separately treated at some extra computation as in the Stability of the basic algorithm and its modification section.
where tr denotes trace operation of matrix, Λ is a symmetric Lagrangian multiplier matrix, and P represents the weight matrix that every point has an isotropic weight and is independent of each other. Substituting the expression of E easily obtained from Eq. (7) into Eq.
(8), one can obtain If and only if the following conditions are satisfied, the Lagrangian extremum exists.
By Eqs. (9) and (10), one gets and thus, Obviously, t is the function form of λ and R. Substituting Eq. (16) into Eq. (7), one gets where I n − (1 n P1 n T ) − 1 P1 n T 1 n is the centering matrix.
, and thus, they are the centralized coordinate matrix, and then, Eq. (17) is written as Derivation of Eq. (19) makes use of the properties of trace operation, i.e., Substituting Eq. (19) into Eq. (11), one gets Obviously, λ is the function form of R. Substituting Eq. (19) into Eq. (12), one gets and thus, Substituting Eq. (19) into Eq. (13), one gets the following equation and the derivation process is given in the Appendix.
further substituting Eq. (24) into Eq. (25), one gets and substituting Eq. (27) into Eq. (24), one gets Let then Eq. (28) can be rewritten as Note that D T D is symmetric, non-negative definitive and so has non-negative real eigenvalues. The inverse of the square root of D T D can thus be computed using eigenvalue-eigenvector decomposition. where d i and v i for i = 1, 2, 3 are the eigenvalues and corresponding eigenvectors of the matrix D T D. So, Eq. (30) can be written as Stability of the basic algorithm and its modification Obviously, the construction of the inverse of the square root of D T D , i.e., Eq. (31), fails if one or two of d i for i = 1, 2, 3 equals to 0. Assume that the eigenvalues of the matrix D T D satisfies the following condition.
When the common points are distributed in a plane, i.e., 2D space, the matrix D T D is singular and of rank 2, and thus, One can compute the rotation matrix as follows.
where the sign ± is chosen so that det(R) = + 1 is satisfied.
The above case is the ideal one, but in the case that the common points are distributed in an approximate plane, although   When the common points are distributed in a line, i.e., 1D space, the matrix D T D is singular and of rank 1, and thus, The rotation matrix is impossible to recover in a whole; however, one can recover at most two rotation angles by the following formula.
and the utilization of Eq. (40) make it feasible to compute the translation parameter and scale parameter. The above case is the ideal one, but in the case that the common points are distributed in an approximate line, although

Comparison to classic Procrustes algorithm and improvement of the classic Procrustes algorithm
The classic Procrustes algorithm presented by Grafarend and Awange (2003) is a well-known analytical algorithm of 3D datum transformation. It is also based on the Lagrangian extremum law, similarly to the presented algorithm in this paper. But differently, it does not constrain the orthonormal matrix condition. For the Procrustes algorithm, due to utilization of the singular value decomposition technique, the computed rotation matrix always satisfies the constraint condition R T R = I 3 ; however, in the cases that the common points are   Translation  Angle  Scale  Total  Translation  Angle  Scale  Total  Translation  Angle  Scale  Total   3 D  3  3  1  7  3  3  1  7  3  3  1  7   2D  3  3  1  7 0-3 0 -3 1 1 -7 3 3 1 7 1D 3 0-2 1 4 -6 3 0 -2 1 4 -6 3 0 -2 1 4 -6 distributed in a rigid or approximate plane, this situation that det(R) = − 1 rather than det(R) = + 1 usually happens. This means the computed R is a reflection instead of a rotation. For the presented algorithm in this paper, the constraint condition R T R = I 3 is imposed in the computation, and thus, it is always satisfied. And in the cases that the common points are distributed in a rigid or approximate plane, the sign ± in Eq. (36) is properly chosen so that det(R) = + 1 is satisfied.
To recover the exact rotation matrix by Procrustes algorithm when det(R) = − 1, the computation formula of R, i.e., Eq. (22) in Grafarend and Awange (2003) should be improved as whereṼ ¼ V 1 V 2 −V 3 ½ , V 1 , V 2 , and V 3 are the column matrix of V, and V 3 is the column matrix that corresponds to the singular value that is 0. Therefore, the general improved computation formula of R, i.e., Eq. (22) in Grafarend and Awange (2003) is which is stable for the cases that the common points are distributed in 3D and 2D spaces including approximate 2D space. In the case that the common points are distributed in 1D space, Eq. (43) can recover at most two rotation angles and be used to compute the translation parameter and scale parameter.

Simulative case
The case data is simulated as follows. In order to investigate the stability performance of the presented algorithm (PA) in this paper, the classic Procrustes algorithm (CPA) and improved Procrustes algorithm (IPA) in the cases that the control points are distributed in 3D, 2D, and 1D spaces, six sets of control point in system B is first given in Table 1, of which set 1 is distributed in 3D space, sets 2, 3, and 4 are distributed in 2D space, and sets 5 and 6 are distributed in 1D space. The distribution of six sets of control point is depicted intuitively in Fig. 1. Set 2 has only three control points, and it is necessary for the least point number to solve the seven parameters. Secondly, the theoretical seven parameters are given in Table 2. For the sake of an efficient test of the algorithms, the rotation angles are designed to be big angles. Thirdly, the coordinates of control points in system A are computed by Eq. (1), and the result is listed in Table 3. In this case, the stability of the three algorithms is focused, so the weight matrix is designated to identity matrix for easy demonstration. Next, the transformation seven parameters are recovered by the three algorithms, and the result is listed in Table 4. ME in Table 4 is the mean error and computed by  Fig. 2 The distribution of seven control points in local system It is seen from Table 4 that the results of seven parameters are identical and accurate for the three algorithms in the case of set 1, i.e., the case that the control points are distributed in 3D space. For the cases that the control points are distributed in 2D space, PA and IPA have the identical and accurate results of seven parameters, but CPA has one correct result of seven parameters for set 4 and two wrong results of rotation angles for sets 2 and 3, because reflections rather than rotations are found, and all translation and scale parameters are correct. For the cases that the control points are distributed in 1D space, i.e., sets 5 and 6, the three algorithms all recover the correct translation and scale parameters and at most recover two rotation angles. The number of solvable transformation parameters for the three algorithms is counted and listed in Table 5.

Actual case
The case data is from Grafarend and Awange (2003). The coordinates of control points in system B (local system) and A (WGS-84 system) is listed in Table 6. The distribution of seven control points in local system is depicted in Fig. 2. From this figure, it is seen that the distribution of control points are in an approximate plane. In the process of PA computation, condition number of matrix D T D is 2.5 × 10 11 , so Eq. (32) is seriously ill-conditioned and yields a biased solution if not processed. When the weight matrix is an identity matrix, the computed results of seven parameters with PA and CPA are listed in Table 7. For the situation that the weight matrix is a point-wise matrix, i.e., every point has isotropic weight and is independent of each other, the point-wise matrix is generated by the way introduced in Grafarend and Awange (2003) and is listed in Table 8. The computed results of seven parameters with PA and CPA are listed in Table 9.
It is seen from Tables 7 and 9 that the results of PA and CPA are identical if the bias caused by decimal rounding is ignored. Hence, the PA is comparable with CPA.

Conclusions
The numerical case study shows that the presented new algorithm and improved Procrustes algorithm are both stable and reliable for the cases that the control points are distributed in 3D, and 2D including approximate 2D space, and can recover at most two angles as well as all translation and scale parameters for the cases that the control points are distributed in 1D space. The classic Procrustes algorithm also can compute all translation and scale parameters for all the cases that the control   points are distributed in 3D, 2D, and 1D spaces; however, it may yield a reflection rather than a rotation in the cases that the common points are distributed in 2D space. The numerical case study also shows the presented algorithm in this paper, and improved Procrustes algorithm is both stable and reliable when the rotation angles are big. And the presented algorithm in this paper is comparable with the classic Procrustes algorithm when the point-wise weight matrix is involved.