还原被摄像机透视的纹理
有人问如何还原被透视纹理？给你一张照片，还原照片上四个点所组成的平面的纹理该怎么做？我们可以从数学上推导一下，为了和三维图形的透视纹理映射对照，我们称照片上四个点在照片上的位置为“屏幕坐标”，那么可以发现：
空间中，三维坐标（x,y,z）和纹理坐标（u, v）承线性关系。根据该问题描述，可以理解为已知四个点的屏幕投影坐标（xi,yi），和对应纹理坐标（u，v），求整个纹理坐标系到屏幕坐标系的反向映射过程，即根据（u,v）求解（xi，yi）。
1. 按照纹理隐射的原理，同平面纹理坐标与空间坐标存在线性关系，设 a1-a12为常数，有：
a1 * x + a2 * y + a3 * z + a4 = u ... 线性关系
a5 * x + a6 * y + a7 * z + a8 = v ... 线性关系
a9 * x + a10 * y + a11 * z + a12 = 0 ... 平面方程
2. 求解上面的方程组，可以得到类似下面的关系，其中b1-b9为常数：
x = b1 * u + b2 * v + b3
y = b4 * u + b5 * v + b6
z = b7 * u + b8 * v + b9
常数 b1-b9如果展开，就是9个关于a1-a12的等式，太长了，这里不展开，有兴趣可以自己求解。
3. 屏幕上投影坐标一般是：
x
xi = D1 * --- + C1
z
x
yi = D2 * --- + C2
z
因为同样一个透视投影矩阵下，能隐射成屏幕上同样形状纹理的平面，在空间中存在无穷多个，而且还存在不同透视投影矩阵下，同样屏幕投影的平面存在更多无穷多个。这里我们不用去求解每个平面，直接设置 D1 = D2 = 1 且 C1 = C2 = 0 有：
x
xi = ---
z
x
yi = ---
z
4. 展开等式：
b1 * u + b2 * v + b3
xi = ----------------------
b7 * u + b8 * v + b9
b4 * u + b5 * v + b6
yi = ----------------------
b7 * u + b8 * v + b9
改变一下，用矩阵 M的各个元素 m00 - m22来代替 b1-b9。
m00 * u + m01 * v + m02
xi = -------------------------
m20 * u + m21 * v + m22
m10 * u + m11 * v + m12
yi = -------------------------
m20 * u + m21 * v + m22
5. 于是整个问题很简单了，就是求矩阵 M。我们知道存在矩阵 M，使得：
/ m00 m01 m02 \ / u \
| m10 m11 m12 | . | v | = [ uo, vo, wo ] = [ xi * wo, yi * wo, wo ]
\ m20 m21 m22 / \ 1 /
m00 * u + m01 * v + m02
xi = uo / wo = ------------------------ .......(4)
m20 * u + m21 * v + m22
m10 * u + m11 * v + m12
yi = uo / wo = ------------------------ .......(5)
m20 * u + m21 * v + m22
6. 再此基础上，我们已知屏幕上四个端点的位置（x0, y0）-> (x3, y3), 也知道他们对应的纹理坐标，（u0, v0) -> (u3, v3)。继续设 m22 = 1，可以列出一个方程：
/ u0 v0 1 0 0 0 -u0*x0 -v0*x0 \ /m00\ /x0\
| u1 v1 1 0 0 0 -u1*x1 -v1*x1 | |m01| |x1|
| u2 v2 1 0 0 0 -u2*x2 -v2*x2 | |m02| |x2|
| u3 v3 1 0 0 0 -u3*x3 -v3*x3 |.|m10|=|x3|
| 0 0 0 u0 v0 1 -u0*y0 -v0*y0 | |m11| |y0|
| 0 0 0 u1 v1 1 -u1*y1 -v1*y1 | |m12| |y1|
| 0 0 0 u2 v2 1 -u2*y2 -v2*y2 | |m20| |y2|
\ 0 0 0 u3 v3 1 -u3*y3 -v3*y3 / \m21/ \y3/
对于这个方程组，四个端点的 x, y是已知的，同时纹理坐标 u, v也是已知的。那么我们解个方程就能求出矩阵 M的各个元素。现在你已经有了 m00 - m22，有了式子4，和式子5：
/ m00 m01 m02 \
M = | m10 m11 m12 | ..... 已知
\ m20 m21 m22 /
m00 * u + m01 * v + m02
xi = ------------------------- ... (4)
m20 * u + m21 * v + m22
m10 * u + m11 * v + m12
yi = ------------------------- ... (5)
m20 * u + m21 * v + m22
7. 在已知 M矩阵的各个元素，已知纹理坐标（u, v），你想求对应的屏幕坐标（x, y），就可以用公式（4）和公式（5）来求解了。假设纹理为 256x256大小，最终求解原纹理图片的代码应该是：
for (int v = 0; v < 256; v++) {
for (int u = 0; u < 256; u++) {
int x = CalculateX(M, u, v); // 带入式子（4）
int y = CalculateY(M, u, v); // 带入式子（5）
texture[v * 256 + u] = GetPixel(x, y);
}
}
8. 优化一下，上面每个点的乘除法太多，我们需要找找关系，设：
A = m00 * u + m01 * v + m02
B = m10 * u + m11 * v + m12
C = m20 * u + m21 * v + m22
于是有：x = A / C, y = B / C，而随着 u坐标的均匀递增 A, B, C成线性变化，那么我们可以用差值来完成我们的优化：
for (int v = 0; v < 256; v++) {
float A1 = m00 * 0 + m01 * v + m02;
float A2 = m00 * 256 + m01 * v + m02;
float B1 = m10 * 0 + m11 * v + m12;
float B2 = m10 * 256 + m11 * v + m12;
float C1 = m20 * 0 + m21 * v + m22;
float C2 = m20 * 256 + m21 * v + m22;
float dA = (A2 - A1) * (1 / 256.0); // 计算步长
float dB = (B2 - B1) * (1 / 256.0);
float dC = (C2 - C1) * (1 / 256.0);
float A = A1 + 0.5; // 方便取整时四舍五入
float B = B1 + 0.5;
float C = C1 + 0.5;
uint32_t *pixel = &texture[v * 256];
for (int u = 0; u < 256; u++) {
float i = 1.0 / C; // 减少一次除法
int x = (int)(A * i);
int y = (int)(B * i);
pixel[0] = GetPixel(x, y);
pixel++;
A += dA;
B += dB;
C += dC;
}
}
整个程序基本上就是这个样子了，内循环已经足够紧凑。再次观察等式关系，可以发现其实：dA = m00, dB = m10, dC = m20。还可以继续进一步化解外层循环。是不是和透视纹理绘制方法很像？一个正过程，一个逆过程而已。
注意1：越界处理，简单判断下越界坐标。
注意2：如果还需要进一步提高性能，可以尝试定点数优化计算过程，尝试继续展开循环。