好的，我来这里是为了寻找答案，但没有找到简单明了的东西，所以我继续做了愚蠢但有效（且相对简单）的事情：蒙特卡洛优化。
简单地说，算法如下：随机扰动你的投影矩阵，直到它将你已知的 3D 坐标投影到你已知的 2D 坐标。
这是托马斯坦克引擎的静态照片：
假设我们使用 GIMP 在地平面上找到我们认为是正方形的二维坐标（是否真的是正方形取决于您对深度的判断）：
我在 2D 图像中得到四个点：(318, 247)、(326, 312)、(418, 241) 和 (452, 303)。
按照惯例，我们说这些点应该对应于 3D 点：(0, 0, 0)、(0, 0, 1)、(1, 0, 0) 和 (1, 0, 1)。也就是说，y=0 平面上的一个单位正方形。
将这些 3D 坐标中的每一个投影到 2D 中是通过将 4D 向量 [x, y, z, 1] 与 4x4 投影矩阵相乘，然后将 x 和 y 分量除以 z 来实际获得透视校正来完成的。这或多或少是gluProject() 所做的，除了gluProject() 还考虑了当前视口并考虑了单独的模型视图矩阵（我们可以假设模型视图矩阵是单位矩阵）。查看 gluProject() 文档非常方便，因为我实际上想要一个适用于 OpenGL 的解决方案，但要注意文档中缺少公式中的 z 除法。
请记住，算法是从一些投影矩阵开始并随机扰动它，直到它给出我们想要的投影。所以我们要做的是投影四个 3D 点中的每一个，看看我们离我们想要的 2D 点有多近。如果我们的随机扰动导致投影的 2D 点更接近我们在上面标记的点，那么我们保留该矩阵作为对我们最初（或之前）猜测的改进。
让我们定义我们的观点：
# Known 2D coordinates of our rectangle
i0 = Point2(318, 247)
i1 = Point2(326, 312)
i2 = Point2(418, 241)
i3 = Point2(452, 303)
# 3D coordinates corresponding to i0, i1, i2, i3
r0 = Point3(0, 0, 0)
r1 = Point3(0, 0, 1)
r2 = Point3(1, 0, 0)
r3 = Point3(1, 0, 1)
我们需要从一些矩阵开始，单位矩阵似乎是一个自然的选择：
mat = [
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
]
我们需要实际实现投影（基本上是矩阵乘法）：
def project(p, mat):
x = mat[0][0] * p.x + mat[0][1] * p.y + mat[0][2] * p.z + mat[0][3] * 1
y = mat[1][0] * p.x + mat[1][1] * p.y + mat[1][2] * p.z + mat[1][3] * 1
w = mat[3][0] * p.x + mat[3][1] * p.y + mat[3][2] * p.z + mat[3][3] * 1
return Point(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.)
这基本上是gluProject() 所做的，720 和 576 分别是图像的宽度和高度（即视口），我们从 576 中减去以计算我们从顶部开始计算 y 坐标的事实，而OpenGL 通常从底部开始计算它们。你会注意到我们没有计算 z，那是因为我们在这里并不需要它（尽管它可以很方便地确保它在 OpenGL 用于深度缓冲区的范围内）。
现在我们需要一个函数来评估我们与正确解决方案的距离。该函数返回的值是我们将用来检查一个矩阵是否优于另一个矩阵的值。我选择按距离平方和，即：
# The squared distance between two points a and b
def norm2(a, b):
dx = b.x - a.x
dy = b.y - a.y
return dx * dx + dy * dy
def evaluate(mat):
c0 = project(r0, mat)
c1 = project(r1, mat)
c2 = project(r2, mat)
c3 = project(r3, mat)
return norm2(i0, c0) + norm2(i1, c1) + norm2(i2, c2) + norm2(i3, c3)
为了扰动矩阵，我们只需在某个范围内随机选取一个元素进行扰动：
def perturb(amount):
from copy import deepcopy
from random import randrange, uniform
mat2 = deepcopy(mat)
mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)
（值得注意的是，我们的project() 函数实际上根本不使用mat[2]，因为我们不计算z，并且因为我们所有的y 坐标都是0，所以mat[*][1] 的值也无关紧要。我们可以利用这个事实，并且永远不要试图扰乱这些值，这会带来一点加速，但这留作练习......）
为方便起见，让我们添加一个函数，通过反复调用perturb() 对我们迄今为止发现的最佳矩阵进行大量近似：
def approximate(mat, amount, n=100000):
est = evaluate(mat)
for i in xrange(n):
mat2 = perturb(mat, amount)
est2 = evaluate(mat2)
if est2 < est:
mat = mat2
est = est2
return mat, est
现在剩下要做的就是运行它...：
for i in xrange(100):
mat = approximate(mat, 1)
mat = approximate(mat, .1)
我发现这已经给出了一个相当准确的答案。跑了一会儿，发现的矩阵是：
[
[1.0836000765696232, 0, 0.16272110011060575, -0.44811064935115597],
[0.09339193527789781, 1, -0.7990570384334473, 0.539087345090207 ],
[0, 0, 1, 0 ],
[0.06700844759602216, 0, -0.8333379578853196, 3.875290562060915 ],
]
错误约为2.6e-5。 （请注意，我们所说的未在计算中使用的元素实际上并没有从我们的初始矩阵中改变；那是因为改变这些条目不会改变评估的结果，所以改变永远不会继续。）
我们可以使用glLoadMatrix() 将矩阵传递到OpenGL（但记得先转置它，并记得用单位矩阵加载你的模型视图矩阵）：
def transpose(m):
return [
[m[0][0], m[1][0], m[2][0], m[3][0]],
[m[0][1], m[1][1], m[2][1], m[3][1]],
[m[0][2], m[1][2], m[2][2], m[3][2]],
[m[0][3], m[1][3], m[2][3], m[3][3]],
]
glLoadMatrixf(transpose(mat))
现在我们可以例如沿 z 轴平移以获得沿轨道的不同位置：
glTranslate(0, 0, frame)
frame = frame + 1
glBegin(GL_QUADS)
glVertex3f(0, 0, 0)
glVertex3f(0, 0, 1)
glVertex3f(1, 0, 1)
glVertex3f(1, 0, 0)
glEnd()
从数学的角度来看，这肯定不是很优雅；你没有得到一个封闭形式的方程，你可以将你的数字插入并得到一个直接（和准确）的答案。但是，它确实允许您添加额外的约束，而不必担心使方程式复杂化；例如，如果我们也想合并高度，我们可以使用房子的那个角落并说（在我们的评估函数中）从地面到屋顶的距离应该是某某，然后再次运行该算法。所以，是的，这是一种蛮力，但很有效，而且效果很好。