PCD与深度图像互转时的点云丢失问题求助

我正在开展PCD（点云数据）与深度图像的互转工作，在将PCD转换为深度图像再转回PCD的过程中，重建后的PCD丢失了部分点云图层，但深度图像显示正常。

运行结果说明 #

• 原始与重建PCD对比：红色为原始PCD，绿色为重建PCD
• 掩码前后的深度图像

核心目标是获得完整的重建点云，恳请提供可行的解决技术方案。

实现代码 #

import numpy as np
import pandas as pd
from pypcd4 import PointCloud
import matplotlib.pyplot as plt
import open3d as o3d
import os


# Defined parameters
vfov = [-10, 10]            # Vertical field of view in degrees
hfov = [-60, 60]            # Horizontal field of view in degrees
size = [360, 880]           # Size of the depth image (height, width)
depth_scale = 7.97          # Scale factor for depth values
log_scale = True  

def pcd2range_scala3(pcd):
    """
    Conversion of point cloud to depth image using Scala3 projection.
    """
    fov_up = np.deg2rad(vfov[1])  # Field of view up in radians
    fov_down = np.deg2rad(vfov[0])  # Field of view down in radians
    fov_range = fov_up - fov_down  # Total field of view in radians

    # Get depth (distance) of all points
    depth = np.linalg.norm(pcd, axis=1)
    print(f"Depth range before filtering: min={np.min(depth)}, max={np.max(depth)}")

    # Mask points out of range
    mask = np.logical_and(depth > 0, depth < 200)
    depth, pcd = depth[mask], pcd[mask]

    # Get scan components
    scan_x, scan_y, scan_z = pcd[:, 0], pcd[:, 1], pcd[:, 2]

    # Calculate yaw and pitch angles
    yaw = np.arctan2(scan_y, scan_x)
    pitch = np.arcsin(scan_z / depth)

    # Get projections in image coordinates
    proj_x = (yaw - np.deg2rad(hfov[0])) / np.deg2rad(hfov[1] - hfov[0])
    proj_y = (pitch - np.deg2rad(vfov[0])) / (fov_range)
    
    # Scale to image size
    proj_x *= size[1]  # Width
    proj_y *= size[0]  # Height

    # Round and clamp for use as index
    proj_x = np.clip(np.floor(proj_x).astype(int), 0, size[1] - 1)
    proj_y = np.clip(np.floor(proj_y).astype(int), 0, size[0] - 1)

    # Order by decreasing depth
    order = np.argsort(depth)[::-1]
    proj_x, proj_y = proj_x[order], proj_y[order]
    depth = depth[order]

    # Project depth
    proj_range = np.full(size, -1, dtype=np.float32)
    proj_range[proj_y, proj_x] = depth

    # Apply transformation
    if log_scale:
        proj_range = np.log2(proj_range + 1e-3 + 1)  # Added a small epsilon
        proj_range = proj_range / depth_scale
        proj_range = proj_range * 2. - 1.

    print(f"Depth image range after transformation: min={np.min(proj_range)}, max={np.max(proj_range)}")

    return proj_range


def range2pcd_scala3(range_img):
    """
    Conversion of depth image back to point cloud using Scala3 projection.
    """
    fov_up = np.deg2rad(vfov[1])  # Field of view up in radians
    fov_down = np.deg2rad(vfov[0])  # Field of view down in radians
    fov_range = fov_up - fov_down  # Total field of view in radians
    size_range = range_img.shape

    # Inverse transform from depth
    depth = (range_img + 1.) / 2.
    depth = (depth * depth_scale).flatten()
    depth = np.exp2(depth) - 1

    print(f"Reconstructed depth range: min={np.min(depth)}, max={np.max(depth)}")

    # Compute image coordinates
    scan_x, scan_y = np.meshgrid(np.arange(size_range[1]), np.arange(size_range[0]))
    scan_x = scan_x.astype(np.float64) / size_range[1]
    scan_y = scan_y.astype(np.float64) / size_range[0]

    yaw = (scan_x * np.deg2rad(hfov[1] - hfov[0]) + np.deg2rad(hfov[0])).flatten()
    pitch = (scan_y * fov_range + np.deg2rad(vfov[0])).flatten()

    pcd = np.zeros((len(yaw), 3))
    pcd[:, 0] = np.cos(yaw) * np.cos(pitch) * depth
    pcd[:, 1] = np.sin(yaw) * np.cos(pitch) * depth
    pcd[:, 2] = np.sin(pitch) * depth

    # Mask out invalid points
    mask = np.logical_and(depth > 0, depth < 200)
    pcd = pcd[mask, :]

    return pcd


def save_pcd_file(pcd, file_path):
    """
    Save the point cloud to a .pcd file.
    """
    rgb = np.ones((pcd.shape[0], 1)) * [0.7]  # Simple gray color for visualization
    pcd_with_rgb = np.hstack((pcd, rgb))  # Add RGB to point cloud

    pc = PointCloud.from_xyzi_points(pcd_with_rgb)
    pc.save(file_path)

def visualize_point_clouds(pcd1, pcd2, title1='Point Cloud 1', title2='Point Cloud 2'):
    """
    Visualize two point clouds side by side using Open3D.
    """
    # Create point clouds
    point_cloud1 = o3d.geometry.PointCloud()
    point_cloud1.points = o3d.utility.Vector3dVector(pcd1)
    point_cloud1.paint_uniform_color([1, 0, 0])  # Red for the original

    point_cloud2 = o3d.geometry.PointCloud()
    point_cloud2.points = o3d.utility.Vector3dVector(pcd2)
    point_cloud2.paint_uniform_color([0, 1, 0])  # Green for the reconstructed

    # Visualize together
    o3d.visualization.draw_geometries([point_cloud1, point_cloud2],
                                      zoom=0.3412,
                                      front=[0.4257, -0.2125, -0.8795],
                                      lookat=[2.6172, 2.0475, 1.532],
                                      up=[-0.0694, -0.9768, 0.2024])


if __name__ == "__main__":
    input_path = os.path.expanduser("~/Desktop/cutting_fov/updated_code/input/880x360/07_3.pcd")
    output_path = os.path.expanduser("~/Desktop/cutting_fov/updated_code/output/800x320/07_3_reconstructed.pcd")

    # Load the point cloud
    pcd = PointCloud.from_path(input_path)
    df_pcd = pd.DataFrame(pcd.pc_data)
    arr_pcd = df_pcd.to_numpy()

    points = arr_pcd.astype(np.float32) + 1e-8

    # Load xyz coordinates
    xyz = points[:, :3]

    # Convert to depth image before applying the FOV mask
    proj_range_before = pcd2range_scala3(xyz)

    # Define the FOV mask
    mask_y = np.zeros(size[0], dtype=bool)
    mask_y[40:320] = True
    mask_x = np.zeros(size[1], dtype=bool)
    mask_x[80:800] = True

    # Apply the FOV mask
    mask = np.outer(mask_y, mask_x)

    # Convert to depth image after applying the FOV mask
    proj_range_after = np.where(mask, proj_range_before, -1)

    # Normalize proj_range for visualization
    proj_range_display_before = np.clip((proj_range_before + 1) / 2, 0, 1)
    proj_range_display_after = np.clip((proj_range_after + 1) / 2, 0, 1)

    # Display depth images
    plt.figure(figsize=(16, 8))

    plt.subplot(1, 2, 1)
    plt.imshow(proj_range_display_before, cmap='jet', interpolation='nearest')
    plt.title('Depth Image Before Applying Mask')

    plt.subplot(1, 2, 2)
    plt.imshow(proj_range_display_after, cmap='jet', interpolation='nearest')
    plt.title('Depth Image After Applying Mask')

    plt.tight_layout()
    plt.show()

    # Convert depth image back to point cloud
    pcd_reconstructed = range2pcd_scala3(proj_range_after)

    # Save the reconstructed point cloud
    save_pcd_file(pcd_reconstructed, output_path)

    # Visualize the original and reconstructed point clouds
    visualize_point_clouds(xyz, pcd_reconstructed, title1='Original Point Cloud', title2='Reconstructed Point Cloud')

1. 修复逆变换中的无效值处理逻辑 #

当前range2pcd_scala3函数未提前过滤深度图像中的无效值（-1），导致无效值参与深度逆变换计算，生成异常点后被过滤，最终丢失有效区域的点。修改如下：

def range2pcd_scala3(range_img):
    """
    Conversion of depth image back to point cloud using Scala3 projection.
    """
    fov_up = np.deg2rad(vfov[1])
    fov_down = np.deg2rad(vfov[0])
    fov_range = fov_up - fov_down
    size_range = range_img.shape

    # 先过滤深度图像中的无效值(-1)
    valid_mask = range_img != -1
    scan_y, scan_x = np.where(valid_mask)
    range_img_valid = range_img[valid_mask]

    # 仅对有效区域做深度逆变换
    depth = (range_img_valid + 1.) / 2.
    depth = (depth * depth_scale)
    depth = np.exp2(depth) - 1

    print(f"Reconstructed depth range: min={np.min(depth)}, max={np.max(depth)}")

    # 计算有效区域的归一化坐标
    scan_x_norm = scan_x.astype(np.float64) / size_range[1]
    scan_y_norm = scan_y.astype(np.float64) / size_range[0]

    yaw = (scan_x_norm * np.deg2rad(hfov[1] - hfov[0]) + np.deg2rad(hfov[0]))
    pitch = (scan_y_norm * fov_range + np.deg2rad(vfov[0]))

    pcd = np.zeros((len(yaw), 3))
    pcd[:, 0] = np.cos(yaw) * np.cos(pitch) * depth
    pcd[:, 1] = np.sin(yaw) * np.cos(pitch) * depth
    pcd[:, 2] = np.sin(pitch) * depth

    # 过滤异常深度值
    mask = np.logical_and(depth > 0, depth < 200)
    pcd = pcd[mask, :]

    return pcd

2. 解决多点投影覆盖问题 #

原pcd2range_scala3函数中，多个点投影到同一像素时仅保留最远点，导致近点丢失。若需保留所有点，可修改为每个像素存储多个深度值：

def pcd2range_scala3(pcd):
    """
    Conversion of point cloud to depth image using Scala3 projection, preserving all points per pixel.
    """
    fov_up = np.deg2rad(vfov[1])
    fov_down = np.deg2rad(vfov[0])
    fov_range = fov_up - fov_down

    depth = np.linalg.norm(pcd, axis=1)
    print(f"Depth range before filtering: min={np.min(depth)}, max={np.max(depth)}")

    mask = np.logical_and(depth > 0, depth < 200)
    depth, pcd = depth[mask], pcd[mask]

    scan_x, scan_y, scan_z = pcd[:, 0], pcd[:, 1], pcd[:, 2]
    yaw = np.arctan2(scan_y, scan_x)
    pitch = np.arcsin(scan_z / depth)

    proj_x = (yaw - np.deg2rad(hfov[0])) / np.deg2rad(hfov[1] - hfov[0])
    proj_y = (pitch - np.deg2rad(vfov[0])) / fov_range
    
    proj_x *= size[1]
    proj_y *= size[0]

    # 用round替代floor提升投影精度
    proj_x = np.clip(np.round(proj_x).astype(int), 0, size[1] - 1)
    proj_y = np.clip(np.round(proj_y).astype(int), 0, size[0] - 1)

    # 存储每个像素的所有深度值
    proj_range = [[[] for _ in range(size[1])] for _ in range(size[0])]
    for x, y, d in zip(proj_x, proj_y, depth):
        proj_range[y][x].append(d)

    # 生成平均深度图像（若需保留所有点，可直接存储列表格式）
    proj_range_avg = np.full(size, -1, dtype=np.float32)
    for y in range(size[0]):
        for x in range(size[1]):
            if proj_range[y][x]:
                proj_range_avg[y][x] = np.mean(proj_range[y][x])

    if log_scale:
        proj_range_avg = np.log2(proj_range_avg + 1e-3 + 1)
        proj_range_avg = proj_range_avg / depth_scale
        proj_range_avg = proj_range_avg * 2. - 1.

    print(f"Depth image range after transformation: min={np.min(proj_range_avg)}, max={np.max(proj_range_avg)}")

    return proj_range_avg

3. 提升投影坐标精度 #

将原投影函数中的floor取整改为round，减少坐标偏移导致的点丢失：

proj_x = np.clip(np.round(proj_x).astype(int), 0, size[1] - 1)
proj_y = np.clip(np.round(proj_y).astype(int), 0, size[0] - 1)

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