open3D体素化

体素化：
定义：体素（voxel）是像素（pixel）、体积（volume）和元素（element）的组合词，相当于3D空间中的像素。
将物体的几何形式表示转换成最接近该物体的体素表示形式，产生体数据，包含模型的表面信息和内部属性。 表示3D模型的体素跟表示2D图像的像素相似，只不过从二维的点扩展到三维的立方体单元。 体素化能够对模 型进行简化，得到均匀一致的网格。
体素化的优点：
1、点云的数据将在内存中有序存储
2、数据有序存储和降采样，能够处理大规模的数据
3、可以将二维的技术用到三维上
体素化的缺点：
1、丢失信息，与分辨率有关
2、占用内存，与分辨率有关
3、稀疏的点云体素化会构建很多空体素而造成

体素化的方式：
比较普遍的算法主要分为以下两种：
• 对三角面片的顶点，边线以及三 角面分别体素化；
• 对每个已经划分好的体素判断与 三角面片是否相交。

Open3D 体素函数解析：
create_from_triangle_mesh(input, voxel_size)
Open3d中从三角网构建体素：Open3D 提供了从三角形网格创建体素网格的方法create_from_triangle_mesh。它返回一个体素网格，其中所有与三角形相交的体素都设置为1 ，所有其他体素都设置为 0。参数voxel_size定义体素网格的分辨率。

Open3d中从点云构建体素：
create_from_point_cloud(input, voxel_size)
体素网格也可以使用create_from_point_cloud方法从点云创建。如果点云的至少一个点在体素内，则体素被占用。体素的颜色是体素内所有点的平均值。参数voxel_size定义体素网格的分辨率。从给定的 PointCloud 创建 VoxelGrid。给定体素的颜色值是落入其中的点的平均颜色值（如果PointCloud 有颜色）。创建的 VoxelGrid 的边界是根据 PointCloud 计算的。

从体素中构建八叉树：
octree = o3d.geometry.Octree(max_depth=X)
octree.create_from_voxel_grid(voxel_grid)
Open3D从 VoxelGrid几何构造八叉树create_from_voxel_grid。输入的每个体素都VoxelGrid被视为 3D 空间中的 一个点，其坐标对应于体素的原点。

体素下采样：体素下采样使用常规体素网格从输入点云创建均匀下采样点云。voxel_size为体素的尺寸大小，体素的尺寸越大，下采样的倍数越大，点云也就越稀疏。voxel_size越小，会生成更细的网格和更多的点。空间八叉树的叶子结点里，选择一个点代表叶子里的所有点。
它通常用作许多点云处理任务的预处理步骤。
该算法分两步运行：

1. 点被存储到体素中。
2. 每个占用的体素通过平均内部所有点生成一个点。 downpcd = pcd.voxel_down_sample(voxel_size=0.01)
   判断点是否在占用的体素内：
   queries = np.asarray(pcd.points)
   output = voxel_grid.check_if_included(o3d.utility.Vector3dVector(queries))

import open3d as o3d
import numpy as np


def xyz_spherical(xyz):
    x = xyz[0]
    y = xyz[1]
    z = xyz[2]
    r = np.sqrt(x * x + y * y + z * z)
    r_x = np.arccos(y / r)
    r_y = np.arctan2(z, x)
    return [r, r_x, r_y]


def get_rotation_matrix(r_x, r_y):
    rot_x = np.asarray([[1, 0, 0], [0, np.cos(r_x), -np.sin(r_x)],
                        [0, np.sin(r_x), np.cos(r_x)]])
    rot_y = np.asarray([[np.cos(r_y), 0, np.sin(r_y)], [0, 1, 0],
                        [-np.sin(r_y), 0, np.cos(r_y)]])
    return rot_y.dot(rot_x)


def get_extrinsic(xyz):
    rvec = xyz_spherical(xyz)
    r = get_rotation_matrix(rvec[1], rvec[2])
    t = np.asarray([0, 0, 2]).transpose()
    trans = np.eye(4)
    trans[:3, :3] = r
    trans[:3, 3] = t
    return trans


def preprocess(model):
    min_bound = model.get_min_bound()
    max_bound = model.get_max_bound()
    center = min_bound + (max_bound - min_bound) / 2.0
    scale = np.linalg.norm(max_bound - min_bound) / 2.0
    vertices = np.asarray(model.vertices)
    vertices -= center
    model.vertices = o3d.utility.Vector3dVector(vertices / scale)
    return model


def voxel_carving(mesh, cubic_size, voxel_resolution, w=300, h=300):
    mesh.compute_vertex_normals()
    camera_sphere = o3d.geometry.TriangleMesh.create_sphere(radius=1.0, resolution=10)

    # Setup dense voxel grid.
    voxel_carving = o3d.geometry.VoxelGrid.create_dense(
        width=cubic_size,
        height=cubic_size,
        depth=cubic_size,
        voxel_size=cubic_size / voxel_resolution,
        origin=[-cubic_size / 2.0, -cubic_size / 2.0, -cubic_size / 2.0],
        color=[1.0, 0.7, 0.0])
    # o3d.visualization.draw_geometries([voxel_carving])

    # Rescale geometry.
    camera_sphere = preprocess(camera_sphere)
    mesh = preprocess(mesh)

    # Setup visualizer to render depthmaps.
    vis = o3d.visualization.Visualizer()
    vis.create_window(width=w, height=h, visible=False)
    vis.add_geometry(mesh)
    vis.get_render_option().mesh_show_back_face = True
    ctr = vis.get_view_control()
    param = ctr.convert_to_pinhole_camera_parameters()

    # Carve voxel grid.
    centers_pts = np.zeros((len(camera_sphere.vertices), 3))
    for cid, xyz in enumerate(camera_sphere.vertices):
        # Get new camera pose.
        trans = get_extrinsic(xyz)
        param.extrinsic = trans
        c = np.linalg.inv(trans).dot(np.asarray([0, 0, 0, 1]).transpose())
        centers_pts[cid, :] = c[:3]
        ctr.convert_from_pinhole_camera_parameters(param)

        # Capture depth image and make a point cloud.
        vis.poll_events()
        vis.update_renderer()
        depth = vis.capture_depth_float_buffer(False)

        # Depth map carving method.
        voxel_carving.carve_depth_map(o3d.geometry.Image(depth), param)
        print("Carve view %03d/%03d" % (cid + 1, len(camera_sphere.vertices)))
        # o3d.visualization.draw_geometries([voxel_carving])
        
    vis.destroy_window()

    return voxel_carving


if __name__ == "__main__":

    mesh = o3d.io.read_triangle_mesh("data//轮子.STL")
    cubic_size = 2.0
    voxel_resolution = 128.0

    carved_voxels = voxel_carving(mesh, cubic_size, voxel_resolution)
    print("Carved voxels ...")
    print(carved_voxels)
    o3d.visualization.draw([carved_voxels])

# 从三角网构建体素

import open3d as o3d
import numpy as np

mesh = o3d.io.read_triangle_mesh("data//轮子.STL")
# 缩放到单位尺寸
mesh.scale(1 / np.max(mesh.get_max_bound() - mesh.get_min_bound()), center=mesh.get_center())
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh], window_name="mesh",
                                  mesh_show_back_face=False)  # 显示mesh


voxel_grid = o3d.geometry.VoxelGrid.create_from_triangle_mesh(mesh, voxel_size=0.01)

o3d.visualization.draw_geometries([voxel_grid], window_name="体素",
                                  mesh_show_back_face=False)  # 显示体素

# 从点云构建体素

import open3d as o3d
import numpy as np

mesh = o3d.io.read_triangle_mesh("data//轮子.STL")

N = 2000
pcd = mesh.sample_points_poisson_disk(N)
# fit to unit cube
pcd.scale(1 / np.max(pcd.get_max_bound() - pcd.get_min_bound()), center=pcd.get_center())
pcd.colors = o3d.utility.Vector3dVector(np.random.uniform(0, 1, size=(N, 3)))
o3d.visualization.draw_geometries([pcd])

voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=0.05)
o3d.visualization.draw_geometries([voxel_grid])

# 测试体素内外

queries = np.asarray(pcd.points)
output = voxel_grid.check_if_included(o3d.utility.Vector3dVector(queries))
print(output[:10])

#  体素的表示

voxels = voxel_grid.get_voxels()

# 构建八叉树分布
voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd,
                                                            voxel_size=0.05)
o3d.visualization.draw_geometries([voxel_grid])

octree = o3d.geometry.Octree(max_depth=4)
octree.create_from_voxel_grid(voxel_grid)
o3d.visualization.draw_geometries([octree])

# 体素下采样
import open3d as o3d
pcd = o3d.io.read_point_cloud("data//bunny.pcd")
o3d.visualization.draw_geometries([pcd])
downpcd = pcd.voxel_down_sample(voxel_size=0.01)
o3d.visualization.draw_geometries([downpcd])