python osgb转换为点云数据 python 点云重建

1.简介

Open3D:一个用于3D数据处理的现代库

Open3D是一个开源库，支持处理3D数据的软件的快速开发。Open3D前端在c++和Python中公开了一组精心选择的数据结构和算法。后端经过高度优化，并设置为并行化。我们欢迎来自开源社区的贡献。

• Open3D的核心功能包括:
• 三维数据结构
• 三维数据处理算法
• 现场重建
• 表面对齐
• 三维可视化
• 物理渲染(PBR)
• 3D机器学习支持PyTorch和TensorFlow
• GPU加速核心3D操作
• c++和Python版本可用

官方：Open3D: A Modern Library for 3D Data Processing — Open3D 0.15.1 documentation

2. 从python开始，深度图转点云

2.1 安装

安装系统ubuntu,mac win10都支持

conda create -n open3d python=3.7
activate open3d -i https://pypi.tuna.tsinghua.edu.cn/simple
# 安装
pip install open3d
# 验证
python -c "import open3d as o3d; print(o3d.__version__)"

测试可视化一个球：test3d.py

import open3d as o3d
mesh = o3d.geometry.TriangleMesh.create_sphere()

mesh.compute_vertex_normals()
o3d.visualization.draw(mesh, raw_mode=True)

2.2可视化人脸点云

OPEN3D支持各种格式的3d文件，pcd,ply等

import pandas as pd
import numpy as np
import open3d as o3d
#ply_point_cloud = o3d.data.PLYPointCloud()
pcd = o3d.io.read_point_cloud("./face.ply")
#pcd = o3d.io.read_point_cloud("./my_points.txt", format='xyz')
#pcd = o3d.io.read_point_cloud("./face.pcd")
print(pcd)
print(np.asarray(pcd.points))
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

2.3从深度图到点云

通常使用TOF等3d摄像头采集的格式一般只是深度图，需要经过转化，python这里的方式，先将深度图转化为3D坐标，存储为numpy格式，然后直接使用open3d转化为可视点云。

原本的csv可视的深度图如下：

data_path = "./face.csv"
w = 320
h = 240
data = pd.read_csv(data_path, header=None)
points = np.zeros((w*h, 3), dtype=np.float32)
n=0
for i in range(h):
    for j in range(w):
        deep = data.iloc[i, j]
        points[n][0] = j
        points[n][1] = i
        points[n][2] = deep
        #points.append([j,i,deep])
        n=n+1

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)
#o3d.io.write_point_cloud("../../test_data/sync.ply", pcd)

print("==========")
print(pcd)
print(np.asarray(pcd.points))
print("==========")




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

点云结果：

这里只是简单的转化，没有根据相机内参进行映射，所以点的距离并不正常

查看相机内参，经过处理后可视化点云：

import pandas as pd
import numpy as np
import open3d as o3d
#ply_point_cloud = o3d.data.PLYPointCloud()
#pcd = o3d.io.read_point_cloud("./face.ply")
#pcd = o3d.io.read_point_cloud("./my_points.txt", format='xyz')
#pcd = o3d.io.read_point_cloud("./face.pcd")

data_path = "./face.csv"

camera_factor = 1;
camera_cx = 180.8664;
camera_cy = 179.088;

camera_fx = 216.75;
camera_fy = 214.62;

w = 320
h = 240
data = pd.read_csv(data_path, header=None)
points = np.zeros((w*h, 3), dtype=np.float32)
n=0
for i in range(h):
    for j in range(w):
        deep = data.iloc[i, j]
        points[n][2] = deep/camera_factor
        points[n][0] = (j-camera_cx)*points[n][2]/camera_fx
        points[n][1] = (i-camera_cy)*points[n][2]/camera_fy
        
        #points.append([j,i,deep])
        n=n+1

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)
#o3d.io.write_point_cloud("../../test_data/sync.ply", pcd)

print("==========")
print(pcd)
print(np.asarray(pcd.points))
print("==========")


#knot_data = o3d.data.KnotMesh()
#mesh = o3d.io.read_triangle_mesh("./face.ply")
#print(mesh)
#o3d.io.write_triangle_mesh("copy_of_knot.ply", mesh)

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

侧面：

2.4 点云分割聚类

保留前景信息，利用聚类和分割函数，也可先进行numpy预处理

import pandas as pd
import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt
#ply_point_cloud = o3d.data.PLYPointCloud()
#pcd = o3d.io.read_point_cloud("./face.ply")
#pcd = o3d.io.read_point_cloud("./my_points.txt", format='xyz')
#pcd = o3d.io.read_point_cloud("./face.pcd")

data_path = "./face.csv"

camera_factor = 10;
camera_cx = 180.8664;
camera_cy = 179.088;

camera_fx = 216.75;
camera_fy = 214.62;

w = 320
h = 240
data = pd.read_csv(data_path, header=None)
points = np.zeros((w*h, 3), dtype=np.float32)

n=0
for i in range(h):
    for j in range(w):
        deep = data.iloc[i, j]
        points[n][2] = deep/camera_factor
        points[n][0] = (j-camera_cx)*points[n][2]/camera_fx
        points[n][1] = (i-camera_cy)*points[n][2]/camera_fy
        
        #points.append([j,i,deep])
        n=n+1

points = points[points[:,2]<100]

pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(points)
#o3d.io.write_point_cloud("../../test_data/sync.ply", pcd)
#pcd.paint_uniform_color([1, 0.706, 0])
print("==========")
print(pcd)
print(np.asarray(pcd.points))
print("==========")

##聚类啦
#with o3d.utility.VerbosityContextManager(
#        o3d.utility.VerbosityLevel.Debug) as cm:
#    labels = np.array(
#        pcd.cluster_dbscan(eps=0.02, min_points=10, print_progress=True))

#max_label = labels.max()
#print(f"point cloud has {max_label + 1} clusters")
#colors = plt.get_cmap("tab20")(labels / (max_label if max_label > 0 else 1))
#colors[labels < 0] = 0
#pcd.colors = o3d.utility.Vector3dVector(colors[:, :3])


plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
                                         ransac_n=3,
                                         num_iterations=1000)
[a, b, c, d] = plane_model
print(f"Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0")

inlier_cloud = pcd.select_by_index(inliers)
inlier_cloud.paint_uniform_color([1.0, 0, 0])
outlier_cloud = pcd.select_by_index(inliers, invert=True)

#画个框玩玩
aabb = pcd.get_axis_aligned_bounding_box()
aabb.color = (1, 0, 0)
obb = pcd.get_oriented_bounding_box()
obb.color = (0, 1, 0)
#knot_data = o3d.data.KnotMesh()
#mesh = o3d.io.read_triangle_mesh("./face.ply")
#print(mesh)
#o3d.io.write_triangle_mesh("copy_of_knot.ply", mesh)

o3d.visualization.draw_geometries([outlier_cloud, aabb],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

3,表面重建

3.1泊松重建

在很多情况下，我们想要生成一个密集的3D几何体，例如，一个三角形网格。然而，从多视角立体视觉方法，或深度传感器，我们只能获得非结构化点云。为了从这个非结构化输入中得到一个三角形网格，我们需要执行表面重建。在文献中有两个方法和Open3D目前实现以下:α形状(Edelsbrunner1983)，球旋转(Bernardini1999)，泊松表面重建[Kazhdan2006]

泊松重建需要法线估计，直接调用：

pcd.normals = o3d.utility.Vector3dVector(np.zeros(
    (1, 3)))  # invalidate existing normals

pcd.estimate_normals()

泊松表面重建还将在点密度低的区域创建三角形，甚至可以外推到某些区域(见上图底部的老鹰输出)。create_from_point_cloud_poisson函数有第二个密度返回值，表示每个顶点的密度。低密度值意味着只支持来自输入点云的少量点。

3.2Alpha shapes重建

alpha形状[Edelsbrunner1983]是凸包的泛化。

tetra_mesh, pt_map = o3d.geometry.TetraMesh.create_from_point_cloud(pcd)
for alpha in np.logspace(np.log10(2.5), np.log10(0.1), num=2):
    print(f"alpha={alpha:.3f}")
    mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(
        pcd, alpha, tetra_mesh, pt_map)
    mesh.compute_vertex_normals()
    o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)

3.3Ball pivoting

球体旋转算法(BPA) [Bernardini1999]是一种与alpha形状相关的表面重建方法。

radii = [0.005, 0.01, 0.02, 0.04]
rec_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(
    pcd, o3d.utility.DoubleVector(radii))
o3d.visualization.draw_geometries([pcd, rec_mesh])