基于多深度图的3D物体重建：旋转对齐与融合问题求助

问题背景 #

基于转盘旋转拍摄的多幅570×570深度图重建3D物体，高光谱相机保持静止。当前已完成高光谱数据立方体到二维深度图的转换，使用旋转矩阵对齐点云时出现重叠与旋转错误，无法正确合并为3D模型。

现有代码 #

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

# Settings
num_pictures = 10  # This is a choice of the user.
rotation_angle = 360 / num_pictures  # Angle difference per picture
image_size = 570  # 2D-array is 570x570
camera_distance = 162  # Depth distance from camera to rotation-axis / center of object

# Simulate a object_distance_map (dummy data for tests)
object_distance_map = np.random.uniform(100, 300, (image_size, image_size))
object_distance_map -= camera_distance  # Correctie toepassen

# Assenstelsel define
x_range = np.linspace(-image_size // 2, image_size // 2, image_size)
y_range = np.linspace(0, image_size, image_size)

X, Y = np.meshgrid(x_range, y_range)
Z = object_distance_map  # Depth information

# 3D-puntenwolk initiëren
points_3d = np.column_stack((X.flatten(), Y.flatten(), Z.flatten()))

# Function to create a rotation matrix around the y-axis
def rotation_matrix_y(theta):
    theta_rad = np.radians(theta)
    return np.array([
        [np.cos(theta_rad), 0, np.sin(theta_rad)],
        [0, 1, 0],
        [-np.sin(theta_rad), 0, np.cos(theta_rad)]
    ])

# Transform 3D points per rotation angle
all_points = []
for i in range(num_pictures):
    theta = i * rotation_angle
    R = rotation_matrix_y(theta)
    rotated_points = points_3d @ R.T  
    all_points.append(rotated_points)

# Combine all points
all_points = np.vstack(all_points)

# Visualisation
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(all_points[:, 0], all_points[:, 1], all_points[:, 2], s=0.5, alpha=0.5)
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis") # y-axis points upwarts
ax.set_zlabel("Z-axis") # z-axis should be depth-axis
ax.set_title("3D-reconstruction of object")
plt.show()

核心问题分析 #

1. 深度值与坐标系映射错误
   • 直接将像素索引作为3D点的X/Y坐标，不符合相机成像的透视投影逻辑，导致点云空间位置失真。
   • object_distance_map -= camera_distance计算逻辑错误：深度值是相机到物体的距离，需转换为相机坐标系下的Z值后再映射到世界坐标系，而非简单减法。
2. 旋转矩阵应用对象错误
   • 代码重复旋转同一张深度图的点云，而非对每张不同角度的深度图点云应用对应旋转，导致点云重叠混乱。

修正方案与代码实现 #

1. 正确的深度图转点云流程 #

先将像素坐标通过相机内参转换为相机坐标系下的3D点，再映射到世界坐标系（考虑相机与旋转轴的位置关系）。

2. 修正后的完整代码 #

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

# 配置参数
num_pictures = 10
rotation_angle = 360 / num_pictures
image_size = 570
camera_distance = 162  # 相机到世界坐标系旋转轴（Y轴）的距离
fx = fy = 500  # 相机焦距（根据实际标定值调整，单位：像素）
cx = cy = image_size // 2  # 相机主点（图像中心）

def depth_to_world_point_cloud(depth_map, rotation_angle_deg):
    """将单张深度图转换为世界坐标系下的点云"""
    # 生成像素坐标网格
    x_pixel = np.linspace(0, image_size - 1, image_size)
    y_pixel = np.linspace(0, image_size - 1, image_size)
    X_pixel, Y_pixel = np.meshgrid(x_pixel, y_pixel)
    
    # 像素坐标转相机坐标系：透视投影逆变换
    Z_cam = depth_map  # Z_cam为相机到物体的距离（深度值）
    X_cam = (X_pixel - cx) / fx * Z_cam
    Y_cam = (Y_pixel - cy) / fy * Z_cam
    
    # 相机坐标系转世界坐标系：
    # 1. 旋转矩阵：物体绕世界Y轴旋转rotation_angle_deg，相机看到的点云需反向旋转得到世界坐标
    theta_rad = np.radians(rotation_angle_deg)
    rotation_mat = np.array([
        [np.cos(theta_rad), 0, np.sin(theta_rad)],
        [0, 1, 0],
        [-np.sin(theta_rad), 0, np.cos(theta_rad)]
    ])
    # 2. 平移：相机位于世界坐标系(0, 0, camera_distance)位置
    points_cam = np.column_stack((X_cam.flatten(), Y_cam.flatten(), Z_cam.flatten()))
    points_world = (points_cam @ rotation_mat.T) + np.array([0, 0, camera_distance])
    
    return points_world

# 处理所有深度图（实际应用中替换为读取真实深度图）
all_world_points = []
for idx in range(num_pictures):
    current_theta = idx * rotation_angle
    # 生成模拟深度图（替换为真实深度图加载：np.load("depth_{}.npy".format(idx))）
    depth_map = np.random.uniform(100, 300, (image_size, image_size))
    world_points = depth_to_world_point_cloud(depth_map, current_theta)
    all_world_points.append(world_points)

# 合并所有点云
merged_points = np.vstack(all_world_points)

# 可视化结果
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(merged_points[:, 0], merged_points[:, 1], merged_points[:, 2], s=0.3, alpha=0.6)
ax.set_xlabel("世界X轴")
ax.set_ylabel("世界Y轴")
ax.set_zlabel("世界Z轴")
ax.set_title("修正后的3D物体重建点云")
plt.show()

3. 旋转矩阵的替代方案：ICP配准 #

如果转盘旋转精度不足，旋转角度存在误差，可使用ICP（迭代最近点）算法自动对齐多视角点云，推荐使用open3d库实现：

import open3d as o3d

# 假设已生成每张深度图对应的Open3D点云列表point_clouds
point_clouds = []
for idx in range(num_pictures):
    current_theta = idx * rotation_angle
    depth_map = np.random.uniform(100, 300, (image_size, image_size))
    points = depth_to_world_point_cloud(depth_map, current_theta)
    # 转换为Open3D点云对象
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(points)
    point_clouds.append(pcd)

# 以第一张点云为基准，ICP配准并合并
target_pcd = point_clouds[0]
merged_pcd = target_pcd

for source_pcd in point_clouds[1:]:
    # ICP配准
    reg_result = o3d.pipelines.registration.registration_icp(
        source_pcd, target_pcd, max_correspondence_distance=1.0,
        estimation_method=o3d.pipelines.registration.TransformationEstimationPointToPoint(),
        criteria=o3d.pipelines.registration.ICPConvergenceCriteria(max_iteration=2000)
    )
    # 将源点云转换到基准坐标系
    source_pcd.transform(reg_result.transformation)
    merged_pcd += source_pcd

# 体素下采样优化点云（减少冗余）
merged_pcd_down = merged_pcd.voxel_down_sample(voxel_size=0.5)

# 可视化合并结果
o3d.visualization.draw_geometries([merged_pcd_down], window_name="ICP配准后3D重建")

# 可选：泊松重建生成网格模型
mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(merged_pcd_down, depth=9)
o3d.visualization.draw_geometries([mesh], window_name="网格模型重建")

关键优化点 #

• 相机内参标定：如果有相机的实际内参（fx, fy, cx, cy），替换代码中的模拟值，可大幅提升点云精度。
• 无效深度值过滤：实际深度图中可能存在无效值（如0或NaN），转换点云时需过滤这些点：

  valid_mask = (depth_map > 0) & ~np.isnan(depth_map)
  X_cam = X_cam[valid_mask]
  Y_cam = Y_cam[valid_mask]
  Z_cam = Z_cam[valid_mask]
  

• 点云去重：合并后可使用体素下采样或统计滤波去除噪声与重叠点。

内容的提问来源于stack exchange，提问作者KrokanteKorneel