图像坐标空间变换：透视变换（Perspective Transformation），或称为单应性（Homography）变换

真实的世界是三维的，而图像是二维的（至少目前是），如果要以二维图像描述三维世界，并且看起来足够真实，那么三维世界向二维图像转化的过程就需要满足一定的几何投影关系，即透视关系，用非常简单的话来讲就是近处的物体在图像中显得比较大，远处的物体在图像中显得比较小。比如下图中的铁轨，越往远处看，图像中铁轨的间距越小，但实际上铁轨的间距是不变的。 透视变换就是对图像中的物体进行空间坐标变换时，使变换结果满足一定的透视关系。透视变换包含以下三个过程：

说到这里可能仍然非常抽象，继续向下看应该能使各位对透视变换的概念逐渐变得清晰一些。

在继续向下了解透视变换的公式推导之前，建议先详细了解一下仿射变换和图像插值的相关知识，这部分知识在其他地方已经写过了，所以这里不会再写，不过请务必详细了解，否则再往下理解起来可能会有点费劲。 图像坐标空间变换：仿射变换（Affine Transformation） 图像插值：最邻近（nearest）与双线性（bilinear）

其中特别注意以下内容：

下面我们先讲投影的概念，然后再讲二维坐标向齐次坐标的变换（请注意，这个讲述顺序与计算顺序相反）

在仿射变换（Affine Transformation）中曾经简单讲过， ( x , y , 1 ) (x,y,1) (x,y,1)就是 ( x , y ) (x,y) (x,y)对应的齐次坐标。现在我们要扩展一下齐次坐标的概念。当我们使用相机进行拍摄时，可以认为三维空间中的物体都投影在了一个二维图像上面，如下图所示，红色的平面就可以认为是投影平面，我们将其定义为 z = 1 z=1 z=1的平面，那么 ( x , y , 1 ) (x,y,1) (x,y,1)就是二维图像上的点 ( x , y ) (x,y) (x,y)在三维空间中的坐标。

接下来，对 ( x , y , 1 ) (x,y,1) (x,y,1)乘以一个系数，即 α ( x , y , 1 ) \alpha(x,y,1) α(x,y,1)，改变 α \alpha α就相当于在红色虚线上移动 ( x , y , 1 ) (x,y,1) (x,y,1)点，这一点可以通过相似三角形的原理进行证明。红色虚线通过连接观察者（或相机）与点 ( x , y , 1 ) (x,y,1) (x,y,1)得到。

我们将 α ( x , y , 1 ) \alpha(x,y,1) α(x,y,1)记为 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)，即 ( X = α x , Y = α y , Z = α ) (X=\alpha x, Y=\alpha y, Z=\alpha) (X=αx,Y=αy,Z=α)。投影就是将 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)投影至 ( x , y , 1 ) (x,y,1) (x,y,1)的过程。若已知 ( X , Y , Z ) (X,Y,Z) (X,Y,Z)，可以通过投影计算出 ( x , y ) (x,y) (x,y)，其计算过程比较简单，根据 α ( x , y , 1 ) = ( X , Y , Z ) \alpha(x,y,1)=(X,Y,Z) α(x,y,1)=(X,Y,Z)的定义， ( X , Y ) (X,Y) (X,Y)同除以 Z Z Z即可： { x = X / Z y = Y / Z \begin{cases} x = X/Z \\[2ex] y = Y/Z \end{cases} ⎩⎨⎧​x=X/Zy=Y/Z​ 上述定义和计算过程表明，所有由 α ( x , y , 1 ) \alpha(x,y,1) α(x,y,1)定义的点有着一个共同的投影坐标，即 ( x , y , 1 ) (x,y,1) (x,y,1)。

此处我们直接推导后向映射过程，后向映射是编程时真正使用的计算方法。 记原始图像中的坐标为 ( u , v ) (u,v) (u,v)，对应的齐次坐标为 ( U , V , W ) (U,V,W) (U,V,W)，变换后图像的坐标为 ( x , y ) (x,y) (x,y)。 那么有： [ U V W ] = [ x y 1 ] [ t 11 t 12 t 13 t 21 t 22 t 23 t 31 t 32 t 33 ] \begin{bmatrix} U & V & W \end{bmatrix} = \begin{bmatrix} x & y & 1 \end{bmatrix} \begin{bmatrix} t_{11} & t_{12} & t_{13} \\ t_{21} & t_{22} & t_{23} \\ t_{31} & t_{32} & t_{33} \end{bmatrix} [U​V​W​]=[x​y​1​]⎣⎡​t11​t21​t31​​t12​t22​t32​​t13​t23​t33​​⎦⎤​ 请注意，这个式子不好求解，因为如果将此式化为方程组形式，你会发现没有常系数项。接下来，我们要运用上面投影部分讲过的概念变换此式，使其容易求解。

根据投影的概念，当我们求出 ( U , V , W ) (U,V,W) (U,V,W)后，我们接着通过 ( u = U / W , v = V / W ) (u = U/W,v=V/W) (u=U/W,v=V/W)来求解 ( u , v ) (u,v) (u,v)，请特别注意，在这个过程中，我们可以对 ( U , V , W ) (U,V,W) (U,V,W)乘以任意一个常系数而不会影响计算结果，也就是说，我们可以对上面的矩阵方程两边同乘以一个系数。

那么我们两边同除以 t 33 t_{33} t33​，得到的结果除了将 t 33 t_{33} t33​改写为1以外，其他符号保持不变，所以可得： [ U V W ] = [ x y 1 ] [ t 11 t 12 t 13 t 21 t 22 t 23 t 31 t 32 1 ] \begin{bmatrix} U & V & W \end{bmatrix} = \begin{bmatrix} x & y & 1 \end{bmatrix} \begin{bmatrix} t_{11} & t_{12} & t_{13} \\ t_{21} & t_{22} & t_{23} \\ t_{31} & t_{32} & 1 \end{bmatrix} [U​V​W​]=[x​y​1​]⎣⎡​t11​t21​t31​​t12​t22​t32​​t13​t23​1​⎦⎤​ 这个式子就容易求解了。 现在可以去跟仿射变换（Affine Transformation）对比一下，然后会发现，这里的变换矩阵多了两个参数： t 13 t_{13} t13​和 t 23 t_{23} t23​。

将上述矩阵方程式展开得： { U = t 11 ∗ x + t 21 ∗ y + t 31 V = t 12 ∗ x + t 22 ∗ y + t 32 W = t 13 ∗ x + t 23 ∗ y + 1 \begin{cases} U = t_{11}*x + t_{21}*y + t_{31} \\ V = t_{12}*x + t_{22}*y + t_{32} \\ W = t_{13}*x + t_{23}*y + 1 \end{cases} ⎩⎪⎨⎪⎧​U=t11​∗x+t21​∗y+t31​V=t12​∗x+t22​∗y+t32​W=t13​∗x+t23​∗y+1​ 由投影过程可得： { u = U W = t 11 ∗ x + t 21 ∗ y + t 31 t 13 ∗ x + t 23 ∗ y + 1 v = V W = t 12 ∗ x + t 22 ∗ y + t 32 t 13 ∗ x + t 23 ∗ y + 1 \begin{cases} u = \frac {U}{W} = \frac {t_{11}*x + t_{21}*y + t_{31}} {t_{13}*x + t_{23}*y + 1} \\[2ex] v = \frac {V}{W} = \frac {t_{12}*x + t_{22}*y + t_{32}} {t_{13}*x + t_{23}*y + 1} \end{cases} ⎩⎨⎧​u=WU​=t13​∗x+t23​∗y+1t11​∗x+t21​∗y+t31​​v=WV​=t13​∗x+t23​∗y+1t12​∗x+t22​∗y+t32​​​ 进一步可得： { u ( t 13 ∗ x + t 23 ∗ y + 1 ) = t 11 ∗ x + t 21 ∗ y + t 31 v ( t 13 ∗ x + t 23 ∗ y + 1 ) = t 12 ∗ x + t 22 ∗ y + t 32 \begin{cases} u(t_{13}*x + t_{23}*y + 1) = t_{11}*x + t_{21}*y + t_{31} \\[2ex] v(t_{13}*x + t_{23}*y + 1) = t_{12}*x + t_{22}*y + t_{32} \end{cases} ⎩⎨⎧​u(t13​∗x+t23​∗y+1)=t11​∗x+t21​∗y+t31​v(t13​∗x+t23​∗y+1)=t12​∗x+t22​∗y+t32​​ 整理得： { x ∗ t 11 + y ∗ t 21 + 1 ∗ t 31 + 0 ∗ t 12 + 0 ∗ t 22 + 0 ∗ t 32 − u x ∗ t 13 − u y ∗ t 23 = u 0 ∗ t 11 + 0 ∗ t 21 + 0 ∗ t 31 + x ∗ t 12 + y ∗ t 22 + 1 ∗ t 32 − v x ∗ t 13 − v y ∗ t 23 = v \begin{cases} x*t_{11} + y*t_{21} + 1*t_{31} + 0*t_{12} + 0*t_{22} + 0*t_{32} - ux*t_{13} - uy*t_{23} = u \\[2ex] 0*t_{11} + 0*t_{21} + 0*t_{31} + x*t_{12} + y*t_{22} + 1*t_{32} - vx*t_{13} - vy*t_{23} = v \end{cases} ⎩⎨⎧​x∗t11​+y∗t21​+1∗t31​+0∗t12​+0∗t22​+0∗t32​−ux∗t13​−uy∗t23​=u0∗t11​+0∗t21​+0∗t31​+x∗t12​+y∗t22​+1∗t32​−vx∗t13​−vy∗t23​=v​ 假设我们现在有n组对应的点： ( u 1 , v 1 ) (u_1,v_1) (u1​,v1​)对应 ( x 1 , y 1 ) (x_1,y_1) (x1​,y1​)， ( u 2 , v 2 ) (u_2,v_2) (u2​,v2​)对应 ( x 2 , y 2 ) (x_2,y_2) (x2​,y2​), …, ( u n , v n ) (u_n,v_n) (un​,vn​)对应 ( x n , y n ) (x_n,y_n) (xn​,yn​)，那么可以得到矩阵形式的线性方程组： [ x 1 y 1 1 0 0 0 − u 1 x 1 − u 1 y 1 0 0 0 x 1 y 1 1 − v 1 x 1 − v 1 y 1 x 2 y 2 1 0 0 0 − u 2 x 2 − u 2 y 2 0 0 0 x 2 y 2 1 − v 2 x 2 − v 2 y 2 . . . x n y n 1 0 0 0 − u n x n − u n y n 0 0 0 x n y n 1 − v n x n − v n y n ] [ t 11 t 21 t 31 t 12 t 22 t 32 t 13 t 23 ] = [ u 1 v 1 u 2 v 2 . . . u n v n ] \begin{bmatrix} x_1 & y_1 & 1 & 0 & 0 & 0 & -u_1x_1 & -u_1y_1 \\ 0 & 0 & 0 & x_1 & y_1 & 1 & -v_1x_1 & -v_1y_1 \\ x_2 & y_2 & 1 & 0 & 0 & 0 & -u_2x_2 & -u_2y_2 \\ 0 & 0 & 0 & x_2 & y_2 & 1 & -v_2x_2 & -v_2y_2 \\ ... \\ x_n & y_n & 1 & 0 & 0 & 0 & -u_nx_n & -u_ny_n \\ 0 & 0 & 0 & x_n & y_n & 1 & -v_nx_n & -v_ny_n \\ \end{bmatrix} \begin{bmatrix} t_{11}\\ t_{21} \\ t_{31} \\ t_{12} \\ t_{22} \\ t_{32} \\ t_{13} \\ t_{23} \end{bmatrix} = \begin{bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ ... \\ u_n \\ v_n \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​x1​0x2​0...xn​0​y1​0y2​0yn​0​101010​0x1​0x2​0xn​​0y1​0y2​0yn​​010101​−u1​x1​−v1​x1​−u2​x2​−v2​x2​−un​xn​−vn​xn​​−u1​y1​−v1​y1​−u2​y2​−v2​y2​−un​yn​−vn​yn​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​t11​t21​t31​t12​t22​t32​t13​t23​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​u1​v1​u2​v2​...un​vn​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​ 记： A = [ x 1 y 1 1 0 0 0 − u 1 x 1 − u 1 y 1 0 0 0 x 1 y 1 1 − v 1 x 1 − v 1 y 1 x 2 y 2 1 0 0 0 − u 2 x 2 − u 2 y 2 0 0 0 x 2 y 2 1 − v 2 x 2 − v 2 y 2 . . . x n y n 1 0 0 0 − u n x n − u n y n 0 0 0 x n y n 1 − v n x n − v n y n ] A = \begin{bmatrix} x_1 & y_1 & 1 & 0 & 0 & 0 & -u_1x_1 & -u_1y_1 \\ 0 & 0 & 0 & x_1 & y_1 & 1 & -v_1x_1 & -v_1y_1 \\ x_2 & y_2 & 1 & 0 & 0 & 0 & -u_2x_2 & -u_2y_2 \\ 0 & 0 & 0 & x_2 & y_2 & 1 & -v_2x_2 & -v_2y_2 \\ ... \\ x_n & y_n & 1 & 0 & 0 & 0 & -u_nx_n & -u_ny_n \\ 0 & 0 & 0 & x_n & y_n & 1 & -v_nx_n & -v_ny_n \\ \end{bmatrix} A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​x1​0x2​0...xn​0​y1​0y2​0yn​0​101010​0x1​0x2​0xn​​0y1​0y2​0yn​​010101​−u1​x1​−v1​x1​−u2​x2​−v2​x2​−un​xn​−vn​xn​​−u1​y1​−v1​y1​−u2​y2​−v2​y2​−un​yn​−vn​yn​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​

T = [ t 11 t 21 t 31 t 12 t 22 t 32 t 13 t 23 ] T = \begin{bmatrix} t_{11}\\ t_{21} \\ t_{31} \\ t_{12} \\ t_{22} \\ t_{32} \\ t_{13} \\ t_{23} \end{bmatrix} T=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡​t11​t21​t31​t12​t22​t32​t13​t23​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​

B = [ u 1 v 1 u 2 v 2 . . . u n v n ] B = \begin{bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ ... \\ u_n \\ v_n \end{bmatrix} B=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​u1​v1​u2​v2​...un​vn​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​ 那么可得最小二乘解为： T = ( A T A ) − 1 A T B T = (A^TA)^{-1}A^TB T=(ATA)−1ATB 最小二乘解的推导过程是： A T = B A T A T = A T B ( A T A ) − 1 A T A T = ( A T A ) − 1 A T B T = ( A T A ) − 1 A T B AT = B \\[2ex] A^TAT = A^TB \\[2ex] (A^TA)^{-1}A^TAT = (A^TA)^{-1}A^TB \\[2ex] T = (A^TA)^{-1}A^TB AT=BATAT=ATB(ATA)−1ATAT=(ATA)−1ATBT=(ATA)−1ATB 注意：由于上述方线性程组中有8个未知数，所以最少需要4组对应的点才能得到方程组的解，如果刚好有4组点，那么得到的解就是精确解，如果多于4组点，那么得到的解就称为最小二乘解。

有的时候我们拍摄角度不太精准的情况下，A4纸在拍摄的图像中不是呈现出标准的长方形，而是不规则的四边形，这时我们可以通过透视变换进行视角校正。

比如下图中有一张A4纸，但是很明显它不是长方形，现在我们就用透视变换来校正它，步骤如下：

下图在程序中的名字是side.jpg