射影变换

射影变换（projective transformation）是在射影空间中由可逆线性映射诱导的变换。对实射影平面，取可逆 $3\times 3$ 矩阵 $H$，它把齐次点 $[p]$ 送到 $[Hp]$。平面上的这种变换也称为单应变换（homography）。它保持点、线及相交关系，但一般不保持平行、长度、角度或中点。

一个具体例子是 $$H=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 1& 1\end{array}\right),\left(\begin{array}{c}{X}^{\prime }\\ Y^{\prime }\\ W^{\prime }\end{array}\right)=H\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)=\left(\begin{array}{c}x\\ y\\ y+1\end{array}\right).$$ 只要 $y\ne -1$，普通坐标中的关系就是 $(x^{\prime },y^{\prime })=(x/(y+1),y/(y+1))$。例如 $(1,1)$ 被送到 $(1/2,1/2)$，而 $(1,0)$ 保持不动。同一幅图中不同位置的缩放比例不同，这正是它与一般仿射变换的区别之一。

可逆矩阵在齐次点上的作用

齐次坐标把非零向量的非零倍数视为同一点。若将代表 $p$ 换为 $\lambda p$，则 $H(\lambda p)=\lambda Hp$，所以输出点不变。可逆性保证 $Hp\ne 0$，输出总是合法的齐次点；逆矩阵 $H^{-1}$ 则给出逆变换。这证明了该规则在整个射影平面上定义良好且为双射。

矩阵 $H$ 与 $\mu H$（$\mu \ne 0$）表示同一射影变换。九个矩阵元素模掉一个共同尺度，使一般平面单应变换具有八个自由度。可以把某个非零元素归一化为一，但不能预先断定右下角 $h_{33}$ 非零。两次变换的复合按矩阵乘法表示：先 $H_1$、后 $H_2$，结果是 $H_2{H}_1$。

写 $H=(h_{ij})$，有限点的像在 $W^{\prime }\ne 0$ 时为 $$x^{\prime }=\frac{h_{11}x+h_{12}y+h_{13}}{h_{31}x+h_{32}y+h_{33}},y^{\prime }=\frac{h_{21}x+h_{22}y+h_{23}}{h_{31}x+h_{32}y+h_{33}}.$$ 分母为零表示像落在目标坐标片的无穷远直线上，并不表示射影变换本身在那里没有定义。若 $h_{31}=h_{32}=0$，可逆性迫使 $h_{33}\ne 0$，归一化后便得到可逆仿射变换；一般情形的分母随位置变化。

直线按逆转置变换

设源直线的非零系数向量为 $\ell$，其点满足 ${\ell }^{𝖳}p=0$。由于 $p^{\prime }=Hp$，可以代入 $p=H^{-1}{p}^{\prime }$： $$0={\ell }^{𝖳}{H}^{-1}{p}^{\prime }=(H^{-𝖳}\ell {)}^{𝖳}{p}^{\prime }.$$ 因此像直线的坐标可取 $${\ell }^{\prime }=H^{-𝖳}\ell ,H^{-𝖳}:=(H^{-1}{)}^{𝖳}.$$ 可逆性保证 ${\ell }^{\prime }\ne 0$。反过来，若一个目标点满足 $\ell {\text{'}}^{𝖳}{p}^{\prime }=0$，其逆像必满足 ${\ell }^{𝖳}{H}^{-1}{p}^{\prime }=0$，故像恰好是整条直线，而不仅是它的一部分。这个双向推导证明了直线与点线关联的保持。

不同点的像仍不同，否则用 $H^{-1}$ 会得到两个原点相同。于是两点的连线被送到像点的连线，两条不同线的交点被送到像线的交点。这个结论同时适用于有限交点与无穷远交点。直接把线系数乘以 $H$ 一般是错的；点和线虽然都写成三元列向量，变换法则却不同。

在开篇算例中， $$H^{-1}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& -1& 1\end{array}\right),H^{-𝖳}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& -1\\ 0& 0& 1\end{array}\right).$$ 例如直线 $x=1$ 的系数为 $\ell =(1,0,-1{)}^{𝖳}$，因此 ${\ell }^{\prime }=(1,1,-1{)}^{𝖳}$，像直线为 $x^{\prime }+y^{\prime }=1$。把像点 $(1,0)$ 和 $(1/2,1/2)$ 代入，均满足该式。

正方形、消失点与消失线

取单位正方形顶点 $A=(0,0)$、$B=(1,0)$、$C=(1,1)$、$D=(0,1)$。按上述矩阵逐点计算， $$A^{\prime }=(0,0),B^{\prime }=(1,0),C^{\prime }=(1/2,1/2),D^{\prime }=(0,1/2).$$ 四条边仍在直线上，像围成一个梯形。正方形内部也映到该梯形内部：由 $y^{\prime }=y/(y+1)$ 得 $y=y^{\prime }/(1-y^{\prime })$；再由 $x=x^{\prime }/(1-y^{\prime })$，条件 $0\le x,y\le 1$ 恰好变为 $0\le y^{\prime }\le 1/2$、$0\le x^{\prime }\le 1-y^{\prime }$。

虚线表示原正方形，着色梯形表示像；A′=A、B′=B。图中C到C′的箭头只标示这一个点的对应。

源平面中的竖直平行线 $x=a$，齐次系数为 $(1,0,-a{)}^{𝖳}$，像线为 $$x^{\prime }+a{y}^{\prime }-a=0.$$ 所有这些线都经过 $V=(0,1)$。原来竖直方向的无穷远点 $v_{\mathrm{\infty }}=(0,1,0{)}^{𝖳}$ 被送到 $H{v}_{\mathrm{\infty }}=(0,1,1{)}^{𝖳}$，这就是该方向的消失点。在普通坐标中，也可固定 $x=a$ 并令 $y\to +\mathrm{\infty }$，得到 $x^{\prime }\to 0$、$y^{\prime }\to 1$。没有一个源有限点被送到这个消失点，因为它的逆像是无穷远点。

所有源无穷远点组成源无穷远线 ${\ell }_{\mathrm{\infty }}=(0,0,1{)}^{𝖳}$。其像为 $${\ell }_{\mathrm{v}}=H^{-𝖳}{\ell }_{\mathrm{\infty }}=(0,-1,1{)}^{𝖳},$$ 即目标平面的直线 $y^{\prime }=1$。这是源仿射平面的消失线；当源平面表示地面时，通常称为地平线。它不必是图像中的水平线，这取决于坐标与成像姿态；本例恰好水平。某些方向的消失点也可能仍在目标无穷远处，例如 $H(1,0,0{)}^{𝖳}=(1,0,0{)}^{𝖳}$，所以水平平行线在本例中仍保持平行。

三条像线分别为x′=0、x′+y′=1、x′+2y′=2。延长线交于V；金色虚线是整个源无穷远直线的像。

另有一条容易混淆的线：哪些源有限点被送到目标无穷远处由 $W^{\prime }=0$ 决定。它的系数是 $H^{𝖳}{\ell }_{\mathrm{\infty }}=(0,1,1{)}^{𝖳}$，即源直线 $y=-1$。例如 $$H(1,-1,1{)}^{𝖳}=(1,-1,0{)}^{𝖳}.$$ 源中的 $y=-1$ 与目标中的消失线 $y^{\prime }=1$ 属于不同平面、承担相反的像与逆像关系，不能混成同一条“分母为零的地平线”。

这里t是源纵坐标y，u是像纵坐标y′。竖虚线t=−1对应像落在无穷远处；横虚线u=1是源无穷远点的像，两者不是同一个坐标值。

精确有理数实现

下面代码实现已知矩阵的点变换、线变换和逆变换。输入为整数或 Fraction，矩阵必须为可逆的 $3\times 3$ 矩阵，点和线必须是非零三元组。它不修改输入；输出按“第一个非零分量归一化为一”选取代表，因此也适用于末坐标为零的点。affine 单独进行有限坐标的除法，像在无穷远处时明确拒绝这一步。

Python 3

from fractions import Fraction


def rational(value):
    if isinstance(value, bool) or not isinstance(value, (int, Fraction)):
        raise TypeError("use int or Fraction, not floating-point input")
    return Fraction(value)


def vector3(values):
    v = tuple(rational(x) for x in values)
    if len(v) != 3 or not any(v):
        raise ValueError("a nonzero homogeneous triple is required")
    return v


def canonical(values):
    v = vector3(values)
    scale = next(x for x in v if x)
    return tuple(x / scale for x in v)


def affine(values):
    x, y, w = vector3(values)
    if w == 0:
        raise ValueError("point is at infinity in the W=1 chart")
    return x / w, y / w


def cross(u, v):
    return (u[1]*v[2] - u[2]*v[1],
            u[2]*v[0] - u[0]*v[2],
            u[0]*v[1] - u[1]*v[0])


def matvec(matrix, vector):
    return tuple(sum((a*b for a, b in zip(row, vector)), Fraction(0))
                 for row in matrix)


class Homography:
    def __init__(self, matrix):
        rows = tuple(tuple(rational(x) for x in row) for row in matrix)
        if len(rows) != 3 or any(len(row) != 3 for row in rows):
            raise ValueError("matrix must be 3 by 3")
        self.H = rows
        self.C = (cross(rows[1], rows[2]), cross(rows[2], rows[0]),
                  cross(rows[0], rows[1]))
        self.det = sum((rows[0][j] * self.C[0][j] for j in range(3)),
                       Fraction(0))
        if self.det == 0:
            raise ValueError("singular matrix is not a projectivity")

    def point(self, values):
        return canonical(matvec(self.H, vector3(values)))

    def line(self, values):
        # C = det(H) * H^(-T); the common nonzero scale cancels.
        return canonical(matvec(self.C, vector3(values)))

    def inverse(self):
        # C^T = det(H) * H^(-1), equally valid projectively.
        return Homography(tuple(zip(*self.C)))


H = Homography(((1, 0, 0), (0, 1, 0), (0, 1, 1)))
q = H.point((1, 1, 1))
print(tuple(map(str, affine(q))))                # ('1/2', '1/2')
print(tuple(map(str, H.line((1, 0, -1)))))        # ('1', '1', '-1')
print(H.inverse().point(q) == canonical((1, 1, 1)))  # True
J = Homography(((1, 0, 0), (0, 0, 1), (0, 1, 0)))
print(tuple(map(str, affine(J.point((2, 3, 1))))))  # ('2/3', '1/3')

代码中的 $C$ 是代数余子式矩阵，三行依次由另两行的叉积构成。记 $H$ 的行为 $r_i$、$C$ 的行为 $c_i$，则 $r_i\cdot c_i=detH$，因为三重积的循环排列不变；当 $j\ne i$ 时，$r_j$ 正是构成 $c_i$ 的一个叉积因子，故 $r_j\cdot c_i=0$。逐项相乘便得到 $H{C}^{𝖳}=(detH)I$，于是 $C^{𝖳}=(detH)H^{-1}$、$C=(detH)H^{-𝖳}$。共同的非零因子不改变齐次结果。矩阵非奇异保证两种乘法都不会把非零输入变为零；归一化只去掉一个公共尺度。这给出了点变换、线变换及逆变换的正确性理由，也说明代码无需固定 $h_{33}=1$。

对固定的三维齐次表示，建立矩阵、变换一个点或一条线以及构造逆变换，各需要常数次有理数运算和常数个有理数存储。处理 $N$ 个对象需要 $O(N)$ 次这样的运算；若保留所有结果，输出占用 $O(N)$ 个三元组。精确有理数的分子、分母可能很长，实际位运算成本还随其长度增长，不能把任意大整数运算视为固定时间。实现拒绝普通浮点数，包括 NaN 和无穷；用精确整数或分数可避免数值溢出和“接近零”判定的歧义，但不能消除输入观测误差。

边界、变式与矩阵估计的条件

代码最后的矩阵 $J$ 有 $j_{33}=0$，但 $detJ=-1$，所以它是合法射影变换。它把 $[X:Y:W]$ 送到 $[X:W:Y]$，因此 $(2,3)$ 的像为 $(2/3,1/3)$。把它强行除以右下角来“归一化”会错误地丢掉这个变换。若将主例点改为 $(2,-1)$，$H$ 的像为 $[2:-1:0]$；point 可以返回它，但 affine 会拒绝除以零。

可逆性也有实际作用。奇异矩阵 $\mathrm{diag}(1,1,0)$ 把 $[0:0:1]$ 送到零向量，其余点的像又都落在同一条无穷远直线上，不能定义整个射影平面的双射。代码因此在构造时就拒绝奇异矩阵。

射影变换虽然保持共线，一般不保持中点。主例中 $A=(0,0)$、$D=(0,1)$ 的中点为 $M=(0,1/2)$，其像是 $M^{\prime }=(0,1/3)$；像端点的普通中点却是 $(0,1/4)$。两者不同，所以不能在透视变换后直接沿用原来的仿射比例。

若 $H$ 未知，点对应 $p_i\leftrightarrow {p_i}^{\prime }$ 给出约束 ${p_i}^{\prime }\times H{p}_i=0$，对每一对非零齐次点只有两条独立线性约束。直接线性变换法（DLT）把这些约束堆叠为 $Ah=0$，其中 $h$ 排列矩阵九个元素。在无噪声情况下，唯一到尺度的候选要求 $\mathrm{rank}A=8$，且重排后的候选矩阵仍须可逆；至少需要四对点，四对点的常用非退化配置要求源、目标两边都无三点共线。仅有“四对”这个数量并不能保证可解。带噪声时常求最小奇异值对应的方向，但代数残差不等于几何重投影误差。本篇的精确代码只应用已知变换，没有实现或证明完整 DLT 估计流程。

单应变换用于平面图像校正与拼接时，还需检查场景模型。理想透视相机拍摄同一平面、投影中心不在该平面上，或两幅图像具有相同投影中心，是产生单应关系的典型条件；相机平移且场景具有一般深度变化时，不能把全部点强行套进同一个平面单应模型。