如何将2D点投影到3D？

我在屏幕空间中有4个2D点，​​我需要将它们反向投影到3D空间。 我知道4个点中的每一个都是3D旋转的刚性矩形的一个angular，我知道矩形的大小。 我怎样才能从这个3D坐标？

我没有使用任何特定的API，也没有现有的投影matrix。 我只是寻找基本的math来做到这一点。 当然，没有足够的数据将单个2D点转换为3D，而没有其他参考，但是我想如果你有4个点，你就知道它们在同一个平面上彼此成直angular，你知道他们之间的距离，你应该能够从那里弄清楚。 不幸的是我不能很好的解决问题。

这可能属于摄影测量的范围之内，但谷歌search没有让我得到任何有用的信息。

Answer 1

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

简而言之，algorithm如下：随机扰动您的投影matrix，直到它将已知的三维坐标投影到您已知的二维坐标上。

这是来自Thomas the Tank Engine的照片：

假设我们使用GIMP来查找我们认为是地平面上的正方形的二维坐标（不pipe它是否真的是正方形取决于您对深度的判断）：

我在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]与4×4投影matrix相乘，然后将x和y分量除以z以实际获得透视校正，将每个3D坐标投影到2D中。 除了gluProject()也考虑了当前的视口并考虑了一个单独的模型视图matrix（我们可以假设模型视图matrix是单位matrix）之外，这或多或less是gluProject（）所做的。 查看gluProject()文档非常方便，因为我实际上需要一个适用于OpenGL的解决scheme，但要小心文档在公式中缺lessz除法。

请记住，该algorithm是从一些投影matrix开始，随机扰动它，直到它给出我们想要的投影。 所以我们要做的就是投射四个3D点中的每一个，看看我们到达我们想要的2D点有多接近。 如果我们的随机扰动导致投影的二维点更接近我们上面标记的点，那么我们保持这个matrix作为我们最初（或之前）猜测的改进。

让我们来定义我们的观点：

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

我们需要从一些matrix开始，单位matrix似乎是一个自然的select：

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

我们需要实际实现投影（基本上是一个matrix乘法）：

 def project(p, mat): x = mat[0][0] * px + mat[0][1] * py + mat[0][2] * pz + mat[0][3] * 1 y = mat[1][0] * px + mat[1][1] * py + mat[1][2] * pz + mat[1][3] * 1 w = mat[3][0] * px + mat[3][1] * py + mat[3][2] * pz + mat[3][3] * 1 return Point(720 * (x / w + 1) / 2., 576 - 576 * (y / w + 1) / 2.) 

这基本上就是gluProject()所做的工作，720和576分别是图像的宽度和高度（即视口），我们从576中减去以计算从顶部计算y坐标的事实，而OpenGL通常计数他们从底部。 你会注意到我们没有计算z，这是因为我们并不需要它（尽pipe确保它落在OpenGL用于深度缓冲区的范围内可能会很方便）。

现在我们需要一个函数来评估我们离正确的解决scheme有多近。 这个函数返回的值是我们将用来检查一个matrix是否比另一个好。 我select走平方距离，即：

 # The squared distance between two points a and b def norm2(a, b): dx = bx - ax dy = by - ay 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) 

为了干扰matrix，我们简单地选取一个元素在一定范围内随机扰动：

 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) 

我发现这已经给出了相当准确的答案。 经过一段时间，我发现matrix是：

 [ [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左右。 （请注意，我们所说的元素在计算中没有被使用，实际上并没有从我们的初始matrix中改变;那是因为改变这些元素不会改变评估的结果，所以变化永远不会被执行。

我们可以使用glLoadMatrix()将matrix传递给OpenGL（但要记住要先将其转置，并记住将模型视图matrix加载到单位matrix中）：

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

从mathangular度来看，这当然不是很优雅; 你不会得到一个封闭的forms方程，你可以插入你的数字，并得到一个直接的（和准确的）答案。 然而，它允许你添加额外的约束，而不必担心复杂的方程; 例如，如果我们想要结合高度，我们可以使用房子的这个angular落，并且（在我们的评估函数中）说从地面到屋顶的距离应该是这样的，然后再运行algorithm。 所以是的，这是一种蛮横的行为，但是工作，并且运作良好。

Answer 2

D. DeMenthondevise了一种algorithm，当知道物体的模型时，从2D图像中的特征点计算物体的姿态 （它在空间中的位置和方向） – 这是你确切的问题 ：

我们描述了一种从单个图像中找出对象姿态的方法。 我们假设我们可以在图像中检测和匹配对象的四个或更多非共面特征点，并且我们知道它们在对象上的相对几何关系。

该algorithm被称为Posit，并在其中描述了经典文章“25行代码中的基于模型的对象姿势”（可在其网站上获得 ，第4节）。

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

这个想法是通过缩放的正交投影反复逼近透视投影，直到收敛到准确的姿态。

Answer 3

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

你有一个方形的标记（二维条码），你想find它的姿势（平移和相对于相机旋转），find标记的四个边缘。 综述图片

我不知道这个领域的最新贡献，但至less在某一点上（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

（PS – 我知道这是一个老话题，但不pipe怎样，这个post可能对某人有帮助）

Answer 4

对于我的OpenGL引擎，下面的剪辑将鼠标/屏幕坐标转换成3D世界坐标。 阅读评论以了解正在发生的事情的实际描述。

 / *function：YCamera :: CalculateWorldCoordinates
     自variables：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; 

}

Answer 5

从2-D空间将会有2个有效的矩形可以build立。 不知道原始matrix投影，你不会知道哪一个是正确的。 这与“盒子”问题是一样的：你看到两个正方形，一个在另一个里面，4个顶点连接到4个外部顶点。 你是从上到下还是从下往上看盒子？

这就是说，你正在寻找一个matrix变换T，其中…

{{x1，y1，z1}，{x2，y2，z2}，{x3，y3，z3}，{x4，y4，z4}} x T = {{x1，y1}，{x2，y2} x3，y3}，{x4，y4}}

（4×3）×T =（4×2）

所以T必须是（3 x 2）matrix。 所以我们有6个未知数。

现在在T上构build一个约束系统，并用Simplex来解决。 为了build立约束条件，你知道通过前两个点的直线必须与通过后两个点的直线平行。 你知道一条通过点1和3的线必须平行于通过点2和4的线。你知道一条通过1和2的线必须与通过点2和3的线正交。你知道长度从1到2的行必须等于从3到4的行的长度。您知道从1到3的行的长度必须等于从2到4的行的长度。

为了使这更容易，你知道矩形，所以你知道所有方面的长度。

这应该给你很多的限制来解决这个问题。

当然，要找回来，你可以findT-inverse。

@Rob：是的，有无数的投影，但不是无数的项目，其中的点必须满足矩形的要求。

@ nlucaroni：是的，只有在投影中有四个点才能解决。 如果矩形投影到2点（即矩形的平面与投影面正交），则无法解决。

嗯…我应该回家写这个小gem。 这听起来很有趣。

更新：

1. 除非你修正了其中一个要点，否则有无数的预测。 如果你固定原始矩形的点，那么有两个可能的原始矩形。

Answer 6

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

find最大距离的两个点：这些最有可能定义一个对angular线（例外：矩形几乎平行于YZ平面，留给学生的特殊情况）。 称他们为A，C。计算BAD，BCDangular度。 这些，与直angular相比，给你在三维空间的方向。 要了解z距离，需要将投影边与已知边相关，然后基于三维投影方法（是1 / z？），您可以在正确的轨道上了解距离。

Answer 7

要跟进Rons的方法：如果你知道如何旋转你的矩形，你可以find你的Z值。

诀窍是find投影的投影matrix。 幸运的是，这是可能的，甚至很便宜。 有关的math可以在Paul Heckbert的论文“图像变形的投影映射”中find。

~dyer/cs766/readings/heckbert-proj.pdf

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

现在你还剩下四条线，而不是点（正如罗恩解释）。 既然你知道你的原始矩形的大小，但没有任何东西丢失。 现在，您可以将来自Ron方法和2D方法的数据插入线性方程求解器，并求解z。 您可以获得每个顶点的精确z值。

注意：这只是因为：

1. 原来的形状是一个矩形
2. 您知道三维空间中矩形的确切大小。

这真是一个特例。

希望它有帮助，尼尔斯

Answer 8

您在2D表面上的投影具有无限多的3D矩形，它们将投影到相同的2D形状。

用这种方式考虑一下：你有四个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投影。 你怎么知道你find了“正确的”？

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

想象一下，在你的桌子上面拿着一张纸。 假装纸张的每个angular落都附有一个重量轻的激光指针，指向桌面。 纸是3D物体，桌子上的激光指针点是2D投影。

现在，你怎么能通过查看激光指针点来判断纸张的高度？

你不能。 直接上下移动纸张。 不pipe纸张的高度如何，激光指示器仍会照在桌子上的相同位置上。

在反向投影中findz坐标就像试图根据桌面上的激光指针点来找出纸张的高度。

Answer 9

当你从3D投影到2D时，你会失去信息。

在单点的简单情况下，反向投影会给你一个无限的光线通过三维空间。

立体重build通常从两个二维图像开始并投影到3D。 然后寻找产生的两条3D射线的交点。

投影可以采取不同的forms。 正交或透视。 我猜你正在假设正交投影？

在你的情况下，假设你有原始的matrix，你将在3D空间中有4条光线。 然后，您可以通过您的三维矩形尺寸限制问题并尝试解决。

解决scheme将不会是唯一的，因为围绕任何一个平行于2D投影平面的轴的旋转在方向上都是不明确的。 换句话说，如果二维图像垂直于z轴，则围绕x轴顺时针或逆时针旋转三维矩形将产生相同的图像。 对于y轴同样如此。

在矩形平面与z轴平行的情况下，您甚至有更多的解决scheme。

由于您没有原始投影matrix，任何投影中存在的任意比例因子都会导致进一步的模糊。 您无法区分投影中的缩放比例和z轴方向上的三维平移。 如果您只关心三维空间中4点的相对位置，而不是二维投影的平面，则这不是问题。

在透视投影事情变得更难…

Answer 10

如果你知道这个形状是一个平面上的矩形，那么你可以进一步限制这个问题。 你当然无法弄清楚“哪个”平面，所以你可以select它躺在z = 0的平面上，其中一个angular落在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的vector，因此其中一个是多余的），并且具有focal_length / image_width w的因素 这些数字变成了一个3×3的投影matrix。

这总共有7个未知数：w / h，cx，cy，cz，nx，ny和w。

你总共有8个知识：4 x + y对。

所以这可以解决。

下一步是使用Matlab或Mathmatica。

Answer 11

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

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

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

Answer 12

如果没有人回答我回到家，我会把我的线性代数书拿出来。 但@ DG，并不是所有的matrix都是可逆的。 奇异matrix是不可逆的 （当行列式为0时）。 这实际上会一直发生，因为投影matrix必须有0和1的特征值，并且是正方形的（因为它是幂等的，所以p ^ 2 = p）。

一个简单的例子是行列式= 0的[[0 1] [0 1]]，这就是x = y！

Answer 13

是的，蒙特卡罗的工作，但我find了这个问题更好的解决scheme。 此代码完美工作（并使用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]）和内部值matrix。 这是你需要的一切。

Answer 14

感谢@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] * px + mat[0][1] * py + mat[0][2] * pz + mat[0][3] * 1 y = mat[1][0] * px + mat[1][1] * py + mat[1][2] * pz + mat[1][3] * 1 w = mat[3][0] * px + mat[3][1] * py + mat[3][2] * pz + 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 = bx - ax dy = by - ay 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左右的误差。