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

[Charles R Harris]

2010/3/1 Charles R Harris <charlesr.harris at gmail.com>

>
>
> 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?
>
>
Are you sure the data you gave is correct? Lengths don't seem to be
preserved under rotation.

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