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

[Massimo Di Stefano]

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


my results :

[[ -9.91564629e+00]
 [ -2.36172028e+01]
 [  1.57283853e+01]
 [ -2.68063331e-07]
 [  3.04330048e-06]
 [ -2.76148185e-06]
 [ -2.31598150e-06]]



anyone here can give it a try ?
the results seems wrong but
i haveb't yet find a solution.

thanks to ALL!

regards,

Massimo.