3d 三角测量与 DLT 一起使用,但不与投影矩阵一起使用 cv2.triangulatePoints

3d triangulation working with DLT but not with projection matrix using cv2.triangulatePoints

大纲

我有一个校准的立体相机设置,每个相机都有 11 个 DLT 系数(使用 easyWand package 估计的系数,该 easyWand package 使用在场景中移动的 'wand' 校准对象)。使用 DLT 方法的 3D 投影和获得的系数工作正常并产生合理的结果(见下图)。

使用 DLT 三角测量获得抛出物体的抛物线轨迹(绘图已旋转以与重力对齐)

但是,当我将 DLT 系数转换为投影矩阵时(P,其中 x = P X,其中 x 是二维像素坐标,P 是3x4 矩阵和 X 是包含对象 3d 坐标的齐次 4x1 矩阵)- 轨迹的 'shape' 有意义但值没有意义(x、y、z 的值接近 ~0.01) .

使用基于投影矩阵的方法获得轨迹cv2.triangulatePoints(绘图已旋转以与重力对齐)

任何人都可以解释为什么基于投影矩阵的方法会产生奇怪的缩放坐标或者这里做错了什么?

差异没有意义,因为 DLT 系数和投影矩阵是可相互转换的(例如,参见 this link)。

复制结果的数据和代码:

# The DLT coefficients for the two cameras 
import numpy as np 
import cv2 
camera1_dlt = np.array([-26.618, -169.460,  57.650,  306.560,  24.733,  11.723,  169.830,
        95.812,  0.220, -0.161,  0.110])
camera2_dlt = np.array([-75.615, -80.218,  43.704,  412.710,  1.655,  21.247,  110.420,
        235.580,  0.061, -0.145,  0.112])

# The projection matrix was derived from the DLT coefficients using the function belo

def extract_P_from_dlt(dltcoefs):
    '''function which recreates the Projection matrix
    '''
    
    # dltcoefs contains basically the correct matrix, it just needs to be
    # reshaped correctly
    all_coefs = np.concatenate((dltcoefs, np.array([1]))).flatten()
    P = all_coefs.reshape((4,3)).T
    return P

# The subset of common *undistorted* 2D points are given below here
cam1_pts = np.array([ [ 98.221,  358.293],
                      [ 142.927,  348.913],
                      [ 186.275,  344.950],
                      [ 225.152,  341.567]],dtype='float32')

cam2_pts = np.array([ [ 314.633,  195.970],
                      [ 335.168,  200.769],
                      [ 354.778,  207.820],
                      [ 372.543,  215.284]], dtype='float32')

# generate projection matrix 
Pcam1 = extract_P_from_dlt(camera1_dlt)
Pcam2 = extract_P_from_dlt(camera2_dlt)

# Now get 3D positions using cv2.triangulatePositions
xyz_cv2 = []
for pt1, pt2 in zip(cam1_pts, cam2_pts):
    pt1_homog, pt2_homog = (X.reshape(1,1,2) for X in [pt1, pt2])
    position = cv2.triangulatePoints(Pcam1, Pcam2, pt1_homog, pt2_homog)
    final_xyz = cv2.convertPointsFromHomogeneous(position.T).flatten()
    xyz_cv2.append(final_xyz)

系数的堆叠应该是下面的格式。

[ L1  L2  L3  L4 ]
[ L5  L6  L7  L8 ]
[ L9 L10 L11 L12 ]

但是你的extract_P_from_dlt给了,

[ L1  L4  L7  L10 ]
[ L2  L5  L8  L11]
[ L3  L6  L9  L12 ]

此外,L12 不是堆叠 1,而是堆叠 0。这是更新后的功能。

def extract_P_from_dlt(dltcoefs):
    '''function which recreates the Projection matrix'''
    dltcoefs = np.append(dltcoefs, 0)
    norm = np.linalg.norm(dltcoefs)
    dltcoefs = dltcoefs/norm

    P = dltcoefs.reshape(3,4)
    return P

整体代码,

# The DLT coefficients for the two cameras 
import numpy as np 
import cv2 

camera1_dlt = np.array([-26.618, -169.460,  57.650,  306.560,  24.733,  11.723,  169.830,
        95.812,  0.220, -0.161,  0.110])
camera2_dlt = np.array([-75.615, -80.218,  43.704,  412.710,  1.655,  21.247,  110.420,
        235.580,  0.061, -0.145,  0.112])

# The projection matrix was derived from the DLT coefficients using the function belo

def extract_P_from_dlt(dltcoefs):
    '''function which recreates the Projection matrix'''
    dltcoefs = np.append(dltcoefs, 0)
    norm = np.linalg.norm(dltcoefs)
    dltcoefs = dltcoefs/norm

    P = dltcoefs.reshape(3,4)
    return P


# The subset of common *undistorted* 2D points are given below here
cam1_pts = np.array([ [ 98.221,  358.293],
                      [ 142.927,  348.913],
                      [ 186.275,  344.950],
                      [ 225.152,  341.567]],dtype='float32')

cam2_pts = np.array([ [ 314.633,  195.970],
                      [ 335.168,  200.769],
                      [ 354.778,  207.820],
                      [ 372.543,  215.284]], dtype='float32')

# generate projection matrix 
Pcam1 = extract_P_from_dlt(camera1_dlt)
Pcam2 = extract_P_from_dlt(camera2_dlt)

# Now get 3D positions using cv2.triangulatePositions
xyz_cv2 = []
for pt1, pt2 in zip(cam1_pts, cam2_pts):
    pt1_homog, pt2_homog = (X.reshape(1,1,2) for X in [pt1, pt2])
    position = cv2.triangulatePoints(Pcam1, Pcam2, pt1_homog, pt2_homog)
    final_xyz = cv2.convertPointsFromHomogeneous(position.T).flatten()
    xyz_cv2.append(final_xyz)

生成下图。我不确定规范化部分。

更新的答案 - 包含了解图像来源约定的课程!

在与 easyWand 包的开发人员来回交流之后 - 事实证明图像的原始约定起着 的作用。

来自 easyWand 包的 DLT 系数是在假设 bottom-left 来源的情况下生成的。我的输入像素数据假定为 top-left 原点。这种起源约定的差异导致了错误的三角测量! 不需要任何类型的规范化!

这是提供验证结果的最终工作代码:


# The DLT coefficients for the two cameras 
import numpy as np 
import cv2 

camera1_dlt = np.array([-26.618, -169.460,  57.650,  306.560,  24.733,  11.723,  169.830,
        95.812,  0.220, -0.161,  0.110])
camera2_dlt = np.array([-75.615, -80.218,  43.704,  412.710,  1.655,  21.247,  110.420,
        235.580,  0.061, -0.145,  0.112])

# The projection matrix was derived from the DLT coefficients using the function below

def extract_P_from_dlt_v2(dltcoefs):
    '''Change in order of placing dlt coefficients 
    into 3x4 matrix - NO normalisation
    '''
    dltcoefs = np.append(dltcoefs, 1)
    dltcoefs = dltcoefs

    P = dltcoefs.reshape(3,4)
    return P


# The subset of common *undistorted* 2D points are given below here with 
# top-left image origin (upper row is 0 and left column is 0)
cam1_pts = np.array([ [ 98.221,  358.293],
                      [ 142.927,  348.913],
                      [ 186.275,  344.950],
                      [ 225.152,  341.567]],dtype='float32')

cam2_pts = np.array([ [ 314.633,  195.970],
                      [ 335.168,  200.769],
                      [ 354.778,  207.820],
                      [ 372.543,  215.284]], dtype='float32')

# move the origin from upper left to bottom left
# Now bottom row is 0, left column is still as before
cam1_pts[:,1] = 511 - cam1_pts[:,1]
cam2_pts[:,1] = 511 - cam2_pts[:,1]

# generate projection matrix 
Pcam11 = extract_P_from_dlt_v2(camera1_dlt)
Pcam22 = extract_P_from_dlt_v2(camera2_dlt)

# Now get 3D positions using cv2.triangulatePositions
xyz_cv2 = []
for pt1, pt2 in zip(cam1_pts, cam2_pts):
    pt1_homog, pt2_homog = (X.reshape(1,1,2) for X in [pt1, pt2])
    position = cv2.triangulatePoints(Pcam11, Pcam22, pt1_homog, pt2_homog)
    final_xyz = cv2.convertPointsFromHomogeneous(position.T).flatten()
    xyz_cv2.append(final_xyz)


# Now get 3D positions using cv2.triangulatePositions
xyz_cv2 = []
for pt1, pt2 in zip(cam1_pts, cam2_pts):
    pt1_homog, pt2_homog = (X.reshape(1,1,2) for X in [pt1, pt2])
    position = cv2.triangulatePoints(Pcam11, Pcam22, pt1_homog, pt2_homog)
    final_xyz = cv2.convertPointsFromHomogeneous(position.T).flatten()
    xyz_cv2.append(final_xyz)

xyz_cv2 = np.array(xyz_cv2).reshape(-1,3)

plt.figure()
a0 = plt.subplot(111, projection='3d')
plt.plot(xyz_cv2[:,0], xyz_cv2[:,2], xyz_cv2[:,1], '*') # the y and z need to be swapped 

输出的 xyz 坐标与直接使用 DLT 重建方法生成的坐标相匹配 - 该方法采用 11 个 DLT 系数对对象位置进行三角测量。