PERFORMING 3D SIMILARITY TRANSFORMATION WITH LARGE ROTATION ANGLES USING CONSTRAINED MULTIVARIATE TOTAL LEAST SQUARES

Wuyong Tao; Xianghong Hua; Shaoquan Feng

resúmenes

secciones

referencias

imágenes

Abstract: 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.

Keywords: 3D similarity transformation3D similarity transformation,constraintconstraint,large rotation anglelarge rotation angle,multivariate total least squaresmultivariate total least squares.

Resumo: A transformação da similaridade 3D é frequentemente encontrada no campo do processamento geodésico de dados, e há muitas aplicações que envolvem grandes ângulos de rotação. Em estudos anteriores, os erros da matriz de coeficientes foram geralmente negligenciados e um algoritmo de menos quadrados foi aplicado para calcular os parâmetros de transformação. No entanto, a matriz de coeficiente é composta pelas coordenadas de ponto no sistema de coordenadas de origem, ou seja, a matriz de coeficiente também está contaminada por erros. Portanto, um algoritmo total de quadrados mínimos deve ser aplicado. Neste artigo, propõe-se um novo método para abordar o problema de transformação da similaridade 3D com grandes ângulos de rotação. Em primeiro lugar, o fator de escala e a matriz de rotação são montados como a matriz de parâmetros para evitar o problema de defeito de classificação. Em seguida, o vetor de tradução é removido e o modelo multivariado é construído. Finalmente, as restrições são introduzidas de acordo com as propriedades da matriz de parâmetros e o algoritmo multivariado multivariado restrito de menos quadrados é derivado para obter os parâmetros de transformação. Os resultados experimentais mostram que o método proposto tem alta eficiência computacional.

Palavras-chave: transformação de similaridade 3D, restrição, grande ângulo de rotação, multivariar total menos quadrados.

Carátula del artículo

ORIGINAL ARTICLE

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

Wuyong Tao wuyong_tao@whu.edu.cnxhhua@sgg.whu.edu.cn2017202140047@whu.edu.cn

Wuhan University, China

Xianghong Hua wuyong_tao@whu.edu.cnxhhua@sgg.whu.edu.cn2017202140047@whu.edu.cn

Wuhan University, China

Shaoquan Feng wuyong_tao@whu.edu.cnxhhua@sgg.whu.edu.cn2017202140047@whu.edu.cn

Wuhan University, China

Boletim de Ciências Geodésicas, vol. 26, no. 4, 2020
Universidade Federal do Paraná

Received: 20 March 2019

Accepted: 08 November 2020

DOI: 10.1590/s1982-21702020000400021

Funding
Funding source: National Natural Science Foundation of China
Contract number: 41674005; 41374011
Funding
Funding source: China Scholarship Council
Contract number: scholarship (201906270179)

1. 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:

1. 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.
2. 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.

2. 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 (X, Y, Z) denote the coordinates in the target system and (x, y, z) denote the coordinates of the corresponding points in the source system. The mathematical model of the transformation from the source coordinate system to the target coordinate system reads

[\begin{matrix} X \\ Y \\ Z \end{matrix}] = k R [\begin{matrix} x \\ y \\ z \end{matrix}] + [\begin{matrix} X_{0} \\ Y_{0} \\ Z_{0} \end{matrix}]

(1)

where k is the scale factor, X₀, Y₀, Z₀ are the three translation parameters, and R is the rotation matrix.

The rotation matrix is the product of three individual rotation matrices as follow

R = R_{z} R_{y} R_{x} = [\begin{matrix} r_{1} & r_{4} & r_{7} \\ r_{2} & r_{5} & r_{8} \\ r_{3} & r_{6} & r_{9} \end{matrix}]

(2)

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

[\begin{matrix} r_{1}^{2} + r_{2}^{2} + r_{3}^{2} \\ r_{4}^{2} + r_{5}^{2} + r_{6}^{2} \\ r_{7}^{2} + r_{8}^{2} + r_{9}^{2} \\ r_{1} r_{4} + r_{2} r_{5} + r_{3} r_{6} \\ r_{1} r_{7} + r_{2} r_{8} + r_{3} r_{9} \\ r_{4} r_{7} + r_{5} r_{8} + r_{6} r_{9} \end{matrix}] = [\begin{matrix} 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \end{matrix}]

(3)

3. The Constrained Multivariate Total Least Squares Algorithm

Transposing equation (1), we can obtain

[\begin{matrix} X & Y & Z \end{matrix}] = k [\begin{matrix} x & y & z \end{matrix}] R^{T} + [\begin{matrix} X_{0} & Y_{0} & Z_{0} \end{matrix}]

(4)

For n number of common points, we have

[\begin{matrix} X_{1} & Y_{1} & Z_{1} \\ X_{2} & Y_{2} & Z_{2} \\ ⋮ & ⋮ & ⋮ \\ X_{n} & Y_{n} & Z_{n} \end{matrix}] = k [\begin{matrix} x_{1} & y_{1} & z_{1} \\ x_{2} & y_{2} & z_{2} \\ ⋮ & ⋮ & ⋮ \\ x_{n} & y_{n} & z_{n} \end{matrix}] R^{T} + 1_{n} \otimes [\begin{matrix} X_{0} & Y_{0} & Z_{0} \end{matrix}]

(5)

where 1_n denotes the n x 1 vector of ones and ⊗ denotes the Kronecker product operator.

Designing a matrix G so that $G (1_{n} \otimes [\begin{matrix} X_{0} & Y_{0} & Z_{0} \end{matrix}]) = 0$ , thus we can eliminate the translation vector. Through deduction, we have

\begin{matrix} G (1_{n} \otimes [\begin{matrix} X_{0} & Y_{0} & Z_{0} \end{matrix}]) = (G \otimes 1) (1_{n} \otimes [\begin{matrix} X_{0} & Y_{0} & Z_{0} \end{matrix}]) \\ = (G \cdot 1_{n} \otimes [\begin{matrix} X_{0} & Y_{0} & Z_{0} \end{matrix}]) \end{matrix}

(6)

We need to make $G \cdot 1_{n} = 0$ , so we can obtain

G = [\begin{matrix} - 1 & 1 & 0 & \dots & 0 \\ - 1 & 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ - 1 & 0 & 0 & \dots & 1 \end{matrix}]

(7)

where G is an (n - 1) x n matrix. Then, multiplying equation (5) by matrix G , we obtain

G [\begin{matrix} X_{1} & Y_{1} & Z_{1} \\ X_{2} & Y_{2} & Z_{2} \\ ⋮ & ⋮ & ⋮ \\ X_{n} & Y_{n} & Z_{n} \end{matrix}] = k G [\begin{matrix} x_{1} & y_{1} & z_{1} \\ x_{2} & y_{2} & z_{2} \\ ⋮ & ⋮ & ⋮ \\ x_{n} & y_{n} & z_{n} \end{matrix}] R^{T}

(8)

Let $T = [\begin{matrix} X_{1} & Y_{1} & Z_{1} \\ X_{2} & Y_{2} & Z_{2} \\ ⋮ & ⋮ & ⋮ \\ X_{n} & Y_{n} & Z_{n} \end{matrix}]$ , $S = [\begin{matrix} x_{1} & y_{1} & z_{1} \\ x_{2} & y_{2} & z_{2} \\ ⋮ & ⋮ & ⋮ \\ x_{n} & y_{n} & z_{n} \end{matrix}]$ . Equation (8) becomes

G T = G S k R^{T}

(9)

The rotation matrix and scale factor are put together as the parameter matrix. That is ξ = kR^T . 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

Y = A ξ

(10)

According to the Kronecker product and the property of the vectorized operator, we obtain

v e c (A) = v e c (G S I_{3}) = (I_{3} \otimes G) v e c (S)

(11)

v e c (Y) = v e c (G T I_{3}) = (I_{3} \otimes G) v e c (T)

(12)

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 I₃ denotes a 3x3 identity matrix. According to the cofactor propagation law, the cofactor matrices of the coefficient matrix and observation matrix are

Q_{A} = (I_{3} \otimes G) Q_{S} (I_{3} \otimes G^{T})

(13)

Q_{Y} = (I_{3} \otimes G) Q_{T} (I_{3} \otimes G^{T})

(14)

where Q_T is the cofactor matrix of the coordinates in the target system, i.e., the cofactor matrix of T . Q_S 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

ξ ξ^{T} = k^{2} R^{T} R = k^{2} I_{3}

(15)

with

ξ = [\begin{matrix} ξ_{1} & ξ_{4} & ξ_{7} \\ ξ_{2} & ξ_{5} & ξ_{8} \\ ξ_{3} & ξ_{6} & ξ_{9} \end{matrix}]

(16)

Then, according to equation (15), we can write

[\begin{matrix} ξ_{1}^{2} + ξ_{4}^{2} + ξ_{7}^{2} & ξ_{1} ξ_{2} + ξ_{4} ξ_{5} + ξ_{7} ξ_{8} & ξ_{1} ξ_{3} + ξ_{4} ξ_{6} + ξ_{7} ξ_{9} \\ ξ_{1} ξ_{2} + ξ_{4} ξ_{5} + ξ_{7} ξ_{8} & ξ_{2}^{2} + ξ_{5}^{2} + ξ_{8}^{2} & ξ_{2} ξ_{3} + ξ_{5} ξ_{6} + ξ_{8} ξ_{9} \\ ξ_{1} ξ_{3} + ξ_{4} ξ_{6} + ξ_{7} ξ_{9} & ξ_{2} ξ_{3} + ξ_{5} ξ_{6} + ξ_{8} ξ_{9} & ξ_{3}^{2} + ξ_{6}^{2} + ξ_{9}^{2} \end{matrix}] = [\begin{matrix} k^{2} & 0 & 0 \\ 0 & k^{2} & 0 \\ 0 & 0 & k^{2} \end{matrix}]

(17)

As we can see, the diagonal elements of ξξ^T have the same value and the off-diagonal elements are equal to 0. Thus, the new constraints for the reduced observation equations are as follows

[\begin{matrix} ξ_{1}^{2} + ξ_{4}^{2} + ξ_{7}^{2} - ξ_{2}^{2} - ξ_{5}^{2} - ξ_{8}^{2} \\ ξ_{1}^{2} + ξ_{4}^{2} + ξ_{7}^{2} - ξ_{3}^{2} - ξ_{6}^{2} - ξ_{9}^{2} \\ ξ_{1} ξ_{2} + ξ_{4} ξ_{5} + ξ_{7} ξ_{8} \\ ξ_{1} ξ_{3} + ξ_{4} ξ_{6} + ξ_{7} ξ_{9} \\ ξ_{2} ξ_{3} + ξ_{5} ξ_{6} + ξ_{8} ξ_{9} \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}]

(18)

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

Y - E_{Y} = (A - E_{A}) ξ_{0} + (A - E_{A}) δ ξ

(19)

C v e c (δ ξ) + W_{0} = 0

(20)

with

C = [\begin{matrix} 2 ξ_{1}^{0} & - 2 ξ_{2}^{0} & 0 & 2 ξ_{4}^{0} & - 2 ξ_{5}^{0} & 0 & 2 ξ_{7}^{0} & - 2 ξ_{8}^{0} & 0 \\ 2 ξ_{1}^{0} & 0 & - 2 ξ_{3}^{0} & 2 ξ_{4}^{0} & 0 & - 2 ξ_{6}^{0} & 2 ξ_{7}^{0} & 0 & - 2 ξ_{9}^{0} \\ ξ_{2}^{0} & ξ_{1}^{0} & 0 & ξ_{5}^{0} & ξ_{4}^{0} & 0 & ξ_{8}^{0} & ξ_{7}^{0} & 0 \\ ξ_{3}^{0} & 0 & ξ_{1}^{0} & ξ_{6}^{0} & 0 & ξ_{4}^{0} & ξ_{9}^{0} & 0 & ξ_{7}^{0} \\ 0 & ξ_{3}^{0} & ξ_{2}^{0} & 0 & ξ_{6}^{0} & ξ_{5}^{0} & 0 & ξ_{9}^{0} & ξ_{8}^{0} \end{matrix}]

(21)

W_{0} = [\begin{matrix} {(ξ_{1}^{0})}^{2} + {(ξ_{4}^{0})}^{2} + {(ξ_{7}^{0})}^{2} - {(ξ_{2}^{0})}^{2} - {(ξ_{5}^{0})}^{2} - {(ξ_{8}^{0})}^{2} \\ {(ξ_{1}^{0})}^{2} + {(ξ_{4}^{0})}^{2} + {(ξ_{7}^{0})}^{2} - {(ξ_{3}^{0})}^{2} - {(ξ_{6}^{0})}^{2} - {(ξ_{9}^{0})}^{2} \\ ξ_{1}^{0} ξ_{2}^{0} + ξ_{4}^{0} ξ_{5}^{0} + ξ_{7}^{0} ξ_{8}^{0} \\ ξ_{1}^{0} ξ_{3}^{0} + ξ_{4}^{0} ξ_{6}^{0} + ξ_{7}^{0} ξ_{9}^{0} \\ ξ_{2}^{0} ξ_{3}^{0} + ξ_{5}^{0} ξ_{6}^{0} + ξ_{8}^{0} ξ_{9}^{0} \end{matrix}]

(22)

where $ξ_{i}^{0} (i = 1,2, \dots, 9)$ are the elements of ξ₀, which is the initial value of the parameter matrix. E_A is the random error matrix of the coefficient matrix, and E_Y 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 $E_{A_{(j)}}$ of the j th iteration instead of E_A in equation (19), as suggested by Li et al. (2012). Noting $A_{(j)} = A - E_{A_{(j)}}$ , we obtain

\{\begin{matrix} Y - E_{Y} = A ξ_{0} - E_{A} ξ_{0} + A_{(j)} δ ξ \\ C v e c (δ ξ) + W_{0} = 0 \end{matrix}

(23)

Next, the CMTLS algorithm will be derived. We start with the objective function

e_{A}^{T} Q_{A}^{- 1} e_{A} + e_{Y}^{T} Q_{Y}^{- 1} e_{Y} = m i n

(24)

Where e_A= vec(E_A) whose cofactor matrix is Q_A , and e_Y= vec(E_Y) whose cofactor matrix is Q_Y . The corresponding Lagrange target function reads

\begin{matrix} Φ (e_{A}, e_{Y}, δ ξ, λ, u) = e_{A}^{T} Q_{A}^{- 1} e_{A} + e_{Y}^{T} Q_{Y}^{- 1} e_{Y} + \\ 2 v e c^{T} (λ) v e c (Y - E_{Y} - A ξ_{0} + E_{A} ξ_{0} - A_{(j)} δ ξ) + 2 u^{T} (C v e c (δ ξ) + W_{0}) \end{matrix}

(25)

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

{\frac{1}{2} \frac{\partial Φ}{\partial e_{A}}|}_{{\hat{e}}_{A}, {\hat{e}}_{Y}, δ \hat{ξ}, \hat{λ}, \hat{u}} = Q_{A}^{- 1} {\hat{e}}_{A} + (ξ_{0} \otimes I_{n - 1}) v e c (\hat{λ}) = 0

(26)

{\frac{1}{2} \frac{\partial Φ}{\partial e_{Y}}|}_{{\hat{e}}_{A}, {\hat{e}}_{Y}, δ \hat{ξ}, \hat{λ}, \hat{u}} = Q_{Y}^{- 1} {\hat{e}}_{Y} - v e c (\hat{λ}) = 0

(27)

{\frac{1}{2} \frac{\partial Φ}{\partial (v e c (δ ξ))}|}_{{\hat{e}}_{A}, {\hat{e}}_{Y}, δ \hat{ξ}, \hat{λ}, \hat{u}} = - (I_{3} \otimes A_{(j)}^{T}) v e c (\hat{λ}) + C^{T} \hat{u} = 0

(28)

{\frac{1}{2} \frac{\partial Φ}{\partial (v e c (λ))}|}_{{\hat{e}}_{A}, {\hat{e}}_{Y}, δ \hat{ξ}, \hat{λ}, \hat{u}} = v e c (Y - A ξ_{0}) - {\hat{e}}_{Y} + (ξ_{0}^{T} \otimes I_{n - 1}) {\hat{e}}_{A} - (I_{3} \otimes A_{(j)}) v e c (δ \hat{ξ}) = 0

(29)

{\frac{1}{2} \frac{\partial Φ}{\partial u}|}_{{\hat{e}}_{A}, {\hat{e}}_{Y}, δ \hat{ξ}, \hat{λ}, \hat{u}} = C v e c (δ \hat{ξ}) + W_{0} = 0

(30)

where hats indicate “estimated” quantities.

From equations (26) and (27), we can obtain

{\hat{e}}_{A} = - Q_{A} (ξ_{0} \otimes I_{n - 1}) v e c (\hat{λ})

(31)

{\hat{e}}_{Y} = Q_{Y} v e c (\hat{λ})

(32)

Inserting equations (31) and (32) into equation (29), we derive

\begin{matrix} v e c (\hat{λ}) = (^{Q_{Y}} (v e c (Y - A ξ_{0}) - (I_{3} \otimes A_{(j)}) v e c (δ \hat{ξ})) \\ = Q_{y}^{- 1} (v e c (Y - A ξ_{0}) - (I_{3} \otimes A_{(j)}) v e c (δ \hat{ξ})) \end{matrix}

(33)

where $Q_{y} = Q_{Y} + (ξ_{0}^{T} \otimes I_{n - 1}) Q_{A} (ξ_{0} \otimes I_{n - 1})$ .

Inserting equation (33) into equation (28), we derive

(I_{3} \otimes A_{(j)}^{T}) Q_{y}^{- 1} (I_{3} \otimes A_{(j)}) v e c (δ \hat{ξ}) + C^{T} \hat{u} = (I_{3} \otimes A_{(j)}^{T}) Q_{y}^{- 1} v e c (Y - A ξ_{0})

(34)

Combining equation (34) and equation (30), we can obtain

[\begin{matrix} N & C^{T} \\ C & 0 \end{matrix}] [\begin{matrix} v e c (δ \hat{ξ}) \\ \hat{u} \end{matrix}] = [\begin{matrix} W \\ - W_{0} \end{matrix}]

(35)

where $N = (I_{3} \otimes A_{(j)}^{T}) Q_{y}^{- 1} (I_{3} \otimes A_{(j)})$ , $W = (I_{3} \otimes A_{(j)}^{T}) Q_{y}^{- 1} v e c (Y - A ξ_{0})$ .

Thus, we get the CMTLS solution

v e c (δ \hat{ξ}) = N^{- 1} (W + C^{T} {(C N^{- 1} C^{T})}^{- 1} (- W_{0} - C N^{- 1} W))

(36)

Then, the updated parameter matrix is calculated as

ξ = ξ_{0} + δ \hat{ξ}

(37)

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

\hat{k} = \sqrt{ξ ξ^{T} (1,1)}

(38)

where ξξ^T (1,1) denotes the first diagonal element of ξξ^T .

Thus, the rotation matrix is given by

\hat{R} = \frac{ξ^{T}}{\hat{k}}

(39)

Because the observation matrix Y is an (n - 1) x 3 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

\hat{σ} = \sqrt{\frac{e_{A}^{T} Q_{A}^{- 1} e_{A} + e_{Y}^{T} Q_{Y}^{- 1} e_{Y}}{3 (n - 1) - 3 \cdot 3 + 5}} = \sqrt{\frac{v e c {(Y - A ξ)}^{T} Q_{y}^{- 1} v e c (Y - A ξ)}{3 n - 7}}

(40)

In addition, we can calculate the estimated random error vectors of S and T with

{\hat{e}}_{S} = v e c ({\hat{E}}_{S}) = Q_{S} (Ι_{3} \otimes G^{T}) Q_{A}^{- 1} {\hat{e}}_{A}

(41)

{\hat{e}}_{T} = v e c ({\hat{E}}_{T}) = Q_{T} (I_{3} \otimes G^{T}) Q_{Y}^{- 1} {\hat{e}}_{Y}

(42)

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

Δ = (T - {\hat{E}}_{T}) - \hat{k} (S - {\hat{E}}_{S}) {\hat{R}}^{T}

(43)

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 Q_S of the coordinates in the source system, and the cofactor matrix Q_T of the coordinates in the target system.

Step 2:Compute the initial value of the parameter matrix ξ₀ and set the initial random error matrix $E_{A_{0}}$ of the coefficient matrix as zero matrix.

v e c (ξ_{(0)}) = {((I_{3} \otimes A^{T}) Q_{Y}^{- 1} (I_{3} \otimes A))}^{- 1} (I_{3} \otimes A^{T}) Q_{Y}^{- 1} v e c (Y)

E_{A_{(0)}} = 0_{(n - 1) \times 3}

Step 3: Using the initial value, we start the following iterative process:

Q_{y} = Q_{Y} + (ξ_{(j)}^{T} \otimes I_{n - 1}) Q_{A} (ξ_{(j)}^{} \otimes I_{n - 1})

A_{(j)} = A - E_{A_{(j)}}

N = (I_{3} \otimes A_{(j)}^{T}) Q_{y}^{- 1} (I_{3} \otimes A_{(j)})

W = (I_{3} \otimes A_{(j)}^{T}) Q_{y}^{- 1} v e c (Y - A ξ_{0})

v e c (δ {\hat{ξ}}_{(j + 1)}) = N^{- 1} (W + C^{T} {(C N^{- 1} C^{T})}^{- 1} (- W_{0} - C N^{- 1} W))

ξ_{(j + 1)} = ξ_{(j)} + δ {\hat{ξ}}_{(j + 1)}

{\hat{e}}_{A_{(j + 1)}} = - Q_{A} (ξ_{(j)} \otimes I_{n - 1}) Q_{y}^{- 1} (v e c (Y - A ξ_{(j)}) - (I_{3} \otimes A_{(j)}) v e c (δ {\hat{ξ}}_{(j + 1)}))

Step 4: Repeat the step (3) until $‖δ {\hat{ξ}}_{(j + 1)}‖ \leq ε$ (ɛ is a preset threshold).

Step 5: Compute the scale factor and rotation matrix by equations (38) and (39), and compute the translation vector using equation (43).

Step 6: Compute the standard deviation using equation (40).

4. Experiments and Analysis

4.1 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

[\begin{matrix} k & X_{0} & Y_{0} & Z_{0} & α_{1} & α_{2} & α_{3} \end{matrix}] = [\begin{matrix} 1.01 & 6 & 7 & 8 & 30^{0} & 45^{0} & 60^{0} \end{matrix}]

(44)

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.

Figure 1:
The ten points in the source coordinate system (a) and target coordinate system (b).

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 Q_S andQ_T . The matrix ξξ^T is calculated according to the calculated parameter matrix, as shown in Table 1.

Table 1:
The calculated results of the two methods

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:

\hat{k} = 1.00997767

\hat{Δ} = [\begin{matrix} 6.05115997 & 6.94359596 & 8.04445605 \end{matrix}]

\hat{R} = [\begin{matrix} 0.35358869 & 0.92676295 & 0.12682848 \\ -0.61234939 & 0.12683902 & 0.78034613 \\ 0.70710908 & -0.35358491 & 0.61235157 \end{matrix}]

The results of the CWTLS are as follows:

\hat{k} = 1.00997271

\hat{Δ} = [\begin{matrix} 6.03809733 & 6.96173005 & 8.02453490 \end{matrix}]

\hat{R} = [\begin{matrix} 0.35357755 & 0.92676721 & 0.12682840 \\ -0.61235068 & 0.12683161 & 0.78034632 \\ 0.70711353 & -0.35357640 & 0.61235134 \end{matrix}]

The results of the CMTLS are as follows:

\hat{k} = 1.00997271

\hat{Δ} = [\begin{matrix} 6.03809733 & 6.96173005 & 8.02453490 \end{matrix}]

\hat{R} = [\begin{matrix} 0.35357755 & 0.92676721 & 0.12682840 \\ -0.61235068 & 0.12683161 & 0.78034632 \\ 0.70711353 & -0.35357640 & 0.61235134 \end{matrix}]

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.

Figure 2:
The standard deviations calculated by the three algorithms for 1000 simulations

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.

Figure 3:
Computation time under different number of common points

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 Q_S and Q_T 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

ρ_{i j} = \frac{1}{1 + (d_{i j} / d_{0})^{2}}

(45)

where d_ij is the distance between points i and j, and 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.

Then, the CWLS, CWTLS and CMTLS are performed to calculate the transformation parameters.The results of the CWLS are as follows:

\hat{k} = 1.01000135

\hat{Δ} = [\begin{matrix} 6.00053850 & 6.98611275 & 8.00795856 \end{matrix}]

\hat{R} = [\begin{matrix} 0.35355997 & 0.92677411 & 0.12682696 \\ -0.61237080 & 0.12683006 & 0.78033078 \\ 0.70710489 & -0.35355886 & 0.61237145 \end{matrix}]

The results of the CWTLS are as follows:

\hat{k} = 1.01000293

\hat{Δ} = [\begin{matrix} 5.99430186 & 6.99458201 & 8.00556967 \end{matrix}]

\hat{R} = [\begin{matrix} 0.35355541 & 0.92677543 & 0.12683003 \\ -0.61237264 & 0.12682506 & 0.78033015 \\ 0.70710559 & -0.35355719 & 0.61237161 \end{matrix}]

The results of the CMTLS are as follows:

\hat{k} = 1.01000293

\hat{Δ} = [\begin{matrix} 5.99430186 & 6.99458201 & 8.00556967 \end{matrix}]

\hat{R} = [\begin{matrix} 0.35355541 & 0.92677543 & 0.12683003 \\ -0.61237264 & 0.12682506 & 0.78033015 \\ 0.70710559 & -0.35355719 & 0.61237161 \end{matrix}]

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.

Table 2:
The estimated variance components and the computation time

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.

4.2 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.

Figure 4:
(a) the original points and the registration result, and (b) the computation time under different sampling rates

5. 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:

1. 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.
2. 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.
3. 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.

Supplementary material

Appendices

6. APPENDIX

The CWLS algorithm and CWTLS algorithm for the 3D similarity transformation are presented in this section. The observation equation is

e_{y} = B X - y

(A1)

where y = vec(T) , $B = [\begin{matrix} I_{3} \otimes S & I_{3} \otimes 1_{n} \end{matrix}]$ , $X = {[\begin{matrix} v e c^{T} (ξ) & Δ \end{matrix}]}^{T}$ and e_y denotes the residual vector of y .

We linearize equation (A1)

e_{y} = B δ X - L

(A2)

where $L = y - B X_{0}$ , $δ X = X - X_{0}$ and X₀ denotes the initiate value of X .

Considering the constraints, we obtain

\{\begin{matrix} e_{y} = B δ X - L \\ [\begin{matrix} C & 0_{5 \times 3} \end{matrix}] δ X + W_{0} = 0 \end{matrix}

(A3)

where 0_5x3 is a 5 x 3 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 Q_T 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 matrix B is Q_B = MQ_SM^T , and the cofactor matrix of L is still Q_T.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.

ACKNOWLEDGMENT

This research was partly supported by the National Natural Science Foundation of China (41674005; 41374011), and partly supported by the China Scholarship Council (CSC) scholarship (201906270179).

REFERENCES

Amiri-Simkooei, A.R. 2018. Parameter estimation in 3D affine and similarity transformation: implementation of variance component estimation. J Geod, 92(11), pp. 1285-1297.

Besl, P. J. and McKay, D. N. 1992. A method for registration of 3-D shapes. IEEE Trans. Pattern Anal. Mach. Intell., 14(2), pp. 239-256.

Chen Y. Shen, Y. Liu, D. 2004. A simplified model of three dimensional-datum transformation adapted to big rotation angle. Geomatics and Information Science of Wuhan University, 29(12), pp. 1101-1105 (In Chinese with English Abstract).

Chang, G. 2016. Closed form least-squares solution to 3D symmetric helmert transformation with rotational covariance structure. Acta Geodaetica et Geophysica, 51(2), pp. 237-244.

De Moor, B. 1993. Structured total least squares and L2 approximation problem. Linear Algebra and Its Application, 188, pp. 163-205.

Felus, Y.A. Burtch, R.C. 2009. On symmetrical three-dimensional datum conversion. GPS Solution, 13(1), pp. 65-74.

Fang, X. 2011. Weighted total least squares solution for applications in geodesy. Dissertation, Leibniz University Hanover, pp. 39-41.

Fang, X. 2014. A total least squares solution for geodetic datum transformation. Acta Geod Geophys, 49(2), pp. 189-207.

Fang, X. 2015. Weighted total least-squares with constraints: a universal formula for geodetic symmetrical transformation. J Geod , 89(5), pp. 459-469.

Golub, G.H. van Loan, C.F. 1980. An analysis of the total least squares problem. SIAM J Numer Anal, 17(6), pp. 883-893.

Ge, X. and Wunderlich, T. 2016. Surface-based matching of 3D point clouds with variable coordinates in source and target system. ISPRS Journal of Photogrammetry and Remote Sensing, 111, pp. 1-12.

Kwuashi, O. Eyoh, E. 2012. 3D coordinate transformation using total least squares. Academic Research International, 3(1), pp. 2223-9944.

Li, B. Shen, Y. Li, W. 2012. The seamless model for three-dimensional datum transformation. Science China Earth Science, 42(7), pp. 1047-1054.

Li, B. Shen, Y. Zhang, X. Li, C. Lou, L. 2013. Seamless multivariate affine error-in-variables transformation and its application to map rectification. International Journal of Geographical Information Science, 27(8), pp. 1572-1592

Lu, J. Chen, Y. Li, B. Fang, X. 2014. Robust total least squares with reweighting iteration for three-dimensional similar transformation, Survey Review, 46(334), pp. 28-36.

Lu, J. Chen, Y. Fang, X. Zheng, B. 2015. Performing 3D similarity transformation using the weighted total least-squares Method. Munich: International Association of Geodesy Symposia, 140, pp. 71-77.

Laari, P.B. Ziggah, Y.Y. Annan, R.F. 2016. Determination of 3D transformation parameters for the Ghana Geodetic reference network using ordinary least squares and total least squares technology. International Journal of Geomatics and Geoscience, 7(3), pp. 245-261.

Mercan, H. Akyilmaz, O. Aydin, C. 2018. Solution of the weighted symmetric similarity transformations based on quaternions. J Geod , 92(10), pp. 1113-1130.

Neitzel, F. 2010. Generalization of total least- squares on example of unweighted and weighted 2D similarity transformation. J Geod , 84(12), pp. 75-762.

Pope, A. 1972. Some pitfalls to be avoided in the iterative adjustment of nonlinear problems. Proceedingsof the38thAnnual Meeting of the American Society Photogrammetry, May 12-17, WashingtonDC, pp. 449-473.

Rapp, R. 1993. Geometric geodesy part II. Ohio: The Ohio State University, pp. 60-75.

Schaffrin, B. Felus, Y.A. 2008. On the multivariate total least squares approach to empirical coordinate transformations: three algorithm. J Geod , 82(6) , pp. 373-383.

Shi, Y. Xu, P. Liu, J. Shi, C. 2015. Alternative formulae for parameter estimation in partial error-in-variables model. J Geod , 89(1), pp. 13-16.

Tao, W., Hua, X., Yu, K., He, X. and Chen, X. 2018. An improved point-to-plane registration method for terrestrial laser scanning data. IEEE Access. 6, pp. 48062-48073.

Wang, L. Xu, G. 2016. Variance component estimation for partial errors-in-variables models. Stud. Geophys. Geod, 60(1), pp. 35-55.

Xu, P. Liu, J. Shi, C. 2012. Total least squares adjustment in partial errors-in-variables: algorithm and statisitical analysis. J Geod , 86(8), pp. 661-675.

Xu, P. Liu, J. 2014 Variance components in error-in-variables models: estimability, stability and bias Analysis. J Geod . 88(8), pp. 719-734.

Yang, Y. 1999. Robust estimation of geodetic datum transformation. J Geod , 73( 5),pp. 268-274.

Yan, S., and Fan, J. 2001. The solution set of the mixed LS-TLS problem. International Journal of Computer Mathematics, 77(4), pp. 545-561.

You, R. Hwang, H. 2006. Coordinate transformation between two geodetic datums of taiwan by least-squares collocation. Journal of Surveying Engineering, 132(2), pp. 64-70.

Yang, J. Wang, Y. Wang, Q. Tao, Y. 2017. Solution for rank-defect EIV model based on TLS estimation,” Survey review, 49(354) , pp. 171-175.

Zeng, W. Liu, J. Yao, Y. 2015. On partial error-in-variables models with inequality constraints of parameters and variables. J Geod , 89(2), pp. 111-119.

Zhao, J. Gui, Q. 2017. Outlier detection in partial errors-in-variables model. Boletim de Ciências Geodésicas, 23(1), pp. 1-20.

Notes

Author notes

AUTHOR´S CONTRIBUTION Wuyong Tao derived the algorithm and wrote the paper; Xianghong Hua performed the experiments and analyzed the results; Shaoquan Feng checked the paper.

Figure 1:
The ten points in the source coordinate system (a) and target coordinate system (b).

Table 1:
The calculated results of the two methods

Figure 2:
The standard deviations calculated by the three algorithms for 1000 simulations

Figure 3:
Computation time under different number of common points

Table 2:
The estimated variance components and the computation time

Figure 4:
(a) the original points and the registration result, and (b) the computation time under different sampling rates