我在屏幕空间中有 4 个 2D 点，我需要将它们反向投影回 3D 空间。我知道这 4 个点中的每一个都是 3D 旋转刚性矩形的一个角，并且我知道矩形的大小。如何从中获得 3D 坐标？
我没有使用任何特定的 API，也没有现有的投影矩阵。我只是在寻找基本的数学来做到这一点。当然没有足够的数据将单个 2D 点转换为 3D 没有其他参考，但我想如果你有 4 个点，你就会知道它们在同一平面上彼此成直角，而且你知道它们之间的距离，你应该能够从那里弄清楚。不幸的是，我无法完全弄清楚如何。
这可能属于摄影测量的范畴，但谷歌搜索并没有让我找到任何有用的信息。
好吧，我来这里是为了寻找答案，但没有找到简单明了的东西，所以我继续做了愚蠢但有效（且相对简单）的事情：蒙特卡洛优化。
简单地说，算法如下：随机扰动您的投影矩阵，直到它将您已知的 3D 坐标投影到您已知的 2D 坐标。
这是托马斯坦克引擎的静态照片：
假设我们使用 GIMP 在地平面上找到我们认为是正方形的 2D 坐标（它是否真的是正方形取决于您对深度的判断）：
我在 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. （注意我们所说的未在计算中使用的元素实际上并没有从我们的初始矩阵中改变；那是因为改变这些条目不会改变评估的结果，所以改变永远不会继续。）
我们可以使用以下方法将矩阵传递给 OpenGL glLoadMatrix()（但请记住先转置它，并记住使用单位矩阵加载模型视图矩阵）：
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()
从数学的角度来看，这肯定不是很优雅。你没有得到一个封闭形式的方程，你可以将你的数字插入并得到一个直接（和准确）的答案。但是，它确实允许您添加额外的约束，而不必担心使方程式复杂化；例如，如果我们也想合并高度，我们可以使用房子的那个角落并说（在我们的评估函数中）从地面到屋顶的距离应该是某某，然后再次运行该算法。所以，是的，这是一种蛮力，但有效，而且运作良好。
这是基于标记的增强现实的经典问题。
您有一个方形标记（2D 条形码），并且您想在找到标记的四个边缘后找到它的姿势（相对于相机的平移和旋转）。 概览-图片
我不知道对该领域的最新贡献，但至少在一定程度上（2009 年）RPP 应该优于上面提到的 POSIT（并且确实是一个经典的方法）请查看链接，他们也提供来源。
http://www.emt.tugraz.at/publications/EMT_TR/TR-EMT-2005-01.pdf
http://www.emt.tugraz.at/system/files/rpp_MATLAB_ref_implementation.tar.gz
（PS - 我知道这是一个有点老的话题，但无论如何，这篇文章可能对某人有帮助）
D. DeMenthon设计了一种算法，在知道对象模型的情况下，根据 2D 图像中的特征点计算对象的姿势（其在空间中的位置和方向）——这正是您的问题：
我们描述了一种从单个图像中找到对象姿势的方法。我们假设我们可以在图像中检测和匹配对象的四个或更多非共面特征点，并且我们知道它们在对象上的相对几何形状。
该算法被称为Posit，在其经典文章“25 行代码中基于模型的对象姿势”中进行了描述（可在其网站上获得，第 4 节）。
文章直接链接：http ://www.cfar.umd.edu/~daniel/daniel_papersfordownload/Pose25Lines.pdf OpenCV 实现： http: //opencv.willowgarage.com/wiki/Posit
这个想法是通过缩放的正交投影反复逼近透视投影，直到收敛到准确的姿势。
对于我的 OpenGL 引擎，以下片段会将鼠标/屏幕坐标转换为 3D 世界坐标。阅读评论以了解正在发生的事情的实际描述。
/* 功能：YCamera :: 计算世界坐标
参数：x 鼠标 x 坐标
y 鼠标 y 坐标
vec 在哪里存储坐标
返回：不适用
描述：将鼠标坐标转换为世界坐标
*/
void YCamera :: CalculateWorldCoordinates(float x, float y, YVector3 *vec)
{
// START
GLint viewport[4];
GLdouble mvmatrix[16], projmatrix[16];
GLint real_y;
GLdouble mx, my, mz;
glGetIntegerv(GL_VIEWPORT, viewport);
glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);
real_y = viewport[3] - (GLint) y - 1; // viewport[3] is height of window in pixels
gluUnProject((GLdouble) x, (GLdouble) real_y, 1.0, mvmatrix, projmatrix, viewport, &mx, &my, &mz);
/* 'mouse' is the point where mouse projection reaches FAR_PLANE.
World coordinates is intersection of line(camera->mouse) with plane(z=0) (see LaMothe 306)
Equation of line in 3D:
(x-x0)/a = (y-y0)/b = (z-z0)/c
Intersection of line with plane:
z = 0
x-x0 = a(z-z0)/c <=> x = x0+a(0-z0)/c <=> x = x0 -a*z0/c
y = y0 - b*z0/c
*/
double lx = fPosition.x - mx;
double ly = fPosition.y - my;
double lz = fPosition.z - mz;
double sum = lx*lx + ly*ly + lz*lz;
double normal = sqrt(sum);
double z0_c = fPosition.z / (lz/normal);
vec->x = (float) (fPosition.x - (lx/normal)*z0_c);
vec->y = (float) (fPosition.y - (ly/normal)*z0_c);
vec->z = 0.0f;
}
从二维空间将有 2 个可以构建的有效矩形。在不知道原始矩阵投影的情况下，您将不知道哪个是正确的。这与“盒子”问题相同：您看到两个正方形，一个在另一个内部，其中 4 个内部顶点连接到 4 个相应的外部顶点。你是从自上而下还是自下而上看一个盒子？
话虽如此，您正在寻找矩阵变换 T ，其中...
{{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}, {x4, y4, z4}} x T = {{x1, y1}, {x2, y2}, { x3, y3}, {x4, y4}}
(4 x 3) x T = (4 x 2)
所以 T 必须是 (3 x 2) 矩阵。所以我们有6个未知数。
现在在 T 上建立一个约束系统并用 Simplex 求解。要构建约束，您知道通过前两个点的线必须平行于通过后两个点的线。您知道通过点 1 和 3 的线必须平行于通过点 2 和 4 的线。您知道通过 1 和 2 的线必须与通过点 2 和 3 的线正交。您知道长度从 1 到 2 的线的长度必须等于从 3 到 4 的线的长度。你知道从 1 到 3 的线的长度必须等于从 2 到 4 的线的长度。
为了让这更容易，你知道矩形，所以你知道所有边的长度。
这应该会给你很多限制来解决这个问题。
当然，要回来，你可以找到T-inverse。
@Rob：是的，有无限数量的投影，但不是无限数量的点必须满足矩形要求的项目。
@nlucaroni：是的，这只有在投影中有四个点时才能解决。如果矩形仅投影到 2 个点（即矩形的平面与投影面正交），则无法解决此问题。
嗯...我应该回家写这个小宝石。这听起来很有趣。
更新：
To follow up on Rons approach: You can find your z-values if you know how you've rotated your rectangle.
The trick is to find the projective matrix that did the projection. Fortunately this is possible and even cheap to do. The relevant math can be found in the paper "Projective Mappings for Image Warping" by Paul Heckbert.
http://pages.cs.wisc.edu/~dyer/cs766/readings/heckbert-proj.pdf
This way you can recover the homogenous part of each vertex back that was lost during projection.
Now you're still left with four lines instead of points (as Ron explained). Since you know the size of your original rectangle however nothing is lost. You can now plug the data from Ron's method and from the 2D approach into a linear equation solver and solve for z. You get the exact z-values of each vertex that way.
Note: This just works because:
It's a special case really.
Hope it helps, Nils
假设这些点确实是矩形的一部分，我给出一个通用的想法：
找到最大间距的两个点：这些最有可能定义对角线（例外：矩形几乎平行于 YZ 平面的特殊情况，留给学生）。称它们为 A、C。计算 BAD、BCD 角度。与直角相比，这些为您提供 3d 空间中的方向。要了解 z 距离，您需要将投影边与已知边相关联，然后根据 3d 投影方法（是 1/z 吗？），您在正确的轨道上了解距离。
感谢@Vegard 的出色回答。我稍微清理了代码：
import pandas as pd
import numpy as np
class Point2:
def __init__(self,x,y):
self.x = x
self.y = y
class Point3:
def __init__(self,x,y,z):
self.x = x
self.y = y
self.z = z
# 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):
#print 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 Point2(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.)
# 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(mat, amount):
from copy import deepcopy
from random import randrange, uniform
mat2 = deepcopy(mat)
mat2[randrange(4)][randrange(4)] += uniform(-amount, amount)
return mat2
def approximate(mat, amount, n=1000):
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(1000):
mat,est = approximate(mat, 1)
print mat
print est
.1 的近似调用对我不起作用，所以我把它拿出来了。我也运行了一段时间，最后我检查它是在
[[0.7576315397559887, 0, 0.11439449272592839, -0.314856490473439],
[0.06440497208710227, 1, -0.5607502645413118, 0.38338196981556827],
[0, 0, 1, 0],
[0.05421620936883742, 0, -0.5673977598434641, 2.693116299312736]]
误差在 0.02 左右。
如果没有人回答，我回家后会拿出我的线性代数书。但是@DG，并非所有矩阵都是可逆的。奇异矩阵是不可逆的（当行列式 = 0 时）。这实际上会一直发生，因为投影矩阵必须具有 0 和 1 的特征值，并且是正方形的（因为它是幂等的，所以 p^2 = p）。
一个简单的例子是，[[0 1][0 1]] 因为行列式 = 0，那是在 x = y 线上的投影！
您在 2D 表面上的投影具有无限多个 3D 矩形，它们将投影到相同的 2D 形状。
这样想：你有四个 3D 点组成 3D 矩形。称它们为 (x0,y0,z0)、(x1,y1,z1)、(x2,y2,z2) 和 (x3,y3,z3)。当您将这些点投影到 xy 平面上时，您会删除 z 坐标：(x0,y0)、(x1,y1)、(x2,y2)、(x3,y3)。
现在，您想投影回 3D 空间，您需要对 z0,..,z3 进行逆向工程。但是任何一组 z 坐标 a）在点之间保持相同的 xy 距离，并且 b）保持矩形的形状都可以。所以，这个（无限）集合的任何成员都可以：{(z0+i, z1+i, z2+i, z3+i) | 我<-R}。
编辑@Jarrett：想象一下你解决了这个问题并最终在 3D 空间中得到了一个矩形。现在，想象一下在 z 轴上上下滑动该矩形。那些无限量的平移矩形都具有相同的 xy 投影。你怎么知道你找到了“正确的”？
编辑#2：好的，这是我对这个问题的评论——一种更直观的推理方法。
想象一下在你的桌子上方放一张纸。假设纸的每个角落都附有一个无重量的激光笔，指向桌子。纸是3D物体，桌上的激光笔点是2D投影。
现在，您如何仅通过激光笔点来判断纸张离桌面有多高？
你不能。上下移动纸张。无论纸张的高度如何，激光笔仍会照在桌子上的相同点上。
在反向投影中找到 z 坐标就像试图仅根据桌子上的激光指针点来找到纸张的高度。
当您从 3D 投影到 2D 时，您会丢失信息。
在单点的简单情况下，逆投影将为您提供穿过 3d 空间的无限光线。
立体重建通常从两个 2D 图像开始，然后将两者都投影回 3D。然后寻找产生的两条 3D 射线的交点。
投影可以采取不同的形式。正交或透视。我猜你假设正交投影？
在您的情况下，假设您拥有原始矩阵，您将在 3D 空间中拥有 4 条光线。然后，您将能够通过 3d 矩形尺寸限制问题并尝试解决。
该解决方案将不会是唯一的，因为围绕平行于 2d 投影平面的任一轴的旋转在方向上将是不明确的。换句话说，如果 2d 图像垂直于 z 轴，那么围绕 x 轴顺时针或逆时针旋转 3d 矩形将产生相同的图像。对于 y 轴也是如此。
在矩形平面平行于 z 轴的情况下，您有更多的解决方案。
由于您没有原始投影矩阵，因此任何投影中都存在任意比例因子会引入进一步的歧义。您无法区分投影中的缩放和在 z 轴方向上的 3d 平移。如果您只对 3d 空间中 4 个点彼此相关而不是 2d 投影平面的相对位置感兴趣，这不是问题。
在透视投影中，事情变得更加困难......
如果您知道该形状是平面中的矩形，则可以大大进一步限制该问题。您当然无法确定“哪个”平面，因此您可以选择它位于 z=0 且其中一个角位于 x=y=0 的平面上，并且边缘平行于 x/y 轴。
因此 3d 中的点是 {0,0,0}、{w,0,0}、{w,h,0} 和 {0,h,0}。我很确定绝对尺寸不会被发现，所以只有 w/h 的比率是相关的，所以这是一个未知数。
相对于该平面，相机必须位于空间中的某个点 cx,cy,cz，必须指向 nx,ny,nz 方向（长度为 1 的向量，因此其中一个是多余的），并且具有 focus_length/image_width w 的因数。这些数字变成了一个 3x3 的投影矩阵。
这给出了总共 7 个未知数：w/h、cx、cy、cz、nx、ny 和 w。
你总共有 8 个已知：4 个 x+y 对。
所以这个可以解决。
下一步是使用 Matlab 或 Mathmatica。
是的，蒙特卡洛有效，但我为这个问题找到了更好的解决方案。此代码完美运行（并使用 OpenCV）：
Cv2.CalibrateCamera(new List<List<Point3f>>() { points3d }, new List<List<Point2f>>() { points2d }, new Size(height, width), cameraMatrix, distCoefs, out rvecs, out tvecs, CalibrationFlags.ZeroTangentDist | CalibrationFlags.FixK1 | CalibrationFlags.FixK2 | CalibrationFlags.FixK3);
此函数采用已知的 3d 和 2d 点、屏幕大小并返回旋转 (rvecs[0])、平移 (tvecs[0]) 和相机的内在值矩阵。这是你需要的一切。