Point-wise weighted solution for the similarity transformation parameters

In the paper, we consider solutions to the similarity transformation under different adjustment scenarios. We consider asymmetric cases, i.e., random errors included either in a target system (Gauss–Markov model, the most common case) or in a source system and a symmetric case when both systems are considered to be contaminated by random errors (Gauss–Helmert or errors-in-variables model). The presented solutions are based on Procrustes approach and all cases use point-wise weighting scheme. The asymmetric cases have closed-form solutions but the symmetric one, in the most general case considered herein, does not. In the latter case, weight matrices for a source and target systems are combined into one equivalent weight matrix which depends on a scale factor. This case is solved through a simple iterative process but there are some special instances which lead to closed-form solutions. The paper discusses the behavior of transformation parameters (rotations, translations, and a scale) under different adjustment scenarios.


Introduction
Similarity transformation known also as Helmert transformation or 7-parameter transformation (4 parameters in 2D) is among the most frequently used operation in geodesy, surveying, photogrammetry, and GIS. Coordinates of source and target systems are related by a single-scale factor, three rotation angles (one rotation angle in 2D), and three translations (two translations in 2D). Usually, the problem of finding transformation parameters, in the least squares sense, is carried out by means of linearization and iterative refinement of successive approximations. In some instances, this nonlinear problem may be simplified and solved as a linear one (at least within the Gauss-Markov framework), e.g., 2D case when it may be solved through a simple substitution or 3D case (Bursa-Wolf or Molodensky-Badekas transformations) when it is simplified by the use of small rotation angles and assumption on a scale factor not much different from unity (Deakin 2006;Li et al. 2012). In this contribution, we use the Procrustes approach to solve for similarity transformation parameters in the least squares framework under Gauss-Markov and Gauss-Helmert (errors-in-variables) models. The Procrustes analysis matches one configuration of points to another configuration by using only a central (isotropic) dilation (a scale factor), rotations, and translations, i.e., the same operations as in its geodetic counterpart. The term Procrustes analysis itself first appeared in the paper by Hurley and Cattell (1962). To the best of the author's knowledge, the first closed-form solution to the "similarity transformation" comes from Procrustean framework, from the field of Psychometrics (Schönemann 1966;Schönemann and Carroll 1970). The problem was solved as a matrix fitting algorithm (unweighted case) using the multivariate least squares framework under the Gauss-Markov model. The Procrustes analysis becomes more and more popular within the geodetic community due to Crosilla (2003) and Awange and Grafarend (2005). This contribution gives a complete solution for the similarity transformation parameters accessible within the Procrustean framework with a row-wise weighting scheme, i.e., all possibilities of error-contamination. In asymmetric cases, the error-contamination concerns either a source system or a target system, on the other hand, in symmetric case, both systems are considered to be the subject of random errors. Within asymmetric cases framework, the most often approach assumes a target system to be the source of random errors (Soler 1998;Shen et al. 2006;Sjöberg 2013). On the other hand, a less crowded path is to consider a source system to be erroneous. This instance is demonstrated in Li et al. (2013a) for the affine transformation model. Here, we also cover this case for the similarity transformation. The most general model, i.e., symmetric one, the term attributed to Teunissen (1988), solved in least squares framework and very often called a total least squares solution due to Golub and van Loan (1980) has been attracting attention for the last decades in the geodetic community. The model and its solution has been applied to geodetic problems (mainly coordinate transformations) under different assumptions by Teunissen (1988), Schaffrin and Wieser (2008), Felus (2008a, 2008b), Mercan et al. (2018), Chang (2016), Neitzel (2010, and Li et al. (2013b) to mention only a few. Interesting studies on error analysis in the similarity transformation may be found in Chang et al. (2017aChang et al. ( , 2017b. Procrustes-based approach presented herein yields attractive closed-form solutions in asymmetric cases and in certain instances of symmetric one. The derived algorithms use a modification to ensure the resulting orthogonal matrix to be a rotation one. This modification is adapted from Markley (1988), see also Umeyama (1991) and Sjöberg (2013) or for a broader exposition of the problem Myronenko and Song (2009). Markley and Mortari (1999) present a review of rotation matrix recovery algorithms in Wahba's problem (Wahba 1965) or equivalently in orthogonal Procrustes problem (Schönemann 1966) with constraint on resulting orthogonal matrix to have determinant equal to unity (a rotation matrix not a reflection one). The solution for the symmetric case does not require linearization but it does require iteration. In the derivation process, the weight matrices in source and target systems are combined into one equivalent weight matrix-dependent on a scale factor (unfortunately). There are some special cases when it has closed-form solution. It also extends the discussion on the behavior of a rotation matrix under different adjustment scenarios.

Asymmetric cases
The first, as far as the author's knowledge goes, closed-form solution of Helmert's transformation parameters in unweighted case was presented by Schönemann and Carroll (1970) in the field of Psychometrics (the reason it remained unknown to geodetic society for some time). Here, we present a point-wise weighted solution of the Helmert's transformation parameters when random errors are included either in a target system (1a) or in a source system (1b): where P is a matrix of coordinates in a target system (m × 2 or m × 3), Q is a matrix of coordinates in a source system (m × 2 or m × 3), E P and E Q are matrices of disturbances (m × 2 or m × 3), R is a proper orthogonal matrix (2 × 2 or 3 × 3), i.e., R T R = RR T = I and det(R) = 1, t is a translation vector (2 × 1 or 3 × 1), s is a scale factor (scalar), u is a vector of ones (m × 1), m is the number of corresponding points for both systems.
Error-affected P (target system) model In this section, a target system is considered to be contaminated with random errors (the most usual case used in geodetic/ surveying practice). As a loss function, a squared weighted Frobenius norm is used, i.e.: where W P (m × m) is a point-wise weight matrix (a positivedefinite matrix). Inserting E P , obtained from (1a), into (2) one obtains: Using properties of trace operator and taking into account that R T R = RR T = I, the above relation may be rewritten in the following, simpler form: where α = u T W P u; which is minimized subject to the constraints By combining (4) with (5) in a Lagrange function, one obtains the following expression: where Λ and λ are Lagrange multipliers, a matrix and a scalar, respectively. Taking partial derivatives of (6) with respect to the transformation parameters R, t, s and equating them to zero the following formulas emerge (suitable expressions for derivatives of traces and determinants may be found in, e.g., Petersen and Pedersen 2012): Equation (7a) gives: where S = (Λ T + Λ) is a symmetric matrix and det(R) was set to unity in (7a). Eq. (7b) yields: Inserting (9) to (7c) and rearranging terms, the scale factor may be expressed as: By inserting expression for t into (8), the following formula emerges: where the right hand side of (11) may be perceived as a polar decomposition of sQ T W P I − 1 α uu T W P À Á P. This formula includes unknown scale factor s (in fact on both sides) but this does not pose a problem since uniform scaling does not change the orthogonal polar factor R (e.g., Gander 1989Gander , 1990; hence, it is enough to find the orthogonal polar factor (restricted to det(R) = 1) of Q T W P I − 1 α uu T W P À Á P in order to obtain the solution. The polar decomposition of a matrix may be obtained, e.g., by means of singular value decomposition (SVD) which is in fact the basic tool in errors-in-variables literature. A firm review of the polar decomposition algorithms and its properties may be found in, e.g., Gander (1989Gander ( , 1990, Markley and Mortari (1999), Higham (1986), and Shoemake and Duff (1992). A contribution by Higham and Noferini (2015) is especially interesting in this respect because it is dedicated to computing polar decomposition of 3 × 3 matrices (in references therein one may also find an effective algorithm for 2 × 2 matrices). Below, we summarize the derivation in a step-by-step computational algorithm.
Algorithm I Weighted solution of Helmert transformation parameters in error-affected target system model 1. Prepare coordinate matrices P, Q, and a weight matrix for a target system W P 2. Compute the weighted centering matrix C = W P − W P uu T W P /α, where α = u T W P u and u T = [1 1 … 1] 3. Compute singular value decomposition (SVD) of Q T CP=UΣV T and restore the rotation matrix as R = Udiag[1 1 det(U)det(V)]V T see, e.g., Markley (1988), Umeyama (1991) 4. Compute the scale factor as s ¼ For W P = I, one obtains the solution presented in Schönemann and Carroll (1970).
Error-affected Q (source system) model A corresponding derivation for the error-affected source system (Q) model will be presented in a similar manner. The general form of the cost function reads: and takes the explicit form (after inserting the expression for E Q obtained from (1b) and some manipulations): where W Q (m × m) is a point-wise weight matrix. The above expression may be rewritten in a form (by using properties of a trace operator and taking into account R T R = RR T = I): where α = u T W Q u. The optimization process is subject to the same constraints as in the previously considered problem; namely; Combining (14) and (15) Lagrangian reads: where Λ and λ are Lagrange multipliers. Taking partial derivatives with respect to R, t, s, and equating them to zero one obtains the following formulas: Equation (17) gives: where S = (Λ T + Λ) is a symmetric matrix and |R| is set to unity in (17). Equation (18) yields: Inserting (21) into (19) and rearranging terms the scale factor may be expressed as: By inserting expression for t into (20) the following relation emerges: where the right hand side of (23) is again a polar decomposition of 1 s Q T W Q I− 1 α uu T W Q À Á P. Previously, this formula includes unknown scale factor s but uniform scaling does not affect the orthogonal polar factor R (e.g., Gander 1989Gander , 1990; hence, it is enough to find the orthogonal polar factor (rotation matrix) of Q T W Q I− 1 α uu T W Q À Á P. Algorithm II Weighted solution of Helmert transformation parameters in error-affected source system model 1. Prepare coordinate matrices P, Q, and a weight matrix for a source system W Q 2. Compute the weighted centering matrix C = W Q − W Q uu T W Q /α, where α = u T W Q u and u T = [1 1 … 1] 3. Compute singular value decomposition (SVD) of Q T CP=UΣV T and restore the rotation matrix as R = Udiag[1 1 det(U)det(V)]V T 4. Compute the scale factor as s ¼

Symmetric case
Error-affected Q (source system) and P (target system) model Now, a solution for Helmert transformation parameters under errors-in-variables model will be considered in a Procrustean framework. The functional relationship may be expressed as: In this optimization problem, we seek transformation parameters R, t, and s that minimize the following cost function: subject to (24). where W P and W Q are positive-definite weight matrices that weigh points (point-wise weighting schema) in a target and in a source system, respectively.
Combination of these two expressions gives the Lagrange function of the form: where Λ is a Lagrange multipliers matrix. Partial derivatives of (26) with respect to disturbance matrices E P and E Q read: Substituting these results into (24) and using the orthogonality relation R T R = RR T = I, one obtains: and this further gives the expression for Λ: Plugging (27), (28), and then (30) into (25), one obtains the cost function expressed explicitly as a function of transformation parameters R, t, s: where W s ð Þ ¼ W −1 P þ s 2 W −1 Q −1 will be called an equivalent weight matrix.
Combining the final form of the loss function (31) with the orthogonality constraint RR T = R T R = I and |R| = 1, we obtain the Lagrange function to be optimized with respect to the transformation parameters: where Λ R and λ are Lagrange multipliers. Expanding the expression under the first trace operator and using its properties, we obtain the following formula: where α = u T W(s)u Partial derivatives with respect to the transformation parameters R, t, and s are given respectively: From (34) we have: where S ¼ Λ T R þΛ R À Á is a symmetric matrix and |R| in (34) was set to unity and from (40), expression for t may be given by: By inserting the above expression for t into (42) we obtain: As in all previous cases, the solution for a rotation matrix has a form of the polar decomposition of sQ T W s ð Þ I− 1 α uu T W s ð Þ À Á P. This expression, in fact, is dependent on scale factor as a multiplier and through the equivalent weight matrix being a function of scale factor too. Relying on the property of the polar decomposition that the uniform scaling does not affect the orthogonal polar factor R, the left hand side of (39) may be reduced to Q T W s ð Þ I− 1 α uu T W s ð Þ À Á P. This expression is still dependent on the unknown scale factor s, thus the entire solution for transformation parameters will require an iterative approach.
An expression for a scale factor is obtained by collecting and rearranging terms in (36) what gives a quadratic equation: where Dependence of the equivalent weight matrix W(s) on the scale factor makes the solution a bit troublesome but fortunately, it does not require linearization but only iteration. Hence, the above formulas may be summarized and cast into the following iterative scheme.
Algorithm III Weighted solution of Helmert transformation parameters for the symmetric case 1. Prepare coordinate matrices P, Q, and weight matrices W P and W Q 2. Find an approximate solution for R, t, and s using, e.g., Algorithm I 3. Update s by solving (40 It is easily verifiable that in case of W P = W Q = W (equally weighted points in both systems, not necessarily W = I) the equivalent weight matrix takes the form: what gives the objective function, counterpart of (31), of the form: If the weight matrix W is taken to be a scalar multiple of an identity matrix (particular case of 44) i.e., W = kI, then the cost function is given by: The latter expression, being a particular case of (31) or (45), gives the same loss function as given in Chang (2015) but in a more compact and comfortable matrix form.
Also, in case of W P ¼ 1 k P W, W Q ¼ 1 k Q W, i.e., scalar multiples of the same weight matrix W (compare Chang 2016) the equivalent weight matrix takes the form: and the loss function will be given as: If W is taken as a diagonal matrix, the loss function (48) is equivalent to the one presented in Chang (2016). Hence, the solution proposed here absorbs solutions presented by Chang (2015Chang ( , 2016.

Discussion
Comparing final formulas for the transformation parameters (rotation, translation, and scale) in the asymmetric cases considered here, we notice that, in general, the two adjustment scenarios produce different results. If W P ≠ W Q , then all transformation parameters for both cases are different (compare formulas (9, 10, 11) to their counterparts (21,22,23)). On the other hand, if W P = W Q = W (not necessarily W = I) we obtain the same rotation parameter R (a rotation matrix) for these two cases but scale factors s and translations t (scaledependent) are different. It is easy to notice that this weighting scheme leads to the full equivalence when s = 1 (a rigid body transformation-a particular case of the similarity transformation) then all transformation parameters R, s, and t in the asymmetric cases would be the same. In the most general case, if W P ≠ W Q , the solution for transformation parameters in the symmetric case leads to an iterative procedure (without linearization) due to the nonlinear involvement of a scale factor in the equivalent weight matrix. In this instance, it is hard to infer anything about similarities or dissimilarities to the previously mentioned asymmetric cases, but there are instances when the solution may be presented in a closed-form, i.e., when a scale factor may be unraveled from the equivalent weight matrix. Among these cases, we may list, e.g., In all these cases, the cost function of type (3) is multiplied by some scalar function involving the scale factor and the equivalent weight matrix is independent of the scale factor. Hence, in these cases, one may expect the equivalence for the rotation parameter between symmetric (Gauss-Helmert) and asymmetric cases (Gauss-Markov) under the condition of the same weight matrix W. The same cannot be stated about a scale factor and translation parameter. Moreover, considering the same weighting scheme, i.e., equal weight matrices for both systems, then in the symmetric case, we will obtain the same rotation and translation parameters as in asymmetric cases for a rigid-body transformation, i.e., when s = 1.

Conclusions
We have presented closed-form and iterative solutions to the similarity transformation under various adjustment scenarios. In general, closed-form expressions, besides their elegance, may potentially give better understanding of how different factors affect the final solution. They may also serve as a reference for testing iterative methods that may turn out to be more efficient, e.g., with respect to time of execution. In this study, we have considered asymmetric cases when only one system is subject to random errors (Gauss-Markov model) and the symmetric case when both systems are contaminated by random errors (Gauss-Helmert model). To solve the stated problems, we employed multivariate Procrustes approach with point-wise weighting scheme (e.g., positional errors may be used to obtain weight matrices). The presented solutions are general and may be applied to both 2D and 3D similarity transformations without modifications. Both considered asymmetric cases lead to closed-form solutions, the symmetric one, in general, does not. The most general symmetric case considered in the paper is solved through a simple iterative procedure but there are also special instances when it has closed-form solutions. The presented solutions if not satisfactory due to the use of special weighting may be used as an initial estimate for more general weighting schemes.
Step-bystep algorithms have been given for the ease of implementation. Although all derivations for restoring a rotation parameter are presented in the polar decomposition form, to obtain the latter mentioned, another factorization has been used, i.e., singular value decomposition.