计算仿射变换六参数(Python)
做保密处理或者坐标系转换的时候经常会用到关于空间配准的问题，那么如何根据已知点对求得坐标转换的参数是一个值得研究的问题。这里用到的编程技巧不多，关键是要用到线性代数和数值分析的知识。纵观当前地图坐标保密处理或者坐标系转换的实例，其无外乎采用旋转、平移、拉伸等方式，于是数值的计算无外乎于解n个n元一次方程组，最后通过误差分析进行拟合。下面就是一个形如 x' = Ax + by +C; y' = Dx + Ey +F 的六参数仿射变换的数值解法，采用Python2语言。
__author__ = 'Lee'
# -*- coding: utf-8 -*-
def affine_fit(from_pts, to_pts):
q = from_pts
p = to_pts
if len(q) != len(p) or len(q) < 1:
print "原始点和目标点的个数必须相同."
return False
dim = len(q)# 维度
if len(q) < dim:
print "点数小于维度."
return False
# 新建一个空的 维度 x (维度+1) 矩阵 并填满
c = [ for i in range(dim+1)]
for j in range(dim):
for k in range(dim+1):
for i in range(len(q)):
qt = list(q) +
c += qt * p
# 新建一个空的 (维度+1) x (维度+1) 矩阵 并填满
Q = [ + for i in range(dim+1)]
for qi in q:
qt = list(qi) +
for i in range(dim+1):
for j in range(dim+1):
Q += qt * qt
# 判断原始点和目标点是否共线，共线则无解. 耗时计算，如果追求效率可以不用。
# 其实就是解n个三元一次方程组
def gauss_jordan(m, eps=1.0/(10**10)):
(h, w) = (len(m), len(m))
for y in range(0, h):
maxrow = y
for y2 in range(y+1, h):
if abs(m) > abs(m):
maxrow = y2
(m, m) = (m, m)
if abs(m) <= eps:
return False
for y2 in range(y+1, h):
c = m / m
for x in range(y, w):
m -= m * c
for y in range(h-1, 0-1, -1):
c = m
for y2 in range(0, y):
for x in range(w-1, y-1, -1):
m -= m * m / c
m /= c
for x in range(h, w):
m /= c
return True
M = + c for i in range(dim+1)]
if not gauss_jordan(M):
print "错误，原始点和目标点也许是共线的."
return False
class transformation:
"""对象化仿射变换."""
def To_Str(self):
res = ""
for j in range(dim):
str = "x%d' = " % j
for i in range(dim):
str +="x%d * %f + " % (i, M)
str += "%f" % M
res += str + "\n"
return res
def transform(self, pt):
res =
for j in range(dim):
for i in range(dim):
res += pt * M
res += M
return res
return transformation()
def test():
from_pt = ((38671803.6437, 2578831.9242), (38407102.8445, 2504239.2774), (38122268.3963, 2358570.38514),
(38126455.4595, 2346827.2602), (38177232.2601, 2398763.77833), (38423567.3485, 2571733.9203),
(38636876.4495, 2543442.3694), (38754169.8762, 2662401.86536), (38410773.8815, 2558886.6518),
(38668962.0430, 2578747.6349))# 输入点坐标对
to_pt = ((38671804.6165, 2578831.1944), (38407104.0875, 2504239.1898), (38122269.2925, 2358571.57626),
(38126456.5675, 2346826.27022), (38177232.3973, 2398762.11714), (38423565.7744, 2571735.2278),
(38636873.6217, 2543440.7216), (38754168.8662, 2662401.86101), (38410774.5621, 2558886.0921),
(38668962.5493, 2578746.94)) # 输出点坐标对
trn = affine_fit(from_pt, to_pt)
if trn:
print "转换公式:"
print trn.To_Str()
err = 0.0
for i in range(len(from_pt)):
fp = from_pt
tp = to_pt
t = trn.Transform(fp)
print ("%s => %s ~= %s" % (fp, tuple(t), tp))
err += ((tp - t)**2 + (tp - t)**2)**0.5
print "拟合误差 = %f" % err
if __name__ == "__main__":
test()
页:
[1]