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

Massimo Di Stefano massimodisasha@yahoo...
Tue Mar 2 17:29:42 CST 2010


Hi,

from here :

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cholesky.html

seems it is applicable only on a square matrix.
is that true ?

reading this pdf :

http://geomatica.como.polimi.it/corsi/misure_geodetiche/MG_03SistRif2.pdf

scroll to the page : 15

i can see it isn't a common system :  Y = AX  
but it is :  Y = AX + B + V


the points are :



p1:
4402553.569 727053.737 4542823.332 
4399375.518 703845.639 4549214.880 
4412911.336 701094.214 4536517.955


p2 :
4402553.334 727053.937 4542823.474 
4399375.347 703845.876 4549215.105 
4412911.150 701094.435 4536518.139




tring to rewrite the code i did :

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[1::3, 3] = L[:,2]
A[2::3, 3] = -L[:,1]
A[ ::3, 4] = -L[:,2]
A[2::3, 4] = L[:,0]
A[ ::3, 5] = L[:,1]
A[1::3, 5] = -L[:,0]
A[ ::3, 6] = L[:,0]
A[1::3, 6] = L[:,1]
A[2::3, 6] = L[:,2]

G = np.loadtxt(p2)

V = np.zeros((3*G.shape[0],1),float)
V[ ::3, 0] = G[:,0] - L[:,0]
V[1::3, 0] = G[:,1] - L[:,1]
V[2::3, 0] = G[:,2] - L[:,2]


i used the code in the previouse mail :

N = np.dot(A.T.conj(), A)
T = np.dot(A.T.conj(), Y)
C = np.dot(linalg.inv(N), T)

to solve a system like : 

 Y = AX  

how can i change the code to solve : 

 Y = AX + B + V

instead ?

the results in the dpf are :

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 

any suggestion, also on how to apply a decomposition can give me a great help!



Il giorno 02/mar/2010, alle ore 19.00, scipy-user-request@scipy.org ha scritto:

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100303/548a759d/attachment.html 


More information about the SciPy-User mailing list