空间直角坐标转换之仿射变换

 

空间直角坐标转换之仿射变换

□文/ 3Echo

一、引言

工作开发中常常会遇到坐标系转换的问题，关于如何实现不同坐标系之间的转换的论述非常之多，基于实际应用项目，大都提出了一种较好的解决方法。两年前，我也从网上下载了一篇文章——《坐标系转换公式》（青岛海洋地质研究所戴勤奋译），文中对各种变换模型都有详细的描述，如莫洛金斯基-巴德卡斯转换模型、赫尔黙特转换模型、布尔莎模型以及多项式转换，算是一篇比较全面介绍坐标系转换方面的文章。

我想大家对常用转换模型的理解方面一般不会有大太困难，如果基于当前流行GIS平台（如超图、ArcGIS、MapInfo）的基础上作二次开发，我想也不会有什么困难，只要找准了它们提供的接口，理顺一下思路，我们也能实现用户提出的需求。但是对于内核算法、参数求解的过程我们却一无所知，很多时候我们自己觉得解决了这个问题，也就不会太去关注底层实现的算法问题了。不过，说实话要去真正弄清楚各个模型之间的关系确实是一件头痛的事情，没有一定的数学功底还真的是不知道它在说些什么。

二、仿射变换

仿射变换是空间直角坐标变换的一种，它是一种二维坐标到二维坐标之间的线性变换，保持二维图形的“平直线”和“平行性”，其可以通过一系列的原子变换的复合来实现，包括平移(Translation)、缩放（Scale）、翻转（Flip）、旋转（Rotation）和剪切(Shear)。

此类变换可以用一个3×3的矩阵来表示，其最后一行为(0, 0, 1)。该变换矩阵将原坐标(x, y)变换为新坐标(x', y')，这里原坐标和新坐标皆视为最末一行为(1)的三维列向量，原列向量左乘变换矩阵得到新的列向量：

[x']   [m00 m01 m02] [x]   [m00*x+m01*y+m02]

[y'] = [m10 m11 m12] [y] = [m10*x+m11*y+m12]

[1 ]   [ 0   0   1 ] [1]   [     1         ]

用代数式表示如下：

x’ = m00*x+m01*y+m02；

             y’ = m10*x+m11*y+m12；

如果将它写成按旋转、缩放、平移三个分量的复合形式，则其代数式如下：

    其示意图如下：

几种典型的仿射变换：

1.public static AffineTransform getTranslateInstance(double tx, double ty)

平移变换，将每一点移动到(x+tx, y+ty)，变换矩阵为：

[   1    0    tx  ]

[   0    1    ty  ]

[   0    0    1   ]

（译注：平移变换是一种“刚体变换”，rigid-body transformation，中学学过的物理，都知道啥叫“刚体”吧，就是不会产生形变的理想物体，平移当然不会改变二维图形的形状。同理，下面的“旋转变换”也是刚体变换，而“缩放”、“错切”都是会改变图形形状的。） 

2.public static AffineTransform getScaleInstance(double sx, double sy)

缩放变换，将每一点的横坐标放大（缩小）至sx倍，纵坐标放大（缩小）至sy倍，变换矩阵为：

[   sx   0    0   ]

[   0    sy   0   ]

[   0    0    1   ]

3.public static AffineTransform getShearInstance(double shx, double shy)

剪切变换，变换矩阵为：

[   1   shx   0   ]

[  shy   1    0   ]

[   0    0    1   ]

相当于一个横向剪切与一个纵向剪切的复合

[   1      0    0   ][   1   shx   0   ]

[  shy    1    0   ][   0     1   0  ]

[   0      0    1   ][   0     0    1  ]

（译注：“剪切变换”又称“错切变换”，指的是类似于四边形不稳定性那种性质，街边小商店那种铁拉门都见过吧？想象一下上面铁条构成的菱形拉动的过程，那就是“错切”的过程。） 

4.public static AffineTransform getRotateInstance(double theta)

旋转变换，目标图形围绕原点顺时针旋转theta弧度，变换矩阵为：

[   cos(theta)    -sin(theta)    0   ]

[   sin(theta)     cos(theta)    0   ]

[       0           0            1   ] 

5.public static AffineTransform getRotateInstance(double theta, double x, double y)

旋转变换，目标图形以(x, y)为轴心顺时针旋转theta弧度，变换矩阵为：

[   cos(theta)    -sin(theta)    x-x*cos+y*sin]

[   sin(theta)     cos(theta)    y-x*sin-y*cos ]

[       0              0          1            ]

相当于两次平移变换与一次原点旋转变换的复合：

[1  0  -x][cos(theta)  -sin(theta)  0][1  0  x]

[0  1  -y][sin(theta)   cos(theta)  0][0  1  y]

[0  0  1 ][   0           0        1 ][0  0  1]

三、仿射变换四参数求解

A、C#自定义函数实现求解：

1、求解旋转参数Rotaion:

 1/**////<summary>
 2
 3        ///获取旋转角度
 4
 5        ///</summary>
 6
 7        ///<param name="fromCoordPoint1">源点1</param>
 8
 9        ///<param name="toCoordPoint1">目标点1</param>
10
11        ///<param name="fromCoordPoint2">源点2</param>
12
13        ///<param name="toCoordPoint2">目标点2</param>
14
15        ///<returns>返回旋转角度</returns>
16
17        private double GetRotation(CoordPoint fromPoint1, CoordPoint toPoint1,CoordPoint fromPoint2,CoordPoint toPoint2)
18
19        {
20
21            double a = (toPoint2.Y - toPoint1.Y) * (fromPoint2.X - fromPoint1.X) - (toPoint2.X - toPoint1.X) * (fromPoint2.Y - fromPoint1.Y);
22
23            double b = (toPoint2.X - toPoint1.X) * (fromPoint2.X - fromPoint1.X) + (toPoint2.Y - toPoint1.Y) * (fromPoint2.Y - fromPoint1.Y);
24
25            
26
27            if (Math.Abs(b) > 0)
28
29                return Math.Tan(a / b);
30
31            else
32
33                return Math.Tan(0);            
34
35        }
36
37

        2、求解缩放比例参数(Scale):

 1 /**////<summary>
 2
 3        ///获取缩放比例因子
 4
 5        ///</summary>
 6
 7        ///<param name="fromCoordPoint1">源点1</param>
 8
 9        ///<param name="toCoordPoint1">目标点1</param>
10
11        ///<param name="fromCoordPoint2">源点2</param>
12
13        ///<param name="toCoordPoint2">目标点2</param>
14
15        ///<param name="rotation">旋转角度</param>
16
17        ///<returns>返回旋转因子</returns>
18
19        private double GetScale(CoordPoint fromPoint1, CoordPoint toPoint1, CoordPoint fromPoint2, CoordPoint toPoint2, double rotation)
20
21        {
22
23            double a = toPoint2.X - toPoint1.X;
24
25            double b = (fromPoint2.X - fromPoint1.X) * Math.Cos(rotation) - (fromPoint2.Y - fromPoint1.Y)*Math.Sin(rotation);
26
27            if (Math.Abs(b) > 0)
28
29                return a / b;
30
31            else
32
33                return 0;
34
35        }
36
37

       3、求解X方向偏移距离参数（XTranslate）:

 1/**////<summary>
 2
 3        ///得到X方向偏移量
 4
 5        ///</summary>
 6
 7        ///<param name="fromCoordPoint1">源点1</param>
 8
 9        ///<param name="toCoordPoint1">目标点1</param>
10
11        ///<param name="rotation">旋转角度</param>
12
13        ///<param name="scale">缩放因子</param>
14
15        ///<returns>返回X方向偏移量</returns>
16
17        private double GetXTranslation(CoordPoint fromPoint1,CoordPoint toPoint1,double rotation,double scale)
18
19        {
20
21            return (toPoint1.X - scale * (fromPoint1.X * Math.Cos(rotation) - fromPoint1.Y * Math.Sin(rotation)));
22
23        }
24
25

        4、求解Y方向偏移距离参数(YTranslate):

 1 /**////<summary>
 2
 3        ///得到Y方向偏移量
 4
 5        ///</summary>
 6
 7        ///<param name="fromCoordPoint1">源点1</param>
 8
 9        ///<param name="toCoordPoint1">目标点1</param>
10
11        ///<param name="rotation">旋转角度</param>
12
13        ///<param name="scale">缩放因子</param>
14
15        ///<returns>返回Y方向偏移量</returns>
16
17        private double GetYTranslation(CoordPoint fromPoint1, CoordPoint toPoint1, double rotation, double scale)
18
19        {
20
21            return (toPoint1.Y - scale * (fromPoint1.X * Math.Sin(rotation) + fromPoint1.Y * Math.Cos(rotation)));
22
23        }

       B、C#+AE求解：

 1/**////<summary>
 2
 3        ///从控制点定义仿射变换程式
 4
 5        ///</summary>
 6
 7        ///<param name="pFromPoints">源控制点</param>
 8
 9        ///<param name="pToPoints">目标控制点</param>
10
11        ///<returns>返回变换定义</returns>
12
13        private ITransformation GetAffineTransformation(IPoint[] pFromPoints, IPoint[] pToPoints)
14
15        {
16
17            //实例化仿射变换对象
18
19            IAffineTransformation2D3GEN tAffineTransformation = new AffineTransformation2DClass();
20
21            //从源控制点定义参数
22
23            tAffineTransformation.DefineFromControlPoints(ref pFromPoints, ref pToPoints);
24
25            //查询引用接口
26
27            ITransformation tTransformation = tAffineTransformation as ITransformation;
28
29            return tTransformation;
30
31        }

四、空间对象转换

求出参数后，再利用公式对相应坐标点进行转换是一件相对简单的事件了。

示例代码：

 1 /**////<summary>
 2
 3        ///转换空间点
 4
 5        ///</summary>
 6
 7        ///<param name="pPoint">点</param>
 8
 9        ///<returns>返回转换后的点</returns>
10
11        private IGeometry TransformPoint(IPoint pPoint)
12
13        {
14
15            //********************************************
16
17            //说明：采用相似变换模型（四参数变换模型）
18
19            // X= ax - by + c
20
21            // Y= bx + ay + d
22
23            //*********************************************
24
25            double A = this.m_Scale * Math.Cos(this.m_RotationAngle);
26
27            double B = this.m_Scale * Math.Sin(this.m_RotationAngle);
28
29            IPoint tPoint = new PointClass();
30
31            tPoint.X = A * pPoint.X - B * pPoint.Y + this.m_DX;
32
33            tPoint.Y = B * pPoint.X + A * pPoint.Y + this.m_DY;
34
35            return tPoint;
36
37        }

       五 总结：

本文主要介绍了如何利用仿射变换方程来进行空间直角坐标转换，对仿射变换的几种典型情况作了详细的讲解，对于具体如何求解参数给出了关键的实现代码，对于空间对象的变换给出了参考示例。如果是ArcGIS用户，完全可以利用它自身提供的接口进行空间转换。

写这篇文章的时候，说实话，对于坐标变换的各个模型我也不是完全的理解，心中存在着许多问题，比如说如何利用最小二乘法公式来求解参数就一直没有弄清楚，还希望各位朋友能够多多指点，不胜感激！

六 备注

希望基于AE开发的朋友注意一下，9.2版本中提供的关于仿射变换模型，其代数形式有误：

其提供的错误代数形式：

X=ax+by+c

Y=-bx+ay+f

正确形式应该如下:

X=ax-by+c

Y=bx+ay+d

参考资料：

1、《坐标系转换公式》（青岛海洋地质研究所戴勤奋译）

2、Java文档帮助之AffineTransform

3、ESRI开发文档