空間直角座標轉換之仿射變換)

一、引言
    工作開發中常常會遇到座標系轉換的問題，關於如何?不同座標系之間的轉換的論述非常之多，基於實際應用項目，大都提出了一種較好的解決方案。兩年前，我也從網上下載了一篇文章——《座標系轉換公式》（青島海洋地質研究所戴勤奮譯），文中對各種變換模型都有詳細的描述，如莫洛金斯基-巴德卡斯轉換模型、赫爾黙特轉換模型、布爾莎模型以及多項式轉換，算是一篇比較全面介紹座標系轉換方面的文章。

    我想大家對常用轉換模型的理解方面一般不會有大太困難，如果基於當前流行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
///目標點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 IGeometry TransformPoint(IPoint pPoint)
{
    //********************************************
    //說明：採用相似變換模型（四參數變換模型）
    // X= ax - by + c
    // Y= bx + ay + d
    //*********************************************

    double A = this.m_Scale * Math.Cos(this.m_RotationAngle);
    double B = this.m_Scale * Math.Sin(this.m_RotationAngle);
   
    IPoint tPoint = new PointClass();
   
    tPoint.X = A * pPoint.X - B * pPoint.Y + this.m_DX;
    tPoint.Y = B * pPoint.X + A * pPoint.Y + this.m_DY;
   
    return tPoint;
}五 總結：
    本文主要介紹了如何利用仿射變換方程來進行空間直角座標轉換，對仿射變換的幾種典型情況作了詳細的講解，對於具體如何求解參數給出了關鍵的實現代碼，對於空間對象的變換給出了參考樣本。如果是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開發文檔

 

本文來自CSDN部落格，轉載請標明出處：http://blog.csdn.net/sjzwl/archive/2008/12/31/3668930.aspx