我们这里所说的透视，简单地说就是将三维物理空间中的点经过与光学中心的连线投影到某个面上，这个面可以是平面，也可以是曲面，如图 1 是一个简单的示例。平面投影透视的应用最为广泛，我们绝大多数的相机镜头都是基于这种方式设计的，不过这不代表镜头就一定是平面的。其特点是空间中的直线在平面上的投影依然是直线，这样在某种意义上我们可以认为是没有畸变的，当然近大远小作为一种物理性质始终不可避免。这种方式实际上就类似于小孔成像，为了方便我们可以把该投影平面当作是最终的成像平面。但这种方式的缺点也很明显，由于成像点到成像平面中心的距离与入射角是正切的关系，其有效视角通常比较小。
图1 透视示例为了扩大镜头的视角，我们可以采用球面投影透视。而球面是一个曲面，我们很难设计一个球面的感光器件，所以，我们还需要一个从球面到平面的映射过程，这通常是基于折射的原理的来实现的。如果折射角小于入射角，那么我们就可以把更大的视角压缩到一张平面图像上了，就是鱼眼镜头的基本设计原理。如图 2 所示，空间中的一点经过球面投影透视投影到了 P P P 点，然后经过折射投影到了 P 1 P_1 P1 点，但通常成像平面转换到最终的图像需要旋转180°，所以其在鱼眼图像上的坐标应该为 P 2 P_2 P2。这种方式明显的缺点是物理空间中的直线在鱼眼图像上通常不再是直线，所以在我们的日常摄影中一般比较少用。关于球面投影透视和鱼眼图像的相关内容这里不多说，接下来我们主要讨论平面投影透视，以及给定两幅图像中一些相匹配的特征点的坐标，我们应该怎样求解相关的变换参数。
图2 鱼眼成像示例一般来说，我们会以镜头的光学中心作为原点建立坐标参考系，并且令成像平面为 z = 1 z=1 z=1，那么对于物理空间中的一点 X = ( X , Y , Z ) {\bf{X}}=(X,Y,Z) X=(X,Y,Z)，其经过平面投影透视在成像平面上的坐标可表示为 x = ( x , y , 1 ) = 1 Z ( X , Y , Z ) = 1 Z X . (1) {\bf{x}} = \left( {x,y,1} \right) = \frac{1}{Z}\left( {X,Y,Z} \right) = \frac{1}{Z}{\bf{X}}. \tag{1} x=(x,y,1)=Z1(X,Y,Z)=Z1X.(1)另外，我们可以使用一个 3x3 的矩阵 H {\bf{H}} H 来描述三维空间中的线性坐标变换，即 ( X 2 Y 2 Z 2 ) = ( h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ) ( X 1 Y 1 Z 1 ) . (2) \left( {\begin{array}{c} {{X_2}}\\ {{Y_2}}\\ {{Z_2}} \end{array}} \right) = \left( {\begin{array}{c} {{h_{11}}}&{{h_{12}}}&{{h_{13}}}\\ {{h_{21}}}&{{h_{22}}}&{{h_{23}}}\\ {{h_{31}}}&{{h_{32}}}&{{h_{33}}} \end{array}} \right)\left( {\begin{array}{c} {{X_1}}\\ {{Y_1}}\\ {{Z_1}} \end{array}} \right). \tag{2} ⎝⎛X2Y2Z2⎠⎞=⎝⎛h11h21h31h12h22h32h13h23h33⎠⎞⎝⎛X1Y1Z1⎠⎞.(2) 那么对于成像平面上的点，有 x 2 = 1 Z 2 X 2 = 1 Z 2 H X 1 = Z 1 Z 2 H x 1 . (3) {{\bf{x}}_2} = \frac{1}{{{Z_2}}}{{\bf{X}}_2} = \frac{1}{{{Z_2}}}{\bf{H}}{{\bf{X}}_1} = \frac{{{Z_1}}}{{{Z_2}}}{\bf{H}}{{\bf{x}}_1}. \tag{3} x2=Z21X2=Z21HX1=Z2Z1Hx1.(3)这就是二维平面图像的透视变换，我们一般称 H \bf{H} H 为单应性（Homography）矩阵。然而透视本身是不可逆的，如果仅有一张平面图像，我们无法知道其原来的 Z 1 {Z_1} Z1 值，所以我们通常令 Z 1 = 1 {Z_1=1} Z1=1。那么可得 x 2 = h 11 x 1 + h 12 y 1 + h 13 h 31 x 1 + h 32 y 1 + h 33 , y 2 = h 21 x 1 + h 22 y 1 + h 23 h 31 x 1 + h 32 y 1 + h 33 . (4) {x_2} = \frac{{{h_{11}}{x_1} + {h_{12}}{y_1} + {h_{13}}}}{{{h_{31}}{x_1} + {h_{32}}{y_1} + {h_{33}}}},{\rm{ }}{y_2} = \frac{{{h_{21}}{x_1} + {h_{22}}{y_1} + {h_{23}}}}{{{h_{31}}{x_1} + {h_{32}}{y_1} + {h_{33}}}}.\tag{4} x2=h31x1+h32y1+h33h11x1+h12y1+h13,y2=h31x1+h32y1+h33h21x1+h22y1+h23.(4)虽然式（2）表示的是一个线性变换，但在透视变换中，因为我们需要对 Z 2 {Z_2} Z2 进行归一化，而各个点的 Z 2 {Z_2} Z2 并不一定相等，这就引入了非线性的过程，所以透视变换不属于线性变换，也就是说第一个图像平面中的平行线在第二个图像平面中不一定平行。但由于采用的是平面投影透视，其中的直线依然还会保持直线。
考虑一种特殊的情况，假如令 h 31 = h 32 = 0 , h 33 = 1 h_{31}=h_{32}=0, h_{33}=1 h31=h32=0,h33=1，即令 A = ( a 11 a 12 a 13 a 21 a 22 a 23 0 0 1 ) . {\bf{A}}= \left( {\begin{array}{c} {{a_{11}}}&{{a_{12}}}&{{a_{13}}}\\ {{a_{21}}}&{{a_{22}}}&{{a_{23}}}\\ 0&0&1 \end{array}} \right). A=⎝⎛a11a210a12a220a13a231⎠⎞. 那么对于式（3），我们始终可以保证 Z 2 = 1 Z_2=1 Z2=1，这就意味着不再需要对变换后的坐标的 z z z 轴进行归一化，这时式（3）就代表了一个平面上的线性变换， z = 1 z=1 z=1 平面上的点经过变换后还是落在 z = 1 z=1 z=1 平面上，我们称这种变换为仿射（Affine）变换。由于仿射变换属于线性变换，平面上的平行线经过变换后依然保持平行。
从透视变换到仿射变换，实际反映了变换参数自由度的缩减。对于透视变换，其总共有 9 个可变化的参数，因此灵活性最高，通过变换可改变物体在平面上的基本形状。而仿射变换只剩下了 6 个可变参数，使得原本的非线性变换转换为线性变换，因此无法改变物体中线条的平行性，在一定程度上保留了物体的基本形状。同样，我们还可以继续缩减变换参数的自由度，例如只允许物体在平面上平移、旋转以及宽高等比例缩放，这时物体的形状除了大小以外没有任何改变，因此称为刚体（Rigid）或欧氏（Euclidean）变换，即 ( x 2 y 2 1 ) = ( s ⋅ cos θ − s ⋅ sin θ t x s ⋅ sin θ s ⋅ cos θ t y 0 0 1 ) ( x 1 y 1 1 ) . (5) \left( {\begin{array}{c} {{x_2}}\\ {{y_2}}\\ 1 \end{array}} \right) = \left( {\begin{array}{c} {s \cdot \cos \theta }&{ - s \cdot \sin \theta }&{tx}\\ {s \cdot \sin \theta }&{s \cdot \cos \theta }&{ty}\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{c} {{x_1}}\\ {{y_1}}\\ 1 \end{array}} \right). \tag{5} ⎝⎛x2y21⎠⎞=⎝⎛s⋅cosθs⋅sinθ0−s⋅sinθs⋅cosθ0txty1⎠⎞⎝⎛x1y11⎠⎞.(5)那么这时变换参数的自由度只剩下 4 了，由于。实际上，式（5）还可以继续分解为自由度更低的变换的组合，例如 ( x 2 y 2 1 ) = ( 1 0 t x 0 1 t y 0 0 1 ) ( s 0 0 0 s 0 0 0 1 ) ( cos θ − sin θ 0 sin θ cos θ 0 0 0 1 ) ( x 1 y 1 1 ) . (6) \left( {\begin{array}{c} {{x_2}}\\ {{y_2}}\\ 1 \end{array}} \right) = \left( {\begin{array}{c} 1&0&{tx}\\ 0&1&{ty}\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{c} s&0&0\\ 0&s&0\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{c} {\cos \theta }&{ - \sin \theta }&0\\ {\sin \theta }&{\cos \theta }&0\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{c} {{x_1}}\\ {{y_1}}\\ 1 \end{array}} \right).\tag{6} ⎝⎛x2y21⎠⎞=⎝⎛100010txty1⎠⎞⎝⎛s000s0001⎠⎞⎝⎛cosθsinθ0−sinθcosθ0001⎠⎞⎝⎛x1y11⎠⎞.(6) 也就是先进行只有 1 自由度的旋转，然后再进行 1 自由度的宽高等比例缩放，最后进行 2 自由度的水平或垂直的平移。实际上，式（5）可以表示任意次数与顺序的平移、旋转、宽高等比例缩放的变换组合的最终形式，但要注意，由于矩阵的乘法不满足交换律，相同的变换参数采用不同的变换顺序有可能导致不同的变换结果。
综上，虽然我们只使用了一个简单的 3x3 矩阵，但这已经足以描述我们日常生活中经常遇到的大部分简单的运动了，其具体形式如图3所示，包括平移、旋转、缩放、仿射以及透视变换等等。
图3 常见的变换与运动模型假设我们得到物体上的一些点（通常为特征点或显著点）在两幅平面图像上对应坐标，例如分别为 ( x i 1 , y i 1 ) (x_{i1}, y_{i1}) (xi1,yi1) 和 ( x i 2 , y i 2 ) (x_{i2}, y_{i2}) (xi2,yi2)，并且该物体的运动形式为前面图 3 所提到的几种之一，那么我们应该怎样根据这些匹配点对来估计相应的变换参数呢？虽然所有这些参数都是来自于一个 3x3 的矩阵，但由于透视变换是非线性变换，而仿射以及其他更低自由度的变换为线性变换，两者的参数优化过程还是有很大不同的。接下来我们先讨论仿射变换的参数优化过程，而其他更低自由度的变换的参数优化过程基本类似，不再赘述。
由于参数优化过程中坐标是已知量，而参数才是未知数，将坐标变换的式子重排可得 ( x 11 y 11 1 0 0 0 x 21 y 21 1 0 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 0 0 x 11 y 11 1 0 0 0 x 21 y 21 1 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ) ( a 11 a 12 a 13 a 21 a 22 a 23 ) = ( x 12 x 22 ⋮ y 12 y 22 ⋮ ) ⇔ A a = b . (7) \left( {\begin{array}{c} {{x_{11}}}&{{y_{11}}}&1&0&0&0\\ {{x_{21}}}&{{y_{21}}}&1&0&0&0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&0&0&{{x_{11}}}&{{y_{11}}}&1\\ 0&0&0&{{x_{21}}}&{{y_{21}}}&1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{array}} \right)\left( {\begin{array}{c} {{a_{11}}}\\ {{a_{12}}}\\ {{a_{13}}}\\ {{a_{21}}}\\ {{a_{22}}}\\ {{a_{23}}} \end{array}} \right) = \left( {\begin{array}{c} {{x_{12}}}\\ {{x_{22}}}\\ \vdots \\ {{y_{12}}}\\ {{y_{22}}}\\ \vdots \end{array}} \right) \Leftrightarrow {\bf{Aa}} = {\bf{b}}.\tag{7} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛x11x21⋮00⋮y11y21⋮00⋮11⋮00⋮00⋮x11x21⋮00⋮y11y21⋮00⋮11⋮⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎛a11a12a13a21a22a23⎠⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛x12x22⋮y12y22⋮⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞⇔Aa=b.(7)上述方程组为超定方程组，那么根据线性最小二乘法可得 a = ( A T A ) − 1 A T b , (8) {\bf{a}} = {\left( {{{\bf{A}}^T}{\bf{A}}} \right)^{ - 1}}{{\bf{A}}^T}{\bf{b}},\tag{8} a=(ATA)−1ATb,(8) 不过这种直接求逆的方法在实际应用过程中不一定是最佳的，我们可以将其转换为求解以下线性方程组 A T A a = A T b , (9) {{\bf{A}}^T}{\bf{Aa}} = {{\bf{A}}^T}{\bf{b}},\tag{9} ATAa=ATb,(9) 其中 A T A = ( ∑ i = 1 n x i 1 x i 1 ∑ i = 1 n x i 1 y i 1 ∑ i = 1 n x i 1 0 0 0 ∑ i = 1 n x i 1 y i 1 ∑ i = 1 n y i 1 y i 1 ∑ i = 1 n y i 1 0 0 0 ∑ i = 1 n x i 1 ∑ i = 1 n y i 1 ∑ i = 1 n 1 0 0 0 0 0 0 ∑ i = 1 n x i 1 x i 1 ∑ i = 1 n x i 1 y i 1 ∑ i = 1 n x i 1 0 0 0 ∑ i = 1 n x i 1 y i 1 ∑ i = 1 n y i 1 y i 1 ∑ i = 1 n y i 1 0 0 0 ∑ i = 1 n x i 1 ∑ i = 1 n y i 1 ∑ i = 1 n 1 ) , {{\bf{A}}^T}{\bf{A}} = \left( {\begin{array}{c} {\sum\limits_{i = 1}^n {{x_{i1}}{x_{i1}}} }&{\sum\limits_{i = 1}^n {{x_{i1}}{y_{i1}}} }&{\sum\limits_{i = 1}^n {{x_{i1}}} }&0&0&0\\ {\sum\limits_{i = 1}^n {{x_{i1}}{y_{i1}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}{y_{i1}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}} }&0&0&0\\ {\sum\limits_{i = 1}^n {{x_{i1}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}} }&{\sum\limits_{i = 1}^n 1 }&0&0&0\\ 0&0&0&{\sum\limits_{i = 1}^n {{x_{i1}}{x_{i1}}} }&{\sum\limits_{i = 1}^n {{x_{i1}}{y_{i1}}} }&{\sum\limits_{i = 1}^n {{x_{i1}}} }\\ 0&0&0&{\sum\limits_{i = 1}^n {{x_{i1}}{y_{i1}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}{y_{i1}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}} }\\ 0&0&0&{\sum\limits_{i = 1}^n {{x_{i1}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}} }&{\sum\limits_{i = 1}^n 1 } \end{array}} \right), ATA=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛i=1∑nxi1xi1i=1∑nxi1yi1i=1∑nxi1000i=1∑nxi1yi1i=1∑nyi1yi1i=1∑nyi1000i=1∑nxi1i=1∑nyi1i=1∑n1000000i=1∑nxi1xi1i=1∑nxi1yi1i=1∑nxi1000i=1∑nxi1yi1i=1∑nyi1yi1i=1∑nyi1000i=1∑nxi1i=1∑nyi1i=1∑n1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞, A T b = ( ∑ i = 1 n x i 1 x i 2 ∑ i = 1 n y i 1 x i 2 ∑ i = 1 n x i 2 ∑ i = 1 n x i 1 y i 2 ∑ i = 1 n y i 1 y i 2 ∑ i = 1 n y i 2 ) T . {{\bf{A}}^T}{\bf{b}} = {\left( {\begin{array}{c} {\sum\limits_{i = 1}^n {{x_{i1}}{x_{i2}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}{x_{i2}}} }&{\sum\limits_{i = 1}^n {{x_{i2}}} }&{\sum\limits_{i = 1}^n {{x_{i1}}{y_{i2}}} }&{\sum\limits_{i = 1}^n {{y_{i1}}{y_{i2}}} }&{\sum\limits_{i = 1}^n {{y_{i2}}} } \end{array}} \right)^T}. ATb=(i=1∑nxi1xi2i=1∑nyi1xi2i=1∑nxi2i=1∑nxi1yi2i=1∑nyi1yi2i=1∑nyi2)T.
上面我们可以看到仿射变换的参数优化过程还是比较简单的，那为什么透视变换的参数优化还要单独去说呢？实际上，由于透视变换引入了一个 z z z 轴归一化的过程，整个参数优化的过程也发生了很大的改变。将式（4）重排可得 ( x 11 y 11 1 0 0 0 − x 11 x 12 − y 11 x 12 − x 12 x 21 y 21 1 0 0 0 − x 21 x 22 − y 21 x 22 − x 22 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 0 0 x 11 y 11 1 − x 11 y 12 − y 11 y 12 − y 12 0 0 0 x 21 y 21 1 − x 21 y 22 − y 21 y 22 − y 22 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ) ( h 11 h 12 h 13 h 21 h 22 h 23 h 31 h 32 h 33 ) = ( 0 0 ⋮ 0 0 ⋮ ) . (10) \left( {\begin{array}{c} {{x_{11}}}&{{y_{11}}}&1&0&0&0&{ - {x_{11}}{x_{12}}}&{ - {y_{11}}{x_{12}}}&{ - {x_{12}}}\\ {{x_{21}}}&{{y_{21}}}&1&0&0&0&{ - {x_{21}}{x_{22}}}&{ - {y_{21}}{x_{22}}}&{ - {x_{22}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&0&0&{{x_{11}}}&{{y_{11}}}&1&{ - {x_{11}}{y_{12}}}&{ - {y_{11}}{y_{12}}}&{ - {y_{12}}}\\ 0&0&0&{{x_{21}}}&{{y_{21}}}&1&{ - {x_{21}}{y_{22}}}&{ - {y_{21}}{y_{22}}}&{ - {y_{22}}}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{array}} \right)\left( {\begin{array}{c} {{h_{11}}}\\ {{h_{12}}}\\ {{h_{13}}}\\ {{h_{21}}}\\ {{h_{22}}}\\ {{h_{23}}}\\ {{h_{31}}}\\ {{h_{32}}}\\ {{h_{33}}} \end{array}} \right) = \left( {\begin{array}{c} 0\\ 0\\ \vdots \\ 0\\ 0\\ \vdots \end{array}} \right).\tag{10} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛x11x21⋮00⋮y11y21⋮00⋮11⋮00⋮00⋮x11x21⋮00⋮y11y21⋮00⋮11⋮−x11x12−x21x22⋮−x11y12−x21y22⋮−y11x12−y21x22⋮−y11y12−y21y22⋮−x12−x22⋮−y12−y22⋮⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛h11h12h13h21h22h23h31h32h33⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛00⋮00⋮⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞.(10) 也就是说，透视变换的参数优化变成了一个求解线性齐次方程组 A h = 0 {\bf{Ah=0}} Ah=0 的问题。我们知道，对于线性齐次方程组，如果 h {\bf{h}} h 是它的解，那么 k h {k\bf{h}} kh 也是它的解，这从式（4）也可以看出，因为 z z z 轴需要进行归一化，所以 h {\bf{h}} h 的长度并不影响最终的结果。更一般地，对于线性非齐次方程组 A x = b {\bf{Ax=b}} Ax=b，假设 A {\bf{A}} A 的行数不少于列数，且列数为 n n n，秩为 r r r，有唯一解的充要条件是 R ( A , b ) = R ( A ) = r = n {R({\bf{A, b}})=R({\bf{A}})=r=n} R(A,b)=R(A)=r=n，而有无限多解的充要条件是 R ( A , b ) = R ( A ) = r < n {R({\bf{A, b}})=R({\bf{A}})=r<n} R(A,b)=R(A)=r<n，且对应的齐次方程组 A x = 0 {\bf{Ax=0}} Ax=0 的基础解系里的最大无关向量的个数为 n − r n-r n−r，那么 A x = b {\bf{Ax=b}} Ax=b 的通解可以表示为 A x = 0 {\bf{Ax=0}} Ax=0 基础解系的线性组合加上 A x = b {\bf{Ax=b}} Ax=b 的任意一个解。
在仿射变换矩阵的优化过程中，只要有三个点不共线，那么式（7）中左侧的矩阵都是列满秩的，那么从直觉上判断式（10）中的矩阵 A {\bf{A}} A 是否也一定是列满秩的？答案是否定的。很明显，如果式（10）中的矩阵 A {\bf{A}} A 是列满秩的，那么该方程组只有唯一解，而 0 {\bf{0}} 0 是任意线性齐次方程组的解，也就是该方程组只有零解。然而，如果我们所用的匹配点对确实是经过某个单应性矩阵变换而得到的，为什么从匹配点对反推回矩阵参数就只能得到零解了呢？原因在于此时矩阵 A {\bf{A}} A 的秩一般最多只有 8，而不是列满秩时的 9。因为这时的秩 r = 8 r=8 r=8，那么基础解系中的最大无关向量个数为 n − r = 1 n-r=1 n−r=1，刚好就对应了实际所用到的单应性矩阵参数。而之所以 A {\bf{A}} A 的秩为 8，是因为其在矩阵中引入了 x 2 x_2 x2 和 y 2 y_2 y2，这就是透视变换与仿射变换的不同之处，这两个值是与 x 1 x_1 x1 和 y 1 y_1 y1 相关的，也就是根据式（4）可得， h ′ 11 x 1 + h ′ 12 y 1 + h ′ 13 − h ′ 31 x 1 x 2 − h ′ 32 y 1 x 2 − h ′ 33 x 2 = 0 , h ′ 21 x 1 + h ′ 22 y 1 + h ′ 23 − h ′ 31 x 1 y 2 − h ′ 32 y 1 y 2 − h ′ 33 y 2 = 0. (11) \begin{array}{l} {{h'}_{11}}{x_1} + {{h'}_{12}}{y_1} + {{h'}_{13}} - {{h'}_{31}}{x_1}{x_2} - {{h'}_{32}}{y_1}{x_2} - {{h'}_{33}}{x_2} = 0,\\ {{h'}_{21}}{x_1} + {{h'}_{22}}{y_1} + {{h'}_{23}} - {{h'}_{31}}{x_1}{y_2} - {{h'}_{32}}{y_1}{y_2} - {{h'}_{33}}{y_2} = 0. \end{array}\tag{11} h′11x1+h′12y1+h′13−h′31x1x2−h′32y1x2−h′33x2=0,h′21x1+h′22y1+h′23−h′31x1y2−h′32y1y2−h′33y2=0.(11) 其中 h ′ {\bf{h'}} h′ 代表实际的单应性矩阵参数，不是式（10）中的未知数。那就意味着， A {\bf{A}} A 中的任意一个非零参数对应的列向量都可以表示成其他列向量的线性组合，那自然地其秩只有 8 了。极端情况下如果 h ′ {\bf{h'}} h′ 前两行 6 个参数都为零，那么 x 2 x_2 x2 和 y 2 y_2 y2 都为零，这时矩阵 A {\bf{A}} A 的后三列都为零，其秩就只有 6 了，不过这种情况一般不考虑。因为通解中只有一个解向量，而其长度并不影响最终的结果，那么我们可以令某个非零的分量为1，例如 h 33 {h_{33}} h33 通常是非零的，那么对应的列就可以移到右侧，如式（10）可表示为 B = ( a 1 a 2 ⋯ a 8 ) , x = ( h 11 h 12 ⋯ h 32 ) T , h 33 = 1 , B x = h 11 a 1 + h 12 a 2 + ⋯ + h 32 a 8 = − a 9 ⇒ B T B x = − B T a 9 . (12) \begin{array}{l} {\bf{B}} = \left( {\begin{array}{c} {{{\bf{a}}_1}}&{{{\bf{a}}_2}}& \cdots &{{{\bf{a}}_8}} \end{array}} \right),{\rm{ }}{\bf{x}} = {\left( {\begin{array}{c} {{h_{11}}}&{{h_{12}}}& \cdots &{{h_{32}}} \end{array}} \right)^T},{\rm{ }}{h_{33}} = 1,\\ {\bf{Bx}} = {h_{11}}{{\bf{a}}_1} + {h_{12}}{{\bf{a}}_2} + \cdots + {h_{32}}{{\bf{a}}_8} = - {{\bf{a}}_9} \Rightarrow {{\bf{B}}^T}{\bf{Bx}} = - {{\bf{B}}^T}{{\bf{a}}_9}. \end{array}\tag{12} B=(a1a2⋯a8),x=(h11h12⋯h32)T,h33=1,Bx=h11a1+h12a2+⋯+h32a8=−a9⇒BTBx=−BTa9.(12) 这时 B {\bf{B}} B 是列满秩的，且有 R ( B , a 9 ) = R ( B ) = 8 {R({\bf{B, a}}_9)=R({\bf{B}})=8} R(B,a9)=R(B)=8，因此有唯一解，且对应于 h ′ {\bf{h'}} h′。因为这时只有 8 个未知数，我们只需要 4 个三点互不共线的匹配点对即可求出相应的单应性矩阵参数，这种方法常用于将某幅图像投影到另外一幅图像上，因为这时我们只需要知道图像四个角对应的位置即可。
不过要注意的是，我们不能选取那些本来为 0 的分量让其等于 1，虽然一般来说 h 33 h_{33} h33 是非零的，但单应性矩阵本身是三维空间中的自由线性变换，除了不让其第三行全为 0 也就是让式（4）可除以外，确实也存在 h 33 {h_{33}} h33 可能为零的的情况。这时，如果 h 33 = 0 {h_{33}=0} h33=0，那么式（11）表示矩阵 A {\bf{A}} A 的前 8 列是线性相关的，也就是 R ( B ) = 7 R({\bf{B}})=7 R(B)=7，而 R ( B , a 9 ) = 8 ≠ R ( B ) {R({\bf{B, a}}_9)=8≠R({\bf{B}})} R(B,a9)=8=R(B)，因此这时式（12）无解。这本身也是因为 0 是不可能通过缩放得到 1 的。然而，由于我们事前是不能知道参数的哪个分量是非零的，我们想要避免将本身为零的分量设为 1，最好的方法还是通过初等行变换将矩阵 A {\bf{A}} A 转换为行最简形矩阵，这样哪些分量为 0 就十分清楚了。另外，我们也可以通过限制解向量的长度来避免将本来为 0 的分量设为 1，例如 ∣ ∣ h ∣ ∣ = 1 {||{\bf{h}}||=1} ∣∣h∣∣=1。但这时由于多了一个二次型约束条件，我们应该怎么去求解？
前面我们虽然讨论了很多，实际上都是基于所有匹配点都对应同一个单应性矩阵这一假设的，但在优化问题中，我们通常只能最小化匹配误差，而不能让其为零。也就是说，矩阵 A {\bf{A}} A 通常是列满秩的，那么这时候我们应该怎么找到一个非零解使得匹配误差最小？对于线性齐次方程组，我们不能直接使用最小二乘法，因为如果要最小化 ∣ ∣ A h ∣ ∣ ||{\bf{Ah}}|| ∣∣Ah∣∣，那很明显只有 h = 0 \bf{h=0} h=0 这个极值点了，不符合我们要优化的问题。类似于前面的分析，如果 A h ≈ 0 {\bf{Ah}} \approx {\bf{0}} Ah≈0，那么对于较小的 k k k 值，也有 A ⋅ k h ≈ 0 {\bf{A}} \cdot k{\bf{h}} \approx {\bf{0}} A⋅kh≈0。因此，解决这个问题的一种简单的方法是令某个未知数为常量，如 h 33 = 1 h_{33}=1 h33=1，那么其在 A {\bf{A}} A 中对应的那一列就可以移到方程组的右侧，这样就得到了一个非齐次的的线性方程组，接下来就可以使用最小二乘法来求解最优解了。但是，这种方法的局限同样在于，如果实际非零最优解（后面有介绍）中 h 33 h_{33} h33 的相比于其他分量的量级很小或者等于 0，那么 k k k 值就很大或者无解，这样造成的误差就不一定能忽视了。因此，我们需要从其他的角度来分析。
我们知道，函数从不同方向接近极值点时其下降速度并不一定相等。以简单的二次型（只包含 x i x j x_ix_j xixj 项的二次函数） f ( x ) = f ( x , y ) = a x 2 + b y 2 (13) f({\bf{x}}) = f(x,y) = a{x^2} + b{y^2}\tag{13} f(x)=f(x,y)=ax2+by2(13) 为例，假设 a > b > 0 a > b > 0 a>b>0，那么其等高线为一系列的同心椭圆，并且原点为唯一的极小值点，如图 4 所示。
图4 椭圆等高线如果我们使用一系列垂直于 X O Y XOY XOY 平面的二维平面 y = k x {y=kx} y=kx 以及 x = 0 x = 0 x=0 去切开该函数，那么在剖面处我们可以得到一系列的二次曲线。令 z = ∣ ∣ x ∣ ∣ ⋅ s i g n ( y ) z=||{\bf{x}}||·sign(y) z=∣∣x∣∣⋅sign(y)，也就是该点在该方向上离原点的距离，那么这些曲线可以表示为 f ( z ) = c z 2 f(z)=cz^2 f(z)=cz2，其中 b ≤ c ≤ a b \le c \le a b≤c≤a，如图 5 所示。也就是说，虽然该二次型 f ( x ) {f(x)} f(x) 只有原点这么一个极小值点，但是从不同方向趋近原点时，它们的下降速度是不一样的。而且从图 5 可以发现，无论我们选择哪个方向，其下降速度都会处在 a a a 和 b b b 之间，而这两者刚好是组成二次型的两个值，如果从二次型矩阵的角度来讲，这刚好就是该矩阵的两个特征值。这是否存在着某些巧合？
图5 剖面二次曲线回到我们前面所讨论的参数优化问题，我们需要最小化 ∣ ∣ A h ∣ ∣ ||{\bf{Ah}}|| ∣∣Ah∣∣，其实也就对应着最小化二次型 f ( h ) = h T A T A h . (14) f({\bf{h}}) = {{\bf{h}}^T}{{\bf{A}}^T}{\bf{Ah}}.\tag{14} f(h)=hTATAh.(14) 虽然其比式（13）的二次型复杂很多，如果我们采用类似的方法使用一些超平面去切割该函数，得到的剖面曲线跟图 5 也是差不多的。这是因为 f ( k h ) = k 2 f ( h ) {f(k{\bf{h}})=k^2f({\bf{h}})} f(kh)=k2f(h)，所以这些曲线其实都是关于缩放系数 k {k} k 的只含二次项的一元二次函数，而因为 A \bf{A} A 并不一定列满秩， A T A {\bf{A}}^T{\bf{A}} ATA 可能是正定或者半正定的，也就是 f ( k h ) ≥ 0 {f(k{\bf{h}}) \ge 0} f(kh)≥0，那么得到的剖面曲线要么是开口向上的二次曲线，要么为 0。前面已经提到，如果 A h ≈ 0 {\bf{Ah}} \approx {\bf{0}} Ah≈0，那么对于较小的 k {k} k 值，也有 A ⋅ k h ≈ 0 {\bf{A}} \cdot k{\bf{h}} \approx {\bf{0}} A⋅kh≈0。从优化的角度来讲，在保持较低误差的情况下， k {k} k 的取值范围应该越大越好，也就是 k k k 对优化误差的影响并不明显。那么，结合图 5，我们应该选择那些最平缓的曲线，那么即便 k k k 值变得很大，实际的误差也不会特别的明显。而对于那些特别陡峭的曲线，我们稍微增大 k k k 值，也有可能导致误差超出我们的预期。因此，式（14）的问题就可以转化为求解可以得到最平缓的剖面曲线的平面方向了。为了描述这些曲线的平缓程度，我们可以比较在各个方向上离原点同样距离（如 ∣ ∣ h ∣ ∣ = 1 ||{\bf{h}}||=1 ∣∣h∣∣=1）时函数值的大小，如果该值为 0，那就意味着 A T A {\bf{A}}^T{\bf{A}} ATA 是半正定的，对应着有一个精确解使得匹配误差为 0。这样，我们实际上就可以得到一个等式约束的最优化问题，即 m i n i m i z e f ( h ) = h T A T A h s u b j e c t t o h T h = 1 (15) \begin{array}{l} {\rm{minimize }}f({\bf{h}}) = {{\bf{h}}^T}{{\bf{A}}^T}{\bf{Ah}}\\ {\rm{subject to }}{{\bf{h}}^T}{\bf{h}} = 1 \end{array}\tag{15} minimizef(h)=hTATAhsubjecttohTh=1(15) 对于等式约束的最优化问题，引入拉格朗日算子 λ \lambda λ，构造拉格朗日函数得 l ( h , λ ) = h T A T A h + λ ( 1 − h T h ) . (16) l\left( {{\bf{h}},\lambda } \right) = {{\bf{h}}^T}{{\bf{A}}^T}{\bf{Ah}} + \lambda \left( {1 - {{\bf{h}}^T}{\bf{h}}} \right).\tag{16} l(h,λ)=hTATAh+λ(1−hTh).(16) 根据拉格朗日条件，如果 h ∗ \bf{h^*} h∗ 是 f ( h ) f(\bf{h}) f(h) 在约束条件下的一个极小值点，存在 λ ∗ \lambda^* λ∗，使得 ∇ l ( h ∗ , λ ∗ ) = ( D h l ( h ∗ , λ ∗ ) T D λ l ( h ∗ , λ ∗ ) T ) = ( 2 A T A h ∗ − 2 λ ∗ h ∗ 1 − h ∗ T h ∗ ) = 0. (17) \nabla l\left( {{{\bf{h}}^*},{\lambda ^*}} \right) = \left( {\begin{array}{c} {{D_{\bf{h}}}l{{\left( {{{\bf{h}}^*},{\lambda ^*}} \right)}^T}}\\ {{D_\lambda }l{{\left( {{{\bf{h}}^*},{\lambda ^*}} \right)}^T}} \end{array}} \right) = \left( {\begin{array}{c} {2{{\bf{A}}^T}{\bf{A}}{{\bf{h}}^*} - 2{\lambda ^*}{{\bf{h}}^*}}\\ {1 - {{\bf{h}}^{*T}}{{\bf{h}}^*}} \end{array}} \right) = {\bf{0}}.\tag{17} ∇l(h∗,λ∗)=(Dhl(h∗,λ∗)TDλl(h∗,λ∗)T)=(2ATAh∗−2λ∗h∗1−h∗Th∗)=0.(17) 于是有 A T A h ∗ = λ ∗ h ∗ , h ∗ T h ∗ = 1. (18) \begin{array}{c} {{\bf{A}}^T}{\bf{A}}{{\bf{h}}^*} = {\lambda ^*}{{\bf{h}}^*},\\ {{\bf{h}}^{*T}}{{\bf{h}}^*} = 1. \end{array}\tag{18} ATAh∗=λ∗h∗,h∗Th∗=1.(18) f ( h ∗ ) = h ∗ T A T A h ∗ = λ ∗ h ∗ T h ∗ = λ ∗ . (19) f({{\bf{h}}^*}) = {{\bf{h}}^{*T}}{{\bf{A}}^T}{\bf{A}}{{\bf{h}}^*} = {\lambda ^*}{{\bf{h}}^{*T}}{{\bf{h}}^*} = {\lambda ^*}.\tag{19} f(h∗)=h∗TATAh∗=λ∗h∗Th∗=λ∗.(19)也就是说，如果 h ∗ \bf{h^*} h∗ 是 f ( h ) f(\bf{h}) f(h) 在约束条件下的极小值点，那么其必定为式（18）中 A T A {\bf{A}}^T{\bf{A}} ATA 的特征向量，同时 f ( h ∗ ) f(\bf{h^*}) f(h∗) 刚好为对应的特征值。这就是为什么前面我们提到的图 5 里面为什么下降最快和最慢的刚好就是两个特征值对应的曲线。
其实，拉格朗日条件只是一个一阶必要条件，也就是说满足式（18）的特征向量 h ∗ \bf{h^*} h∗ 还不一定是一个极值点。因为 A T A {\bf{A}}^T{\bf{A}} ATA 是对称正定或半正定矩阵，其特征值必定非负，且不同特征值对应的特征向量相互正交，而重根特征值的特征向量也可进行正交化得到相同数量的正交特征向量，所以这些特征向量可以组成一个标准正交基。那么，符合约束条件的任意向量可表示为 y = ∑ i b i h i , ∑ i b i 2 = 1. (20) {\bf{y}} = \sum\limits_i {{b_i}{{\bf{h}}_i}} ,{\rm{ }}\sum\limits_i {b_i^2} = 1.\tag{20} y=i∑bihi,i∑bi2=1.(20) 可得 y T A T A y = ∑ i b i h i T ⋅ ∑ i λ i b i h i = ∑ i λ i b i 2 . (21) {{\bf{y}}^T}{{\bf{A}}^T}{\bf{Ay}} = \sum\limits_i {{b_i}{\bf{h}}_i^T} \cdot \sum\limits_i {{\lambda _i}{b_i}{{\bf{h}}_i}} = \sum\limits_i {{\lambda _i}b_i^2} .\tag{21} yTATAy=i∑bihiT⋅i∑λibihi=i∑λibi2.(21) 如果将特征值从小到大排列，那么有 h 1 T A T A h 1 = λ 1 ≤ y T A T A y ≤ λ 9 = h 9 T A T A h 9 . (22) {\bf{h}}_1^T{{\bf{A}}^T}{\bf{A}}{{\bf{h}}_1} = {\lambda _1} \le {{\bf{y}}^T}{{\bf{A}}^T}{\bf{Ay}} \le {\lambda _9} = {\bf{h}}_9^T{{\bf{A}}^T}{\bf{A}}{{\bf{h}}_9}.\tag{22} h1TATAh1=λ1≤yTATAy≤λ9=h9TATAh9.(22) 因此，最小特征值对应的特征向量一定可以保证目标函数值最小，最大特征值对应的特征向量一定可以保证目标函数值最大。如果最小或者最大特征值是唯一的，没有重根，那么就可以保证对应的特征向量是全局极小值点或者全局极大值点。例如当 A \bf{A} A 的秩只有 8，也就是可以得到精确解时，刚好有一个 0 特征值，且仅对应唯一特征向量，这时可以保证其是全局的最优解，也就是精确解。
综上所述，透视变换的参数优化问题可以转化为求解式（18）的特征向量的问题，而且最优的参数向量为最小特征值对应的特征向量。由于特征向量在某些分量上的值可能很小或者为零，那么在最初提到的直接让 h 33 h_{33} h33 为非零常量的方法中，其首先就排除了那些 h 33 ≈ 0 h_{33} \approx 0 h33≈0 的特征向量方向，并且不同方向上的 ∣ ∣ h ∣ ∣ ||\bf{h}|| ∣∣h∣∣ 不再一致，这样很难比较不同方向的真正平缓程度，求出来的最优解也不一定是使用式（18）求出来的特征向量解。所以，令 h 33 h_{33} h33 为常量这种方法虽然比较简单，但具有比较大的局限性。另外，因为这两种方法都引入了一个约束条件，所以透视变换虽然有 9 个参数，但实际只有 8 自由度。
例如我们有以下一张图片，由于手机摄像头的焦距比较短，所以有明显的透视效应。我们想要更换屏幕上的壁纸，相当于将一张四平八方的图片通过透视变换投影到屏幕包围的区域。前面已经提到，要求解相应的透视变换参数，最少只需要 4 个匹配点对即可。所以，我们只需找到屏幕四个角对应的坐标，然后加上我们想要更换的图片的四个角的坐标，构造矩阵 A T A {\bf{A}}^T{\bf{A}} ATA，求解特征值和特征向量，取最小特征值对应的特征向量作为变换的参数即可。以下是相应的程序脚本以及变换后的效果。
# -*- coding: utf-8 -*- import numpy as np from PIL import Image def solve_homo(ptsrc, ptdst): npts = ptsrc.shape[0] A = np.zeros([2*npts, 9]) A[:npts, [0, 1]] = np.copy(ptsrc) A[:npts, 2] = np.ones(npts) A[npts:, [3, 4]] = np.copy(ptsrc) A[npts:, 5] = np.ones(npts) A[:npts, 6] = -ptsrc[:, 0] * ptdst[:, 0] A[npts:, 6] = -ptsrc[:, 0] * ptdst[:, 1] A[:npts, 7] = -ptsrc[:, 1] * ptdst[:, 0] A[npts:, 7] = -ptsrc[:, 1] * ptdst[:, 1] A[:npts, 8] = -ptdst[:, 0] A[npts:, 8] = -ptdst[:, 1] A = A.T.dot(A) vals, vecs = np.linalg.eig(A) homo_mat = vecs[:, np.argmin(vals)] if np.max(np.abs(homo_mat)) < 1e-7: homo_mat[-1] = 1 if np.abs(homo_mat[-1]) > 1e-7: homo_mat /= homo_mat[-1] return homo_mat if __name__ == '__main__': src = Image.open('src.jpg').convert('RGBA') # 希望更换的壁纸 dst = Image.open('dst.jpg').convert('RGBA') # 目标图 w, h = src.size W, H = dst.size # 给出对应坐标对， 源坐标和目标坐标均为 Nx2 的数组 ptsrc = np.array([[0, 0], [w-1, 0], [0, h-1], [w-1, h-1]]) ptdst = np.array([[610, 210], [3322, 388], [378, 1594], [2988, 2268]]) # 映射到目标图实际是为目标图上的每个像素找出对应的原图坐标然后插值， # 所以这里实际上应该是求目标图到原始图的变换矩阵 homo_mat = solve_homo(ptdst, ptsrc) # 利用PIL变换到目标图，注意其默认h33=1，所以只输入前8个参数，并且需要进行缩放 img = src.transform((W, H), Image.PERSPECTIVE, homo_mat[:-1], Image.BICUBIC) # 利用alpha通道将两张目标图粘贴到一起即可得到最终的图片 dst.paste(img, img.split()[-1]) dst = dst.convert('RGB') dst.save('new_wallpaper.jpg')