Extension of Procrustes algorithm
for 3D affine coordinate transformation

Béla ((Paláncz)^1),  Piroska Zaletnyi(k^2) and Joseph L. Awang(e^3)

1Department of Photogrammetry and Geoinformatics
Budapest University of Technology and Economy,H-1521,Hungary
e-mail:palancz@epito.bme.hu

 2Department of Geodesy and Surveying and
Research Group of Physical Geodesy and Geodynamics of the Hungarian Academy of Sciences
Budapest University of Technology and Economy,H-1521,Hungary
e-mail: zaletnyikp@hotmail.com

3Western Australian Centre for Geodesy
Department of Spatial Sciences,
Division of Resource and Environmental
Curtin University of Technology, Australia
e-mail: J.awange@curtin.edu.au

Abstract

Procrustes method is a very effective method to solve parameter estimation problem for Helmert coordinate transformation, because it  requires  neither initial start values nor iteration. That is why, it is a very tempting idea to extend it to 3D affine transformation, where scale factors are different in the 3 principal directions. In this  study,  it will be illustrated, that the direct extension is restricted for cases where the values of these 3 scale factors are just slightly different. We also demonstrate that this fact does not mean that it is unreasonable to use affine transformation model instead of Helmert one, because it can provide  improvement in the quality of the transformation even in these restricted cases. However, employing iteration for the scale parameters, the extension can be fully  achieved, even for strongly different values of  scale parameters. It goes without saying that in this way the computation time will increase considerable, but it is even competitive with those of other methods.

Introduction

Three-dimensional coordinate transformations play a central role in contemporary Euclidean point positioning. In precise positioning with global positioning system (GPS), coordinates given in the World Geodetic System 1984 (WGS84) often have to be transformed into local geodetic coordinate system.  The transformation between the traditional terrestrial coordinate system and the satellite observations derived network is a difficult task due to the heterogeneity of the data. 
3D affine transformation problem arises when the scale is not uniform within two or more configurations of points such that they vary within given dimensions. It means that beside the  parameters of the rotation matrix (3)  and the translation vector (3), in addition the 3 scales should be also computed.
Recently a comprehensive study  comparing  different symbolic - numeric computation techniques to estimate these 9 parameters  has been published, see Paláncz, Lewis, Zaletnyik and Awange (2008). Although this work did not considered the application of the Procrustes algorithm. 
The first paper trying to employ  Procrustes algorithm to the parameter estimation problem of the 3D affine model was written by Awange, Bae and Claessens (2008). However their suggested algorithm proved to be successful only in case of very mild anisotropy in scaling. In this work we shall give a more generally valid solution of the problem.

Procrustes problems

The Procrustean approach solving least square problem appeared already nearly 25 years ago, see Gower (1984)  is very fast method, which needs neither initial starting value nor iteration. Procrustes being a technique of matching one configuration into another and producing a measure of match, seeks the isotropic dilation and the rigid translation, reflection and rotation needed to best match one configuration to another, Cox and Cox (1994).
The Procrustes problem is concerned with fitting a configuration ℬ  into 𝒜  as close as possible. The simplest Procrustes case is one in which both configurations have the same dimensionality and the same number of points, which can be brought into a 1 - 1  correspondence by substantive considerations, Borg  and Groenen (1997). Let us consider the case where both  𝒜   and  ℬ  are of the  same dimension. The partial Procrustes problem is then formulated as

𝒜 = ℬ  T

The rotation  matrix T is then solved by measuring the distances between corresponding points in both configurations, square these values, and add them to obtain the sum of squares ‖ 𝒜 -ℬ  T‖ which is then minimized. One proceeds via Frobenius norm as  follows

‖ 𝒜 -ℬ  T‖=(Sqrt[tr(((𝒜^T)- (T^T)((ℬ )^T))(𝒜 -ℬ  T))])-> (min_T)

where

(T^T) T = I

This partial Procrustes problem can be extended  for additional translation and scaling, beside the rotation, Awange and Grafarend (2005)

𝒜 = 𝓈 ℬ  T + 𝒟

Then the objective function to be minimized

(Subscript[‖ 𝒜 -𝓈(ℬ  T+𝒟)‖, W])=(Sqrt[tr(((𝒜^T)- 𝓈((T^T)((ℬ )^T)+𝒟))W(𝒜 -𝓈(ℬ  T+𝒟)))])-> (min_(T, s, 𝒟 ))

where W is a weight matrix, s is the scaling parameter and 𝒟 is translation vector. This Procrustes problem is called as general Procrustes problem.

These Procrustes problems arise in applications related to, e.g. rigid body movement and psychometrics, factor analysis, multivariate analysis, multidimensional scaling and global positioning system, see for example, Gower (1984) and Meridith  (1977). 
The solution of the general Procrustes problem provides the solution of the parameter estimation problem of the Helmert transformation with 7 parameter.

The general Procrustes method

The original Procrustes method with weighting can be applied to the parameter estimation problem of Helmert transformation with 7 parameter,  see Awange and Grafarend (2005).

The form of the general Procrustes algorithm for (Subscript[C, 7])(3,3) problem in case of N points,

(({{(Subscript[x, 1]), (Subscript[y, 1]), (Subscript[z, 1])}, {., ., .}, {., ., .}, {(Subscript[x, N]), (Subscript[y, N]), (Subscript[z, N])}}))=(({{(Subscript[X, 1]), (Subscript[Y, 1]), (Subscript[Z, 1])}, {., ., .}, {., ., .}, {(Subscript[X, N]), (Subscript[Y, N]), (Subscript[Z, N])}}))(R^T)s+(({{1}, {.}, {.}, {1}}))((Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0]))+E

where s is a scalar valued unknown for scale factor, R ∈(ℝ^(3×3)) unknown  orthonormal  matrix for the rotation matrix ,  and  E = (Subscript[e, i,j])  is the error matrix , i = 1..N, j = 1..3.

One may use a more compact form

xyz=XYZ   (R^T) s+1  ((Subscript[XYZ, 0])^T)+E

where 1= (({{1}, {.}, {.}, {.}, {1}}))  with length of N.  Let us express the error matrix explicitly,

E=xyz - XYZ   (R^T)s- 1 (((Subscript[XYZ, 0]) )^T)

then introducing W  positive definite weight matrix

W= diag((Subscript[w, 1]),...,(Subscript[w, N]))

Then the square of the Frobenius norm of the weighted error matrix is

((Subscript[(∥E∥), 𝓌])^2)= tr ((E^T) W E)

The problem is to minimize,

((Subscript[(∥E(s,(Subscript[XYZ, 0]),R)∥), w])^2)->min

under the condition

(R^T)R= (Subscript[I, 3])  and ∥R∥= +1

This means, that in case of (Subscript[C, 7])(3,3), one has to solve a constrained minimization problem, where

-  the scale factor, s corresponds to   s∈ℝ.
-  the translation vector,  ((Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])),corresponds to  (Subscript[XYZ, 0]) ∈(ℝ^3)
              -  the rotation matrix   R∈ (ℝ^(3×3))

The algorithm solving  this problem can be found in Awange and Grafarend (2005)  as it follows

1)  Compute the centralization matrix C

C =(Subscript[I, N])-(1/N)*1*(1^T)

where (Subscript[I, N]) is the unit matrix.

2) Now, we can compute the A matrix

A = (xyz^T) (C^T)  W C XYZ

3) Employing singular value decomposition (SVD)  of A,  we get the matrices U, V and diag((Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3])), where

A = U diag((Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3]))V

4) Then the rotation matrix can be computed

R= U (V^T)

5) As we have seen the rotation matrix can be expressed by the skew matrix as

R =((((Subscript[I, 3])-Ξ))^(-1))((Subscript[I, 3])+Ξ),

Consequently the skew matrix is

Ξ= (R-(Subscript[I, 3]))(R+(Subscript[I, 3]))

and then the parameters a, b and c can be computed from

Ξ=(({{0, -c, b}, {c, 0, -a}, {-b, a, 0}}))

7) The scale parameter can be expressed as

s=((tr((xyz^T)* (1^T)*W*  1 * XYZ *  (R^T)))/(tr((XYZ^T)*  (1^T)*  W*   1*  XYZ)))

where tr is the trace of a matrix.

8) While the translation vector satisfies the following normal equation

(1^T) * W*  1 *  (Subscript[XYZ, 0]) = (((xyz-XYZ  * (R^T) s))^T)*W * 1

therefore

(Subscript[XYZ, 0])=((((1^T)  W *1))^(-1)) (((xyz-XYZ  * (R^T) s))^T)W* 1

The implementation of this algorithm  in Mathematica is straightforward:

Procrustes[xyz_,XYZ_,W_]:=Module[
{n,unit,CC,A,U,V,Ξ,ℛ,s,T,Δ,Err},
n=Dimensions[xyz][[1]];
unit=Table[{1},{i,1,n}];
CC=IdentityMatrix[n]-(1/n)*unit.Transpose[unit];
A=Transpose[xyz].Transpose[CC].W.CC.XYZ;
{U,Ξ,V}=SingularValueDecomposition[SetPrecision[A,20]];
ℛ=U.Transpose[V];
s=Tr[Transpose[xyz].Transpose[CC].W.CC.XYZ.Transpose[ℛ]]/Tr[Transpose[XYZ].Transpose[CC].W.CC.XYZ];
T=Flatten[Inverse[Transpose[unit].W.unit]][[1]]*Transpose[xyz-XYZ.Transpose[ℛ] s].W.unit;
Δ=xyz-XYZ.Transpose[ℛ] s-unit.Transpose[T];
Err=Tr[Transpose[Δ].W.Δ];
{ℛ,s,T,Err}];

Definition of the 3D affine transformation problem

The (Subscript[C, 9])(3,3) 3D affine transformation is one possible generalization of the (Subscript[C, 7])(3,3) Helmert transformation, using three different scale ((Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3])) parameters instead of a single one. In this case the scale factors can be modeled by a diagonal matrix (S).

(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))=S R (({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))+(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}}))

where S =(({{(Subscript[s, 1]), 0, 0}, {0, (Subscript[s, 2]), 0}, {0, 0, (Subscript[s, 3])}}))- scale matrix, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]) - 3 translation parameters and R the rotation matrix

The rotation matrix can be expressed with the skew-symmetric matrix (Ξ) (see Awange and Grafarend (2003))

Ξ=(({{0, -c, b}, {c, 0, -a}, {-b, a, 0}}));

The rotation matrix is

R=((((Subscript[I, 3])-Ξ))^(-1))((Subscript[I, 3])+Ξ),

where (Subscript[I, 3]) is a 3×3 identity matrix.

(Subscript[I, 3])=IdentityMatrix[3];

R=Inverse[((Subscript[I, 3])-Ξ)] .((Subscript[I, 3])+Ξ)//Simplify; MatrixForm[R]

for which the next restriction is true:

R.Transpose[R]//Simplify//MatrixForm

Awange - Bae - Claessens's  solution  (ABC method)

Recently,  an extension of general Procrustes method for 3D affine 9 -parameter transformation  has been published, see Awange, Bae and Claessens (2008). The algorithm can be summarized as it follows:

1) Compute the center of gravity of the two systems

(Subscript[a, 0])=(1/N) xyz

(Subscript[b, 0])=(1/N) XYZ

2) Let us translate the systems  into the center of the coordinate system

xyzC = xyz-1 * (Subscript[a, 0])

XYZC = XYZ-1 * (Subscript[b, 0])

3) Compute matrix A

A=(XYZC^T) * xyzC

4) Compute the rotation matrix R via SVD decomposition of A,  getting  the matrices U, V and diag((Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3])), where

A = U diag((Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3]))V

and

R = U * (V^T)

5) First approximation of the scale matrix,

S=((((R^T)*(XYZC^T)*XYZC*R))^(-1))*(R^T)*(XYZC^T)*xyzC

only the diagonal elements will be considered

(Subscript[S, 0])=diag(S)

5) Iteration to improve the scale matrix

w=xyz-XYZ*R*(Subscript[S, 0])-1*((Subscript[a, 0])-(Subscript[b, 0]) R*(Subscript[S, 0]))

dS=diag(((((R^T)*(XYZ^T)*XYZ*R))^(-1))*(R^T)*(XYZ^T)*w)

(Subscript[S, 0])=(Subscript[S, 0])+dS

6) After the iteration the scale matrix is

S=diag((Subscript[S, 0]))

7) Then the translational vector

(Subscript[XYZ, 0])=(Subscript[a, 0])-(Subscript[b, 0]) R*S

The algorithm can be easily implemented in Mathematica,

ProcrustesABC[xyz_,XYZ_,𝒩_]:=Module[{xyzC,XYZC,U,V,Ξ,ℛ,a0,b0,𝒮,n,𝒮0,i,w,d𝒮,T,Err},
(*center of mass of the two systems*)
a0=Mean[xyz];
b0=Mean[XYZ];
(*translation*)
xyzC=Map[#-a0&,xyz];
XYZC=Map[#-b0&,XYZ];
(*singular value decomposition*)
{U,Ξ,V}=SingularValueDecomposition[SetPrecision[Transpose[XYZC].xyzC,20]]; 
(*rotation matrix*)
ℛ=U.Transpose[V];
(*first approximation of the scale matrix*)
𝒮=Inverse[Transpose[ℛ].Transpose[XYZC].XYZC.ℛ].Transpose[ℛ].Transpose[XYZC].xyzC;
(*number of stations*)
n=First[Dimensions[xyz]];
(*initial diagonalized scale matrix*)
𝒮0=DiagonalMatrix[Diagonal[𝒮]];
(*iteration loop for improving the scale matrix*)
Do[
w=xyz-XYZ.ℛ.𝒮0-Table[(a0-b0.ℛ.𝒮0),{n}];
d𝒮=Inverse[Transpose[ℛ].Transpose[XYZ].XYZ.ℛ].Transpose[ℛ].Transpose[XYZ].w;
d𝒮=DiagonalMatrix[Diagonal[d𝒮]];
𝒮=𝒮0+d𝒮;
𝒮0=𝒮,
{i,1,𝒩}];
(* scale matrix after m iterations*)
𝒮=DiagonalMatrix[Diagonal[𝒮0]];
(*translation parameters*)
T=a0-b0.ℛ.𝒮;
(*error of the data fitting*)
Err=Tr[Transpose[w].w];
(*output*)
{Err,ℛ,𝒮,T}];

Paláncz - Zaletnyik's solution (PZ method)

The ABC method is very fast, and the computation effort for iteration to  improve the S scale matrix is negligible.  The main problem with  ABC  method is that after shifting the two systems into the origin of the coordinate system, the rotation matrix is computed via SVD. It means that the physical rotation as well as the deformation (distortion)  caused by different scaling in the different directions of  the 3 different principal axes  will be also involved in the computation of the rotation matrix. This can be approximately allowed only when these scale factors do not differ from each other considerably. In addition, the necessary condition for getting optimal  scale matrix is not restricted for diagonal matrix, and the off -diagonal elements are simply deleted in every iteration step, which will not ensure real global minimum for the trace of the error matrix.

The PZ method can cure these problems, but there is a price for it! This method start with a initially guessed diagonal S scale matrix, and first using the inverse of this matrix to eliminate  the deformation caused by scaling, and only after the elimination of this "distortion" will be  SVD applied to compute the rotation matrix itself. It goes without saying that in this way, the computation of the rotation matrix should be repeated in iterative way in order to decrease the error. That is why that this method is precise and correct for any ratios of scale parameters, but takes  longer time to achieve the results, than in case of the ABC method.

The algorithm can be summarized as it follows:

1) Guess the elements of the diagonal scale matrix, S

S =(({{(Subscript[s, 1]), 0, 0}, {0, (Subscript[s, 2]), 0}, {0, 0, (Subscript[s, 3])}}))

The function for computing the error in case of a given scale matrix S is the following

2) Compute the center of gravity of the two systems

(Subscript[a, 0])=(1/N) xyz

(Subscript[b, 0])=(1/N) XYZ

3) Let us translate the systems  into the origo

xyzC = xyz-1 * (Subscript[a, 0])

XYZC = XYZ-1 * (Subscript[b, 0])

4) Eliminate rotation (distortion) caused by scaling with different values in different directions in the computation of  matrix A

A=(XYZC^T) * xyzC* (S^(-1))

5) Compute the rotation matrix R via SVD decomposition of A,  getting  the matrices U, V and diag((Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3])), where

A = U diag((Subscript[λ, 1]),(Subscript[λ, 2]),(Subscript[λ, 3]))V

and

R = U * (V^T)

6) Then the translational vector

(Subscript[XYZ, 0])=(Subscript[a, 0])-(Subscript[b, 0]) R*S

7) Compute the error matrix

E=xyz - XYZ   (R^T).S- 1 (((Subscript[XYZ, 0]) )^T)

then introducing W  positive definite weight matrix

W= diag((Subscript[w, 1]),...,(Subscript[w, N]))

The square of the Frobenius norm of the weighted error matrix is

((Subscript[(∥E∥), 𝓌])^2)= tr ((E^T) W E)

If this error is to high, then modify the scale matrix S, and repeat  the computation. One iteration step, the computation of the objective function, is implemented in Mathematica as

ProcrustesPZ[xyz_,XYZ_,S_]:=Module[{xyzC,XYZC,U,V,Ξ,ℛ,a0,b0,T,n,w,Err},
(*center of mass of the two systems*)
a0=Mean[xyz];
b0=Mean[XYZ];
(*translation*)
xyzC=Map[#-a0&,xyz];
XYZC=Map[#-b0&,XYZ];
(*singular value decomposition*)
{U,Ξ,V}=SingularValueDecomposition[SetPrecision[Transpose[XYZC].xyzC.Inverse[S],20]]; 
(*rotation matrix*)
ℛ=U.Transpose[V];
(*translation parameters*)
T=a0-b0.ℛ.S;
(*number of stations*)
n=First[Dimensions[xyz]];
(*error of the data fitting*)
w=xyz-XYZ.ℛ.S-Table[T,{n}];
Err=Tr[Transpose[w].w];
(*output*)
{Err,ℛ,T}];

The objective function, the trace of the error matrix  is a function of the scale matrix, which should be minimized,

F[s1_?NumberQ,s2_?NumberQ,s3_?NumberQ]:=First[ProcrustesPZ[xyz,XYZ,W,(({{s1, 0, 0}, {0, s2, 0}, {0, 0, s3}}))]];

Minimization can be carried out with different methods built in Mathematica, but the most effective is the Newton method.

Although this is a local method and needs a good guess for the starting S matrix.  In case of the 7 parameter similarity transformation  (Helmert) the scale value (s) can be estimated by dividing the sum of length in both systems from the centre of gravity, see Albertz and Kreiling (1975).

The center of the gravity of the systems are,

(Subscript[a, 0])=(({{(Subscript[x, s])}, {(Subscript[y, s])}, {(Subscript[z, s])}}))     and    (Subscript[b, 0])=(({{(Subscript[X, s])}, {(Subscript[Y, s])}, {(Subscript[Z, s])}}))

Then the estimated scale parameter according to Albertz and Kreimlig in the Helmert similarity transformation

(Subscript[s, priori])=(((∑_(i=1)^N)(Sqrt[((((Subscript[x, i])-(Subscript[x, s])))^2)+((((Subscript[y, i])-(Subscript[y, s])))^2)+((((Subscript[z, i])-(Subscript[z, s])))^2)]))/((∑_(i=1)^N)(Sqrt[((((Subscript[X, i])-(Subscript[X, s])))^2)+((((Subscript[Y, i])-(Subscript[Y, s])))^2)+((((Subscript[Z, i])-(Subscript[Z, s])))^2)])))

In case of the 9 parameter affine transformation where 3 different scale values ((Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3]) are applied according to the 3 coordinate axes, a good approach for the scale parameters can be given  modifying the Albertz-Kreiling expression. Instead of the quotient of the two lengths in the centre of gravity system we can use the quotients of the sum of the lengths in the corresponding coordinate axes directions.

The estimated scale parameters according to the modified Albertz-Kreimlig expression for the 9 parameter transformation, see Zaletnyik (2008)

(Subscript[s, 1,priori])=(((∑_(i=1)^N)(Sqrt[((((Subscript[x, i])-(Subscript[x, s])))^2)]))/((∑_(i=1)^N)(Sqrt[((((Subscript[X, i])-(Subscript[X, s])))^2)]))),	(Subscript[s, 2,priori])=(((∑_(i=1)^N)(Sqrt[((((Subscript[y, i])-(Subscript[y, s])))^2)]))/((∑_(i=1)^N)(Sqrt[((((Subscript[Y, i])-(Subscript[Y, s])))^2)]))),	(Subscript[s, 3,priori])=(((∑_(i=1)^N)(Sqrt[((((Subscript[z, i])-(Subscript[z, s])))^2)]))/((∑_(i=1)^N)(Sqrt[((((Subscript[Z, i])-(Subscript[Z, s])))^2)]))),

Application to a network with mild anisotropy

Sometimes it is useful to introduce the inverse of the scale matrix  (Subscript[σ, i]) =(1/(Subscript[s, i]))and Ω = (S^(-1))

Ω=(({{(Subscript[σ, 1]), 0, 0}, {0, (Subscript[σ, 2]), 0}, {0, 0, (Subscript[σ, 3])}}));

For the determination of the 9 parameters of the transformation (a, b, c, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]), (Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3])) we need N ≥3 non-collinear points with known coordinates in both coordinate systems. In further,  instead of the scale parameters ((Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3])), we will use ((Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])) to get more simple equations. In case of N stations, we have 3 N nonlinear equations, a system of multivariable polynomial to solve.

First we consider the transformation problem of the Australian Geodetic Datum (AGD 84) to Geocentric Datum Australia (GDA 94), see Awange, Bae and Claessens (2008)
Now we have   N = 82 points in both systems, expressing the rotation matrix with the skew-symmetric matrix and using the inverse of the scale matrix (Ω) the nonlinear  system to be solved for ((Subscript[f, i]) =0 ) determining the 9 parameters leads to:

𝒩=82;

EQm=Table[(({{(Subscript[f, 1+3(j-1)])}, {(Subscript[f, 2+3(j-1)])}, {(Subscript[f, 3+3(j-1)])}}))=(((Subscript[I, 3])-Ξ).Ω.(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))- ((Subscript[I, 3])+Ξ).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3])-Ξ).Ω.(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i->j//Expand,{j,1,𝒩}];

sys=Table[(Subscript[f, i]),{i,3 𝒩}];Short[sys,10]

The number of the equations is

Length[sys]

Direct Numerical Solution via Global Minimization

In this section we solve this  system of  246  equations for the unknown parameters a, b, c, s, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]). First the direct numerical solution via minimization of the least square objective will be presented.

The objective function can be constructed by considering the least square residuum as objective function, namely

obj=Apply[Plus,Map[((#)^2)&,sys]];Short[obj,10]

which should be minimized. 
The numerical values of the local and global coordinates are the following:

dataxyz=ReadList["G:\Awange\WGS84.txt",{Number,Number,Number}];

Short[dataxyz,10]

dataXYZ=ReadList["G:\Awange\AGD84.txt",{Number,Number,Number}];

Short[dataXYZ,10]

dataxyzS=Drop[dataxyz,0];

dataXYZS=Drop[dataXYZ,0];

nxyz=Length[dataxyzS]

dataxyzT=Transpose[dataxyzS];

numericalValuesxyz=Table[{(Subscript[x, i])->dataxyzT[[1,i]],(Subscript[y, i])->dataxyzT[[2,i]],(Subscript[z, i])->dataxyzT[[3,i]]},{i,nxyz}]//Flatten;

dataXYZT=Transpose[dataXYZS];

numericalValuesXYZ=Table[{(Subscript[X, i])->dataXYZT[[1,i]],(Subscript[Y, i])->dataXYZT[[2,i]],(Subscript[Z, i])->dataXYZT[[3,i]]},{i,nxyz}]//Flatten;

numericalValues=Join[numericalValuesxyz,numericalValuesXYZ];

To carry out global minimization via genetic algorithm, the built -in function NMinimize can be employed.

solN=NMinimize[obj/.numericalValues,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])}]//Timing

ρ=SetPrecision[solN[[2,2]],20]

The scales are

Map[1/#&,{(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])}/.ρ]

The rotation matrix

Transpose[R/.ρ]//MatrixForm

The residual error  of the solution in least square sense

ΔN=(Sqrt[solN[[2,1]]])

Symbolic preparation of the local coordinates

xyz=Table[{(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])},{i,nxyz}];Short[xyz,10]

Their numerical values

xyz=xyz/.numericalValues;

Similarly for the global coordinates

XYZ=Table[{(Subscript[X, i]),(Subscript[Y, i]),(Subscript[Z, i])},{i,nxyz}];Short[XYZ,10]

XYZ=XYZ/.numericalValues;

Application of the Helmert transformation with original Procrustes

W=IdentityMatrix[𝒩];

{ℛ,s,T,Err}=Procrustes[xyz,XYZ,W];//Timing

The scale value

SetPrecision[s,20]

The rotation matrix

SetPrecision[ℛ,20]//MatrixForm

The translation parameters

SetPrecision[T,20]//Flatten

The error

SetPrecision[Err,20]//Sqrt

Application of ABC method

{Err,ℛ,𝒮,T}=ProcrustesABC[xyz,XYZ,𝒩];//Timing

The error

Sqrt[Err]

The rotation

ℛ//MatrixForm

The scale matrix

SetPrecision[𝒮, 20] // MatrixForm

The translations

T//SetPrecision[#,20]&

Application of PZ method

First we should estimate the initial starting values for the scales.

(Subscript[x, s])=(((∑_(i=1)^𝒩)(Subscript[x, i]))/3)/.numericalValues;(Subscript[y, s])=(((∑_(i=1)^𝒩)(Subscript[y, i]))/3)/.numericalValues;(Subscript[z, s])=(((∑_(i=1)^𝒩)(Subscript[z, i]))/3)/.numericalValues;(Subscript[X, s])=(((∑_(i=1)^𝒩)(Subscript[X, i]))/3)/.numericalValues;(Subscript[Y, s])=(((∑_(i=1)^𝒩)(Subscript[Y, i]))/3)/.numericalValues;(Subscript[Z, s])=(((∑_(i=1)^𝒩)(Subscript[Z, i]))/3)/.numericalValues;

then

sxpriori=(((∑_(i=1)^𝒩)Abs[(Subscript[x, i])-(Subscript[x, s])])/((∑_(i=1)^𝒩)Abs[(Subscript[X, i])-(Subscript[X, s])]))/.numericalValues;sypriori=(((∑_(i=1)^𝒩)Abs[(Subscript[y, i])-(Subscript[y, s])])/((∑_(i=1)^𝒩)Abs[(Subscript[Y, i])-(Subscript[Y, s])]))/.numericalValues;szpriori=(((∑_(i=1)^𝒩)Abs[(Subscript[z, i])-(Subscript[z, s])])/((∑_(i=1)^𝒩)Abs[(Subscript[Z, i])-(Subscript[Z, s])]))/.numericalValues;

Map[ SetPrecision[#,20]&,{sxpriori,sypriori,szpriori}]

The objective function is

F[s1_?NumberQ,s2_?NumberQ,s3_?NumberQ]:=First[ProcrustesPZ[xyz,XYZ,(({{s1, 0, 0}, {0, s2, 0}, {0, 0, s3}}))]];

Using minimization in order to find the scales

Off[FindMinimum::"lstol"];

solNewton=FindMinimum[F[s1,s2,s3],{{s1,sxpriori},{s2,sypriori},{s3,szpriori}}, Method->"Newton",WorkingPrecision->20]//Timing

With the known scales we can compute the other parameters and the error as well

{Err,ℛ,T}=ProcrustesPZ[xyz,XYZ,(({{s1, 0, 0}, {0, s2, 0}, {0, 0, s3}}))/.solNewton[[2,2]]];

Sqrt[Err]

ℛ//MatrixForm

SetPrecision[T,20]

Table 1. Results of the different methods in case of network with mild anisotropy

({{Method, Time,sec
, Error, ℛ =(R^T), S
s1, s2, s3, T
X0, Y0, Z0}, {Procrustes
for Helmert, 0.062, 6.866908, (({{0.99999999999, -1.7957307265579*((10)^(-6)), 1.8616292227837*((10)^(-6))}, {1.795731807008 *((10)^(-6)), 0.99999999999, -5.80377489174479*((10)^(-7))}, {-1.8616281805787*((10)^(-6)), 5.80380832159346*((10)^(-7)), 0.99999999999}})), 1.00000368981, -115.837707
-48.3732068
144.7595508}, {ABC, 0.282, 6.78838, (({{0.99999999999, 1.79573180951*(10^(-6)), -1.86162818141*(10^(-6))}, {-1.79573072906*((10)^(-6)), 0.99999999999, 5.80380831238*((10)^(-7))}, {1.86162922361*((10)^(-6)), -5.80377488254*(10^(-7)), 0.99999999999}})), ({{1.00000396085}, {1.00000354605}, {1.00000345982}}), ({{-115.0618939996}, {-47.67569050658}, {144.09550540243}})}, {PZ, 0.687, 6.64167, (({{0.99999999999, 2.119134673260*((10)^(-6)), -2.11439131407*(10^(-6))}, {-2.11913384350*((10)^(-6)), 0.99999999999, 3.92440675404*((10)^(-7))}, {2.11439214570*((10)^(-6)), -3.92432894724*(10^(-7)), 0.99999999999}})), 1.00000416838
1.00000310067
1.00000346938, -112.1689410633
-44.04637841973
144.31147376889}, {Global optimization
Genetic algorithm, 3.562, 6.64165, (({{0.99999999999, 2.11915502430*(10^(-6)), -2.11439079495*(10^(-6))}, {-2.11915419458*(10^(-6)), 0.99999999999, 3.92418312412*(10^(-7))}, {2.11439162654*(10^(-6)), -3.9241383169*(10^(-7)), 0.99999999999}})), 1.00000416842
1.00000310066
1.00000346943, -112.168784046180
-44.046246596546
144.31172203392}})

Application to network with strong anisotropy

For this test an artificial network is generated from Hungarian Datum points. Now, the angles of rotations are not small values and the values of the 3 scale parameters are considerably different from each other. We have N = 81 points in both systems, expressing the rotation matrix with the skew-symmetric matrix and using the inverse of the scale matrix (Ω) the nonlinear  system to be solved for ((Subscript[f, i]) =0 ) determining the 9 parameters leads to:

𝒩=81;

EQm=Table[(({{(Subscript[f, 1+3(j-1)])}, {(Subscript[f, 2+3(j-1)])}, {(Subscript[f, 3+3(j-1)])}}))=(((Subscript[I, 3])-Ξ).Ω.(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))- ((Subscript[I, 3])+Ξ).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3])-Ξ).Ω.(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i->j//Expand,{j,1,𝒩}];

sys=Table[(Subscript[f, i]),{i,3 𝒩}];Short[sys,10]

The number of the equations is

Length[sys]

Direct Numerical Solution via Global Minimization

In this section we solve this  system of 243 equations for the unknown parameters a, b, c, s, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]). First the direct numerical solution via minimization of the least square objective will be presented.

The objective function can be constructed by considering the least square residuum as objective function, namely

obj=Apply[Plus,Map[((#)^2)&,sys]];Short[obj,10]

which should be minimized. 
The numerical values of the local and global coordinates are the following:

dataxyz=ReadList["G:\\testdataPIRI\\testxyz.txt",{Number,Number,Number}];

Short[dataxyz,10]

dataXYZ=ReadList["G:\\testdataPIRI\\ETR589XYZGlobal.txt",{Number,Number,Number}];

Short[dataXYZ,10]

dataxyzS=Drop[dataxyz,0];

dataXYZS=Drop[dataXYZ,0];

nxyz=Length[dataxyzS]

dataxyzT=Transpose[dataxyzS];

numericalValuesxyz=Table[{(Subscript[x, i])->dataxyzT[[1,i]],(Subscript[y, i])->dataxyzT[[2,i]],(Subscript[z, i])->dataxyzT[[3,i]]},{i,nxyz}]//Flatten;

dataXYZT=Transpose[dataXYZS];

numericalValuesXYZ=Table[{(Subscript[X, i])->dataXYZT[[1,i]],(Subscript[Y, i])->dataXYZT[[2,i]],(Subscript[Z, i])->dataXYZT[[3,i]]},{i,nxyz}]//Flatten;

numericalValues=Join[numericalValuesxyz,numericalValuesXYZ];

To carry out global minimization, the built -in function NMinimize can be employed. Now we consider here and further, the first three points,

solN=NMinimize[obj/.numericalValues,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3]),(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])}]//Timing

ρ=SetPrecision[solN[[2,2]],20]

The scales are

Map[1/#&,{(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])}/.ρ]

The rotation matrix

Transpose[R/.ρ]//MatrixForm

The residual error  of the solution in least square sense

ΔN=(Sqrt[solN[[2,1]]])

Symbolic preparation of the local coordinates

xyz=Table[{(Subscript[x, i]),(Subscript[y, i]),(Subscript[z, i])},{i,nxyz}];Short[xyz,10]

Their numerical values

xyz=xyz/.numericalValues;

Similarly for the global coordinates

XYZ=Table[{(Subscript[X, i]),(Subscript[Y, i]),(Subscript[Z, i])},{i,nxyz}];Short[XYZ,10]

XYZ=XYZ/.numericalValues;

Application of the Helmert transformation with original Procrustes

W=IdentityMatrix[𝒩];

{ℛ,s,T,Err}=Procrustes[xyz,XYZ,W];//Timing

The e scale value

SetPrecision[s,20]

The rotation matrix

SetPrecision[ℛ,20]//MatrixForm

The translation parameters

SetPrecision[T,20]//Flatten

The error

SetPrecision[Err,20]//Sqrt

Application of ABC method

{Err,ℛ,𝒮,T}=ProcrustesABC[xyz,XYZ,𝒩];//Timing

Err

Sqrt[Err]

ℛ//MatrixForm

SetPrecision[𝒮,20]//MatrixForm

T

Application of PZ method

First we should estimate the initial starting values for the scales.

(Subscript[x, s])=(((∑_(i=1)^𝒩)(Subscript[x, i]))/3)/.numericalValues;(Subscript[y, s])=(((∑_(i=1)^𝒩)(Subscript[y, i]))/3)/.numericalValues;(Subscript[z, s])=(((∑_(i=1)^𝒩)(Subscript[z, i]))/3)/.numericalValues;(Subscript[X, s])=(((∑_(i=1)^𝒩)(Subscript[X, i]))/3)/.numericalValues;(Subscript[Y, s])=(((∑_(i=1)^𝒩)(Subscript[Y, i]))/3)/.numericalValues;(Subscript[Z, s])=(((∑_(i=1)^𝒩)(Subscript[Z, i]))/3)/.numericalValues;

then

sxpriori=(((∑_(i=1)^𝒩)Abs[(Subscript[x, i])-(Subscript[x, s])])/((∑_(i=1)^𝒩)Abs[(Subscript[X, i])-(Subscript[X, s])]))/.numericalValues;sypriori=(((∑_(i=1)^𝒩)Abs[(Subscript[y, i])-(Subscript[y, s])])/((∑_(i=1)^𝒩)Abs[(Subscript[Y, i])-(Subscript[Y, s])]))/.numericalValues;szpriori=(((∑_(i=1)^𝒩)Abs[(Subscript[z, i])-(Subscript[z, s])])/((∑_(i=1)^𝒩)Abs[(Subscript[Z, i])-(Subscript[Z, s])]))/.numericalValues;

Map[ SetPrecision[#,20]&,{sxpriori,sypriori,szpriori}]

The objective function is

First we should estimate the initial starting values for the scales.

F[s1_?NumberQ,s2_?NumberQ,s3_?NumberQ]:=First[ProcrustesPZ[xyz,XYZ,(({{s1, 0, 0}, {0, s2, 0}, {0, 0, s3}}))]];

Using minimization

Off[FindMinimum::"lstol"];

solNewton=FindMinimum[F[s1,s2,s3],{{s1,sxpriori},{s2,sypriori},{s3,szpriori}}, Method->"Newton",WorkingPrecision->20]//Timing

{Err,ℛ,T}=ProcrustesPZ[xyz,XYZ,(({{s1, 0, 0}, {0, s2, 0}, {0, 0, s3}}))/.solNewton[[2,2]]];

Err

Sqrt[Err]

ℛ//MatrixForm

SetPrecision[T,20]

Table 2. Results of the different methods in case of network with strong anisotropy

({{Method, Time,sec
, Error, ℛ =(R^T), S
s1, s2, s3, T
X0, Y0, Z0}, {Procrustes
for Helmert, 0.047, 117285.16, (({{0.901598216943688, 0.291216325804948., 0.319865138439}, {-0.355506948061613, 0.9201085955924943, 0.16436235030178.}, {-0.246445663552596, -0.261903081122529, 0.933097696393387}})), 1.2335164851, -4.82635678889*(10^6)
253120.161800.
6.355619148*(10^6)}, {ABC, 0.281, 93431.63, (({{0.90159821694369, -0.355506948061613, -0.24644566355258619}, {0.29121632580495, 0.92010859559249449, -0.26190308112252379}, {0.31986513843903, 0.16436235030177868, 0.933097696393391705}})), 1.099389091087
1.287089145364
1.136301910900, -4.0744058635024456*(10^6)
219870.67323442164
6.642437966975877`*(10^6)}, {PZ, 0.906, 1.3581, (({{0.94185020502006, -0.2598139255449, 0.213107755364746696}, {0.275416078439841, 0.960195515761870, -0.046589218147696507}, {-0.19252058342028, 0.1025733669318080, 0.97591727587730796466}})), 0.6200004792118
1.3000001184524
1.86999586512189, 1339.036188432248
-236.060661877505
152.7108742780983}, {Global optimization
Genetic algorithm, 4.625, 1.23799, (({{0.941850203067, -0.2598139259119, 0.21310776354879}, {0.275416078816, 0.96019551586879, -0.0465892137160}, {-0.19252059243, 0.10257336500117, 0.975917274301750}})), 0.62000047769
1.30000011779
1.86999590054, 1339.071922447
-236.0467167038
152.46054233354}})

Conclusion

In case of mild anisotropy of the network,  the ABC method using 3D affine transformation gives  better approximation than the general Procrustes method employing Helmert transformation model, while, the PZ method provides precise, geometrically correct solution. The ABC method is about 2 times faster than the PZ method, and the later is roughly 5 times faster than the global optimization method applying to 3D affine model, see Table 1.
In case of strong anisotropy Helmert transformation is not a proper model and  the ABC method using 3D affine transformation completely fails, while PZ method provides again correct solution.
The results of  the computations have been verified with the original MATLAB code of  ABC methods, too.

References

Albertz J. and Kreiling W. (1975)  Photogrammetric Guide, Herbert Wichmann Verlag, Karlsruhe, pp. 58-60.

Awange J.L, Grafarend E.W. (2003). : Closed form solution of the overdetermined nonlinear 7 parameter datum transformation. Allgemeine VermessungsNachrichten (AVN) 110, pp .130 - 148.

Awange and Grafarend (2005). Solving Algebraic Computational Problems in Geodesy and Geoinformatics, Springer.

Awange, Bae and Claessens (2008). Procrustean solution of the 9 -parameter transformation problem.

Borg I. and Groenen P. (1997)  Modern multidimensional scaling. Springer, New York

Cox T.F. and Cox M.A.(1994)   Multidimensional scaling. St.Edmundsbury Press, St. Edmunds, Suffolk

Gower J.C. (1984) Multivariate Analysis. Ordination, Multidimensional Scaling and Allied Topics. Handbook of Applicable Mathematics, VI.: Statistics (B)

Gower J.C. and Dijksterhuis G.B. (2004)  Procrustes Problems.  Oxford University Press.

Meridith W. (1977) On weighted Procrustes and Hyperplane fitting in factor analytic rotation. Psyhometrika, 42, pp.491-522

Paláncz B., Lewis R, Zaletnyik P. and Awange J.  (2008). Computational Study of 3D Affine Coordinate Transformation: Part I. 3-point Problem, MathSource /7090/

Paláncz B., Lewis R, Zaletnyik P. and Awange J.  (2008). Computational Study of 3D Affine Coordinate Transformation: Part II. N-point Problem, MathSource /7121/

Zaletnyik, P. (2008). Coordinate transform using Computer Algebra Systems (CAS) and Artificial Neural Netwotks (ANN). Ph.D.  Dissertation, Budapest University of Technology and Economy, http://www.agt.bme.hu/index_e.html