首页 / 知识

关于几何：如何将2D点反向投影到3D中？

2023-04-16 03:54:00

How do I reverse-project 2D points into 3D?

我在屏幕空间中有4个2D点，我需要将它们反向投影回3D空间。我知道这4个点中的每一个都是3D旋转的刚性矩形的一个角，并且知道矩形的大小。如何从中获得3D坐标？

我没有使用任何特定的API，并且没有现有的投影矩阵。我只是在寻找基本的数学来做到这一点。当然，没有足够的数据将单个2D点转换为没有其他参考的3D，但是我想像一下，如果您有4个点，您就知道它们在同一平面上彼此成直角，并且您知道它们之间的距离，那么您应该能够从那里弄清楚。不幸的是我无法完全解决。

这可能属于摄影测量学的范畴，但是google搜索并没有带给我任何有用的信息。

好吧，我来这里是为了寻找答案，却没有找到简单明了的东西，所以我继续做愚蠢但有效(相对简单)的事情：蒙特卡洛优化。

好。

简而言之，算法如下：随机扰动投影矩阵，直到将已知的3D坐标投影到已知的2D坐标为止。

好。

这是Thomas坦克引擎的静态照片：

好。

Thomas the Tank Engine

好。

假设我们使用GIMP在地平面上找到我们认为是正方形的2D坐标(是否真的是正方形取决于您对深度的判断)：

好。

With an outline of the square

好。

我在2D图像中得到四个点：(318, 247)，(326, 312)，(418, 241)和(452, 303)。

好。

按照惯例，我们说这些点应对应3D点：(0, 0, 0)，(0, 0, 1)，(1, 0, 0)和(1, 0, 1)。换句话说，y = 0平面中的单位正方形。

好。

通过将4D向量[x, y, z, 1]与4x4投影矩阵相乘，然后将x和y分量除以z，以实际获得透视校正，可以将这些3D坐标中的每一个投影到2D中。这几乎是gluProject()所做的，除了gluProject()还将当前视口考虑在内，并将一个单独的modelview矩阵也考虑在内(我们可以假设modelview矩阵是恒等矩阵)。查看gluProject()文档非常方便，因为我实际上想要一个适用于OpenGL的解决方案，但是请注意，该文档缺少公式中的z除法。

好。

请记住，算法是从一些投影矩阵开始，然后随机扰动它，直到给出我们想要的投影为止。因此，我们要做的是投影四个3D点中的每一个，并查看我们离想要的2D点有多近。如果我们的随机扰动导致投影的2D点越来越接近上面标记的点，那么我们保留该矩阵作为对我们最初(或先前)猜测的改进。

好。

让我们定义点：

好。

1
2
3
4
5
6
7
8
9
10
11

# 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)

我们需要从一些矩阵开始，单位矩阵似乎是自然的选择：

好。

1
2
3
4
5
6

mat = [
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
]

我们需要实际实现投影(基本上是矩阵乘法)：

好。

1
2
3
4
5

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，这是因为我们在这里确实不需要z(尽管确保它在OpenGL用于深度缓冲区的范围内可能很方便)。

好。

现在，我们需要一个函数来评估我们与正确解决方案的距离。此函数返回的值是我们将用来检查一个矩阵是否优于另一个矩阵的值。我选择按平方距离的总和进行计算，即：

好。

1
2
3
4
5
6
7
8
9
10
11
12

# 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)

要扰动矩阵，我们只需选择一个要扰动某个范围内的随机量的元素即可：

好。

1
2
3
4
5

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()来进行到目前为止我们发现的最佳矩阵，该近似值是：

好。

1
2
3
4
5
6
7
8
9
10
11

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

现在剩下要做的就是运行它...：

好。

1
2
3

for i in xrange(100):
mat = approximate(mat, 1)
mat = approximate(mat, .1)

我发现这已经给出了非常准确的答案。运行一段时间后，我发现的矩阵是：

好。

1
2
3
4
5
6

[
[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中(但请记住先将其转置，并记住用恒等矩阵加载您的modelview矩阵)：

好。

1
2
3
4
5
6
7
8
9

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轴平移以获取沿轨道的不同位置：

好。

1
2
3
4
5
6
7
8
9

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()

With 3D translation

好。

从数学的角度来看，这肯定不是很优雅。您不会得到一个封闭式方程，您只需将数字插入其中就可以得到直接(准确)的答案。但是，它确实允许您添加其他约束，而不必担心使方程复杂化。例如，如果我们也想合并高度，则可以使用房屋的那个角落并说(在评估函数中)从地面到屋顶的距离应该是某某某点，然后再次运行该算法。是的，这是种蛮力，但是行之有效。

好。

Choo choo!

好。

这是基于标记的增强现实的经典问题。

您有一个方形标记(2D条形码)，并且要在找到标记的四个边缘后找到其姿势(相对于摄像机的平移和旋转)。
概述-图片

我尚不清楚该领域的最新贡献，但至少在某种程度上(2009年)，RPP的表现要优于上面提到的POSIT(这确实是一种经典方法)
请查看链接，它们也提供了源。

http://www.emt.tugraz.at/~vmg/schweighofer
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

(附言：我知道这是一个有点老的话题，但是无论如何，该帖子可能会对某人有所帮助)

D. DeMenthon设计了一种算法，可以在知道对象模型时根据2D图像中的特征点计算对象的姿态(其在空间中的位置和方向)-这是您的确切问题：

We describe a method for finding the pose of an object from a single image. We assume that we can detect and match in the image four or more noncoplanar feature points of the object, and that we know their relative geometry on the object.

该算法称为Posit，并在其经典文章"代码行25行中基于模型的对象姿势"中进行了描述(可在其网站上的第4节中找到)。

文章的直接链接：http://www.cfar.umd.edu/~daniel/daniel_papersfordownload/Pose25Lines.pdf
OpenCV实施：http://opencv.willowgarage.com/wiki/Posit

这个想法是用缩放的正投影法反复逼近透视投影，直到收敛到准确的姿势为止。

对于我的OpenGL引擎，以下片段将鼠标/屏幕坐标转换为3D世界坐标。阅读评论，以了解发生了什么。

1
2
3
4
5
6
7

/* FUNCTION: YCamera :: CalculateWorldCoordinates
ARGUMENTS: x mouse x coordinate
y mouse y coordinate
vec where to store coordinates
RETURN: n/a
DESCRIPTION: Convert mouse coordinates into world coordinates
*/

void YCamera :: CalculateWorldCoordinates(float x, float y, YVector3 *vec)
{
// START
GLint viewport[4];
GLdouble mvmatrix[16], projmatrix[16];

}

在二维空间中，将可以构建2个有效的矩形。不知道原始矩阵投影，就不会知道哪一个是正确的。它与"盒子"问题相同：您看到两个正方形，一个正方形在另一个正方形内，其中四个内部顶点分别连接到四个外部顶点。您是从上到下还是自下而上查看一个框？

话虽如此，您正在寻找一个矩阵变换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逆。

@Rob：是的，有无限数量的投影，但是没有无限多个项目的点必须满足矩形的要求。

@nlucaroni：是的，只有在投影中有四个点时，这才可以解决。如果矩形仅投影到2个点(即矩形的平面与投影表面正交)，则无法解决。

嗯...我应该回家写这个小宝石。这听起来很有趣。

更新：

除非您固定其中一个点，否则投影的数量是无限的。如果固定原始矩形的点，则有两个可能的原始矩形。

遵循Rons方法：如果知道如何旋转矩形，则可以找到z值。

诀窍是找到进行投影的投影矩阵。幸运的是，这样做是可能的，甚至是便宜的。相关数学可以在Paul Heckbert的论文"图像变形的投影映射"中找到。

http://pages.cs.wisc.edu/~dyer/cs766/readings/heckbert-proj.pdf

这样，您可以恢复投影过程中丢失的每个顶点的同质部分。

现在您仍然剩下四行而不是点(如罗恩所解释)。由于您知道原始矩形的大小，因此不会丢失任何内容。现在，您可以将Ron方法和2D方法中的数据插入线性方程求解器，并求解z。您可以通过这种方式获得每个顶点的确切z值。

注意：这样做是因为：

原始形状是矩形

您知道3D空间中矩形的确切大小。

真的是特例。

希望能帮助到你，
尼尔斯

假设这些点确实是矩形的一部分，我给出一个通用的想法：

找到两个最大距离的点：最可能定义对角线的一个点(例外：特殊情况下，矩形与YZ平面几乎平行，留给学生)。称它们为A，C。计算BAD，BCD角度。与直角相比，这些使您可以在3d空间中定位。要了解z距离，您需要将投影的边与已知的边相关联，然后基于3d投影方法(是1 / z？)，您将在正确的轨道上了解距离。

感谢@Vegard的出色回答。我整理了一下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

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进行的近似通话对我不起作用，因此我将其取出。我也跑了一段时间，最后我检查了一下

1
2
3
4

[[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。

http://library.wolfram.com/infocenter/Articles/2794/

http://campar.in.tum.de/Students/SepPoseEstimation

http://opencvlibrary.sourceforge.net/Posit可以工作，但是可以
可以发散或循环。

如果您知道形状是平面上的矩形，则可以极大地限制该问题。当然，您不能确定"哪个"平面，因此可以选择它位于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。

从3D投影到2D时，您会丢失信息。

在单点的简单情况下，逆投影会给您无限的光线穿过3d空间。

立体重建通常从两个2d图像开始，然后再投影回3D。然后寻找产生的两条3D射线的交点。

投影可以采用不同的形式。正交或透视。我猜您正在假设正交投影？

在您的情况下，假设您拥有原始矩阵，那么在3D空间中将有4条射线。然后，您将能够通过3d矩形尺寸来约束问题并尝试解决。

该解决方案将不是唯一的，因为围绕平行于2d投影平面的任一轴的旋转在方向上将是模棱两可的。换句话说，如果2d图像垂直于z轴，则围绕x轴顺时针或逆时针旋转3d矩形将产生相同的图像。对于y轴也是如此。

在矩形平面平行于z轴的情况下，您还有更多解决方案。

由于您没有原始的投影矩阵，因此任何投影中都存在一个不确定的比例因子，这进一步引入了歧义。您无法区分投影的缩放比例和沿z轴方向的3d平移。如果您仅对3d空间中4个点的相对位置感兴趣，而彼此相关而不与2d投影的平面相关，那么这不是问题。

从透视图的角度看，事情变得越来越难...

您在2D曲面上的投影具有无限多个3D矩形，这些矩形将投影为相同的2D形状。

这样考虑：您有四个组成3D矩形的3D点。称它们为(x0，y0，z0)，(x1，y1，z1)，(x2，y2，z2)和(x3，y3，z3)。将这些点投影到x-y平面上时，将放下z坐标：(x0，y0)，(x1，y1)，(x2，y2)，(x3，y3)。

现在，要投影回3D空间，需要对z0，..，z3进行反向工程。但是，任何a组z坐标(a)在点之间保持相同的x-y距离，b)保持矩形的形状都可以。因此，该(无限)集合的任何成员都将做：{(z0 + i，z1 + i，z2 + i，z3 + i)|我<-R}。

编辑@Jarrett：假设您解决了这个问题，最后在3D空间中得到一个矩形。现在，想象一下将矩形沿z轴上下滑动。那些无限数量的平移矩形都具有相同的x-y投影。您怎么知道您找到了"正确"的？

编辑2：好的，这是我对这个问题的评论-一种更直观的推理方法。

想象一下，在办公桌上方拿着一张纸。假装纸的每个角上都装有一个失重的激光指示器，该指示器向下指向桌子。纸是3D物体，桌子上的激光指示器点是2D投影。

现在，您如何仅看激光笔的点数就能知道纸离办公桌有多远？

你不能垂直向上和向下移动纸张。无论纸张的高度如何，激光笔仍会照在桌子上的相同位置上。

在反投影中找到z坐标就像试图仅在桌子上基于激光指示器点找到纸张的高度一样。

如果没有人回答，我回家时会拿出线性代数书。但是@ D G，并非所有矩阵都是可逆的。奇异矩阵是不可逆的(行列式= 0时)。实际上这将一直发生，因为投影矩阵必须具有0和1的特征值，并且必须是平方的(因为它是幂等的，所以p ^ 2 = p)。

一个简单的例子是[[0 1] [0 1]]，因为行列式= 0，并且这是x = y线上的投影！

是的，蒙特卡洛工作，但我找到了解决此问题的更好解决方案。这段代码运行完美(并使用OpenCV)：

1	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])和相机的内在值矩阵。这就是您需要的一切。

3d 反向投影 2D 空间