admin管理员组

文章数量:1276086

Recently I write a simple code which generates 3D points (cube of 4x4x4 or just some random sampled 11 3D points). Then I rotate this points along Y-axis then add some Z-offset and then calculate the two 2-d perspective projections for original array of 3d points and rotated one. By that I simulate different view on the same 3-d object.

And I getting totally unsatisfactory results while calculating cv2.findFundamentalMat(points1, points2, cv2.USAC_MAGSAC, 0.2, 0.99999, 50000) it is work is very unstable despite the fact that data provided to it may be always the same. I do not see any significant problem in my data its seems (to me) to be more than enough, there 64-points without any noise whats so ever, so why do findFundamentalMatrix works so poorly for me?

Code:

import cv2
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def make_camera(cam_x, cam_y, f):
    X = np.linspace(-cam_x, cam_x, 2)
    Y = np.linspace(-cam_y, cam_y, 2)
    X, Y = np.meshgrid(X, Y)
    # Define the Z values for the surface
    Z = np.zeros_like(X)  # For a flat square surface at Z=0
    O_cam = np.array([0, 0, -f])
    return np.concatenate([np.stack([X, Y, Z], axis=2).reshape(-1, 3), [O_cam]], axis=0)

def show_camera(ax, cam_xyz, color='green'):
    # Plot the surface
    O_cam = cam_xyz[4]
    X, Y, Z = [e[..., 0] for e in np.split(cam_xyz[:4].reshape(2, 2, 3), 3, axis=2)]
    ax.plot_surface(X, Y, Z, color=color, edgecolor='cyan', alpha=0.3)
    ax.scatter([O_cam[0]], [O_cam[1]], [O_cam[2]], s=4, color=color)
    for i in range(4):
        ax.plot([O_cam[0], cam_xyz[i, 0]],
                [O_cam[1], cam_xyz[i, 1]],
                [O_cam[2], cam_xyz[i, 2]], color='cyan')


t = np.pi / 6
R_z_man = np.array([[np.cos(t), -np.sin(t), 0],
                    [np.sin(t), np.cos(t),  0],
                    [0,         0,          1]])

R_x_man = np.array([[1,         0,          0],
                    [0, np.cos(t), -np.sin(t)],
                    [0, np.sin(t),  np.cos(t)]])

R_y_man = np.array([[np.cos(t), 0, -np.sin(t)],
                    [0,         1,          0],
                    [np.sin(t), 0, np.cos(t)]])

R_man = R_y_man

points1 = []
for x in np.linspace(start=-1, stop=1, num=4):
    for y in np.linspace(start=-1, stop=1, num=4):
        for z in np.linspace(start=-1, stop=1, num=4):
            points1.append([x, y, z])

#np.random.seed(42)
#points1 = np.random.rand(11, 3) - 0.5

points1 = np.array(points1)
points2 = points1 @ R_man

points1 += np.array([0, 0, 4])
points2 += np.array([0, 0, 4])

orign_points = np.copy(points1)

points1 = points1[:, :2] / points1[:, [2]]
points2 = points2[:, :2] / points2[:, [2]]

F, mask = cv2.findFundamentalMat(points1, points2, cv2.USAC_MAGSAC, 0.2, 0.99999, 50000)

print(F)

# Camera intrinsic parameters (example values, replace with actual calibration data)
fx = 1 # Focal length in x
fy = 1 # Focal length in y
cx = 0.5  # Principal point x
cy = 0.5  # Principal point y

# Camera calibration matrix
K = np.array([[fx, 0, cx],
              [0, fy, cy],
              [0, 0, 1]])

# Compute the essential matrix
E = K.T @ F @ K
print("Essential Matrix:\n", E)

# Decompose the essential matrix to get rotation and translation
U, S, Vt = np.linalg.svd(E)

# Ensure a proper rotation matrix
W = np.array([[0, -1, 0],
              [1, 0, 0],
              [0, 0, 1]])

print(R_man.round(2))   // Manual rotation

R_cal = U @ W @ Vt

print(R_cal.round(2).T) // Calculated rotation is rarely the same as manual why?!

t = U[:, 2]

# Normalize translation vector
print(t)
t = t / np.linalg.norm(t)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
cam_1_xyz = make_camera(cam_x=cx, cam_y=cy, f=(fx + fy)/2)
cam_2_xyz = make_camera(cam_x=cx, cam_y=cy, f=(fx + fy)/2)
show_camera(ax, cam_xyz=cam_1_xyz, color='red')
cam_2_xyz = cam_2_xyz @ R_cal# + t
show_camera(ax, cam_xyz=cam_2_xyz, color='blue')

ax.scatter(orign_points[:, 0], orign_points[:, 1], orign_points[:, 2], s=1.0, c='g')

ax.set_xlabel('X axis')
ax.set_ylabel('Y axis')
ax.set_zlabel('Z axis')
ax.set_title('3D Square Surface')
ax.set_aspect('equal')

plt.show()

本文标签: pythonRotation matrix reconstruction from set of 2D projectionsStack Overflow