最近换了新的工作,可恶的是这个地方不能上网,因此上论坛的时间大量减少。
工作中遇到了将图形进行旋转,查了几本书,发现了仿射变换,才后悔大学时没有好好学习数学和图形学。
将找到的一篇文章转帖在这儿。
空间直角坐标转换之仿射变换
作者:3echo 文章来源:http://www.cnblogs.com/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;
如果将它写成按旋转、缩放、平移三个分量的复合形式,则其代数式如下:
X= (x-x0)cosθ+(y-y0)sinθ+x0
Y=-(x-x0)sinθ+(y-y0)cosθ+y0
其示意图如下:
几种典型的仿射变换:
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
- ///目标点1
- ///源点2
- ///目标点2
- ///返回旋转角度
- private double GetRotation(CoordPoint fromPoint1, CoordPoint toPoint1,CoordPoint fromPoint2,CoordPoint toPoint2)
- {
- double a = (toPoint2.Y - toPoint1.Y) * (fromPoint2.X - fromPoint1.X) - (toPoint2.X - toPoint1.X) * (fromPoint2.Y - fromPoint1.Y);
- double b = (toPoint2.X - toPoint1.X) * (fromPoint2.X - fromPoint1.X) + (toPoint2.Y - toPoint1.Y) * (fromPoint2.Y - fromPoint1.Y);
- if (Math.Abs(b) > 0)
- return Math.Tan(a / b);
- else
- return Math.Tan(0);
- }
2、求解缩放比例参数(Scale):
- /**////
- ///获取缩放比例因子
- ///
- ///源点1
- ///目标点1
- ///源点2
- ///目标点2
- ///旋转角度
- ///返回旋转因子
- private double GetScale(CoordPoint fromPoint1, CoordPoint toPoint1, CoordPoint fromPoint2, CoordPoint toPoint2, double rotation)
- {
- double a = toPoint2.X - toPoint1.X;
- double b = (fromPoint2.X - fromPoint1.X) * Math.Cos(rotation) - (fromPoint2.Y - fromPoint1.Y)*Math.Sin(rotation);
- if (Math.Abs(b) > 0)
- return a / b;
- else
- return 0;
- }
3、求解X方向偏移距离参数(XTranslate):
- /**////
- ///得到X方向偏移量
- ///
- ///源点1
- ///目标点1
- ///旋转角度
- ///缩放因子
- ///返回X方向偏移量
- private double GetXTranslation(CoordPoint fromPoint1,CoordPoint toPoint1,double rotation,double scale)
- {
- return (toPoint1.X - scale * (fromPoint1.X * Math.Cos(rotation) - fromPoint1.Y * Math.Sin(rotation)));
- }
4、求解Y方向偏移距离参数(YTranslate):
- /**////
- ///得到Y方向偏移量
- ///
- ///源点1
- ///目标点1
- ///旋转角度
- ///缩放因子
- ///返回Y方向偏移量
- private double GetYTranslation(CoordPoint fromPoint1, CoordPoint toPoint1, double rotation, double scale)
- {
- return (toPoint1.Y - scale * (fromPoint1.X * Math.Sin(rotation) + fromPoint1.Y * Math.Cos(rotation)));
- }
B、C#+AE求解:
- /**////
- ///从控制点定义仿射变换程式
- ///
- ///源控制点
- ///目标控制点
- ///返回变换定义
- private ITransformation GetAffineTransformation(IPoint[] pFromPoints, IPoint[] pToPoints)
- {
- //实例化仿射变换对象
- IAffineTransformation2D3GEN tAffineTransformation = new AffineTransformation2DClass();
- //从源控制点定义参数
- tAffineTransformation.DefineFromControlPoints(ref pFromPoints, ref pToPoints);
- //查询引用接口
- ITransformation tTransformation = tAffineTransformation as ITransformation;
- return tTransformation;
- }
四、空间对象转换
求出参数后,再利用公式对相应坐标点进行转换是一件相对简单的事件了。
示例代码:
- /**////
- ///从控制点定义仿射变换程式
- ///
- ///源控制点
- ///目标控制点
- ///返回变换定义
- private ITransformation GetAffineTransformation(IPoint[] pFromPoints, IPoint[] pToPoints)
- {
- //实例化仿射变换对象
- IAffineTransformation2D3GEN tAffineTransformation = new AffineTransformation2DClass();
- //从源控制点定义参数
- tAffineTransformation.DefineFromControlPoints(ref pFromPoints, ref pToPoints);
- //查询引用接口
- ITransformation tTransformation = tAffineTransformation as ITransformation;
- return tTransformation;
- }
五 总结:
本文主要介绍了如何利用仿射变换方程来进行空间直角坐标转换,对仿射变换的几种典型情况作了详细的讲解,对于具体如何求解参数给出了关键的实现代码,对于空间对象的变换给出了参考示例。如果是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开发文档
因为现在工作的地方保密条例比较严格,因此不将自己的代码贴出来了。
关于备注,作者所讲的代数式子
X=ax-by+c
Y=bx+ay+d
是正确的。
但是在实际编程应用中,因为C,C++等语言所提供的sin函数得到的结果是有符号的,即sin(负数)得到的结果是个负数,因此在编程中的式子应该是
X=ax+by+c
Y=-bx+ay+f