[SciPy-User] helmert implementation (7 parameters - geometric transformation)
Massimo Di Stefano
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.
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Charles R Harris
Charles R Harris
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
Charles R Harris
Are you sure the data you gave is correct? Lengths don't seem to be preserved under rotation.
Chuck
Ewald Zietsman
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.
>>> ? 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?
>>
Massimo Di Stefano
Charles R Harris
I've attached a script for the computation. Two things of note: the angles in the pdf results are in degrees, the value of μ is pretty much in the noise (the condition number of the matrix is 10^10), and the fit you get with my version differs from the data by about 7mm. Note that the Helmert parameters result from a linearization of the rotation matrix, see Gauss-Newton for more about linearization and least squares, but in this case I think the limiting factor is the accuracy of the coordinates, not the linearization.
<snip>
Chuck
Pauli Virtanen
> 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 ?
Cholesky decomposition works only on Hermitian matrices -- so being
square is one requirement.
Also, if you are going to do least-squares fitting, use linalg.lstsq; it
does the same thing as (A^T A)^-1 A^T b, but should be numerically more
stable.
--
Pauli Virtanen