图像坐标空间变换：仿射变换（Affine Transformation）

仿射变换（Affine Transformation）简介

一般来说，仿射变换是指将一个向量空间进行一次线性变换（linear transformation）和一次平移（translation），变换到另一个向量空间的操作。
通常仿射变换对点、线、面具有一定的保持性，这种保持性体现在：

• 经过变换后，点还是点，线还是线，面还是面（如果不是仿射变换而是3D投影的话，在一定的视角下，面有可能变成线，线有可能变成点）
• 经过变换后，平行线和平行面依然平行
• 经过变换后，图形间的某些比例关系保持不变，比如两条平行线的长度比不变，点在线段中的位置比例保持不变

注意：仿射变换不能保证角度不变

仿射变换的基础类型

仿射变换的基础类型包括恒等（identity），尺度（scaling），旋转（rotation），剪切（shear），镜像（reflection），平移（translation），除了平移以外，其他都属于线性变换。线性变换通过矩阵乘法实现，平移通过向量加法实现。
以上基础变换类型主要从变换效果上进行分类，从数学上讲并不是很严格，因为存在包含关系，比如恒等变换可以理解为宽高方向放缩倍数均为1的尺度变换。不过也不必对此过分纠结，因为从效果进行分类更符合我们对图像的日常认知。

接下来我们使用公式对各种变换进行描述，公式采用两种形式进行表达：线性方程组和矩阵。
开始前我们需要进行符号约定：(u,v)用来表示原始图像中的坐标，(x,y)用来表示变换后图像的坐标。

在矩阵乘法的表达式中，我们将坐标以行向量的形式放在左侧，将变换矩阵放在右侧。这样做的好处是使用numpy或者tensorflow这种具有向量化的矩阵乘法算子时，不必再对坐标矩阵调整通道并reshape。当然如果你不在乎这类写代码的便利性的话，那么将变换矩阵放在左侧，将坐标以列向量形式放在右侧也没问题，而且确实有一些教科书是这么做的。

恒等

$$\begin{gathered} \{\begin{array}{l}x=u\\ y=v\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right] \end{gathered}$$

尺度

$$\begin{gathered} \{\begin{array}{l}x=\alpha \ast u\\ y=\beta \ast v\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}\alpha & 0\\ 0& \beta \end{array}\right] \end{gathered}$$

旋转

$$\begin{gathered} \{\begin{array}{l}x=cos(\theta )\ast u-sin(\theta )\ast v\\ y=sin(\theta )\ast u+cos(\theta )\ast v\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}cos(\theta )& sin(\theta )\\ -sin(\theta )& cos(\theta )\end{array}\right] \end{gathered}$$

剪切

分为水平剪切和垂直剪切，两种分开来写

水平剪切

$$\begin{gathered} \{\begin{array}{l}x=u+s_v\ast v\\ y=v\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}1& 0\\ s_v& 1\end{array}\right] \end{gathered}$$

垂直剪切

$$\begin{gathered} \{\begin{array}{l}x=u\\ y=s_u\ast u+v\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}1& s_u\\ 0& 1\end{array}\right] \end{gathered}$$

镜像

$$\begin{gathered} \{\begin{array}{l}x=-u\\ y=v\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right] \end{gathered}$$

平移

$$\begin{gathered} \{\begin{array}{l}x=u+t_x\\ y=v+t_y\end{array} \\ \left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]+\left[\begin{array}{cc}{t}_x& t_y\end{array}\right] \end{gathered}$$

我们可以乘以多个线性变换矩阵来实现各种基础线性变换类型的叠加，比如如果我们想同时实现旋转和尺度变换，就可以使用下式。实际计算时一般先把后面两个2×2矩阵乘算出来，然后再与左侧的坐标值做计算，这样我们只需要一次坐标变换就可以实现两种变换效果。（公式中我们只写了一个点的坐标，但实际操作中，我们需要对图像中$height\times width$个点做坐标变换，所以上述方式可以节省很多计算量）
$$\left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}\alpha & 0\\ 0& \beta \end{array}\right]\left[\begin{array}{cc}cos(\theta )& sin(\theta )\\ -sin(\theta )& cos(\theta )\end{array}\right]$$

仿射变换通式

我们可以将各种基础变换融合起来用一个式子表达，也就是仿射变换的通式，其线性方程组形式为：
$$\{\begin{array}{l}x=t_{11}\ast u+t_{21}\ast v+t_{31}\\ y=t_{12}\ast u+t_{22}\ast v+t_{32}\end{array}$$
可以有三种方法将其写为矩阵形式：

• 将线性变换和平移分开来写，这种写法比较符合仿射变换的基础概念，但是通常我们在实际应用时不使用这种方法，因为根据这种写法我们需要做两次矩阵运算，一次乘法和一次加法，也就是说它不够简单。
  $$\left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}{t}_{11}& t_{12}\\ t_{21}& t_{22}\end{array}\right]+\left[\begin{array}{cc}{t}_{31}& t_{32}\end{array}\right]$$
• 常用的是直接根据线性方程组的形式获得一种融合式的写法。此时我们使用一个矩阵乘法就完成了线性变换和平移两个运算。
  $$\left[\begin{array}{cc}x& y\end{array}\right]=\left[\begin{array}{ccc}u& v& 1\end{array}\right]\left[\begin{array}{cc}{t}_{11}& t_{12}\\ t_{21}& t_{22}\\ t_{31}& t_{32}\end{array}\right]$$
  当我们根据现行方程组的形式将其改写为矩阵形式时，我们很自然地就会把等号右侧两个矩阵的第一个写成$\left[\begin{array}{ccc}u& v& 1\end{array}\right]$，这种形式的坐标称为齐次坐标，它是大有来头的。在2D的图像仿射变换中，我们对其最直观的感受是使用这种坐标形式可以将仿射变换的通式用一个矩阵乘法完成，也就是将线性变换和平移合并了起来，这已经是很大的便利了，但是给人的感觉还是不够冲击，可能因为在仿射变换中它并非不可或缺。当我们进一步推导透视变换（Projective Transformation）时，就会发现必需要引入齐次坐标才行。另外齐次坐标也是图形学的基石之一。不过我们这里不讨论透视变换和图形学。
• 第三种写法是第二种写法的延伸，简单来说就是将等号左侧也写为齐次坐标形式，这种写法可能是教科书中最常见的形式了，但是我个人比较推荐第二种，因为在仿射变换中，变换矩阵的第三列固定为${\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T$，没有什么实际用处，并且多少浪费了点计算量。
  $$\left[\begin{array}{ccc}x& y& 1\end{array}\right]=\left[\begin{array}{ccc}u& v& 1\end{array}\right]\left[\begin{array}{ccc}{t}_{11}& t_{12}& 0\\ t_{21}& t_{22}& 0\\ t_{31}& t_{32}& 1\end{array}\right]$$

例子：人脸位置对齐

接下来我们以人脸识别中的人脸位置对齐为例，看看如何求解仿射变换矩阵以及如何应用这个知识点。这一部分对知识背景做简要描述。

人脸识别是目前（2020年）为止深度学习最成功的应用之一。跟很多其他的识别任务比起来，人脸识别任务本身其实并不算难。比如给定相同数量的训练集（1000万）和测试集（100万），要求大家在测试集上的准确率达到95%才算合格，那么相比较行人重识别，或者ImageNet那种通类识别任务，人脸识别可能是最容易合格的，这就是为什么说任务本身不算难。那么其难点在哪呢，难点在于大家对人脸识别的精度要求非常高，绝不是95%的准确率就可以让人满意的，试想如果每100个人脸识别支付的订单中都有5个可能出错，那将是一场大灾难，比如某个跟你长得有点相似的路人甲在超市买了点东西，刷了一下脸，然后被识别成了你，并且从你的账户里扣了钱，那就真的悲剧了。人脸识别的准确率要达到99.999…%，小数点后要n个9才行（2019年看到过一个关于人脸识别准确率的要求：要求在某个特定的拒识率下，误识率不高于千万分之一。这个意思也就是说对于少部分非常难以识别的人脸图片，可以拒绝识别，但是一旦同意识别，就必需得识别对）。

那么在这么高的准确率的要求下，人脸识别在数据预处理上与其他的图像识别任务就有着一些区别。比如人脸识别要求进行人脸位置对齐，这样可以减少由于位置偏差而带来的误识别率。

具体而言就是在做完人脸检测后，还要进行人脸特征点回归，一般五点就够用：两个眼睛中心，鼻尖，两个嘴角。回归出五个特征点后，将其配准到我们预先定义好的五点模板上。回归出的五个特征点坐标以$(u,v)$表示，五点模板坐标以$(x,y)$表示，这些信息都是已知的，所以在仿射变换通式中，只有变换矩阵是未知的，此时我们可以通过最小二乘法反推出变换矩阵。然后我们再正向运用变换矩阵，将原图中的人脸通过坐标变换变到五点模板上，完成位置对齐的任务。

这里多讲一点。并不是所有的图像识别任务都可以做位置对齐，人脸识别之所以可以做位置对齐，是因为人脸图片大体上存在一种固定的模式：眉毛、眼睛、鼻子和嘴巴的相对位置关系不会发生非常显著的变化（或者可以反过来理解，眼睛不会跑到嘴巴和鼻子中间去）。

一个需要小心的坑：图像索引与坐标的关系

在推导公式的时候，我们通常不需要引入图像索引这个概念，也就是我们常用的(i, j)，但是写代码的时候就必然得用到了，这个时候需要特别小心图像索引与坐标的对应关系。图像索引中我们常常使用i来表示行，用j来表示列，此时需注意以下对应关系：

i --> 行 --> y
j --> 列 --> x

所以像素索引为(i, j)的点的坐标为(y, x)。这个对应关系如果搞错了的话，可能会发生莫名其妙的错误，比如明明我们想要顺时针旋转图片，结果却是逆时针旋转。

数据准备

我们利用opencv写一个非常简单的交互程序来获取图像中某些像素的坐标：

import cv2

WIN_NAME = 'pick_points'
IMAGE_FILE = 'test.png'

def pick_points(event, x, y, flags, param):
    if event == cv2.EVENT_LBUTTONDOWN:
        print('x = %d, y = %d' % (x, y))


cv2.namedWindow(WIN_NAME, 0)
cv2.setMouseCallback(WIN_NAME, pick_points)
image = cv2.imread(IMAGE_FILE, cv2.IMREAD_UNCHANGED)
while True:
    cv2.imshow(WIN_NAME, image)
    key = cv2.waitKey(30)
    if key == 27:  # ESC
        break
cv2.destroyAllWindows()

你可能需要修改程序中图片文件的路径，注意路径中不要包含中文字符。运行程序后，左键点击图片可以获得当前鼠标位置的坐标，通过print函数将坐标打印出来。

本来这里应该用一张人脸图片的，不过我怕随意使用别人的图片侵犯别人肖像权，但是又不想用自己的大头照，所以思来想去搞了一张大猩猩的脸。图片在我本地是550 x 333（宽 x 高）的分辨率，不知道上传后会不会发生变化，如果变了的话，可以使用上述代码自己再拾取一下坐标。

拾取坐标的顺序为：左眼，右眼，鼻尖，左嘴角，右嘴角。注意这里的左右以我们看到的图片为基准，不用去想着图片里的对象的真实左右是什么，因为有些相机在拍摄的时候做过镜像，所以真实的左右很难讲明白。我拾取到的坐标为（还记得我们之前的符号约定，原图用u, v表示坐标）：

u1 = 288, v1 = 115
u2 = 334, v2 = 119
u3 = 278, v3 = 160
u4 = 250, v4 = 193
u5 = 318, v5 = 192

做人脸对齐时，必须有一个已知的人脸五点模板。然而现在我们正在做大猩猩脸对齐，假设有一个大猩猩脸的五点模板（模板尺寸是100 x 100）为：

x1 = 31.3, y1 = 45.2
x2 = 67.8, y2 = 45.2
x3 = 49.5, y3 = 68.7
x4 = 35.2, y4 = 84.8
x5 = 62.6, y5 = 84.8

这个模板随便搞的，毕竟是大猩猩脸，没那么数据来算模板。

求解仿射变换矩阵

前向映射

将五组对应的点代入仿射变换通式可得：
$$\{\begin{array}{l}{x}_1=t_{11}\ast u_1+t_{21}\ast v_1+t_{31}\ast 1\\ y_1=t_{12}\ast u_1+t_{22}\ast v_1+t_{32}\ast 1\\ x_2=t_{11}\ast u_2+t_{21}\ast v_2+t_{31}\ast 1\\ y_2=t_{12}\ast u_2+t_{22}\ast v_2+t_{32}\ast 1\\ ......\\ x_5=t_{11}\ast u_5+t_{21}\ast v_5+t_{31}\ast 1\\ y_5=t_{12}\ast u_5+t_{22}\ast v_5+t_{32}\ast 1\end{array}$$
上式中，$u,v,x,y$都是已知量，待求解的是$t_{ij}$，另外注意到$t_{11},t_{21},t_{31}$只与$x$相关，$t_{12},t_{22},t_{32}$只与$y$相关。所以上述线性方程组可以写成两组：
$$\{\begin{array}{l}{x}_1=t_{11}\ast u_1+t_{21}\ast v_1+t_{31}\ast 1\\ x_2=t_{11}\ast u_2+t_{21}\ast v_2+t_{31}\ast 1\\ ......\\ x_5=t_{11}\ast u_5+t_{21}\ast v_5+t_{31}\ast 1\end{array}$$
$$\{\begin{array}{l}{y}_1=t_{12}\ast u_1+t_{22}\ast v_1+t_{32}\ast 1\\ y_2=t_{12}\ast u_2+t_{22}\ast v_2+t_{32}\ast 1\\ ......\\ y_5=t_{12}\ast u_5+t_{22}\ast v_5+t_{32}\ast 1\end{array}$$
记：
$$\begin{gathered} T_1={\left[\begin{array}{c}{t}_{11},t_{21},t_{31}\end{array}\right]}^T \\ {T}_2={\left[\begin{array}{c}{t}_{12},t_{22},t_{32}\end{array}\right]}^T \end{gathered}$$
$$\begin{gathered} X={\left[\begin{array}{c}{x}_1,x_2,...,x_5\end{array}\right]}^T \\ Y={\left[\begin{array}{c}{y}_1,y_2,...,y_5\end{array}\right]}^T \end{gathered}$$
$$A=\left[\begin{array}{ccc}{u}_1& v_1& 1\\ u_2& v_2& 1\\ ......& & \\ u_5& v_5& 1\end{array}\right]$$
那么上面两个线性方程组可写为：$A{T}_1=X$，$A{T}_2=Y$。现在$A,X,Y$是已知的，分别求出$T_1,T_2$即可。这两个线性方程组是超定方程组，即方程数目多于未知量数目，所以方程组没有精确解，只有最小二乘解。
下面以$A{T}_1=X$为例推导最小二乘解：
$$\begin{gathered} A{T}_1=X \\ {A}^{T}A{T}_1=A^{T}X \\ (A^{T}A{)}^{-1}{A}^{T}A{T}_1=(A^{T}A{)}^{-1}{A}^{T}X \\ {T}_1=(A^{T}A{)}^{-1}{A}^{T}X \end{gathered}$$
简单解释上面的推导过程：

• 两边同乘以$A^T$可以将方程组变为正定，即现在的系数矩阵$A^{T}A$是一个方阵
• 两边再同乘以$A^{T}A$的逆，此时左侧的系数矩阵$(A^{T}A{)}^{-1}{A}^{T}A$变为单位阵，等式右侧就是$T_1$的解

同理可得$T_2=(A^{T}A{)}^{-1}{A}^{T}Y$

以上求解过程从原始图像的坐标$(u,v)$出发，经过矩阵变换，得到变换后的坐标$(x,y)$，该过程称为前向映射，前向映射中$(u,v)$是整数，$(x,y)$一般是浮点数。前向映射是比较容易理解的求解方法，但是实际上在真正的计算过程中不使用前向映射，而是使用后向映射。原因是坐标变换后的值通常是浮点数，而图像只记录整数坐标位置上的像素值，所以坐标变换后一般要接上一个插值的流程，而前向映射的计算结果对插值不友好，后向映射则顺滑得多。

后向映射

后向映射从变换后图像的坐标$(x,y)$出发，经矩阵变换，得到其在原始图像中对应的坐标$(u,v)$，后向映射中$(x,y)$是整数，$(u,v)$一般是浮点数。

后向映射与前向映射的推导基本雷同，只需将$(x,y)$与$(u,v)$的位置调换一下就可以，下面罗列一下推导流程。
将五组对应的点代入通式：
$$\{\begin{array}{l}u=t_{11}\ast x+t_{21}\ast y+t_{31}\\ v=t_{12}\ast x+t_{22}\ast y+t_{32}\end{array}$$
可得：
$$\{\begin{array}{l}{u}_1=t_{11}\ast x_1+t_{21}\ast y_1+t_{31}\ast 1\\ u_2=t_{11}\ast x_2+t_{21}\ast y_2+t_{31}\ast 1\\ ......\\ u_5=t_{11}\ast x_5+t_{21}\ast y_5+t_{31}\ast 1\end{array}$$
$$\{\begin{array}{l}{v}_1=t_{12}\ast x_1+t_{22}\ast y_1+t_{32}\ast 1\\ v_2=t_{12}\ast x_2+t_{22}\ast y_2+t_{32}\ast 1\\ ......\\ v_5=t_{12}\ast x_5+t_{22}\ast y_5+t_{32}\ast 1\end{array}$$
记：
$$\begin{gathered} T_1={\left[\begin{array}{c}{t}_{11},t_{21},t_{31}\end{array}\right]}^T \\ {T}_2={\left[\begin{array}{c}{t}_{12},t_{22},t_{32}\end{array}\right]}^T \end{gathered}$$
$$\begin{gathered} U={\left[\begin{array}{c}{u}_1,u_2,...,u_5\end{array}\right]}^T \\ V={\left[\begin{array}{c}{v}_1,v_2,...,v_5\end{array}\right]}^T \end{gathered}$$
$$A=\left[\begin{array}{ccc}{x}_1& y_1& 1\\ x_2& y_2& 1\\ ......& & \\ x_5& y_5& 1\end{array}\right]$$
那么方程组的解为：
$$\begin{gathered} T_1=(A^{T}A{)}^{-1}{A}^{T}U \\ {T}_2=(A^{T}A{)}^{-1}{A}^{T}V \end{gathered}$$
记：$\bar{A}=(A^{T}A{)}^{-1}{A}^T$，先把$\bar{A}$计算出来可以节省一些计算量，此时有：
$$\begin{gathered} T_1=\bar{A}U \\ {T}_2=\bar{A}V \end{gathered}$$

后向映射代码实现

请注意，本文档的主要内容是介绍仿射变换，但是如果想要看到仿射变换后的图片，坐标变换后还需要接一个插值的步骤，工程中一般用双线性插值会多一些，不过插值不是这里的重点，所以图个简单省事使用了最邻近插值方法，插值函数为nearest_sampler。
如果对插值有额外兴趣的话，可以参考下文：
图像插值：最邻近（nearest）与双线性（bilinear）

# -*- coding: utf-8 -*-

import cv2
import numpy as np


def nearest_sampler(image, coords):
    """
    Nearest Neighbour sampler
    
    Parameters
    ----------
    image: ndarray
        source image, whose shape is [height, width, channels]
    coords: ndarray
        coordinates to be interpolated, the length of last axis should be 2,
        meaning 2D coordinate
    
    Returns
    -------
    output: ndarray
        the interpolated image, same shape as coords except the last axis
    """
    height, width, channels = image.shape[0:3]
    output_shape = list(coords.shape)
    coords = np.reshape(coords, (-1, output_shape[-1]))
    output_shape[-1] = channels
    coords = np.round(coords).astype(np.int32)
    idx = (coords[:, 0] >= 0) & (coords[:, 0] < width) & \
          (coords[:, 1] >= 0) & (coords[:, 1] < height)
    output = np.zeros((coords.shape[0], channels), dtype=np.uint8)
    output[idx] = image[coords[idx, 1], coords[idx, 0], :]
    output = np.reshape(output, output_shape)
    return output


if __name__ == '__main__':
    # feature points
    uv = np.array([[288, 115],
                   [334, 119],
                   [278, 160],
                   [250, 193],
                   [318, 192]])
    # model points
    xy = np.array([[31.3, 45.2],
                   [67.8, 45.2],
                   [49.5, 68.7],
                   [35.2, 84.8],
                   [62.6, 84.8]])

    # compute transform matrix
    U = np.reshape(uv[:, 0], [uv.shape[0], 1])
    V = np.reshape(uv[:, 1], [uv.shape[0], 1])
    A = np.zeros((xy.shape[0], 3))
    A[:, 0:2] = xy
    A[:, 2] = 1
    A_bar = np.linalg.inv(A.T @ A) @ A.T
    T1 = A_bar @ U
    T2 = A_bar @ V
    T = np.concatenate((T1, T2), axis=1)

    # coordinates mapping
    height = 100
    width = 100
    coords = np.meshgrid(np.arange(0, width), np.arange(0, height))
    # shape = [height, width, 2] after transpose
    coords = np.array(coords).transpose([1, 2, 0])
    ones = np.ones([height, width, 1])
    # homogeneous  coordinates
    coords = np.concatenate((coords, ones), axis=2)
    # transformed coordinates
    coords = coords @ T

    # sampler
    image = cv2.imread('2.jpg')
    img = nearest_sampler(image, coords)

    cv2.imshow('affine', img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

大猩猩脸对齐的结果为：

这个脸看起来有一些变形，追溯回基础变换类型的话有两个因素导致了变形：

• 宽高方向的缩放尺度不一样;
• 剪切变换。

在实际的人脸位置对齐中，一般不希望引入这种变形，我们通常希望宽高方向的缩放尺度是一样的，并且不希望有剪切变换存在。所以就需要对仿射变换加一些限制，下面一部分就介绍这种特殊的仿射变换矩阵的求解方法。

一种特殊的仿射变换及变换矩阵求解方法

这一节介绍的放射变换为了解决之前的问题，要求：

• 宽高方向的缩放尺度一样
• 不引入剪切变换

请注意，人脸位置对齐时镜像变换也不会被引入，所以事实上就剩下了三种基础变换类型：

• 旋转
• 尺度（宽高方向缩放因子相同）
• 平移

所以变换过程可以表示为（直接推导后向映射过程）：
$$\left[\begin{array}{cc}u& c\end{array}\right]=\left[\begin{array}{cc}x& y\end{array}\right]\left[\begin{array}{cc}s& 0\\ 0& s\end{array}\right]\left[\begin{array}{cc}cos(\theta )& sin(\theta )\\ -sin(\theta )& cos(\theta )\end{array}\right]+\left[\begin{array}{cc}{t}_u& t_v\end{array}\right]$$
将上述式子合并，变为：
$$\left[\begin{array}{cc}u& v\end{array}\right]=\left[\begin{array}{ccc}x& y& 1\end{array}\right]\left[\begin{array}{cc}s\ast cos(\theta )& s\ast sin(\theta )\\ -s\ast sin(\theta )& s\ast cos(\theta )\\ t_u& t_v\end{array}\right]$$
现在记：
$$\begin{gathered} a=s\ast cos(\theta ) \\ b=s\ast sin(\theta ) \\ c=t_u \\ d=t_v \end{gathered}$$
那么仿射变换可以写为：
$$\left[\begin{array}{cc}u& v\end{array}\right]=\left[\begin{array}{ccc}x& y& 1\end{array}\right]\left[\begin{array}{cc}a& b\\ -b& a\\ c& d\end{array}\right]$$
接下来我们求解这个变换矩阵。类似之前的过程，我们将五组点代入上述线性方程组可得（此处不再推导前向映射过程，直接写后向映射过程）：
$$\{\begin{array}{l}{u}_1=a\ast x_1-b\ast y_1+c\ast 1+d\ast 0\\ u_2=a\ast x_2-b\ast y_2+c\ast 1+d\ast 0\\ ...\\ u_5=a\ast x_5-b\ast y_5+c\ast 1+d\ast 0\\ v_1=a\ast y_1+b\ast x_1+c\ast 0+d\ast 1\\ v_2=a\ast y_2+b\ast x_2+c\ast 0+d\ast 1\\ ...\\ v_5=a\ast y_5+b\ast x_5+c\ast 0+d\ast 1\end{array}$$
因为变换矩阵的第一列和第二列存在共同变量，也就是变量耦合，所以没法再像之前一样分开两个方程组分别求解第一列和第二列。这里通过一次最小二乘求解得到变换矩阵中的所有4个参数。
记：
$$A=\left[\begin{array}{cccc}{x}_1& -y_1& 1& 0\\ x_2& -y_2& 1& 0\\ ...& & & \\ x_5& -y_5& 1& 0\\ y_1& x_1& 0& 1\\ y_2& x_2& 0& 1\\ ...& & & \\ y_5& x_5& 0& 1\end{array}\right]$$
$$B={\left[\begin{array}{cccccccc}{u}_1& u_2& ...& u_5& v_1& v_2& ...& v_5\end{array}\right]}^T$$
$$\bar{T}={\left[\begin{array}{cccc}a& b& c& d\end{array}\right]}^T$$
那么有：
$$\bar{T}=(A^{T}A{)}^{-1}{A}^{T}B$$
有了$\bar{T}$以后我们重新排列一下元素的形状和顺序就可以得到真正的变换矩阵。下面是代码示例，除了compute transform matrix部分之外，其它与上面的代码一样，还是用上面大猩猩脸的数据。

# -*- coding: utf-8 -*-

import cv2
import numpy as np


def nearest_sampler(image, coords):
    """
    Nearest Neighbour sampler
    
    Parameters
    ----------
    image: ndarray
        source image, whose shape is [height, width, channels]
    coords: ndarray
        coordinates to be interpolated, the length of last axis should be 2,
        meaning 2D coordinate
    
    Returns
    -------
    output: ndarray
        the interpolated image, same shape as coords except the last axis
    """
    height, width, channels = image.shape[0:3]
    output_shape = list(coords.shape)
    coords = np.reshape(coords, (-1, output_shape[-1]))
    output_shape[-1] = channels
    coords = np.round(coords).astype(np.int32)
    idx = (coords[:, 0] >= 0) & (coords[:, 0] < width) & \
          (coords[:, 1] >= 0) & (coords[:, 1] < height)
    output = np.zeros((coords.shape[0], channels), dtype=np.uint8)
    output[idx] = image[coords[idx, 1], coords[idx, 0], :]
    output = np.reshape(output, output_shape)
    return output


if __name__ == '__main__':
    # feature points
    uv = np.array([[288, 115],
                   [334, 119],
                   [278, 160],
                   [250, 193],
                   [318, 192]])
    # model points
    xy = np.array([[31.3, 45.2],
                   [67.8, 45.2],
                   [49.5, 68.7],
                   [35.2, 84.8],
                   [62.6, 84.8]])

    # compute transform matrix
    num = xy.shape[0]
    A = np.zeros((2 * num, 4))
    A[0:num, 0] = xy[:, 0]
    A[0:num, 1] = -xy[:, 1]
    A[0:num, 2] = 1
    A[num:2 * num, 0] = xy[:, 1]
    A[num:2 * num, 1] = xy[:, 0]
    A[num:2 * num, 3] = 1
    B = np.reshape(uv.T, [2 * num, 1])
    T_bar = np.linalg.inv(A.T @ A) @ A.T @ B
    T = np.zeros((3, 2))
    T[0, 0] = T_bar[0]
    T[0, 1] = T_bar[1]
    T[1, 0] = -T_bar[1]
    T[1, 1] = T_bar[0]
    T[2, 0] = T_bar[2]
    T[2, 1] = T_bar[3]

    # coordinates mapping
    height = 100
    width = 100
    coords = np.meshgrid(np.arange(0, width), np.arange(0, height))
    # shape = [height, width, 2] after transpose
    coords = np.array(coords).transpose([1, 2, 0])
    ones = np.ones([height, width, 1])
    # homogeneous  coordinates
    coords = np.concatenate((coords, ones), axis=2)
    # transformed coordinates
    coords = coords @ T

    # sampler
    image = cv2.imread('2.jpg')
    img = nearest_sampler(image, coords)

    cv2.imshow('affine', img)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

得到的结果是：

再把之前得到的带有变形的结果放在这里，各位可以对比一下：

下面是我直接用软件旋转原图并裁减后得到的结果，因为是随便手动裁减的，所以位置跟程序算出来的肯定不一样，各位主要看一下这两种仿射变换并跟原图比较，观察一下第一种更加通用的仿射变换带来的变形效果：

从图片拍摄的角度来看，导致这种变形效果的原因是脸比较侧，假如脸比较正的话，即使使用第一种仿射变换也不会导致太大的变形，甚至歪个脖子也不会。

又一个坑：前向映射与后向映射的变换矩阵不互逆

有些同学写程序时图个概念上的方便，因为教科书上通常教的都是前向映射，所以写程序时也按照前向映射进行求解，求出前向映射的变换矩阵后，通过求逆来得到后向映射的变换矩阵，请注意这种做法是错误的，两者一般不互逆。为什么说一般呢，有一种情况互逆关系可以成立：求解的线性方程组为正定，也就是只有三组不线性相关的特征点，这种情况在实践中不是太多见。

下面我们用程序进行验证。请注意，上面我们求解坐标变换时使用的是3 x 2的矩阵，如果要验证互逆性的话，需要拼个${\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T$组成一个3 x 3的矩阵。

# -*- coding: utf-8 -*-

import numpy as np

# feature points
uv = np.array([[288, 115],
               [334, 119],
               [278, 160],
               [250, 193],
               [318, 192]])
# model points
xy = np.array([[31.3, 45.2],
               [67.8, 45.2],
               [49.5, 68.7],
               [35.2, 84.8],
               [62.6, 84.8]])

# compute backward transform matrix
U = np.reshape(uv[:, 0], [uv.shape[0], 1])
V = np.reshape(uv[:, 1], [uv.shape[0], 1])
A = np.zeros((xy.shape[0], 3))
A[:, 0:2] = xy
A[:, 2] = 1
A_bar = np.linalg.inv(A.T @ A) @ A.T
T1 = A_bar @ U
T2 = A_bar @ V
T = np.concatenate((T1, T2), axis=1)
T_backward = np.concatenate((T, np.array([[0], [0], [1]])), axis=1)

# compute forward transform matrix
X = np.reshape(xy[:, 0], [xy.shape[0], 1])
Y = np.reshape(xy[:, 1], [xy.shape[0], 1])
A = np.zeros((xy.shape[0], 3))
A[:, 0:2] = uv
A[:, 2] = 1
A_bar = np.linalg.inv(A.T @ A) @ A.T
T1 = A_bar @ X
T2 = A_bar @ Y
T_forward = np.concatenate((T1, T2, np.array([[0], [0], [1]])), axis=1)

# verification
T_forward_inv = np.linalg.inv(T_forward)
T_mul = T_forward @ T_backward

上面求到的几个矩阵为：
$$\begin{gathered} T_{forward}=\left[\begin{array}{ccc}0.488518& -0.0123897& 0\\ 0.177793& 0.51996& 0\\ -121.849& -11.6322& 1\end{array}\right] \\ {T}_{backward}=\left[\begin{array}{ccc}1.69603& 0.0564699& 0\\ -0.68687& 1.90407& 0\\ 255.175& 27.8433& 1\end{array}\right] \\ {T}_{forward\_inv}=\left[\begin{array}{ccc}2.02941& 0.048357& 0\\ -0.69393& 1.90669& 0\\ 239.21& 28.0712& 1\end{array}\right] \\ {T}_{mul}=\left[\begin{array}{ccc}0.837049& 0.00399568& 0\\ -0.0556026& 1.00008& 0\\ 56.5053& -1.18607& 1\end{array}\right] \end{gathered}$$
可以看到，$T_{forward\_inv}$和$T_{backward}$虽然有点相近，但是毕竟还是差了点意思，另外$T_{forward}$和$T_{backward}$的乘积$T_{mul}$也不是单位阵，所以两者不互逆。