[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116 IC™ Value: 3.00 CODEN : IJESS7 http: // www.ijesrt.com © International Journal of… [615885]

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[355]

IJESRT
INTERNATIONAL JOURNA L OF ENGINEERING SCI ENCES & RESEARCH
TECHNOLOGY
THREE -DIMENSIONAL TRANSFORMATION OF COORDINATE SYSTEMS
USING NONLINEAR ANALYSIS – PROCRUSTES ALGORITHM
Corina Daniela P ăun *, Valeria Ersilia O niga, Petre Iuliu D ragomir
* Department of Topography and Cadastre , Technical University of Civil Engineering, Bucharest,
Faculty of Geodesy , Romania

DOI : 10.5281/zenodo.291839

ABSTRACT
The coordinates transformation between two projection systems is increasingly approached in geodesy,
photogrammetry and computer vision. This issue led to widespread development of algorithms that offer high
accuracy and fast processing time. One of this algorithms was pr esented in this paper, namely the Procrustes
algorithm, which involves the usage of seven parameters of transformation (three of translation, three of rotation
and one scalar). The main advantages of this algorithm is that you don’t have to know the initia l values of the
parameters, like in the case of iterative numerical methods and the process of equations linearization in order to
obtain the equation correction system, isn’t necessary. So, to the standard equation of 3D – transformation, the
Lagrange func tion is applied, together with some constraints on rotation matrix (R must be orthogonal) and
some derivation conditions. In order to obtain more accurate results in the 3D transformation of coordinates the
weights matrix is introduced, which is calculated separately based on variance -covariance matrix, using the
error matrix. This paper also investigates the stability of using this algorithm for the registration of two TLS
point clouds, by comparing the results with those obtained by applying the 3D confo rmal transformation. When
using the last algorithm, the initial values of the 7 – parameters were calculated and an iterative optimization of
these parameters was applied.

KEYWORDS : 2 Procrustes algorithm, 7 – parameters transformation, TLS, registration, comparison.

INTRODUCTION
The tridimensional coordinates transformation is a process that have frequently been used in geodesy, survey
engineering, but also in photogrammetry. The determination of transformation parameters requires knowledge
of a minimum n umber of common points in the two reference systems.
Until now, researchers have been developed and presented a large number of algorithms to determine these
parameters that can be divided into two categories: numerical iterative algorithm and analytical algorithm [1].
In the case of algorithms belonging to the first category, knowledge of initial value (approximate) parameters,
correction equations linearization and iterative calculation, are necessary. In the case when rotation angles are
large, the appr oximate values of the parameters are difficult to determinate or impossible in some cases, leading
to the impossibility of method implementation (method failure) [1].
The algorithms belonging to the second category does not involve knowledge of initial va lues for the
parameters, equation linearization or iterative calculation, so that the solutions are calculated quickly and
accurately. But these algorithms have a small drawback, namely the high degree complexity of mathematical
derivation, therefore their use and implementation was done fewer times.
One of these algorithms has been implemented by [2], called Procrustes algorithm, and further was improved by
[1] introducing the technique of singular values decomposition (SVD).
The singular values decomposi tion is one of the most complex methods used in linear algebra, specifically
matrix calculation. Using this method many problems were solved, such as calculating the rank of a matrix,
orthogonal bases for the linear subspace, the least squares, or in this case, determining the values of rotation
matrix.

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[356]

MATERIALS AND METHODS
The Classic Procrustes Algorithm
The coordinates transformation involves the calculation of seven parameters, three of translation, three
of rotation and one scalar and it is given by the following mathematical formulation:
1xyz n XYZ R xyz T E  
(1)
where : (X,Y,Z) and (x,y,z) are the coordinates of a point P in the reference coordinate system, respectively
arbitrary coordinate system, Txyz is the vector of the three translation parameters, R is the rotation matrix, λ is a
scalar, 1n=[1,1,1,..1]T is a vector with n components and all elements have a value of 1 – unity vector, E is the
error transformation matrix (we added this term because we have supposed that the points have errors).
The rotation matrix is orthogonal and for this reason the following conditions must be met [3]:
3, det( ) 1TR R I R 
(2)
In other words, R is calculated based on three rotation angles (
,,   ) around coordinate axes X,Y
and Z and has the following form:

cos cos cos sin sin sin cos sin sin cos sin cos
cos sin cos cos sin sin sin sin cos cos sin sin
sin sin cos cos cosk k k k k
R k k k k k      
      
     
   
  
(3)

If the elements of this matrix are known, the values of the three rotation angles can be determinate,
using the following formulas:
1 1 1 32 21
31
33 11tan , sin ( ), tanr rrrr        
 
(4)
For the first equation, the minimum condition (least squares) was apply, which will resolve based on
Lagrangian multiplier matrix
 . Taking into account the constraints given by Eq. 2, Lagrange function will
write, depending on the terms of equation 1 ( tr-trance ):
3 ( , , , ) ( ) ( ( )) minTTL T R tr EE tr R R I     
(5)
The method of Lagrangian multiplier, solves the constrained optimization problem (Eq.1) by
transforming it into a non -constrained optimization (see Eq. 6). So, we compute the derivates of Lagrangian
function (Eq. 8) and the result is the optimum values of ( λ, T, R ). In other words, by introducing this function the
linearization of correction equations will not be necessary attending to solve the system.
So, the expression of E will be replaced (written on eq. 1) in eq. 5 and will get:
3
3( , , , ) ( ) ( ( ))
(( 1 )( 1 ) ) ( ( )) minTT
TT
nnL T R tr EE tr R R I
tr A RB T A RB T tr R R I
     
        
(6)
For an easier writing of the above equation, will make the following notations:
[ ]; [ ]i i i i i i A X Y Z B x y z
(7)
The Lagrangian function exists if the followings condition are satisfied:
0; 0; 0; 0.L L L L
TR          
(8)
Based on the above conditions, the following equations results:
3
3( , , , ) (( 1 )( 1 )) ( ( ))
(( )( ) ( )1 1 ( )
1 1 ) ( ( )) minT T T T T T
nn
T T T T T T T T
nn
T T T
nnL T R tr A RB T A B R T tr R R I
tr A RB A B R A RB T T A B R
T T tr R R I  
            
      
    
(9)
( )1 ( )1 2 1 1 2( )1 2 1 1 0T T T T T
n n n n n n nLA RB A RB T A RB TT          
(10)

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[357]

Eq. 10 will give the translation vector, which is given by the following equation:
1(1 1 ) ( )1TT
n n n T A RB 
(11)
It is noted that the vector T depends on the scalar and the rotation matrix R (as expected). Next, Eq.11
is replaced in the translation vector in Eq. 1, obtaining:
11(1 1 ) ( )1 1 ( )( (1 1 ) 1 1 )T T T T
n n n n n n n n n E A RB A RB A RB I        
(12)
where:
1 1(1 1 ) 1 1 1 1T T T
n n n n n n n nIIn   , is called the centering matrix .
Further, for an easier writing of the equation, the following notations are done:
11( 1 1 ), ( 1 1 )TT
n n n n n n A A I B B Inn     
(13)
Eq. 13 are the centralized coordinate matrix of two coordinate systems, reduces the coordinates to the
centering matrix . Eq. 12 will be rewritten as:
E A R B  
(14)
Replacing Eq. 14 into Eq. 6, we obtain:
2( , , ) (( )( ) ( ( ))
( 2 ) ( ( ))TT
T T T T T TL R tr A R B A R B tr R R
tr A A A B R R B B R tr R R  
         
          
(15)
The derivation of Lagrange function (Eq.15) shall be conform matrix operations, e.g. the matrix trance:
( ) (( ) ) ( )T T T T T Ttr A B R tr A B R tr R B A       
.
So, the function L is derived based on scalar λ, resulting:
2 ( ) 2 ( ) 2 ( ) 2 ( ) 0T T T T T T T Ltr A B R tr R B B R tr A B R tr B B             
(16)
From Eq. 16, the value of scalar can be determined:
()
()TT
Ttr A B R
tr B B
(17)
Next, the function L will be derived based on the rotation matrix, using the Eq. 15:
22( )( ) ( ) 2 2 2 0T T T T LA R B B R A B R B B RR                  
(18)
The rotation matrix is given by the following equation:
21()TTR A B B B     
(19)
By derivation of Lagrange function (
( , , )LR ), based on Lagrangian multiplier matrix (
 ), we
get the following equation, that represent the constraint required at the beginning of the algorithm:
30TR R I
(20)
Replacing the Eq. 19 in Eq. 20, one gets:

2 2 1 2 1
3 ( ) ( )T T T TB B B A A B B B I            
(21)
1
2 2()T T TB B B A A B      
(22)
From Eq. 19 and 22, we obtained the formula for the calculation of the rotation matrix coefficients:
1
2()T T TR AP B BP A AP B     
(23)
After some notations, the rotation matrix becomes:
1
2()TR D D D
, where:
TD A B  (24)

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[358]

The Weight Procrustes Algorithm
The determination of the weight matrix
For determining the weight matrix of each point, the determination of the variance -covariance matrix is
necessary.
In order to determine the dispersion for each set of coordinates, the generation of variance -covariance
matrix has to be done. So, the following steps are necessary:
-compute the errors matrix:
( 1 )xyz n E A R B T  
(25)
-compute the mean for each set of coordinates, for both systems (reference and arbitrary system);
1 1 1;;n n n
i i i
i i i
AX AY AZX Y Z
M M Mn n n       
              
(26)
1 1 1;;n n n
ii
i i i
Bx By Bzx y z
M M Mn n n       
              
(27)
– compute the dispersion;
1
1
1 1 1
1…
… .. .. ..
…TAX
i AX
AY
T T T T T T
AX i AX AY i AY AZ i AZ A
i AY
AZ
i AZXM
XM
YM
D M X M X M Y M Y M Z M Z M
YM
ZM
ZM  



         
(28)
1
1
1 1 1
1…
… .. .. ..
…TAx
i Ax
Ay
TT T T T T
Ax i Ax Ay i Ay Az i Az B
i Ay
Az
i AzxM
xM
yM
D M x M X M y M y M z M z M
yM
zM
zM  



         

(29)
-compute the variance -covariance matrix (errors dispersion matrix)
   
     [ ][ ]
[ ( )] [ ( )]T T T T T T
vec
T T T T T T T T TD E M vecE M vecE vecE M vecE
M vecA M vecA vecRB M vecRB vecA M vecA vecRB M vecR B    
       
(30)
       
   2[ ][ ] [ ][ ]
2 [ ][ ]T T T T T T T T T T T
vec
T T T T TD E M vecA M vecA vecA M vecA M vecRB M vecRB vecRB M v ecRB
M vecA M vecA vecRB M vecRB
      
  
(31)
2 ( , )T T T T T
vec vec vec vecD E D A D RB D A vec RB    
(32)
where:
1 1 1[ , , ,…, , , ]T
n n n vecA X Y Z X Y Z and
1 1 1[ , , ,…, , , ]T
n n n vecB x y z x y z transforms the transpose matrix of the points
coordinates in both systems, into a vector with 3
1n dimension (n is the number of common points).
In Eq. 32, the following developments are done:
– it is known that [4]:
( ) , ,n m m q
q vecAB I A vecB for A B   
(33)
So based on the above formula, vecRBT becomes:
3 3 3( ) , ,T T T n
n vecRB I R vecB for R B    
(34)
Further, the Eq. 34 is replaced into Eq. 32 and the variance -covariance matrix becomes:
( ) ( ) 2 ( ,( ) )T T T T T T
vec vec n vec n vec nD E D A I R D B I R D A I R vecB         
(35)

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[359]

Eq. 35 can be interpreted in detail:
T
vecDA and
T
vecDB are the dispersion matrix (variance -covariance)
of the coordinates sets from the reference and arbitrary system, respectively;
( ,( ) )TT
vec nD A I R vecB  is the
dispersion computed between vector
TvecA and the term
()T
nI R vecB given by Kronecker -Zehfuss
determination.
After performing the specific operations, the dispersion matrix of errors E is ob tained, having the size
of
33nn . Within this matrix, only the elements on the main diagonal are taken into consideration.
1
1
12
2
2
2
2
2. . . . . .
n
n
nx
y
y
T
vec
x
x
xDE















(36)
From the above matrix, using only the elements of the main diagonal the dispersion of each point is
calculated like a geometric mean:
2 2 2 2, 1,
i i i i x y y in       
(37)
Finally, after performing all operations, the error matrix E result (size
nn ),with elements only on the
principal diagonal:
2
1
2
2
2. . . .
nE







(38)
The weight matrix P is equal to the inverse error matrix E.
1
2 1
. . . .
np
pPE
p




(39)
The transformation parameters calculation by using the weight Procrustes algorithm
The determination of seven transformation parameters will be achieved by applying weight Procrustes
algorithm, using the weight matrix calculated in previous section.
The weight Procrustes algorithm is based on Lagrange function, that solves the problem of the
minimum condition, without requiring the linearization of the equation system.
The equations used to calculate the weighted transformation parameters, are the same as in case of the
classical algorithm, noting that the weights matrix will be introduce.
3 ( , , , ) ( ) ( ( )) minTTL T R tr EPE tr R R I     
(40)
By applying a equations from section 2 (Eq. 9 – Eq. 22 ) and a derivation conditions, (Eq. 8) the
translation vector, the scalar and the rotation matrix are given by the following equations:
1(1 1 ) ( ) 1TT
n n n T P A RB P 
(41)
()
()TT
Ttr AP B R
tr BP B
(42)
21()TTR AP B BP B     
(43)
If the Lagrange function is derived (
( , , )LR ) according to Lagrangian multiplier matrix (
 ),the constraint
equation is obtained:
30TR R I
(44)
Substituting the Eq. 44 into the rotation matrix equation (Eq. 43), we obtain:

2 2 1 2 1
3 ( ) ( )T T T TBP B BP A AP B BP B I            
(45)
1
2 2()T T TBP B BP A AP B      
(46)

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[360]

Based on Eq.43 and Eq. 46, the equation for calculating the coefficients of the rotation matrix is
obtained:
1
2()T T TR AP B BP A AP B     
(47)
After some notations the R matrix becomes:
1
2()TR D D D
, where:
TD A B  (48)
The rotation angles can be computed later as in the case of the classic algorithm, using the Eq. 4 or more simply,
by applying the singular values decomposition algorithm.

THE CASE STUDY
For this case study, we considered the three dimensional coordinates of five points, manually chose and
measured from two point clouds acquired with ScanStation 2 terrestrial laser scanner, representing natural tie
points (Fig. 1). The points coordinates were measured after the two point clouds were represented in the same
coordinate system, namely the one of the left point cloud and can be found in [5]. The first point cloud, colored
in red, is considered as reference and the second one, colored in blue, as arbitrary. The tie points distribution is
random, but were chosen so as to cover approximately the entire point cloud.
The ScanStation2, produced by Leica Geosystems, is a terrestrial scanner system using laser -pulsed
technology for distance measuremen t with a precision of 6 mm at a distance of 50 m, which has an integrated 1
MP resolution digital camera [5].

Figure 1. Graphical representation of the TLS point clouds in the left point cloud coordinate system and the
five tie points pairs [ 5]
In this section, we present the classical and the weighted Procustes algorithms in the process of TLS
point clouds registration.
In the first stage, the seven transformation parameters were computed (Table 1), the place of the weight
matrix being taken by the iden tity matrix. So, the dispersion matrix of errors was calculated using the scalar and
the rotation matrix values. The resulting dispersion matrix with a size of 3n x3n, will be the errors matrix for
these 5 points and only the elements from the main diagona l will be considered (section 3.1). Using the Eq. 42,
the root mean square error (RMS) for each of the five points was calculated, the result being a diagonal matrix
with a size of n x n , called error matrix E.
By applying the inverse of this matrix, the weights specify to each given point will be calculated. In
this case, every point has isotropic weight and is independent of each other [2], weight matrix is generated and is
given in Table 2.
In the second stage, the weight Procrustes algorithm for the tr ansformation parameters determination
using the weights of each point (section 3.2), was computed. The new set for the seven transformation
parameters is presented in Table 3. For a comparison purpose, we have also introduced the results of the 3D
conforma l transformation obtained in [5] (Table 5).

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[361]

Table 1.The transformation parameters calculated by using the Classic Procrustes Algorithm
Values
Rotation 0.38268 0.92388 -0.0022594
matrix -0.92388 0.38268 0.00089502
R 0.0016915 0.0017449 1.0000
Rotation ω -0.051281
angle φ -0.12945
k -67.5
Translation
Tx
Ty
Tz -19.896
21.22
-3.8812
Scale 1.0007
Mean error 0.044124

Table 2.Weights matrix
0.040119 0 0 0 0
0 0.089652 0 0 0
0 0 0.081375 0 0
0 0 0 0.093759 0
0 0 0 0 0.000458

Table 3. The transformation parameters calculated by using the Weight Procrustes Algorithm
Values
Rotation 0.38287 0.9238 -0.0020468
matrix -0.9238 0.38287 0.00011697
R 0.00089171 0.0018461 1.00000000
Rotation ω -0.041807
angle φ -0.14975
k -67.564
Translation
Tx -19.879
Ty 21.243
Tz -3.8737
Scale 1.002
Mean error 0.012523

In the success of this algorithm, the Eq. 2 carries a very important role, because it is a necessary
condition to be met by the rotation matrix.
Both Procustes algorithms, classical and weighted, were implemented into Matlab programming
language, the implementation and processing times being much faster than any other 3D coordinate
transformation method, i.e. 3D conformal transformation.
In order t o establish the success and the stability of the Procustes algorithm in the process of TLS point
clouds registration, the 7 – transformation parameters obtained by using both methods (classical and weighted),
will be compared with those obtained by using th e 3D conformal transformation.
The 3D conformal transformation is one of the iterative methods, which involves performing a
preliminary determination of the 7 – transformation parameters (initial values) and an iteratively solve of the
unknowns (the paramet ers corrections).
In Table 4, the residual of tie points coordinates when applying the CP and WCP, respectively are
presented, as well as the differences between them. It can be seen that the differences are quite small, of
centimeters or even millimeters.

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[362]

Table 4. Residual of tie points coordinates when applying the CP and WP transformation and the differences
between them
N
o Residual of tie points
coordinates when applying the
Classic – Procrustes algorithm Residual of tie points
coordinates when applying the
Weight – Procrustes algorithm Differences
X(m) Y(m) Z(m) X(m) Y(m) Z(m) diffX(m) diffY(m) diffZ(m)
1 0.03730 -0.01477 -0.01384 0.02162 0.00376 -0.02000 0.01569 -0.01853 0.00616
2 -0.03019 0.06842 -0.01171 -0.02200 0.05775 -0.01526 -0.00819 0.01067 0.00355
3 0.02944 0.00920 0.01463 -0.01169 -0.00168 0.00259 0.04113 0.01088 0.01204
4 0.02000 -0.04846 0.02371 0.02193 -0.05537 0.02090 -0.00194 0.00691 0.00281
5 -0.05654 -0.01440 -0.01278 -0.10234 -0.02113 -0.02866 0.04580 0.00673 0.01587

The differences between the tie points residuals obtained by applying the CP and WCP algorithms,
shows that weighted Procrustes transformation provides better results, because the residuals of tie points
coordinates are smaller and the root mean squared er ror when applying this algorithm has a values of 1.25 cm ,
compared to the value of 4. 41 cm achieved when applying the classical algorithm.
The values of the transformation parameters, obtained by applying the three methods of transformation
described above, i.e. 3D Conformal, Classical and Weighted Procrustes algorithms, as well as the differences
between them, are given in Table 5.

Table 5. Parameters of 3D conformal transformation, Classical Procrustes and Weighted Procrustes
respectively and the di fferences between them
Parameter Values obtained
by 3D conformal
transformation Values
obtained by
Classical
Procrustes Values
obtained by
Weighted
Procrustes Differences
3D Conformal –
CP Differences
3D Conformal –
WP
λ 1.000675 1.0007 1.002 -0.00002 -0.00133
ω -0.051281o -0.051281o -0.041807o 0.00000o -0.00947o
φ -0.129454o -0.12945o -0.14975o 0.00000o 0.02030o
k -67.500083o -67.5o -67.564o -0.00008o 0.06392o
TX -19.896 m -19.896 m -19.879 m 0.00000 m -0.01700 m
TY 21.220 m 21.220 m 21.243 m 0.00000 m -0.02300 m
TZ -3.881 m -3.8812 m -3.8737 m 0.00020 m -0.00730 m

CONCLUSION
In this paper, the steps for 3D coordinate transformation between two coordinate systems when applying the
classical and the weighted Procrustes algorithms, were presented. Then, an analysis of their potential in the
process of TLS point clouds registration, based on pairs of common points, randomly chosen and evenly
distributed in the point clouds, was performed. So, the 7 – transformation parameters were calculate d by using
both methods, classical and weighted Procustes algorithm and the results were compared to those obtained by
using the 3D conformal transformation method.
We conclude by saying that the results contain quite small differences and the goal was ach ieved. In other
words, Procrustes algorithm has a higher complexity, but if we are looking for a faster method, this is the better
option.

REFERENCES
[1] Huaien Zeng , "Analytical algorithm of weighted 3D datum transformation using the constraint of
orthonormal matrix ", SPRING OPEN JOURNAL EARTH, PLANET AND SPACE , 2015 .
[2] Grafarend, E. W., Awange, J.L. "Nonlinear analysis of the three -dimensional datum transfo rmation
[conformal group C7(3)]", J Geod 77:66 –76, 2003.
[3] Schonemann P. H. – A generalized solu tion of the orthogonal Procrustes problem. Psychometrika, vol.
31, no. 1, USA, 1966.

ISSN: 2277 -9655
[Păun * et al., 6(2): February , 2017 ] Impact Factor: 4.116
IC™ Value: 3.00 CODEN : IJESS7
http: // www.ijesrt.com © International Journal of Engineering Sciences & Research Technology
[363]

[4] Grafarend E, Schaffrin B – Ausgleichungsrechnung in Linearen Modellen. B.I. Wissenschaftsverlag,
Mannheim, 1993.
[5] Oniga, E., Savu, A., Negrila, A., – The Evaluation of Clo udCompare Software in the Process of TLS
Point Clouds Registration, GeoCAD 2016, Alba Iulia, Romania.
[6] S. R. Baig, F. U. Rehman, and M. J. Mughal, “Performance Comparison of DFT, Discrete Wavelet
Packet and Wavelet Transforms in an OFDM Transceiver for Mult ipath Fading Channel,”, 9th IEEE
International Multitopic Conference , pp. 1 -6, Dec. 2005.
[7] N. Ahmed, Joint Detection Strategies for Orthogonal Frequency Division Multiplexing , Dissertation
for Master of Science, Rice University, Houston, Texas. pp. 1 -51, Ap r. 2000.