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

[Charles R Harris]

On Mon, Mar 1, 2010 at 6:17 PM, Massimo Di Stefano
<massimodisasha at yahoo.it>wrote:

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

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100301/5ca0d197/attachment.html>