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

[Charles R Harris]

On Tue, Mar 2, 2010 at 4:29 PM, Massimo Di Stefano
<massimodisasha at yahoo.it>wrote:

> 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 :
>
> *x**0     -9.256 m *
> ***y**0     -23.701 m *
> ***z**0     16.792 m*
> ***R**x     -0.0001990982 rad *
> ***R**y     0.0001778762 rad *
> ***R**z     0.00015 rad *
> **μ    0.00000046
>
> any suggestion, also on how to apply a decomposition can give me a great
> help!
>
>
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100302/113cb4b3/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: helmert.py
Type: text/x-python
Size: 876 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100302/113cb4b3/attachment.py>