Computational Study of 3D Affine Coordinate Transformation  
Part I. 3-point Problem

Bela (Palancz^1),  Robert H. (Lewis^2),  Piroska Zaletnyi(k^3)and Joseph Awang(e^4)

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

(^2)Department of Mathematics
Fordham University, Bronx, NY 10458, USA
e-mail: rlewis @fordham.edu

 3Department 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

4Western 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

In case of  considerable nonlinearity  e.g.  in geodesy, photogrammetry, robotics, it is difficult to find proper initial values to solve the parameter estimation problem of 3D affine transformation with 9 parameters  via linearization and/or iteration. In this paper we develop a symbolic - numeric method to achieve the solution without initial guess. Our method employs explicit analytical expressions developed by the computer algebra technique  Dixon resultant as well as by reduced Grobner basis for solving 3 points problem. Criteria for the proper selection of the 3 points from the N ones, is also given. Numerical illustration is presented with real world geodetic coordinates representing Hungarian Datum.

Keywords:   9-parameter 3D affine transformation, solution of polynomial system, symbolic solution, Dixon resultant, early discovery factors, reduced Groebner basis,  homotopy, numerical Gröbner basis , global minimization, genetic algorithm.

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

Due to the distorsions  between the traditional terrestrial and the GPS derived networks, the 7-parameter similarity transformations in some cases may not offer satisfactory  precision. For example transforming GPS local coordinates to the local Hungarian system with global similarity transformation gives 0.5 meter maximal residuals, see Papp and Szucs (2005). To reduce the remaining  residuals after the transformation, other transformation models with more parameters can be used. The 9-parameter affine transformation is not only a logical extention but even a generalization of  the 7-parameter similarity transformation model.  This transformation is the modification of the Helmert C7(3) transformation, where 3 different scales are used in the corresponding coordinate axes instead of one scale factor, while in case of the 3 scale parameters are equal, we get back the similarity transformation model. The estimation of the model parameters was achieved by Spath (2004) using numerical minimization technique of the residium vector, by Papp and Szucs (2005) using linearized least squares method. Watson (2006) pointed out that the Gauss-Newton method or its variants can be easily implemented for the nine parameter problem using separation of variables and iteration with respect to the rotation parameters alone, while other parameters can be calculated via simple linear least square solution. The method he suggested is analogous to other methods for separated least squares problems, which goes back at least to Golub and Pereyra (1973). The 9-parameter affine transformation is also included in some coordinate-transformation software developed by the request of GPS users (see e.g. Mathes 2002 and Frohlich and Broker 2003). Here we should mention transformation models with more then 9 parameters. Wolfrum (1992) even added one more parameter to the previous nine, one (horizontal) direction  of maximal scale distortion.  Grafarend and Kampmann (1996) applied ten parameter conformal group for geodetic datum transformation employing maximum likelihood estimations for numerical estimation of the parameter values. There are other models with even more parameters, like polynomial transformations (Volgyesi et al 1996, Cai and Grafarend 2004) and models using artificial neural networks (Barsi  2001, Zaletnyik 2004).

In this paper we solve the 3D affine transformation problem in symbolic form via Dixon resultant employing enhanced Dixon KSY and the Dixon EDF algorithms. We also investigate  global numerical techniques,  like global minimization with genetic algorithm, global polynomial solver method employing numerical Gröbner basis using an efficient monomial ordering and eigensystem methods to extract numerical roots  as well as  linear homotopy continuation method. In case of  3 point problem, symbolic solution represented by explicit expressions for the parameters, proved to be very fast and robust, while in addition, it is infinite precision, comparing with the global numerical methods.

Our suggested method employs explicit, analytical expressions developed by computer algebra technique  via Dixon resultant for solving 3 point problem. Then the result of the 3 point problem can be used as a good initial guess for a Newton -Raphson- Krylov numerical  method to solve the N point problem in a form of determined system of 9 polynomials developed by symbolic computation. 

 Numerical illustration is presented with real world geodetic coordinates representing Hungarian Datum. The  computations were carried out with hp workstation xw 4100 with XP operation system, 3 GHz P4 Intel processor and 1 GB RAM.

2. Definition of the 3 points 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 (W).

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

where W=(({{(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 in general is given by using 3 axial rotation  (α,β,γ - Cardan angles).

R=(Subscript[R, 1])(α) (Subscript[R, 2])(β) (Subscript[R, 3])(γ)            with              (Subscript[R, 1])(α)=(({{1, 0, 0}, {0, cos α, sin α}, {0, -sin α, cos α}})), (Subscript[R, 2])(β)=(({{cos β, 0, -sin β}, {0, 1, 0}, {sin β, 0, cos β}})), (Subscript[R, 3])(γ)=(({{cos γ, sin γ, 0}, {-sin γ, cos γ, 0}, {0, 0, 1}}))

leading to

R(α,β,γ)=(({{cos β cos γ, cos β sin γ, -sinβ}, {sin α sin β cos γ-cos α sin γ, sin α sin β sin γ+cos α cos γ, sin αcos β}, {cos α sin β cos γ+sin α sin γ
, cos α sin β sin γ-sin α cos γ, cos α cos β}}))

In the traditional 7 or 9 parameter transformation solution (Papp E., Szucs L. (2005)) the single three by - three rotation matrix is simplified from three separate rotation matrices by assuming that each axial rotation is differentially small (typically less than five arc seconds for most geodetic networks), thus permitting binomial series expansions of the sine and cosine terms for radian measure.

The rotation matrix can be expressed with the skew-symmetric matrix (S) also (see Awange and Grafarend (2003)), and this facilitates the symbolical-numerical solution of the problem without using  simplifications.

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

The rotation matrix is

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

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

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

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

for which the next restriction is true:

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

The axial rotation angles (Cardan angles) can be obtained from the rotation matrix (R) through:

tan α =((Subscript[r, 23])/(Subscript[r, 33])) => (Subscript[Λ, Γ])=arctan ((Subscript[r, 23])/(Subscript[r, 33]))

tan β =((-(Subscript[r, 31]))/(Sqrt[((Subscript[r, 11])^2)+((Subscript[r, 12])^2)]))=>(Subscript[ϕ, Γ])=arctan ((-(Subscript[r, 31]))/(Sqrt[((Subscript[r, 11])^2)+((Subscript[r, 12])^2)])), vagy β=-arcsin (Subscript[r, 13])

tan γ=((Subscript[r, 12])/(Subscript[r, 11]))=>(Subscript[Ξ, Γ])=arctan ((Subscript[r, 12])/(Subscript[r, 11]))

Cardan angles in degree

αR[R_]:=ArcTan[R[[3,3]],R[[2,3]]]*(180/π);βR[R_]:=-ArcSin[R[[1,3]]]*(180/π);γR[R_]:=ArcTan[R[[1,1]],R[[1,2]]]*(180/π);

Cardan[R_]:={αR[R],βR[R],γR[R]}

Cardan angles in seconds

αRs[R_]:=αR[R]*3600;βRs[R_]:=βR[R]*3600;γRs[R_]:=γR[R]*3600;

CardanS[R_]:={αRs[R],βRs[R],γRs[R]}

Instead of the scale matrix (W) we can use the inverse of the scales getting simplier equations in this way.

Let us call (Subscript[σ, i]) =(1/(Subscript[s, i]))and Ω = (W^(-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 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.
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:

For i = 1 point

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

For  i = 2 point

(({{(Subscript[f, 4])}, {(Subscript[f, 5])}, {(Subscript[f, 6])}}))=(((Subscript[I, 3])-S).Ω.(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))- ((Subscript[I, 3])+S).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3])-S).Ω.(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i->2//Expand;

For  i = 3 point

(({{(Subscript[f, 7])}, {(Subscript[f, 8])}, {(Subscript[f, 9])}}))=(((Subscript[I, 3])-S).Ω.(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))- ((Subscript[I, 3])+S).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3])-S).Ω.(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i->3//Expand;

or

sys=Table[(Subscript[f, i]),{i,1,9}];
eq=Table[{(Subscript["f", i]),sys[[i]]},{i,1,9}];TableForm[eq]

Similarly to (Subscript[C, 7])(3,3) problem  (see Awange and Grafarend (2005)  we can reduce the equation system by subtracting the equations. In this way the translation parameters can be eliminated. The elimination of the translation parameters can be done according to the next substractions

r1=(Subscript[f, 1])-(Subscript[f, 7]);

r2=(Subscript[f, 4])-(Subscript[f, 7]);

r3=(Subscript[f, 2])-(Subscript[f, 8]);

r4=(Subscript[f, 5])-(Subscript[f, 8]);

r5=(Subscript[f, 3])-(Subscript[f, 9]);

r6=(Subscript[f, 6])-(Subscript[f, 9]);

In the remaining 6 equations there is  only 6 unknown parameters (a, b, c,  (Subscript[σ, 1]),  (Subscript[σ, 2]),  (Subscript[σ, 3])).

sysR={r1,r2,r3,r4,r5,r6}

Employing Groebner basis and Dixon resultant, the above described reduced equation system can be easily reduced further , if we introduce some new variables, let us call them relative coordinates instead of original  ones, see Awange-Grafarend (2003). Let

SYSR=sysR//FullSimplify;SYSR//TableForm

It is clear that the introduction of following variables can simplify the description of the system

newVarsA={x12->(Subscript[x, 1])-(Subscript[x, 2]),x13->(Subscript[x, 1])-(Subscript[x, 3]),x23->(Subscript[x, 2])-(Subscript[x, 3]),y12->(Subscript[y, 1])-(Subscript[y, 2]),y13->(Subscript[y, 1])-(Subscript[y, 3]),y23->(Subscript[y, 2])-(Subscript[y, 3]),
z12->(Subscript[z, 1])-(Subscript[z, 2]),z13->(Subscript[z, 1])-(Subscript[z, 3]),z23->(Subscript[z, 2])-(Subscript[z, 3]),
X12->(Subscript[X, 1])-(Subscript[X, 2]),X13->(Subscript[X, 1])-(Subscript[X, 3]),X23->(Subscript[X, 2])-(Subscript[X, 3]),
Y12->(Subscript[Y, 1])-(Subscript[Y, 2]),Y13->(Subscript[Y, 1])-(Subscript[Y, 3]),Y23->(Subscript[Y, 2])-(Subscript[Y, 3]),
Z12->(Subscript[Z, 1])-(Subscript[Z, 2]),Z13->(Subscript[Z, 1])-(Subscript[Z, 3]),Z23->(Subscript[Z, 2])-(Subscript[Z, 3])};

then our equation system is

(Subscript[g, 1])=-X13+c Y13-b Z13+x13 (Subscript[σ, 1])+c y13 (Subscript[σ, 2])-b z13 (Subscript[σ, 3]);

(Subscript[g, 2])=-X23+c Y23-b Z23+x23 (Subscript[σ, 1])+c y23 (Subscript[σ, 2])-b z23 (Subscript[σ, 3]);

(Subscript[g, 3])=-c X13-Y13+a Z13-c x13 (Subscript[σ, 1])+y13 (Subscript[σ, 2])+a z13 (Subscript[σ, 3]);

(Subscript[g, 4])=-c X23-Y23+a Z23-c x23 (Subscript[σ, 1])+y23 (Subscript[σ, 2])+a z23 (Subscript[σ, 3]);

(Subscript[g, 5])=b X13-a Y13-Z13+b x13 (Subscript[σ, 1])-a y13 (Subscript[σ, 2])+z13 (Subscript[σ, 3]);

(Subscript[g, 6])=b X23-a Y23-Z23+b x23 (Subscript[σ, 1])-a y23 (Subscript[σ, 2])+z23 (Subscript[σ, 3]);

Let us check these equations

({(Subscript[g, 1]),(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5]),(Subscript[g, 6])}/.newVarsA)-SYSR//FullSimplify

3. Numerical Solutions

Here we shall consider three global methods working without guess of initial values: 

a) a general solver for polynomial systems based on numerical Groebner basis and eigensystem methods, built in the Mathematica system as a function  NSolve
b) Global minimization based on genetic algorithm, built -in in Mathematica as NMinimize
c) A linear homotopy methods written in Mathematica as a function  LinearHomotopyFR, see Palancz (2008) and using the built in function FindRoot employing Newton-Raphson method.

Let us consider the numerical values of  3 Hungarian points in the system of ETRS89 ((Subscript[x, 1]), (Subscript[y, 1]), (Subscript[z, 1]), ... , (Subscript[z, 3]) ) and in the local Hungarian system HD72 (Hungarian Datum 1972) ((Subscript[X, !]), (Subscript[Y, 1]), (Subscript[Z, 1]), ..., (Subscript[Z, 3])).

numericalValues={(Subscript[x, 1])->4171409.677,(Subscript[x, 2])->4146957.889,(Subscript[x, 3])->3955632.880,(Subscript[y, 1])->1470823.777,(Subscript[y, 2])->1277033.850,(Subscript[y, 3])->1611863.197,
(Subscript[z, 1])->4580140.907,(Subscript[z, 2])->4659439.264,(Subscript[z, 3])->4720991.316,(Subscript[X, 1])->4171352.311,(Subscript[X, 2])->4146901.301,(Subscript[X, 3])->3955575.649,
(Subscript[Y, 1])->1470893.887,(Subscript[Y, 2])->1277104.509,(Subscript[Y, 3])->1611933.124,(Subscript[Z, 1])->4580150.178,(Subscript[Z, 2])->4659448.287,(Subscript[Z, 3])-> 4721000.952};

The system in numerical form

Neqs={(Subscript[g, 1]),(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5]),(Subscript[g, 6])}/.newVarsA/.numericalValues

3.1 Global polynomial solver based on numerical Groebner basis and eigensystem methods

The solution providing all roots,

sol=NSolve[Neqs,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])}]

From the eight solutions we need the only one, where the scale variables are positive, (Subscript[σ, i]) > 0, therefore

sol3NSolve=Select[sol,(#[[4,2]]>0)∧(#[[5,2]]>0)∧(#[[6,2]]>0)&]

where

(Subscript[σ, 1])=1.0000054081636032

The computation time

sol=Timing[NSolve[Neqs,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])}];]

Rotation angles in seconds :

CardanS[R/.sol3NSolve[[1]]]

Scale parameters :

SetPrecision[{1/(Subscript[σ, 1]),1/(Subscript[σ, 2]),1/(Subscript[σ, 3])}/.sol3NSolve[[1]],10]

3.2 Global minimization

Using directly the least squares method, we can define the an objective  function, the sum of the squares of the errors, which should be minimized. the objective function is

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

NMin=NMinimize[nobj,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])},Method->"DifferentialEvolution"]//Timing

where for example

(Subscript[σ, 1])=1.0000054081636032

Rotation angles in seconds :

CardanS[R/.NMin[[2,2]]]

Scale parameters :

SetPrecision[{1/(Subscript[σ, 1]),1/(Subscript[σ, 2]),1/(Subscript[σ, 3])}/.NMin[[2,2]],10]

3.3 Homotopy Solution

We can employ the convex  linear or Bézout homotopy, see Drexler(1977), Garcia and Zangwill (1979).  The homotopy function is

H (x, β) = γ (1 - β) G (x) + β F (x)

and while β is changing from 0 to 1, the starting system G(x) will be transformed into the original system F(x) to be solved. The method has been implemented into Mathematica as a function LinearHomotopyFR, see Paláncz(2008) and using the built in function FindRoot employing Newton-Raphson method, however another approach is to set up and solve a system of ODEs that track the path, for example see  Binous (2007).

LinearHomotopyFR[F_,G_,X_,X0_,γ_,n_]:=Module[{H,X0L,λ0,i,β,m,R,RR,j,k,sol},
Off[FindRoot::"lstol"];
λ0=Table[i (1/n),{i,0,n}];
H=Flatten[(1-β) Thread[G  γ]+β F];
m=Length[X];
k=Length[X0];
RR={};
Do[
X0L={X0[[j]]};
Do[AppendTo[X0L,Map[#[[2]]&,FindRoot[H/.β->λ0[[i+1]],MapThread[{#1,#2}&,{X,X0L[[i]]}]]]],{i,1,n}];
R={};
Do[AppendTo[R,Interpolation[MapThread[{#1,#2[[i]]}&,{λ0,X0L}]]],{i,1,m}];AppendTo[RR,{Map[Chop[N[#[1]]]&,R],R}],
{j,1,k}];
sol=Transpose[RR];{sol[[1]],Table[MapThread[#1[λ]->#2[λ]&,{X,sol[[2,i]]}],{i,1,Length[X0]}]}]

This function will compute the homotopy paths using successive Newton-Raphson method, 

Input variables
F   - list of functions of the polynomial system to be solved, F = {(Subscript[f, 1])(x), (Subscript[f, 2])(x),...,(Subscript[f, n])(x)}
G  - list of the starting system, G = {(Subscript[g, 1])(x),(Subscript[g, 2])(x),...,(Subscript[g, n])(x)}, 
X  - list of the independent variables X =  {(Subscript[x, 1]),(Subscript[x, 2]),...,(Subscript[x, n])}
X0  -  list of initial values, X0 = {{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_1,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_2,...,{(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_m}
         where m is the number of the roots of  F 
 γ   -  list of complex numbers, {(Subscript[γ, 1]), (Subscript[γ, 2]),...,(Subscript[γ, n])}, it means (Subscript[γ, i]) can be different for every (Subscript[g, i])(x), 
          If the starting  system is generated for polynomials all  (Subscript[γ, i]) = 1.
 n  -   number of the subintervalls  in [0, 1].

Output variables

sol[1]    -  list of the i-th solutions correspondig to the i-th initial values,  {(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_i, i = 1...m
sol[2]    -  list of the path of  i-th solutions in form of interpolating functions of the variables  
                 correspondig to the i-th  initial values,  {(Subscript[x0, 1]),(Subscript[x0, 2]),...,(Subscript[x0, n])}_i,   i = 1...m,
                 {{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_1,{(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }_2,...,{(Subscript[(Subscript[φ, 1]),(Subscript[φ, 2]),...,(Subscript[φ, n]) }, m])},  (Subscript[φ, i])=(Subscript[φ, i])(λ)

Instead of computing all of the roots bounded by the Bezout bound, which is 64 in this case (six quadratics),  let us consider our starting sytem as the linear part of the nonlinear system,

G1=-X13+c Y13-b Z13+x13 (Subscript[σ, 1]);

G2=-X23+c Y23-b Z23+x23 (Subscript[σ, 1]);

G3=-c X13-Y13+a Z13+y13 (Subscript[σ, 2]);

G4=-c X23-Y23+a Z23+y23 (Subscript[σ, 2]);

G5=b X13-a Y13-Z13+z13 (Subscript[σ, 3]);

G6=b X23-a Y23-Z23+z23 (Subscript[σ, 3]);

In numerical form

Geqs={G1,G2,G3,G4,G5,G6}/.newVarsA/.numericalValues

The starting values can be easily computed by solving  this linear system

solG=NSolve[Geqs,{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])}]//Flatten

The list of the variables

XX={a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])};

The initial values

X0={Map[#[[2]]&,solG]}

Now, the starting system is Geqs and the target system is Neqs. In this case , we can work with real values, there is no need to use complex number to avoid singularity of the homotopy function

γ={1,1,1,1,1,1};

sol= LinearHomotopyFR[Neqs,Geqs,XX,X0,γ,10];

sol[[1]]

where

(Subscript[σ, 1])=1.000005408163603

Rotation angles in seconds :

CardanS[R/.{a->sol[[1,1,1]],b->sol[[1,1,2]],c->sol[[1,1,3]]}]

Scale parameters :

SetPrecision[{1/(Subscript[σ, 1]),1/(Subscript[σ, 2]),1/(Subscript[σ, 3])}/.{(Subscript[σ, 1])->sol[[1,1,4]],(Subscript[σ, 2])->sol[[1,1,5]],(Subscript[σ, 3])->sol[[1,1,6]]},10]

The computation time

sol= Timing[LinearHomotopyFR[Neqs,Geqs,XX,X0,γ,10];]

Remarks : According to our numerical experiences, the generation of the starting system for the homotopy function from the linear part of the original nonlinear one is proved to be reliable in case of many different systems. Homotopy solution is roughly two-three times faster than the global polynomial solver  (NSolve) and  20- 30 times faster than global minimization (NMinimize). In addition  there is no need to select the proper (positive) solution like in case of using NSolve.

Table 1.
Comparing Global Numerical Methods for solving 3 Point Problem

({{Method, Time of computation
(sec), Solution}, {Numerical Groebner basis
&  eigensystem method, 0.125, many}, {Global minimization, 1.110, one}, {Homotopy method, 0.031, one}})

4. Symbolic solutions

4.1 Improved Dixon resultant - Kapur, Saxena  and Yang method

Our intention is to reduce the computation time of the 3 point solution, because, for example, in case of Gauss-Jacobi method for N points, it is very important to use an effective method for  solving one of the many combinations. This reduction can be done by Dixon resultant -  see  Nakos and Williams  (2002) and originally suggested by Kapur, Saxena  and Yang (1994) -  or by reduced Groebner basis as well as by using other multipolynomial resultant. Here we employ the Dixon resultant.

<<Resultant`Dixon`

Unfortunately, even the reduced equation system with 6 nonlinear equations and  with 6 unknown parameters had no direct reduction to a univariate polynomial  using  Dixon resultant via Mathematica. But we can use either method to reduce the equation system to 3 equations with 3 parameters. For the sake  of illustration and completion we will carry out this computation. However , a far superior method  will be presented in the next section. To begin, one has to  three times select   4 - 4 independent equations (fourplexes) from the 6 equations, arbitrarily. From a fourplex, we can eliminate a, b, c parameters using Dixon or Grobner method (both give the same result, here we employed the Dixon resultant). Choosing three independent combinations of  the fourplexes, we get three nonlinear equations for the scale parameters (respectively for the inverse of the 3 scale parameters, (Subscript[σ, i]) ).
Elimination of a, b and c from the fourplexes:

dr1σ123=DixonResultant[{(Subscript[g, 1]),(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4])},{a,b,c},{s1,s2,s3}]

and similarly

dr2σ123=DixonResultant[{(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5])},{a,b,c},{s1,s2,s3}];

dr3σ123=DixonResultant[{(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5]),(Subscript[g, 6])},{a,b,c},{s1,s2,s3}];

The succesfull application of Dixon - KSY, in order to carry out further reduction, needs to compute compact coefficients of the remained three polynomials.

In order to determine the compact form of the polynomial equations,  the following functions, providing to find coefficients in a multivariate polynomial, are developed,

Expon[equ_]:=Union[Table[Exponent[equ[[i]],Reverse[vari]],{i,Length[equ]}]]

Koeff[poly_,{i_,j_,k_}]:=CoefficientList[poly,Reverse[vari]][[i+1,j+1,k+1]]

KoeffList[poly_]:=Table[Koeff[poly,Expon[poly][[i]]],{i,Length[Expon[poly]]}]

ExpVari[{x_,y_,z_},{i_,j_,k_}]:=x^i y^j z^k

ExpVariList[poly_]:=Table[ExpVari[Reverse[vari],Expon[poly][[i]]],{i,Length[Expon[poly]]}]

vari={x,y,z};

pl=15+x^2+2y+3x y+6 x^3+8 x + 9 x+20 z

ExpVariList[pl]

KoeffList[pl]

Let us apply this method for our system. The variables are

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

The first equation is

dr1P=dr1σ123//Expand

The corresponding terms are

d1=ExpVariList[dr1P]

and their coefficients

C1=KoeffList[dr1P]

It is easy to check our result

dr1P-(C1.d1)//Simplify

Now we shall introduce the following new coefficients representing the original ones

HH=Table[(Subscript[H, i]),{i,0,Length[d1]-1}]

Then the proper assignments are

C1s=MapThread[#1->#2&,{HH,C1}]

Consequently, our first equation can be written in a compact form, namely

e1=HH.d1

Similarly, one can carry out these computations for the two other polynomials.

The compact form of the second and third equations,

f1=(Subscript[J, 0])+(Subscript[J, 1]) (Subscript[σ, 1])+(Subscript[J, 2]) (Subscript[σ, 1]^2)+(Subscript[J, 3]) (Subscript[σ, 1]^3)+(Subscript[J, 4]) (Subscript[σ, 2]^2)+(Subscript[J, 5]) (Subscript[σ, 1]) (Subscript[σ, 2]^2)+(Subscript[J, 6]) (Subscript[σ, 3])+(Subscript[J, 7]) (Subscript[σ, 1]) (Subscript[σ, 3])+(Subscript[J, 8]) (Subscript[σ, 1]^2) (Subscript[σ, 3])+(Subscript[J, 9]) (Subscript[σ, 1]^3) (Subscript[σ, 3])+(Subscript[J, 10]) (Subscript[σ, 2]^2) (Subscript[σ, 3])+(Subscript[J, 11]) (Subscript[σ, 1]) (Subscript[σ, 2]^2) (Subscript[σ, 3])+(Subscript[J, 12]) (Subscript[σ, 3]^2)+(Subscript[J, 13]) (Subscript[σ, 1]) (Subscript[σ, 3]^2)+(Subscript[J, 14]) (Subscript[σ, 3]^3)+(Subscript[J, 15]) (Subscript[σ, 1]) (Subscript[σ, 3]^3);

g1=(Subscript[K, 0])+(Subscript[K, 1]) (Subscript[σ, 1])+(Subscript[K, 2]) (Subscript[σ, 1]^2)+(Subscript[K, 3]) (Subscript[σ, 2]^2)+(Subscript[K, 4]) (Subscript[σ, 1]) (Subscript[σ, 2]^2)+(Subscript[K, 5]) (Subscript[σ, 1]^2) (Subscript[σ, 2]^2)+(Subscript[K, 6]) (Subscript[σ, 3]^2)+(Subscript[K, 7]) (Subscript[σ, 1]) (Subscript[σ, 3]^2)+(Subscript[K, 8]) (Subscript[σ, 1]^2) (Subscript[σ, 3]^2);

and their coefficients

C2s={(Subscript[J, 0])->X13 (X23^2) Z13+X13 (Y23^2) Z13-(X13^2) X23 Z23+X23 (Y13^2) Z23-2 X13 Y13 Y23 Z23+X23 (Z13^2) Z23-X13 Z13 (Z23^2),(Subscript[J, 1])->x13 (X23^2) Z13+x13 (Y23^2) Z13+(X13^2) x23 Z23-2 x13 X13 X23 Z23+x23 (Y13^2) Z23-2 x13 Y13 Y23 Z23+x23 (Z13^2) Z23-x13 Z13 (Z23^2),(Subscript[J, 2])->-X13 (x23^2) Z13+2 x13 X13 x23 Z23-(x13^2) X23 Z23,(Subscript[J, 3])->-x13 (x23^2) Z13+(x13^2) x23 Z23,(Subscript[J, 4])->-X13 (y23^2) Z13-X23 (y13^2) Z23+2 X13 y13 y23 Z23,(Subscript[J, 5])->-x13 (y23^2) Z13-x23 (y13^2) Z23+2 x13 y13 y23 Z23,(Subscript[J, 6])->X13 (X23^2) z13+X13 (Y23^2) z13-(X13^2) X23 z23+X23 (Y13^2) z23-2 X13 Y13 Y23 z23+X23 (Z13^2) z23-2 X13 Z13 z23 Z23+X13 z13 (Z23^2),(Subscript[J, 7])->x13 (X23^2) z13+x13 (Y23^2) z13+(X13^2) x23 z23-2 x13 X13 X23 z23+x23 (Y13^2) z23-2 x13 Y13 Y23 z23+x23 (Z13^2) z23-2 x13 Z13 z23 Z23+x13 z13 (Z23^2),(Subscript[J, 8])->-X13 (x23^2) z13+2 x13 X13 x23 z23-(x13^2) X23 z23,(Subscript[J, 9])->-x13 (x23^2) z13+(x13^2) x23 z23,(Subscript[J, 10])->-X13 (y23^2) z13-X23 (y13^2) z23+2 X13 y13 y23 z23,(Subscript[J, 11])->-x13 (y23^2) z13-x23 (y13^2) z23+2 x13 y13 y23 z23,(Subscript[J, 12])->-X13 Z13 (z23^2)-X23 (z13^2) Z23+2 X13 z13 z23 Z23,(Subscript[J, 13])->-x13 Z13 (z23^2)-x23 (z13^2) Z23+2 x13 z13 z23 Z23,(Subscript[J, 14])->-X23 (z13^2) z23+X13 z13 (z23^2),(Subscript[J, 15])->-x23 (z13^2) z23+x13 z13 (z23^2)};

C3s={(Subscript[K, 0])->-(X23^2) (Y13^2)+2 X13 X23 Y13 Y23-(X13^2) (Y23^2)-(X23^2) (Z13^2)+2 X13 X23 Z13 Z23-(X13^2) (Z23^2),(Subscript[K, 1])->-2 x23 X23 (Y13^2)+2 X13 x23 Y13 Y23+2 x13 X23 Y13 Y23-2 x13 X13 (Y23^2)-2 x23 X23 (Z13^2)+2 X13 x23 Z13 Z23+2 x13 X23 Z13 Z23-2 x13 X13 (Z23^2),(Subscript[K, 2])->-(x23^2) (Y13^2)+2 x13 x23 Y13 Y23-(x13^2) (Y23^2)-(x23^2) (Z13^2)+2 x13 x23 Z13 Z23-(x13^2) (Z23^2),(Subscript[K, 3])->(X23^2) (y13^2)-2 X13 X23 y13 y23+(X13^2) (y23^2),(Subscript[K, 4])->2 x23 X23 (y13^2)-2 X13 x23 y13 y23-2 x13 X23 y13 y23+2 x13 X13 (y23^2),(Subscript[K, 5])->(x23^2) (y13^2)-2 x13 x23 y13 y23+(x13^2) (y23^2),(Subscript[K, 6])->(X23^2) (z13^2)-2 X13 X23 z13 z23+(X13^2) (z23^2),(Subscript[K, 7])->2 x23 X23 (z13^2)-2 X13 x23 z13 z23-2 x13 X23 z13 z23+2 x13 X13 (z23^2),(Subscript[K, 8])->(x23^2) (z13^2)-2 x13 x23 z13 z23+(x13^2) (z23^2)};

Now, we can reduce our three polynomial system to two polynomials with two variables via pairwise elimination technique. First, let us consider the first and third polynomials and eliminate the variable (Subscript[σ, 2]) from both, employing reduced Gröbner basis computed via Groebner - walk and then do the same with the first and second polynomails,

eg13=GroebnerBasis[{e1,g1},{(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])},{(Subscript[σ, 2])}]

ef13=GroebnerBasis[{e1,f1},{(Subscript[σ, 1]),(Subscript[σ, 2]),(Subscript[σ, 3])},{(Subscript[σ, 2])}]

The succesfull application of Dixon - KSY, in order to carry out further reduction, needs to compute compact coefficients of the remaining three polynomials, similarly as before.

The compact form of the two equations are,

u1=(Subscript[U, 0])+(Subscript[U, 1]) (Subscript[σ, 1])+(Subscript[U, 2]) (Subscript[σ, 1]^2)+(Subscript[U, 3]) (Subscript[σ, 1]^3)+(Subscript[U, 4]) (Subscript[σ, 1]^4)+(Subscript[U, 5]) (Subscript[σ, 3])+(Subscript[U, 6]) (Subscript[σ, 1]) (Subscript[σ, 3])+(Subscript[U, 7]) (Subscript[σ, 1]^2) (Subscript[σ, 3])+(Subscript[U, 8]) (Subscript[σ, 1]^3) (Subscript[σ, 3])+(Subscript[U, 9]) (Subscript[σ, 1]^4) (Subscript[σ, 3])+(Subscript[U, 10]) (Subscript[σ, 3]^2)+(Subscript[U, 11]) (Subscript[σ, 1]) (Subscript[σ, 3]^2)+(Subscript[U, 12]) (Subscript[σ, 1]^2) (Subscript[σ, 3]^2)+(Subscript[U, 13]) (Subscript[σ, 1]^3) (Subscript[σ, 3]^2)+(Subscript[U, 14]) (Subscript[σ, 1]^4) (Subscript[σ, 3]^2)+(Subscript[U, 15]) (Subscript[σ, 3]^3)+(Subscript[U, 16]) (Subscript[σ, 1]) (Subscript[σ, 3]^3)+(Subscript[U, 17]) (Subscript[σ, 1]^2) (Subscript[σ, 3]^3)+(Subscript[U, 18]) (Subscript[σ, 3]^4)+(Subscript[U, 19]) (Subscript[σ, 1]) (Subscript[σ, 3]^4)+(Subscript[U, 20]) (Subscript[σ, 1]^2) (Subscript[σ, 3]^4);

v1=(Subscript[V, 0])+(Subscript[V, 1]) (Subscript[σ, 1])+(Subscript[V, 2]) (Subscript[σ, 1]^2)+(Subscript[V, 3]) (Subscript[σ, 1]^3)+(Subscript[V, 4]) (Subscript[σ, 3])+(Subscript[V, 5]) (Subscript[σ, 1]) (Subscript[σ, 3])+(Subscript[V, 6]) (Subscript[σ, 1]^2) (Subscript[σ, 3])+(Subscript[V, 7]) (Subscript[σ, 1]^3) (Subscript[σ, 3])+(Subscript[V, 8]) (Subscript[σ, 3]^2)+(Subscript[V, 9]) (Subscript[σ, 1]) (Subscript[σ, 3]^2)+(Subscript[V, 10]) (Subscript[σ, 1]^2) (Subscript[σ, 3]^2)+(Subscript[V, 11]) (Subscript[σ, 1]^3) (Subscript[σ, 3]^2)+(Subscript[V, 12]) (Subscript[σ, 3]^3)+(Subscript[V, 13]) (Subscript[σ, 1]) (Subscript[σ, 3]^3)+(Subscript[V, 14]) (Subscript[σ, 1]^2) (Subscript[σ, 3]^3)+(Subscript[V, 15]) (Subscript[σ, 1]^3) (Subscript[σ, 3]^3)+(Subscript[V, 16]) (Subscript[σ, 3]^4)+(Subscript[V, 17]) (Subscript[σ, 1]) (Subscript[σ, 3]^4)+(Subscript[V, 18]) (Subscript[σ, 3]^5)+(Subscript[V, 19]) (Subscript[σ, 1]) (Subscript[σ, 3]^5);

and their coefficients

C13s={(Subscript[U, 0])->-(Subscript[H, 2]) (Subscript[K, 0])+(Subscript[H, 0]) (Subscript[K, 3]),(Subscript[U, 1])->-(Subscript[H, 2]) (Subscript[K, 1])+(Subscript[H, 0]) (Subscript[K, 4]),(Subscript[U, 2])->-(Subscript[H, 2]) (Subscript[K, 2])+(Subscript[H, 1]) (Subscript[K, 3])+(Subscript[H, 0]) (Subscript[K, 5]),(Subscript[U, 3])->(Subscript[H, 1]) (Subscript[K, 4]),(Subscript[U, 4])->(Subscript[H, 1]) (Subscript[K, 5]),(Subscript[U, 5])->-(Subscript[H, 5]) (Subscript[K, 0])+(Subscript[H, 3]) (Subscript[K, 3]),(Subscript[U, 6])->-(Subscript[H, 5]) (Subscript[K, 1])+(Subscript[H, 3]) (Subscript[K, 4]),(Subscript[U, 7])->-(Subscript[H, 5]) (Subscript[K, 2])+(Subscript[H, 4]) (Subscript[K, 3])+(Subscript[H, 3]) (Subscript[K, 5]),(Subscript[U, 8])->(Subscript[H, 4]) (Subscript[K, 4]),(Subscript[U, 9])->(Subscript[H, 4]) (Subscript[K, 5]),(Subscript[U, 10])->-(Subscript[H, 8]) (Subscript[K, 0])+(Subscript[H, 6]) (Subscript[K, 3])-(Subscript[H, 2]) (Subscript[K, 6]),(Subscript[U, 11])->-(Subscript[H, 8]) (Subscript[K, 1])+(Subscript[H, 6]) (Subscript[K, 4])-(Subscript[H, 2]) (Subscript[K, 7]),(Subscript[U, 12])->-(Subscript[H, 8]) (Subscript[K, 2])+(Subscript[H, 7]) (Subscript[K, 3])+(Subscript[H, 6]) (Subscript[K, 5])-(Subscript[H, 2]) (Subscript[K, 8]),(Subscript[U, 13])->(Subscript[H, 7]) (Subscript[K, 4]),(Subscript[U, 14])->(Subscript[H, 7]) (Subscript[K, 5]),(Subscript[U, 15])->-(Subscript[H, 5]) (Subscript[K, 6]),(Subscript[U, 16])->-(Subscript[H, 5]) (Subscript[K, 7]),(Subscript[U, 17])->-(Subscript[H, 5]) (Subscript[K, 8]),(Subscript[U, 18])->-(Subscript[H, 8]) (Subscript[K, 6]),(Subscript[U, 19])->-(Subscript[H, 8]) (Subscript[K, 7]),(Subscript[U, 20])->-(Subscript[H, 8]) (Subscript[K, 8])};

F13s={(Subscript[V, 0])->-(Subscript[H, 2]) (Subscript[J, 0])+(Subscript[H, 0]) (Subscript[J, 4]),(Subscript[V, 1])->-(Subscript[H, 2]) (Subscript[J, 1])+(Subscript[H, 0]) (Subscript[J, 5]),(Subscript[V, 2])->-(Subscript[H, 2]) (Subscript[J, 2])+(Subscript[H, 1]) (Subscript[J, 4]),(Subscript[V, 3])->-(Subscript[H, 2]) (Subscript[J, 3])+(Subscript[H, 1]) (Subscript[J, 5]),(Subscript[V, 4])->-(Subscript[H, 5]) (Subscript[J, 0])+(Subscript[H, 3]) (Subscript[J, 4])-(Subscript[H, 2]) (Subscript[J, 6])+(Subscript[H, 0]) (Subscript[J, 10]),(Subscript[V, 5])->-(Subscript[H, 5]) (Subscript[J, 1])+(Subscript[H, 3]) (Subscript[J, 5])-(Subscript[H, 2]) (Subscript[J, 7])+(Subscript[H, 0]) (Subscript[J, 11]),(Subscript[V, 6])->-(Subscript[H, 5]) (Subscript[J, 2])+(Subscript[H, 4]) (Subscript[J, 4])-(Subscript[H, 2]) (Subscript[J, 8])+(Subscript[H, 1]) (Subscript[J, 10]),(Subscript[V, 7])->-(Subscript[H, 5]) (Subscript[J, 3])+(Subscript[H, 4]) (Subscript[J, 5])-(Subscript[H, 2]) (Subscript[J, 9])+(Subscript[H, 1]) (Subscript[J, 11]),(Subscript[V, 8])->-(Subscript[H, 8]) (Subscript[J, 0])+(Subscript[H, 6]) (Subscript[J, 4])-(Subscript[H, 5]) (Subscript[J, 6])+(Subscript[H, 3]) (Subscript[J, 10])-(Subscript[H, 2]) (Subscript[J, 12]),(Subscript[V, 9])->-(Subscript[H, 8]) (Subscript[J, 1])+(Subscript[H, 6]) (Subscript[J, 5])-(Subscript[H, 5]) (Subscript[J, 7])+(Subscript[H, 3]) (Subscript[J, 11])-(Subscript[H, 2]) (Subscript[J, 13]),(Subscript[V, 10])->-(Subscript[H, 8]) (Subscript[J, 2])+(Subscript[H, 7]) (Subscript[J, 4])-(Subscript[H, 5]) (Subscript[J, 8])+(Subscript[H, 4]) (Subscript[J, 10]),(Subscript[V, 11])->-(Subscript[H, 8]) (Subscript[J, 3])+(Subscript[H, 7]) (Subscript[J, 5])-(Subscript[H, 5]) (Subscript[J, 9])+(Subscript[H, 4]) (Subscript[J, 11]),(Subscript[V, 12])->-(Subscript[H, 8]) (Subscript[J, 6])+(Subscript[H, 6]) (Subscript[J, 10])-(Subscript[H, 5]) (Subscript[J, 12])-(Subscript[H, 2]) (Subscript[J, 14]),(Subscript[V, 13])->-(Subscript[H, 8]) (Subscript[J, 7])+(Subscript[H, 6]) (Subscript[J, 11])-(Subscript[H, 5]) (Subscript[J, 13])-(Subscript[H, 2]) (Subscript[J, 15]),(Subscript[V, 14])->-(Subscript[H, 8]) (Subscript[J, 8])+(Subscript[H, 7]) (Subscript[J, 10]),(Subscript[V, 15])->-(Subscript[H, 8]) (Subscript[J, 9])+(Subscript[H, 7]) (Subscript[J, 11]),(Subscript[V, 16])->-(Subscript[H, 8]) (Subscript[J, 12])-(Subscript[H, 5]) (Subscript[J, 14]),(Subscript[V, 17])->-(Subscript[H, 8]) (Subscript[J, 13])-(Subscript[H, 5]) (Subscript[J, 15]),(Subscript[V, 18])->-(Subscript[H, 8]) (Subscript[J, 14]),(Subscript[V, 19])->-(Subscript[H, 8]) (Subscript[J, 15])};

The last step is to eliminate (Subscript[σ, 3]) from both polynomials in order to get a monomial containing only one variable, namely (Subscript[σ, 1]). This is a very difficult problem because of the big size of the result. The enhanced Dixon resultant as well as the Buchberger and Groebner - walk algorithms all failed because they employ expanded form of the polynomials during the computation, however classical Dixon method working with factorized form was successful.

uv1=ClassicalDixonResultant[{u1,v1},{(Subscript[σ, 3])},{s1}]

The result is very big in size, in printed form with normal style is about 30 pages! Expanding and count of the number of terms was impossible with our machine!

uv1=ClassicalDixonResultant[{u1,v1},{(Subscript[σ, 3])},{s1}];//Timing

Substituting the numerical values into the monomial containing variable (Subscript[σ, 1])

uv1N=uv1//.F13s/.C13s/.C1s/.C2s/.C3s/.newVarsA/.numericalValues//Expand

which can be simplified as,

uv1Ns=uv1N (10^(-395))//Simplify

Fortunately we do not need to find all of the roots of this polynomial, because for the value of (Subscript[σ, 1]) a very good estimation can be given. In case of the 7 parameter similarity transformation 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). In case of the 9 parameter affine transformation where 3 different scale values ((Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3]) are applied accordingt to the 3 coordinate axis, a good approach for the scale parameters can be given  modifying the Albertz-Kreiling expression. Instead of the quotient of the two lenghtes in the centre of gravity system we can use the quotients of the sum of the lengthes in the corresponding coordinate axes directions.

The center of gravity in the two systems ((Subscript[x, s]),(Subscript[y, s]),(Subscript[z, s])) and ((Subscript[X, s]),(Subscript[Y, s]),(Subscript[Z, s])).

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

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

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

The estimated scale parameters according to the modified Albertz-Kreimlig expression for the 9 parameter transformation

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

and (Subscript[σ, i])=1/(Subscript[s, i]).

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

σ1priori=1/mxpriori; SetPrecision[σ1priori,10]

Therefore Newton - Raphson method can be easily employed

solσ1=FindRoot[SetPrecision[uv1Ns==0,16],{(Subscript[σ, 1]),σ1priori}, WorkingPrecision->16]

4.2 Accelerated Dixon resultant - Lewis EDF method

The univariate polynomial for (Subscript[σ, 1]) also can be computed employing the accelerated Dixon resultant by the Early Discovery Factors (EDF) algorithm, which was suggested and implemented in the computer algebra system Fermat  by Lewis (2007 a, b). Using this method one can get the result via direct elimination in one step in the following factored form in less than two seconds of CPU time

(∏_(i=1)^5) (Subscript[φ, i])((((Subscript[σ, 1])) )^(Subscript[K, i]))

where  (Subscript[φ, i]) ((Subscript[σ, 1]))  are  irreducible polynomials  with low degree, but  their powers, (Subscript[K, i])  are   big positive integer numbers, so expanding this expression would result in millions of terms!

Consequently, it suffices to use (Subscript[K, i])= 1, for i = 1,...,5, namely

(∏_(i=1)^5) (Subscript[φ, i])((Subscript[σ, 1]))

as the determinant of the Dixon matrix  for the resultant. These  polynomials are the following

(Subscript[φ, 11])=y13*z23-y23*z13;

(Subscript[φ, 21])=x13^2*y23*z23*(Subscript[σ, 1])^2-x13*x23*y13*z23*(Subscript[σ, 1])^2-x13*x23*y23*z13*(Subscript[σ, 1])^2+x23^2*y13*z13*(Subscript[σ, 1])^2-Z13^2*y23*z23-Y13^2*y23*z23-X13^2*y23*z23+Z13*Z23*y13*z23+Y13*Y23*y13*z23+X13*X23*y13*z23+Z13*Z23*y23*z13+Y13*Y23*y23*z13+X13*X23*y23*z13-Z23^2*y13*z13-Y23^2*y13*z13-X23^2*y13*z13;

(Subscript[φ, 31])=x13*y23*(Subscript[σ, 1])-x23*y13*(Subscript[σ, 1])+X13*y23-X23*y13;

(Subscript[φ, 41])=Z13*x13*x23*z23*1^(Subscript[σ, 1])^2-Z23*x13^2*z23*(Subscript[σ, 1])^2-Z13*x23^2*z13*(Subscript[σ, 1])^2+Z23*x13*x23*z13*(Subscript[σ, 1])^2+X13*Z13*x23*z23*(Subscript[σ, 1])-2*X13*Z23*x13*z23*(Subscript[σ, 1])+X23*Z13*x13*z23*(Subscript[σ, 1])+X13*Z23*x23*z13*(Subscript[σ, 1])-2*X23*Z13*x23*z13*(Subscript[σ, 1])+X23*Z23*x13*z13*(Subscript[σ, 1])-X13^2*Z23*z23+X13*X23*Z13*z23+X13*X23*Z23*z13-X23^2*Z13*z13;

(Subscript[φ, 51])=Z13*z23-Z23*z13;

Let us compute the roots of these factors

Map[NSolve[SetPrecision[#/.newVarsA/.numericalValues,16],(Subscript[σ, 1]) ,WorkingPrecision->16]&,{(Subscript[φ, 11]),(Subscript[φ, 21]),(Subscript[φ, 31]),(Subscript[φ, 41]),(Subscript[φ, 51])}]

We have here two positive solutions, but the only one is realistic, see modified Albertz and Kreiling estimation,

(Subscript[σ, 1])=1.000005408163603

The numerical solution with global minimization,

(Subscript[σ, 1])=1.0000054081636032

The result computed using  Dixon-KSY method is,

(Subscript[σ, 1])->1.00000540816750316612852439935`16.

It means that Dixon - EDF method provides a bit more precise solution.

The result for (Subscript[σ, 1]) in symbolic form,considering the positive root of (Subscript[φ, 2]) = 0 is

solσ1=(Solve[(Subscript[φ, 21])==0,(Subscript[σ, 1])][[2]]//Simplify)

Similarly we can get simple explicite form for (Subscript[σ, 2]) and (Subscript[σ, 3]), too.

For (Subscript[σ, 2]) we get

(Subscript[φ, 12])=x13*z23-x23*z13;

(Subscript[φ, 22])=x13*y13*y23*z23*(Subscript[σ, 2])^2-x23*y13^2*z23*(Subscript[σ, 2])^2-x13*y23^2*z13*(Subscript[σ, 2])^2+x23*y13*y23*z13*(Subscript[σ, 2])^2+Z13^2*x23*z23+Y13^2*x23*z23+X13^2*x23*z23-Z13*Z23*x13*z23-Y13*Y23*x13*z23-X13*X23*x13*z23-Z13*Z23*x23*z13-Y13*Y23*x23*z13-X13*X23*x23*z13+Z23^2*x13*z13+Y23^2*x13*z13+X23^2*x13*z13;

(Subscript[φ, 32])=x13*y23*(Subscript[σ, 2])-x23*y13*(Subscript[σ, 2])-Y13*x23+Y23*x13;

(Subscript[φ, 42])=Z13*y13*y23*z23*(Subscript[σ, 2])^2-Z23*y13^2*z23*(Subscript[σ, 2])^2-Z13*y23^2*z13*(Subscript[σ, 2])^2+Z23*y13*y23*z13*(Subscript[σ, 2])^2+Y13*Z13*y23*z23*(Subscript[σ, 2])-2*Y13*Z23*y13*z23*(Subscript[σ, 2])+Y23*Z13*y13*z23*(Subscript[σ, 2])+Y13*Z23*y23*z13*(Subscript[σ, 2])-2*Y23*Z13*y23*z13*(Subscript[σ, 2])+Y23*Z23*y13*z13*(Subscript[σ, 2])-Y13^2*Z23*z23+Y13*Y23*Z13*z23+Y13*Y23*Z23*z13-Y23^2*Z13*z13;

(Subscript[φ, 52])=Z13*z23-Z23*z13;

Let us compute the roots of these factors

Map[NSolve[SetPrecision[#/.newVarsA/.numericalValues,16],(Subscript[σ, 2]) ,WorkingPrecision->16]&,{(Subscript[φ, 12]),(Subscript[φ, 22]),(Subscript[φ, 32]),(Subscript[φ, 42]),(Subscript[φ, 52])}]

The only positive real solution is

(Subscript[σ, 2])->0.9999978437466494

The numerical solution  with global minimization

(Subscript[σ, 2])->0.9999978437466492

The result computed using  Dixon-KSY method is

(Subscript[σ, 2])->0.9999978437461572

It means that Dixon - EDF method provides a bit more precise solution again

The result for (Subscript[σ, 2])  in symbolic form,considering the positive root of (Subscript[φ, 2]) = 0 is

solσ2=(Solve[(Subscript[φ, 2])==0,(Subscript[σ, 2])][[2]]//Simplify)

For (Subscript[σ, 3]) we get only four factors

(Subscript[φ, 13])=y13*z23*(Subscript[σ, 3]) - y23*z13*(Subscript[σ, 3]) - Z13*y23 + Z23*y13;

(Subscript[φ, 23])= x13*y23 - x23*y13;

(Subscript[φ, 33])=x13*z23*(Subscript[σ, 3]) - x23*z13*(Subscript[σ, 3]) - Z13*x23 + Z23*x13;

(Subscript[φ, 43])=x13*y13*z23^2*(Subscript[σ, 3])^2 - x13*y23*z13*z23*(Subscript[σ, 3])^2 - x23*y13*z13*z23*(Subscript[σ, 3])^2 +  
x23*y23*z13^2*(Subscript[σ, 3])^2 - Z13^2*x23*y23 - Y13^2*x23*y23 - X13^2*x23*y23 + Z13*Z23*x13*y23 + Y13*Y23*x13*y23 + X13*X23*x13*y23 +  Z13*Z23*x23*y13 + Y13*Y23*x23*y13 + X13*X23*x23*y13- Z23^2*x13*y13 - Y23^2*x13*y13 - X23^2*x13*y13;

Let us compute the roots of these factors

Map[NSolve[SetPrecision[#/.newVarsA/.numericalValues,16],(Subscript[σ, 3]) ,WorkingPrecision->16]&,{(Subscript[φ, 13]),(Subscript[φ, 23]),(Subscript[φ, 33]),(Subscript[φ, 43])}]

The only positive solution is

(Subscript[σ, 3])->0.99998929165675978852219190266`16.000000000000004

The numerical solution  with global minimization

(Subscript[σ, 3])->0.9999892916567604

The result computed using  Dixon-KSY method is

(Subscript[σ, 3])->0.9999892916621562

It means that Dixon - EDF method provides a bit more precise solution again.

The result for (Subscript[σ, 3])  in symbolic form,considering the positive root of (Subscript[φ, 4]) = 0 is

solσ3=(Solve[(Subscript[φ, 43])==0,(Subscript[σ, 3])][[2]]//Simplify)

Remark : The result of Dixon-EDF method is not only faster and more elegant but also a bit more precise than that of the iterative Dixon- KSY method. Although, one should check the solutions of  all polynomials with degree 1 and 2 in order to choose the proper result!

4.3 Reduced Groebner basis

The same result can be computed wia reduced Grobner basis, namely

(Subscript[gbσ, 1])=GroebnerBasis[{(Subscript[g, 1]),(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5]),(Subscript[g, 6])},(Subscript[σ, 1]),{a,b,c,(Subscript[σ, 2]),(Subscript[σ, 3])},MonomialOrder->EliminationOrder]

(Subscript[gbσ, 2])=GroebnerBasis[{(Subscript[g, 1]),(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5]),(Subscript[g, 6])},(Subscript[σ, 2]),{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 3])},MonomialOrder->EliminationOrder]

(Subscript[gbσ, 3])=GroebnerBasis[{(Subscript[g, 1]),(Subscript[g, 2]),(Subscript[g, 3]),(Subscript[g, 4]),(Subscript[g, 5]),(Subscript[g, 6])},(Subscript[σ, 3]),{a,b,c,(Subscript[σ, 1]),(Subscript[σ, 2])},MonomialOrder->EliminationOrder]

These polynomials are the same as the polynomials  selected from the result of the Dixon - EDF algorithm

(Subscript[gbσ, 1])-(Subscript[φ, 21])//Simplify

(Subscript[gbσ, 2])-(Subscript[φ, 22])//Simplify

(Subscript[gbσ, 3])-(Subscript[φ, 43])//Simplify

4.4 Symbolic expressions  for a, b and c

Further unknown parameters can be easily computed  from the system of the reduced  equations, (Subscript[g, 1]),...,(Subscript[g, 6]).

Let us express a from the equations (Subscript[g, 5])and (Subscript[g, 6]),

drba=DixonResultant[{(Subscript[g, 5]),(Subscript[g, 6])},{b},{s1}]

{grba}=GroebnerBasis[{(Subscript[g, 5]),(Subscript[g, 6])},{a,b},{b}]

Then a can be computed form

sola=Solve[grba==0,a]//Simplify//Flatten

The parameter  b is given by one of the two equations, let say from (Subscript[g, 5]),

solb=Solve[(Subscript[g, 5])==0,b]//Simplify//Flatten

and parameter c

solc=Solve[(Subscript[g, 1])==0,c]//Simplify//Flatten

4.5 Symbolic expressions of  (Subscript[X, 0]), (Subscript[Y, 0]) and (Subscript[Z, 0]).

The translation parameters can be similarly computed, but now from the original system of equations, (Subscript[f, i]).

{grbX0}=GroebnerBasis[{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])},{(Subscript[Y, 0]),(Subscript[Z, 0])}]

solX0=Solve[grbX0==0, (Subscript[X, 0])]//Simplify//Flatten

similarly we get

{grbY0}=GroebnerBasis[{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])},{(Subscript[X, 0]),(Subscript[Y, 0]),(Subscript[Z, 0])},{(Subscript[X, 0]),(Subscript[Z, 0])}]

solY0=Solve[grbY0==0, (Subscript[Y, 0])]//Simplify//Flatten

and

solZ0=Solve[(Subscript[f, 1])==0, (Subscript[Z, 0])]//Simplify//Flatten

The numerical result using these symbolic results is the following

{tσ1,σ1}=SetPrecision[((Subscript[σ, 1])/.solσ1/.newVarsA/.numericalValues),16]//Chop//Timing

{tσ2,σ2}=SetPrecision[((Subscript[σ, 2])/.solσ2/.newVarsA/.numericalValues),16]//Chop//Timing

{tσ3,σ3}=SetPrecision[((Subscript[σ, 3])/.solσ3/.newVarsA/.numericalValues),16]//Chop//Timing

sσ={(Subscript[σ, 1])->σ1,(Subscript[σ, 2])->σ2,(Subscript[σ, 3])->σ3}

{ta,aa}=SetPrecision[(a/.sola/.newVarsA/.sσ/.numericalValues),16]//Chop//Timing

{tb,bb}=SetPrecision[(b/.solb/.newVarsA/.sσ/.{a->aa}/.numericalValues),16]//Chop//Timing

{tc,cc}=SetPrecision[(c/.solc/.newVarsA/.sσ/.{b->bb}/.numericalValues),16]//Chop//Timing

sabc={a->aa,b->bb,c->cc}

{tX0,X0}=SetPrecision[((Subscript[X, 0])/.solX0/.newVarsA/.sσ/.sabc/.numericalValues),16]//Chop//Timing

{tY0,Y0}=SetPrecision[((Subscript[Y, 0])/.solY0/.newVarsA/.sσ/.sabc/.numericalValues),16]//Chop//Timing

{tZ0,Z0}=SetPrecision[((Subscript[Z, 0])/.solZ0/.newVarsA/.sσ/.sabc/.{(Subscript[X, 0])->X0,(Subscript[Y, 0])->Y0}/.numericalValues),16]//Chop//Timing

Rotation angles in seconds :

CardanS[R/.sabc]

Scale parameters :

SetPrecision[{1/(Subscript[σ, 1]),1/(Subscript[σ, 2]),1/(Subscript[σ, 3])}/.sσ,10]

The total running time of the evaluation of the analytic expressions developed by computer algebra

Apply[Plus,{tσ1,tσ2,tσ3,ta,tb,tc,tX0,tY0,tZ0}]

Let us summarize our results in case of the 3  points problem:

Table 2.
Results in case of  3points Problem
({{Method, Running Time
(sec)}, {Numerical Groebner basis
&  eigensystem method, 0.11}, {Global minimization
 with genetic algorithm, 1.17}, {Linear Homotopy, 0.03}, {Analytic solution
 computed via computer algebra, ~ 0.00}})

Remarks:  Keep in mind, that built-in functions are quicker than functions written in Mathematica own language, which is an interpreter, but built-in functions have a big overhead, sometimes their code can be many hundred pages, therefore they are not always the faster!

Considering Table 2. it is clear that for 3 points problem the best choice is the application of the analytic form. The advantage of the analytic expression is not only the short computation time. The running time of the homotopy and any iterative algorithm may depend the actual values of the coefficients, while in case of the evaluation of analytic expressions the running time practically is constant. In addition these analytic expressions can be directly transformed into C which is  faster with roughly two magnitudes than any interpreter languages.
For example, in case (Subscript[σ, 1])

CForm[(1/(Sqrt[(x23 y13-x13 y23) (x23 z13-x13 z23)]))(√((X23^2) y13 z13+y13 (Y23^2) z13+(X13^2) y23 z23+(Y13^2) y23 z23+y23 (Z13^2) z23-X13 X23 (y23 z13+y13 z23)-Y13 Y23 (y23 z13+y13 z23)-y23 z13 Z13 Z23-y13 Z13 z23 Z23+y13 z13 (Z23^2)))]

We should emphasize again, there is no guarantee, which polynomials give the proper result for a (Subscript[σ, i]). Therefore in case of  any index i, all corresponding terms should be investigated in order to choose the proper polynomial!

5. The proper selection of the 3 points for initial guess values

Using the result of  3 point solution as initial value for the N point problem, it is important to decide which 3 points from the N points are chosen to calculate symbolically the initial guess values. Let us choose two different triplet from a data set containing 81 points of the Hungarian Datum  (Fig. 1)  and calculate with the symbolical solution the parameters of the coordinate transformation ((Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]), a, b, c, (Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])) between the global WGS84 and the Hungarian System.

Figure 1. Two chosen triplet from the 81 points

The results for the two triplet are quite different  (see Table 3.)

Table 3.
Calculated coordinate transformation parameters from two different triplets

({{., 1. triplet, 2. triplet}, {(Subscript[X, 0]), -77.523, -496.192}, {(Subscript[Y, 0]), 90.366, 124.797}, {(Subscript[Z, 0]), 25.151, 543.061}, {a, 1.520·(10^(-6)), 5.470·(10^(-6))}, {b, 1.526·(10^(-6)), 27.692·(10^(-6))}, {c, -0.543·(10^(-6)), -2.415·(10^(-6))}, {(Subscript[σ, 1]), 0.999998865, 0.999957317}, {(Subscript[σ, 2]), 1.000001182, 0.999989413}, {(Subscript[σ, 3]), 1.000001642, 1.000069623}})

Applying these parameters as initial values to a local method, for example for the Newton-Rapshon method, to solve the N point solution (here 81 points), in the case of the first triplet the method converges rapidly in 4 iteration steps, but in the second case the method is not converging, not even after 100 iteration steps.

There is correlation between the geometry of the chosen triplet and the goodness of the calculated initial values. According to our numerical calculations we get the best initial values when the geometry of the triplet is similar to a equilateral triangle, and the worse when the geometry of the three points is nearly a line. A geometrical index can be introduced to represent the geometry of the selected 3 points to avoid the solutions which provide disadvantageous starting values. This geometrical index is the sine of the minimum angle in the triangle, it's maximal value is ((Sqrt[3])/2)when the triangle is an equilateral triangle, and around 0 when the three points are nearly collinear. In the earlier examples the geometrical index in the first case, which gave good initial values, was 0.529 and in the second case was 0.002.
To check the correlation between the geometry and the goodness of the initial values we calculated the transformation parameters for all 3 points combination from the 81 points, so we examined (({{81}, {3}}))= 85320 combinations. 
Let check the (Subscript[Z, 0]) values for all combinations. The real (Subscript[Z, 0]) value calculated with Newton-Rapshon method from the 81 points is 50.342 m, but the values calculated from the different triplets can be very different from this, for example the maximum value for (Subscript[Z, 0]) was 24 679 629 m! The geometrical index of this extreme triplet was 0.003, so these points were almost collinear. 
In Fig. 2 the calculated (Subscript[Z, 0]) values are represented in function of the geometrical index for all of the 85320 combinations. The figure is only a window of the whole graphics, for the better representation \[RightBracketingBar](Subscript[Z, 0])\[RightBracketingBar]>1000 values are not represented (the minimum (Subscript[Z, 0]) is -1 798 501 and the maximum is 24 679 629). In the figure the adjusted value of (Subscript[Z, 0]) is also represented by a line.

Figure 2. Values of (Subscript[Z, 0]) in function of the geometrical index for all combinations

The triplets which give extremely different solutions from the adjusted value all have geometrical indexes less then 0.1. This fact represents that the bad geometry of a triangle leads to ill-conditioned problem generating fatal round-off errors.  In general the more similar is the geometry of the selected three points to a equilateral triangle the better is the initial value. Accordingly to these results it is very important to examine the geometry of the selected three points for calculating the initial values and we should avoid the nearly collinear triplets.

Conclusions

For the 3 point problem the computer algebra method, namely the accelerated Dixon resultant with the technique of Early Discovery Factors as well as reduced Groebner basis provide a very simple, elegant symbolic solution. It is enormously better than using the Dixon Resultant Application developed for Mathematica, which could only succeed by a long and tedious step-by step reduction and even then, the answer given is difficult to work with it.
Although, numerical methods without initial guess values, like linear homotopy and numerical Groebner basis with eigensystem method are also very efficient. The main advantages of the symbolic solution originated from its iteration-free feature, are the very short computation time and the independence on value of  the actual numerical data.
The solution of the 9 parameter affine coordinate transformation problem for 3 points can provide a good inital value for the N point coordinate transformation. Solving this problem, it is not only the choice of the employed method, but the computation of a good initial guess value on the bases of  the  geometry of the properly selected three points are also important (collinear triplets should be avoided).

Acknowledgement

The author are indebted to professor  Daniel Lichtblau (Wolfram Research) for his valueable remarks and suggestions, especially calling up their attention to the proper implementation of NSolve as well as to the application of the reduced Groebner basis.

References

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

Awange J.L.  (2002).: Gröbner bases, multipolynomial resultants and Gauss-Jacobi combinatorial algorithms - adjustment of nonlinear GPS/LPS observations, Dissertation (D93), Geodätisches Institute der Universität Stuttgart

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 J.L. and Grafarend E.W. (2005): Solving Algebraic Computational Problems in Geodesy and Geoinformatics, Springer, Berlin

Bancroft, S. (1985).: An algebraic solution of the GPS equations, IEEE Transaction on Aerospace and Electronic Systems AES-21, pp. 56-59.

Barsi A. (2001): Performing coordinate transfornations by artificial neural network, Allgemeine VermessungsNachrichten 108, pp. 134-137.

Binous H. (2007) Homotopy Continuation Method tos Solve a System of Nonlinear Algebraic Equations. http://library.wolfram.com/infocenter/MathSource/6710/

Cai, J., Grafarend, E. W. (2004): Systematische Analyse der Transformation zwischen Gauß-Krüger-Koordinaten/DHDN und UTM-Koordinaten/ETRS89 angewandt auf Baden-Wurttemberg, Geodatische Woche 2004 - Stuttgart 12.10.-15.10. poster presentation

Drexler, F.J. (1977): Eine Methode zur Berechnung sämtlicher Lösungen von Ploynomgleichungssystemen, Numer. Math. 29. pp.45-58.

Fröhlich  H, Bröker G. (2003): Trafox version 2.1. - 3d-kartesicsche Helmert - Transformation,  http://www.koordinatentransformation.de/data/trafox.pdf

Garcia C.B. and Zangwill, W.I. (1979): Determinung all solutions to certain systems of nonlinear equations, Math. Operations Res. 4. pp. 1-14

Golub, G. H., Pereyra, V. (1973): The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate, SIAM J. Num. Anal. 10, pp. 413-432.

Grafarend E. and Schaffrin, B. (1993): Ausglechungsrechnung in Linearen Modellen, B.I. Wissenschaftsverlag, Mannheim

Grafarend E. and Kampmann G. (1996): (Subscript[C, 10])(3): The tenparameter conformal group as a datum transformation in threedimensional Euclidean space, eitschrift für Vermessungswesen 121,pp. 68-77.

Kapur D, Saxena T. and Yang L. (1994) : Algebraic and geometric reasoning using Dixon resultants. In : ACM ISSAC 94. Oxford, England, July 1994, pp. 99 - 107.

Kleusberg A. (1994): Die direkte Lösung des räumlichen Hyperbelschnitts. Zeitschrift für Vermessungswesen, 119,  pp. 188 - 192.

Kleusberg  A. (2003): Analytical GPS navigation solution. In Grafarend EW., Krumm FW. and Schwarze VS. (eds) Geodesy -the Challenge of the 3rd Millennium, Springer, Heidelberg, pp. 93 - 96.

Lewis, R.H. and Bridgett, S.(2003): Conic tangency equations and Apollonius problems in biochemistry and pharmacology, Mathematics and Computers in Simulation 61. pp. 101-114

Lewis, R. H. (2007 a): Heuristics to accelerate the Dixon resultant, Mathematics and Computers in Simulation /In Press/, Available online 11 April 2007

Lewis, R. H. (2007 b): Computer algebra system Fermat, http://home.bway.net/lewis/

Lichtenegger H. (1995): Eine direkte Lösung des räumlichen Bogenschnitts., Osterreichische Zeitschrift fur Vermessung und Geoinformation 83,  pp. 224 - 226.

Mathes A. (2002): EasyTrans Pro-Edition, Professionelle Koordinatentransformation für Navigation, Vermessung und GIS, ISBN 978-3-87907-367-2, CD-ROM mit Benutzerhandbuch

Nakos G. and Williams R. (2002): A fast algorithm implemented in Mathematica provides one-step elimination of a block of unknowns from a system of polynomial equations, http://wolfram.com/infocenter/Articles/2597

Palancz, B. (2008) Introduction to Linear Homotopy, MathSource, http://library.wolfram.com/infocenter/MathSource/7119/ 

Papp E., Szucs L. (2005): Transformation Methodes of the Traditional and Satellite Based Networks (in Hungarian with English abstract), Geomatikai Kozlemenyek VIII. 85-92. 

Singer P., Strobel D., Hordt R., Bahndorf J., Linkwitz K. (1993): Direkte Losung des raumlichen Bogenschnitts. Zeitschrift für Vermessungswesen, 124, S. 295 - 297.

Späth,H. (2004): A numerical method for determining the spatial Helmert transformation in case of different scale factors, Zeitschrift für Geodäsie, Geoinformation und Landmanagement 129, pp. 255-257.

Volgyesi L, Toth Gy, Varga J. (1996):  Conversion between Hungarian Map Projection Systems. Periodica Polytechnica Civ.Eng., Vo1.40, Nr.1, pp. 73-83. 

Watson , G.A. (2006): Computing Helmert transformations, Journal of Computational and Applied Mathematics 197, pp. 387-395.

Wolfrum, O. (1992): Merging terrestrial and satellite networks by a ten-parameter transformation model, Manuscripta Geodetica 17, pp. 210-214.

Zaletnyik P. (2004): Coordinate transformation with neural networks and with polynomials in Hungary , International Symposium on Modern technologies, education, and professional practice in geodesy and related fields, 4-5 November 2004, Sofia, Bulgaria, pp. 471-479.