计算仿射变换六参数(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[0])  # 维度
    if len(q) < dim:
        print "点数小于维度."
        return False

    # 新建一个空的 维度 x (维度+1) 矩阵 并填满
    c = [[0.0 for a in range(dim)] 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[i]) + [1]
                c[k][j] += qt[k] * p[i][j]

    # 新建一个空的 (维度+1) x (维度+1) 矩阵 并填满
    Q = [[0.0 for a in range(dim)] + [0] for i in range(dim+1)]
    for qi in q:
        qt = list(qi) + [1]
        for i in range(dim+1):
            for j in range(dim+1):
                Q[i][j] += qt[i] * qt[j]

    # 判断原始点和目标点是否共线，共线则无解. 耗时计算，如果追求效率可以不用。
    # 其实就是解n个三元一次方程组
    def gauss_jordan(m, eps=1.0/(10**10)):
        (h, w) = (len(m), len(m[0]))
        for y in range(0, h):
            maxrow = y
            for y2 in range(y+1, h):    
                if abs(m[y2][y]) > abs(m[maxrow][y]):
                    maxrow = y2
            (m[y], m[maxrow]) = (m[maxrow], m[y])
            if abs(m[y][y]) <= eps:     
                return False
            for y2 in range(y+1, h):    
                c = m[y2][y] / m[y][y]
                for x in range(y, w):
                    m[y2][x] -= m[y][x] * c
        for y in range(h-1, 0-1, -1):  
            c = m[y][y]
            for y2 in range(0, y):
                for x in range(w-1, y-1, -1):
                    m[y2][x] -= m[y][x] * m[y2][y] / c
            m[y][y] /= c
            for x in range(h, w):       
                m[y][x] /= c
        return True

    
    M = [Q[i] + c[i] 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[i][j+dim+1])
                str += "%f" % M[dim][j+dim+1]
                res += str + "\n"
            return res

        def transform(self, pt):
            res = [0.0 for a in range(dim)]
            for j in range(dim):
                for i in range(dim):
                    res[j] += pt[i] * M[i][j+dim+1]
                res[j] += M[dim][j+dim+1]
            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[i]
            tp = to_pt[i]
            t = trn.Transform(fp)
            print ("%s => %s ~= %s" % (fp, tuple(t), tp))
            err += ((tp[0] - t[0])**2 + (tp[1] - t[1])**2)**0.5

        print "拟合误差 = %f" % err

if __name__ == "__main__":
    test()