如何用Open3D的project_to_depth_image实现点云投影与正射视图生成

需求与问题背景 #

目标是使用Open3D（可通过pip install open3d快速安装）生成点云的俯视/侧视正射（或近似正射）视图。希望找到最简方法，不想通过Visualizer的capture_screen_float_buffer进行两次渲染，因此尝试使用open3d.t.geometry.PointCloud.project_to_depth_image。

修改官方示例代码，加载任意点云而非从RGBD图像转换的点云后，得到的深度图几乎全为零；放大后仅能看到微小点云。后续添加键盘控制调整外参，但仍存在以下核心问题：

• 能否用任意点云通过project_to_depth_image生成深度图？若可以，该如何操作？若不行，除可视化渲染外还有哪些Open3D技术可选？
• 如何计算外参以获得俯视图？如何设置参数实现正交投影？

修改后的测试代码 #

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
from math import radians

if __name__ == '__main__':

    def get_intrinsic(width, height):
        return o3d.core.Tensor([[1, 0, width * 0.5], 
                                [0, 1, height * 0.5],
                                [0, 0, 1]])

    def get_extrinsic(x = 0, y = 0, z = 0, rx = 0, ry = 0, rz = 0):
        extrinsic = np.eye(4)
        extrinsic[:3,  3] = (x, y, z)
        extrinsic[:3, :3] = o3d.geometry.get_rotation_matrix_from_axis_angle([radians(rx),radians(ry), radians(rz)])
        return extrinsic

    def compute_show_reprojection(pcd, width, height, intrinsic, extrinsic):
        depth_reproj = pcd.project_to_depth_image(width,
                                                  height,
                                                  intrinsic,
                                                  extrinsic,
                                                  depth_scale=5000.0,
                                                  depth_max=10.0)

        
        fig, axs = plt.subplots(1, 2)
        axs[0].imshow(np.asarray(depth.to_legacy()))
        axs[1].imshow(np.asarray(depth_reproj.to_legacy()))
        plt.show()
        
    
    width, height = 640, 480
    intrinsic = get_intrinsic(width, height)
    extrinsic = get_extrinsic()
    # original example data
    tum_data = o3d.data.SampleTUMRGBDImage()
    depth = o3d.t.io.read_image(tum_data.depth_path)
    pcd = o3d.t.geometry.PointCloud.create_from_depth_image(depth,
                                                            intrinsic,
                                                            extrinsic,
                                                            depth_scale=5000.0,
                                                            depth_max=10.0)

    compute_show_reprojection(pcd, width, height, intrinsic, get_extrinsic(z=1, rz=-45))

    # testing a differen point cloud
    pcd_data = o3d.data.PCDPointCloud()
    pcd = o3d.io.read_point_cloud(pcd_data.path)
    pcd = o3d.t.geometry.PointCloud.from_legacy(pcd)
    c   = pcd.get_center().numpy()
    minb= pcd.get_min_bound().numpy()
    maxb= pcd.get_max_bound().numpy()
    print('point cloud center', c)
    print('min', minb)
    print('max', maxb)
    compute_show_reprojection(pcd, width, height, intrinsic, get_extrinsic(c[0], c[1] + 1, c[2] + 1))

带键盘控制的后续测试代码 #

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
from math import radians

import cv2
from time import time

if __name__ == '__main__':

    def get_intrinsic(width, height):
        return o3d.core.Tensor([[1, 0, width * 0.5], 
                                [0, 1, height * 0.5],
                                [0, 0, 1]])

    def get_extrinsic(x = 0, y = 0, z = 0, rx = 0, ry = 0, rz = 0):
        extrinsic = np.eye(4)
        extrinsic[:3,  3] = (x, y, z)
        extrinsic[:3, :3] = o3d.geometry.get_rotation_matrix_from_axis_angle([radians(rx),radians(ry), radians(rz)])
        return extrinsic

    def compute_show_reprojection(pcd, width, height, intrinsic, extrinsic, window_wait=3000):
        now = time()
        depth_reproj = pcd.project_to_depth_image(width,
                                                  height,
                                                  intrinsic,
                                                  extrinsic,
                                                  depth_scale=5000.0,
                                                  depth_max=10.0)

        
        # fig, axs = plt.subplots(1, 2)
        # axs[0].imshow(np.asarray(depth.to_legacy()))
        # axs[1].imshow(np.asarray(depth_reproj.to_legacy()))
        # plt.show()
        print(time() - now)
        cv2.imshow("depth", np.asarray(depth_reproj.to_legacy()))
        return cv2.waitKey(window_wait)
        
    
    width, height = 640, 480
    intrinsic = get_intrinsic(width, height)
    extrinsic = get_extrinsic()
    # original example data
    tum_data = o3d.data.SampleTUMRGBDImage()
    depth = o3d.t.io.read_image(tum_data.depth_path)
    pcd = o3d.t.geometry.PointCloud.create_from_depth_image(depth,
                                                            intrinsic,
                                                            extrinsic,
                                                            depth_scale=5000.0,
                                                            depth_max=10.0)

    # compute_show_reprojection(pcd, width, height, intrinsic, get_extrinsic(z=1, rz=-45))

    # testing a differen point cloud
    # pcd_data = o3d.data.PCDPointCloud()
    # pcd = o3d.io.read_point_cloud(pcd_data.path)
    # pcd = o3d.t.geometry.PointCloud.from_legacy(pcd)
    # test random points point cloud
    points = np.random.rand(100000, 3)
    point_cloud = o3d.geometry.PointCloud()
    point_cloud.points = o3d.utility.Vector3dVector(points)
    pcd = o3d.t.geometry.PointCloud.from_legacy(point_cloud)

    c   = pcd.get_center().numpy()
    minb= pcd.get_min_bound().numpy()
    maxb= pcd.get_max_bound().numpy()
    print('point cloud center', c)
    print('min', minb)
    print('max', maxb)
    key = ' '
    x, y, z = -0.95, -0.95, -0.95
    while key != ord('q'):
        key = compute_show_reprojection(pcd, width, height, intrinsic, get_extrinsic(x, y, z), 40)
        if key == ord('x'):
            x -= 0.1
            print(f"x, y, z = {x, y, z}")
    
        if key == ord('X'):
            x += 0.1
            print(x, y, z)
    
        if key == ord('y'):
            y -= 0.1
            print(f"x, y, z = {x, y, z}")
    
        if key == ord('Y'):
            y += 0.1
            print(f"x, y, z = {x, y, z}")
    
        if key == ord('z'):
            z -= 0.1
            print(f"x, y, z = {x, y, z}")
    
        if key == ord('Z'):
            z += 0.1
            print(f"x, y, z = {x, y, z}")

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