how to perform coordinates affine transformation using python? -

how to perform coordinates affine transformation using python? -

i perform transformation illustration info set. there 4 known points coordinates x, y, z in 1 coordinate[primary_system] scheme , next 4 known points coordinates x, y, h belong coordinate system[secondary_system]. points correspond; illustration primary_system1 point , secondary_system1 point same point have it's coordinates in 2 different coordinate systems. have here 4 pairs of adjustment points , want transform point coordinates primary scheme secondary scheme according adjustment.

primary_system1 = (3531820.440, 1174966.736, 5162268.086) primary_system2 = (3531746.800, 1175275.159, 5162241.325) primary_system3 = (3532510.182, 1174373.785, 5161954.920) primary_system4 = (3532495.968, 1175507.195, 5161685.049) secondary_system1 = (6089665.610, 3591595.470, 148.810) secondary_system2 = (6089633.900, 3591912.090, 143.120) secondary_system3 = (6089088.170, 3590826.470, 166.350) secondary_system4 = (6088672.490, 3591914.630, 147.440) #transform point x = 3532412.323 y = 1175511.432 z = 5161677.111<br>

@ moment seek average translation x, y , z axis using each of 4 pairs of points like:

#x axis xt1 = secondary_system1[0] - primary_system1[0] xt2 = secondary_system2[0] - primary_system2[0] xt3 = secondary_system3[0] - primary_system3[0] xt4 = secondary_system4[0] - primary_system4[0] xt = (xt1+xt2+xt3+xt4)/4 #averaging

...and on y , z axis

#y axis yt1 = secondary_system1[1] - primary_system1[1] yt2 = secondary_system2[1] - primary_system2[1] yt3 = secondary_system3[1] - primary_system3[1] yt4 = secondary_system4[1] - primary_system4[1] yt = (yt1+yt2+yt3+yt4)/4 #averaging #z axis zt1 = secondary_system1[2] - primary_system1[2] zt2 = secondary_system2[2] - primary_system2[2] zt3 = secondary_system3[2] - primary_system3[2] zt4 = secondary_system4[2] - primary_system4[2] zt = (zt1+zt2+zt3+zt4)/4 #averaging

so above attempted calculate average translation vector every axis

if translation , rotation, transformation known affine transformation.

it takes form:

secondary_system = * primary_system + b

where a 3x3 matrix (since you're in 3d), , b 3x1 translation.

this can equivalently written

secondary_system_coords2 = a2 * primary_system2,

where

secondary_system vector [secondary_system,1], primary_system vector [primary_system,1], and

a2 4x4 matrix:

[ b ] [ 0,0,0,1 ]

(see wiki page more info).

so basically, want solve equation:

y = a2 x

for a2, y consist of points secondary_system 1 stuck on end, , x points primary_system 1 stuck on end, , a2 4x4 matrix.

now if x square matrix solve like:

a2 = y*x^(-1)

but x 4x1. however, lucky , have 4 sets of x 4 corresponding sets of y, can build x 4x4 so:

x = [ primary_system1 | primary_system2 | primary_system3 | primary_system4 ]

where each of primary_systemi 4x1 column vector. same y.

once have a2, transform point system1 scheme 2 do:

transformed = a2 * point_to_transform

you can set (e.g. in numpy) this:

import numpy np def solve_affine( p1, p2, p3, p4, s1, s2, s3, s4 ): x = np.transpose(np.matrix([p1,p2,p3,p4])) y = np.transpose(np.matrix([s1,s2,s3,s4])) # add together ones on bottom of x , y x = np.vstack((x,[1,1,1,1])) y = np.vstack((y,[1,1,1,1])) # solve a2 a2 = y * x.i # homecoming function takes input x , transforms # don't need homecoming 4th row homecoming lambda x: (a2*np.vstack((np.matrix(x).reshape(3,1),1)))[0:3,:]

then utilize this:

transformfn = solve_affine( primary_system1, primary_system2, primary_system3, primary_system4, secondary_system1, secondary_system2, secondary_system3, secondary_system4 ) # test: transform primary_system1 , should secondary_system1 np.matrix(secondary_system1).t - transformfn( primary_system1 ) # np.linalg.norm of above 0.02555 # transform point (x,y,z). transformed = transformfn((x,y,z))

note: there of course of study numerical error here, , may not best way solve transform (you might able sort of to the lowest degree squares thing).

also, error converting primary_systemx secondary_systemx (for example) of order 10^(-2).

you'll have consider whether acceptable or not (it seem large, might acceptable when compared input points of order 10^6).

python coordinate-transformation