[SciPy-User] helmert implementation (7 parameters - geometric transformation)

[Ewald Zietsman]

Hi,

I see your coordinates are very big numbers. These can lead to an
ill-conditioned (AtA) matrix which can lead to a wrong inverse. Try
using cholesky decomposition instead of using linalg.inv. The method
is outlined here http://en.wikipedia.org/wiki/Cholesky_decomposition
and numpy has routines for doing just this. The other option is to
subtract the centre-of-mass from your coordinates i.e. P1 -> P1x -
P1xave, P1y - P1yave etc. Do the same for P2. Calculate the
transformation again. The shifts should now be zero but your rotations
and scale should be statistically independent and correct. The
difference between the two coordinate systems' centres of mass gives
the shifts. Transformations far away from the origin rarely behaves
well.

Good luck.

>>
>>> Hi All,
>>>
>>> i'm tring to implement a function to perform a geometric trasformation
>>> between a set of points in 2 different reference systems
>>>
>>> it is a 7 parameters trasformation (helmert) , based on last square method
>>> it use as input 2 set of x,y,z coordinates in 2 different reference
>>> systems
>>> and give as output 7 parameters needed to performs transformation
>>> from one reference system to the other one.
>>>
>>> http://en.wikipedia.org/wiki/Helmert_transformation
>>>
>>> this function is used in geodesy to "reproject" coordinates between 2
>>> different Datum
>>>
>>> i wrote it as :
>>>
>>> import numpy as np
>>> from scipy import linalg
>>>
>>> def helmert(p1, p2):
>>>    L = np.loadtxt(p1)
>>>    A = np.zeros((3*L.shape[0],7),float)
>>>    A[ ::3, 0] = 1.0
>>>    A[1::3, 1] = 1.0
>>>    A[2::3, 2] = 1.0
>>>    A[ ::3, 3] = L[:,0]
>>>    A[1::3, 3] = L[:,1]
>>>    A[2::3, 3] = L[:,2]
>>>    A[1::3, 4] = L[:,2]
>>>    A[2::3, 4] = -L[:,1]
>>>    A[ ::3, 5] = -L[:,2]
>>>    A[2::3, 5] = L[:,0]
>>>    A[ ::3, 6] = L[:,1]
>>>    A[1::3, 6] = -L[:,0]
>>>    G = np.loadtxt(p2)
>>>    Y = np.zeros((3*G.shape[0],1),float)
>>>    Y[ ::3, 0] = G[:,0] - L[:,0]
>>>    Y[1::3, 0] = G[:,1] - L[:,1]
>>>    Y[2::3, 0] = G[:,2] - L[:,2]
>>>    N = np.dot(A.T.conj(), A)
>>>    T = np.dot(A.T.conj(), Y)
>>>    C = np.dot(linalg.inv(N), T)
>>>    print C
>>>
>>>
>>> from a pdf i find online
>>> it has numerical examples
>>> and i get different results :'(
>>>
>>>
>>> p1 :
>>>
>>> 4402553.334 727053.937 4542823.474
>>> 4399375.347 703845.876 4549215.105
>>> 4412911.150 701094.435 4536518.139
>>>
>>>
>>> p2 :
>>>
>>> 4402553.569 727053.737 4542823.332
>>> 4399375.518 703845.639 4549214.880
>>> 4412911.336 701094.214 4536517.955
>>>
>>>
>>> pdf results :
>>>
>>> x0    -9.256 m
>>> y0    -23.701 m
>>> z0    16.792 m
>>> Rx    -0.0001990982 rad
>>> Ry    0.0001778762 rad
>>> Rz    0.00015 rad
>>> ?      0.00000046
>>>
>>>
>> The angles here look like degrees, not radians. Also, the last factor is
>> pretty definitely wrong. It is given directly as the square root of the
>> ratio of the variance of the two sets of vectors around their means minus 1.
>> For instance
>>
>> sqrt(pt2.var(0).sum()/pt1.var(0).sum()) - 1
>>
>>
>>>
>>> my results :
>>>
>>> [[ -9.91564629e+00]
>>>  [ -2.36172028e+01]
>>>  [  1.57283853e+01]
>>>  [ -2.68063331e-07]
>>>  [  3.04330048e-06]
>>>  [ -2.76148185e-06]
>>>  [ -2.31598150e-06]]
>>>
>>>
>>>
>> And my result agrees with neither ;) BTW, are you going from pt2 -> pt1 or
>> vice versa?
>>