如何使用Python执行坐标仿射变换？

Question

我希望对这个数据集进行转换。
在一个坐标系[primary_system]中，有四个已知点的坐标为x、y、z，在另一个坐标系[secondary_system]中，有下一个四个已知点的坐标为x、y和h。 这些点是对应的；例如，primary_system1点和secondary_system1点实际上是同一个点，但我们在两个不同的坐标系中都有它的坐标。 因此，我这里有四组调整点，并希望根据调整将另一个点的坐标从primary system转换到secondary system。

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 this point
x = 3532412.323 
y = 1175511.432
z = 5161677.111<br>

目前，我正在尝试使用每组四个点来对x、y和z轴进行平均翻译，例如：

#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

...以及y轴和z轴的类似操作

#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

所以，我尝试计算每个轴的平均翻译向量

- daikini

---

你的问题非常模糊不清！那些数字是什么？ - Rik Poggi

你尝试过什么？你有什么代码？ - Marcin

问题不清楚。你到底想做什么？请编辑你的问题以澄清。 - Joel Cornett

这取决于这两个坐标系之间的关系 - 它是旋转和平移吗？还是更复杂的东西，像 GPS -> UTM？ - mathematical.coffee

我已经点赞了这个问题，以抵消负面评价，因为对我来说很清楚 :) - vsnyc

Answer 1

如果只是翻译和旋转，那么这就是一种称为仿射变换的变换。

它基本上采用以下形式：

secondary_system = A * primary_system + b

其中A是一个3x3的矩阵（因为您在三维空间中），b是一个3x1的平移向量。

这也可以等价地写成：

secondary_system_coords2 = A2 * primary_system2,

在哪里

• secondary_system_coords2 is the vector [secondary_system,1],
• primary_system2 is the vector [primary_system,1], and

• A2 is the 4x4 matrix:

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

（更多信息请参见维基页面）。

所以基本上，您想解决方程：

y = A2 x

对于 A2，其中 y 由来自 secondary_system 的点组成，末尾附加了1，而 x 为来自 primary_system 的点组成，末尾附加了1，并且 A2 是一个4x4矩阵。

现在如果 x 是一个方阵，我们可以这样解决：

A2 = y*x^(-1)

但是x是4x1的。不过，你很幸运地有4组x和4组对应的y，因此你可以构建一个4x4的x，如下所示：

x = [ primary_system1 | primary_system2 | primary_system3 | primary_system4 ]

其中每个primary_systemi都是一个4x1列向量。 y同理。

一旦您拥有A2，要将点从系统1转换为系统2，只需执行以下操作：

transformed = A2 * point_to_transform

你可以像这样设置它（例如在numpy中）：

import numpy as 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 ones on the bottom of x and y
    x = np.vstack((x,[1,1,1,1]))
    y = np.vstack((y,[1,1,1,1]))
    # solve for A2
    A2 = y * x.I
    # return function that takes input x and transforms it
    # don't need to return the 4th row as it is 
    return lambda x: (A2*np.vstack((np.matrix(x).reshape(3,1),1)))[0:3,:]

然后像这样使用它：

transformFn = solve_affine( primary_system1, primary_system2, 
                            primary_system3, primary_system4,
                            secondary_system1, secondary_system2,
                            secondary_system3, secondary_system4 )

# test: transform primary_system1 and we should get secondary_system1
np.matrix(secondary_system1).T - transformFn( primary_system1 )
# np.linalg.norm of above is 0.02555

# transform another point (x,y,z).
transformed = transformFn((x,y,z))

注意: 当然，这里存在数值误差，这可能不是解决变换的最佳方法(您可以尝试使用某种最小二乘法)。

另外，将primary_systemx转换为secondary_systemx的误差(对于此示例)约为10^(-2)。

您需要考虑这是否可接受(看起来很大，但与所有输入点都是10^6级别相比可能是可以接受的)。

- mathematical.coffee

---

非常感谢您的完整回答！ - daikini

这可能是我见过的对这个问题最完整的答案之一。干得好。 - b10hazard

Answer 2

你要查找的映射似乎是仿射变换。四个不在同一平面上的三维点是恢复仿射变换所需的确切点数。后者粗略地说，是矩阵乘法和向量加法。

secondary_system = A * primary_system + t

问题现在被简化为找到适当的矩阵A和向量t。我认为这段代码可能会对你有帮助（抱歉代码风格不好——我是数学家，不是程序员）。

import numpy as np
# input data
ins = np.array([[3531820.440, 1174966.736, 5162268.086],
                [3531746.800, 1175275.159, 5162241.325],
                [3532510.182, 1174373.785, 5161954.920],
                [3532495.968, 1175507.195, 5161685.049]]) # <- primary system
out = np.array([[6089665.610, 3591595.470, 148.810],
                [6089633.900, 3591912.090, 143.120],
                [6089088.170, 3590826.470, 166.350],
                [6088672.490, 3591914.630, 147.440]]) # <- secondary system
p = np.array([3532412.323, 1175511.432, 5161677.111]) #<- transform this point
# finding transformation
l = len(ins)
entry = lambda r,d: np.linalg.det(np.delete(np.vstack([r, ins.T, np.ones(l)]), d, axis=0))
M = np.array([[(-1)**i * entry(R, i) for R in out.T] for i in range(l+1)])
A, t = np.hsplit(M[1:].T/(-M[0])[:,None], [l-1])
t = np.transpose(t)[0]
# output transformation
print("Affine transformation matrix:\n", A)
print("Affine transformation translation vector:\n", t)
# unittests
print("TESTING:")
for p, P in zip(np.array(ins), np.array(out)):
  image_p = np.dot(A, p) + t
  result = "[OK]" if np.allclose(image_p, P) else "[ERROR]"
  print(p, " mapped to: ", image_p, " ; expected: ", P, result)
# calculate points
print("CALCULATION:")
P = np.dot(A, p) + t
print(p, " mapped to: ", P)

这段代码演示了如何将仿射变换恢复为矩阵+向量，并测试初始点是否被映射到应该到达的位置。您可以使用 Google colab 测试此代码，因此无需安装任何软件。

关于此代码背后的理论：它基于 "入门指南：仿射映射" 中提出的方程，矩阵恢复在“规范符号恢复”部分中描述，并且在“我们需要多少个点？”一节中讨论了确定精确仿射变换所需的点数。同一作者发表了 "仿射映射工作手册"，其中包含许多此类实用示例。

- guest