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

Question

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

我没有使用任何特定的API,我没有现有的投影矩阵.我只是在寻找基本的数学来做到这一点.当然没有足够的数据将单个2D点转换为3D而没有其他参考,但我想如果你有4个点,你知道它们在同一个平面上都是直角相交的,而且你知道它们之间的距离,你应该能够从中找出它.不幸的是,我无法解释如何.

这可能属于摄影测量的范畴,但谷歌搜索它并没有让我获得任何有用的信息.

Answer 1

好吧,我来到这里寻找答案并没有找到简单明了的东西,所以我继续前进并做了一个愚蠢但有效(而且相对简单)的事情:蒙特卡洛优化.

非常简单地说,算法如下:随机扰动您的投影矩阵,直到它将您已知的3D坐标投影到您已知的2D坐标.

这是Thomas the Tank Engine的静态照片:

假设我们使用GIMP来找到地平面上我们认为是正方形的二维坐标(它是否真的是一个正方形取决于您对深度的判断):

我2D图像中得到4分:(318, 247),(326, 312),(418, 241),和(452, 303).

按照惯例,我们说,这些点应该对应于三维点:(0, 0, 0),(0, 0, 1),(1, 0, 0),和(1, 0, 1).换句话说,y = 0平面中的单位平方.

将这些3D坐标中的每一个投影到2D中是通过将4D矢量[x, y, z, 1]与4×4投影矩阵相乘,然后将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()

从数学的角度来看,这不是很优雅; 你没有得到一个封闭形式的等式,你只需将你的数字插入并得到一个直接(和准确)的答案.但是,它确实允许您添加其他约束,而不必担心使方程复杂化; 例如,如果我们想要合并高度,我们可以使用房子的那个角,并说(在我们的评估函数中)从地面到屋顶的距离应该是某某,并再次运行算法.所以,是的,这是一种蛮力,但有效,并且运作良好.

Answer 2

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

您有一个方形标记(2D条形码),并且在找到标记的四个边缘后,您想要找到它的Pose(相对于相机的平移和旋转). 综述图片

我不知道该领域的最新贡献,但至少有一点(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 - 我知道这有点老话题,但无论如何,帖子可能对某人有帮助)

Answer 3

对于我的 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;
}

Answer 4

D. DeMenthon设计了一种算法,用于在了解物体模型时从2D图像中的特征点计算物体的姿态(其在空间中的位置和方向) - 这是您的确切问题:

我们描述了一种从单个图像中找到对象姿势的方法.我们假设我们可以在图像中检测并匹配对象的四个或更多非共面特征点,并且我们知道它们在对象上的相对几何形状.

该算法被称为Posit,并在其经典文章"25行代码中的基于模型的对象姿势"(可在其网站上获得,第4节)中描述.

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

该想法是通过缩放的正交投影重复近似透视投影,直到收敛到精确的姿势.

Answer 5

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

找到具有最大间距的两个点：这些点最有可能定义对角线（例外：矩形几乎平行于 YZ 平面的特殊情况，留给学生）。称它们为 A、C。计算 BAD、BCD 角。与直角相比，这些角度可以为您提供 3D 空间中的方向。要了解 z 距离，您需要将投影边与已知边相关联，然后根据 3d 投影方法（是 1/z 吗？），您就可以正确了解距离。