Awange J.,Grafarend E., Palancz B., Zaletnyik P. Algebraic Geodesy and Geoinformatics

Подождите немного. Документ загружается.

9-4 The General Procrustes solution 127

Theorem 9.1 (W-LESS of Y

= Y

∗

+ 1x

∗

+ E).

(i) The parameter array {x

, x

, X

} is W-LESS of (9.23) if

trY

∗

WCY

∗

trY

∗

WCY

(9.46)

= (1

∗

W1)

−1

− Y

∗

)

∗

W1 (9.47)

= UV

∗

, (9.48)

subject to the singular value decomposition of the general 3 × 3 matrix

∗

WCY

= UDiag(σ

, σ

∗

, (9.49)

namely



[(Y

∗

WCY

)

∗

WCY

) −σ

I]v

= 0 ∀i ∈ {1, 2, 3}

V =



, v



, VV

∗

= I

(9.50)



U = Y

∗

WCY

VDiag(σ

−1

, σ

−1

, σ

−1

)

∗

U = I

(9.51)

and as well as to the centering matrix

C := I

− (1

∗

W1)

−1

∗

W. (9.52)

(ii) The empirical error matrix of type W-LESS accounts for

= [I

− 11

∗

W(1

∗

W1)

−1

]



− Y

∗

trY

∗

WCY

∗

trY

∗

WCY



, (9.53)

with the related Frobenius matrix W–semi-norm







k E k

:= tr(E

∗

) =

tr{(Y

− Y

∗

trY

∗

WCY

∗

trY

∗

WCY

)

∗

.[I

− 11

∗

W(1

∗

W1)

−1

]

∗

W[I

− 11

∗

W(1

∗

W1)

−1

.(Y

− Y

∗

trY

∗

WCY

∗

trY

∗

WCY

)},

(9.54)

and the representative scalar measure of the error of type W -LES S given by

|k E

tr(E

∗

)/3n. (9.55)

Corollary 9.4 (I-LESS of Y

= Y

∗

+ 1x

∗

+ E ).

(i) The parameter array {x

, x

, X

} is I-LESS of (9.23) if

trY

∗

trY

∗

(9.56)

− Y

∗

)

∗

1 (9.57)

= UV

∗

, (9.58)

128 9 Procrustes solution

subject to the singular value decomposition of the general 3 × 3 matrix

∗

= UDiag(σ

, σ

∗

, (9.59)

namely



[(Y

∗

)

∗

) −σ

I]v

= 0 ∀i ∈ {1, 2, 3}

V =



, v



, VV

∗

= I

(9.60)

U = Y

∗

VDiag(σ

−1

, σ

−1

, σ

−1

)

∗

U = I

(9.61)

and as well as to the centering matrix

C := I

−

∗

. (9.62)

(ii) The empirical error matrix of type I-LESS accounts for

= [I

−

∗

]



− Y

∗

trY

∗

trY

∗



, (9.63)

with the related Frobenius matrix W–semi-norm







k E k

:= tr(E

∗

) =

tr{(Y

− Y

∗

trY

∗

trY

∗

)

∗

.[I

−

∗

.(Y

− Y

∗

trY

∗

trY

∗

)}

(9.64)

and the representative scalar measure of the error of type I-LESS

|k E

tr(E

∗

)/3n.

In the proof of Corollary 9.4, we only sketch the results that the matrix I

−

∗

idemp otent:







−

∗

)(I

−

∗

) =

= I

−

∗

(11

∗

)

= I

−

∗

n11

∗

= I

−

∗

Solution 9.3 (General Procrustes algorithm).









Step 1

Read : Y





. . .





and





. . .





= Y









Step 2

9-4 The General Procrustes solution 129

Compute : Y

∗

subject to C := I

−

∗









Step 3

Compute : SV D Y

∗

= UDiag(σ

, σ

∗

3-1 | (Y

∗

)

∗

) −σ

I |= 0 ⇒ (σ

, σ

)

3-2

((Y

∗

)

∗

) −σ

I)v

= 0 ∀i ∈ {1, 2, 3}

⇒ V =



, v



right eigenvectors

(right eigencoloumns)

3-3

U = Y

∗

VDiag(

) left eigenvectors

(lef t eigencoloumns)









Step 4

Compute : X

= UV

∗

(rotation)









Step 5

Compute : x

trY

∗

trY

∗

(dilatation)









Step 6

Compute : x

− Y

∗

)

∗

1 (translation)









Step 7

Compute : E

= C



− Y

∗

trY

∗

trY

∗



(error matrix)

“optional control”

:= Y

− (Y

∗

+ 1x

∗

)









Step 8

Compute : k E

tr(E

∗

) (error matrix)









Step 9

Compute : |k E

tr(E

∗

)/3n (mean error matrix)

130 9 Procrustes solution

9-5 Extended general Procrustes solution

9-51 Mild anisotropy in scaling

The 3D aﬃne transformation is one possible generalization of the C

(3,3) Helmert

transformation, using three diﬀerent scale (s

, s

) parameters instead of a single

one. In this case the scale factors can be modeled by a diagonal matrix (S).









= S ∗ R

















. (9.65)

Recently, an extension of general Procrustes method for solving the 3D aﬃne

9-parameter transformation has been developed by Awange et al. [42] and is sum-

marized as follows:

1) Compute the center of gravity of the two systems

xyz (9.66)

XYZ (9.67)

2) Translate the systems into the center of the coordinate system

xyzC = xyz − 1 ∗ a

(9.68)

XYZC = XYZ − 1 ∗ b

(9.69)

3) Compute matrix A

A = XYZC

∗ xyzC (9.70)

4) Compute the rotation matrix R via SVD decomposition of A, getting the

matrices U, V and diag (λ

, λ

), where

A = Udiag(λ

, λ

and

R = U ∗ V

5) First approximation of the scale matrix considers only the diagonal elements

of,

S =



∗ XYZC

∗ XYZC ∗ R



−1

∗ R

∗ XYZC

∗ xyzC (9.71)

i.e.,

= diag(S) (9.72)

5) Iterate to improve the scale matrix

9-5 Extended general Procrustes solution 131

w = xyz − XYZ ∗ R ∗ S

− 1 ∗ (a

− b

∗ R ∗ S

) (9.73)

dS = diag





∗ XYZ

∗ XYZ ∗ R



−1

∗ R

∗ XYZ

∗ w



(9.74)

= S

+ dS (9.75)

6) After the iteration the scale matrix is

S = diag (S

) (9.76)

7) Then the translational vector

XYZ

= a

− b

∗ R ∗ S. (9.77)

9-52 Strong anisotropy in scalling

The ABC method is very fast, and the computation eﬀort for iteration to improve

the S scale matrix is negligible. The main problem with ABC method is that after

shifting the two sys tems into the origin of the coordinate system, the rotation

matrix is computed via SVD. This means that the physical rotation as well as

the deformation (distortion) caused by diﬀerent scaling in the diﬀerent directions

of the 3 diﬀerent 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 diﬀer from each other considerably. In addition, the necessary condition

for getting optimal scale matrix is not restricted for diagonal matrix, and the oﬀ

-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 (see Pal´ancz et al. [324]) can cure these problems, but there is

a price for it! This method start with an initially guessed diagonal S scale matrix,

and ﬁrst 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 follows:

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

S =





0 0

0 s

0 0 s





(9.78)

The steps for computing the error in case of a given scale matrix S are the follow-

ing.

2)Using equations (9.66) and (9.67) compute the center of gravity of the two sys-

tem.

132 9 Procrustes solution

3) Translate the system using (9.68) and (9.69)

4) Eliminate rotation (distortion) caused by scaling with diﬀerent values in diﬀer-

ent directions in the computation of matrix A using

A = XYZC

∗ xyzC ∗ S

−1

(9.79)

5) Compute the rotation matrix using SVD decomposition of A, getting the ma-

trices U, V and diag (λ

, λ

), where

A = Udiag(λ

, λ

and

R = U ∗ V

where A = Udiag(. . .)V(see, e .g., 9.59).

6) Compute the translation vector using (9.77).

7) Compute the error matrix

E = xyz − XYZ ∗ R

∗ S − 1 ∗XYZ

(9.80)

8) Compute the square of the norm of the mean error matrix

|k E

tr(E

∗

)/3n (9.81)

If this error is to high, then modify the scale matrix S, and repeat the computation.

To do that reasonably, this computation is implemented a Newton- algorithm to

carry out local minimization of the Frobenius norm of the error matrix. In order

to illustrate the prope r application area of the two extended Procrustes methods,

let us consider two numerical examples. In the 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 [8].

The center of the gravity of the systems are,









and b









(9.82)

Then the estimated scale parameter according to Albertz and Kreimlig [8] in the

Helmert similarity transformation is

priori

i=1

− x

)

+ (y

− y

)

+ (z

− z

)

i=1

− X

)

+ (Y

− Y

)

+ (Z

− Z

)

(9.83)

In case of the 9 parameter aﬃne transformation where 3 diﬀerent scale values

, s

). are applied according to the 3 coordinate axes, a good approach for

the scale parameters can b e given by modifying the Albertz-Kreiling [8] expression.

Instead of the quotient of the two lengths in the c entre of gravity system we can

9-6 Weighted Procrustes transformation 133

use the quotients of the sum of the lengths in the corresponding coordinate axes

directions.

The estimated scale parameters according to the modiﬁed Albertz-Kreiling [8]

expression for the 9 parameter transformation (see also [430]) is

i=1

− x

)

i=1

− X

)

, s

i=1

− y

)

i=1

− Y

)

, s

i=1

− z

)

i=1

− Z

)

(9.84)

where (s

, s

) are the initial elements of S in the iteration.

9-6 Weighted Procrustes transformation

As already stated in Sect. 17-1, other than the problems asso ciated with lineariza-

tion and iterations, the 7-datum transformation problem (conformal group C

(3))

is compounded with the problem of incorporating the weights of the systems in-

volved. This section presents Procrustes algorithm II; a reliable means of solving

the problem. We have already seen that the problem could be solved using Gauss-

Jacobi combinatorial algorithm in Sect. 17-23. Procrustes algorithm II presented

in the preceding section oﬀers therefore an alternative that is not computationally

intensive as the combinatorial method.

To obtain the weight matrix used in Corollaries 9.1, 9.2 and 9.3 in the weighted

Procrustes problem, we proceed via the variance-covariance matrix of Theorem

9.2, whose proof is given in Solution 9.4 where the dispersion matrices of two sets

of coordinates in R

are presented in (9.86) and (9.88). They are used in (9.90) and

(9.91) to obtain the dispersion of the error matrix E in (9.92). In-order to simplify

(9.92), we make use of Corollary 9.5 adopted from [180, Appendix A, p. 419] to

express vec x

∗

as in (9.93) and substitute it in (9.92) to obtain (9.94). We

provide as a summary the following:

Theorem 9.2 (variance-covariance matrix). Let vec E

∗

denote the vector val-

ued form of the transposed error matrix E := Y

− Y

∗

− 1x

∗

. Then



v ecE

∗

= Σ

v ec Y

∗

+ (I

⊗ x

)Σ

v ec Y

∗

⊗ x

)

∗

−2Σ

v ecY

∗

,(I

⊗x

)v ecY

∗

(9.85)

is the exact representation of the dispersion matrix (variance-covariance matrix)

v ec E

∗

of vec E

∗

in terms of dispersion matrices (variance-covariance matrices)

v ec Y

∗

and Σ

v ec Y

∗

of the two coordinates sets vec Y

∗

and vec Y

∗

as well as of

their covariance matrix

v ecY

∗

,(I

⊗x

)v ecY

∗.

Proof. By means of Solution 9.4 we deﬁne the dispersion matrices, also called

variance-covariance matrices, of vec Y

∗

and vec Y

∗

of the two sets of coordinates.

134 9 Procrustes solution

Solution 9.4 (Dispersion matrices of two sets of coordinates in R

vec Y

∗

= E{







− E{x

}

− E{x

}

...

− E{x

}









− E{x

})

∗

. . . (x

− E{x

})

∗



} (9.86)

E{(vecY

∗

− E{Y

∗

})(vecY

∗

− E{Y

∗

})

∗

} (9.87)

vec Y

∗

= E{







− E{X

}

− E{X

}

...

− E{X

}









− E{X

})

∗

. . . (X

− E{X

})

∗



} (9.88)

Next from the transposed error matrix

∗

:= Y

∗

− (x

∗

+ x

∗

) (9.89)

we compute the dispersion matrix (variance-covariance matrix) Σ

vec E

∗







vec E

∗

:= E{[vec E

∗

− E{vec E

∗

}][vec E

∗

− E{vec E

∗

}]

∗

}

= E {[vec Y

∗

− E{vec Y

∗

} −x

(vec X

∗

− E{vec X

∗

})]

× [vec Y

∗

− E{vec Y

∗

} −x

(vec X

∗

− E{vec X

∗

})]

∗

}}

(9.90)







vec E

∗

= E{[vec Y

∗

− E{vec Y

∗

][vec Y

∗

− E{vec Y

∗

]

∗

}

E{[vec X

∗

− E{vec X

∗

}][vec X

∗

− E{vec X

∗

]

∗

−2x

E{[vec Y

∗

− E{vec Y

∗

][vec X

∗

− E{vec X

∗

]

∗

}

(9.91)



vec E

∗

= Σ

vec Y

∗

+ Σ

vec x

∗

− 2Σ

vec Y

∗

,vec x

∗

(9.92)

Corollary 9.5.



vec AB = (I

⊗ A)vec B for all A ∈ R

n×m

, B ∈ R

m×q

vec X

∗

= (I

⊗ x

)vec Y

∗

for all X

∈ R

3×3

, Y

∗

∈ R

3×n

(9.93)

As soon as we implement the Kronecker-Zehfuß decomposition of the vec AB we

arrive at the general representation of the dispersion matrix Σ

v ec E

∗

, namely



v ec E

∗

= Σ

v ec Y

∗

+ (I

⊗ x

)Σ

v ec Y

∗

⊗ x

)

∗

−

−2Σ

v ecY

∗

,(I

⊗x

)v ecY

∗

(9.94)

♣

The results of Theorem 9.2 are interpreted in more detail as follows: The

variance-covariance matrix of the vectorized error matrix E

∗

depends on;

9-7 Concluding remarks 135

(i) the variance-covariance matrix Σ

v ec Y

∗

of the local coordinate set (x

, y

, z

. . ., x

, y

, z

(ii) the variance-covariance matrix Σ

v ec Y

∗

of the global coordinate set (X

, Y

, Z

. . ., X

, Y

, Z

(iii) the covariance matrix between vec Y

∗

and (I

⊗ x

)vecY

∗

of the global

coordinate set vec Y

∗

as well as

(iv) the nonlinearity of the parameter model on the unknowns x

, X

of type

“scale factor” and “rotation matrix” coupled to (I

⊗ x

So as to take advantage of the equivalence theorem between least squares approxi-

mation and best linear uniformly unbiased estimation, e.g., [180, §3, pp. 339–340],

which holds for linear Gauss-Markov model, it is tempting to identify the weight

matrix W of W-LESS with Σ

−1

v ec E

∗

shrunk to a locally isotropic error situation.

Such a shrinking procedure is outlined in Example 17.3, namely by taking in ac-

count isotropic, but inhomogeneous criterion matrices.

9-7 Concluding remarks

The partial Procrustes algorithm presented in this chapter provides a powerful tool

for solving rotation and orientation related problems in general. The approach is

straight forward and does not require linearization, which bog down least squares

and other techniques commonly used. In Chap. 17, it shall be demonstrated how

the general Procrustes approach determines scale and translation parameters of

transformation, in addition to the rotation elements.

The 9-parameter Procrustean algorithm considered in this Chapter can thus

be used for

(i) quicker and eﬀective generation of 9 transformation parameters given coordi-

nates in two systems as matrix conﬁguration,

(ii) quick checking of the transformation parameters obtained from other methods

(iii) generating three scale parameters which could be useful in correcting dis-

tortions following procedures which ﬁrst determine the rotation and translation

parameters independent of scale.

For complete exposition of Procrustes approach, we refer to the works of [79,

86, 87, 93, 110, 111, 116, 120, 121, 130, 151, 153, 157, 191, 197, 198, 293, 294, 333,

355, 356, 378, 397].

Part II

Applications to geodesy and geoinformatics