Open3D 三维重建-Marching Cubes (行进立方体)

目录

一、概述

1.1原理

1.2实现步骤

1.3应用场景

二、代码实现

2.1关键函数

2.1.1步骤

2.1.2函数代码

2.2完整代码

三、实现效果

3.1原始点云

3.2重建后点云

---

Open3D点云算法汇总及实战案例汇总的目录地址：

Open3D点云算法与点云深度学习案例汇总（长期更新）-CSDN博客

---

一、概述

        行进立方体算法（Marching Cubes）是一种用于从三维标量场数据（如体素网格）生成等值面（isosurface）的算法。其基本原理是通过遍历三维标量场中的每个体素单元，根据每个体素的顶点值，构建出相应的三角网格。

        Open3D目前并没有直接的行进立方体算法实现，需要结合其他库来完成，例如scikit-image中的marching_cubes方法。结合 Open3D 和 scikit-image编写行进立方体算法，可以实现从点云数据生成高质量的等值面

1.1原理

1. 遍历三维标量场中的每个体素单元。
2. 根据体素单元的八个顶点值，确定该单元内部是否包含等值面。
3. 使用预定义的查找表，根据顶点值的组合生成对应的三角形面片。
4. 将所有生成的三角形面片合并，构成完整的等值面。

1.2实现步骤

1. 加载点云数据：读取点云数据文件。
2. 估计法向量：计算点云中每个点的法向量。
3. 生成体素网格：从点云数据生成体素网格。
4. 使用行进立方体算法：从体素网格生成等值面。
5. 可视化和保存结果：可视化生成的等值面，并将其保存到文件中。

1.3应用场景

• 医学成像：如CT、MRI数据的三维重建。

• 计算机图形学：生成多边形网格模型。

• 科学可视化：如气象数据、地质数据的三维可视化

二、代码实现

2.1关键函数

• 这段代码使用行进立方体算法（Marching Cubes）从体素网格生成三维等值面。具体来说，它将输入的体素网格转换为三维体积数据，并利用
• scikit-image 库的 marching_cubes 方法生成等值面。生成的等值面被转换为 Open3D 的三角网格对象 TriangleMesh，可以用于可视化和进一步处理。

2.1.1步骤

1. 获取体素索引：从体素网格中提取所有体素的索引，并将其转换为 NumPy 数组。
2. 获取边界和体素大小：获取体素网格的最小和最大边界，以及体素的大小。
3. 计算体素网格形状：根据边界和体素大小计算体素网格的形状。
4. 创建体素数组：创建一个空的体素数组，用于存储体积数据。
5. 填充体素数组：将体素网格中的每个体素标记为1，生成体积数据。
6. 行进立方体算法：使用 scikit-image 的 marching_cubes 方法生成等值面。
7. 转换为三角网格：将生成的顶点、面片和法向量转换为 Open3D 的 TriangleMesh 对象。
8. 返回生成的三角网格：返回生成的 TriangleMesh 对象，可以用于可视化和进一步处理。

2.1.2函数代码

def marching_cubes_from_voxel_grid(voxel_grid, level):
    """
    使用行进立方体算法从体素网格生成等值面。

    参数:
    voxel_grid (open3d.geometry.VoxelGrid): 输入的体素网格。
    level (float): 等值面的阈值。

    返回:
    open3d.geometry.TriangleMesh: 生成的三角网格。
    """
    # 获取体素网格中的所有体素索引
    voxel_indices = np.asarray([v.grid_index for v in voxel_grid.get_voxels()])
    # 获取体素网格的最小和最大边界
    min_bound = voxel_grid.get_min_bound()
    max_bound = voxel_grid.get_max_bound()
    # 获取体素大小
    voxel_size = voxel_grid.voxel_size
    # 计算体素网格的形状
    volume_shape = ((max_bound - min_bound) / voxel_size).astype(int) + 1
    # 创建一个空的体素数组
    volume = np.zeros(volume_shape, dtype=np.float32)

    # 将体素网格中的每个体素标记为1
    for voxel in voxel_indices:
        # 计算体素在体素数组中的索引
        grid_index = ((voxel * voxel_size - min_bound) / voxel_size).astype(int)
        # 确保索引在有效范围内
        if np.all(grid_index >= 0) and np.all(grid_index < volume_shape):
            volume[tuple(grid_index)] = 1

    # 使用 scikit-image 的 marching_cubes 方法生成等值面
    verts, faces, normals, _ = measure.marching_cubes(volume, level=level)

    # 将结果转换为 Open3D 的 TriangleMesh
    mesh = o3d.geometry.TriangleMesh()
    mesh.vertices = o3d.utility.Vector3dVector(verts * voxel_size + min_bound)
    mesh.triangles = o3d.utility.Vector3iVector(faces)
    mesh.vertex_normals = o3d.utility.Vector3dVector(normals)
    return mesh

2.2完整代码

import open3d as o3d
import numpy as np
from skimage import measure

def create_voxel_grid_from_point_cloud(pcd, voxel_size):
    """
    将点云转换为体素网格。

    参数:
    pcd (open3d.geometry.PointCloud): 输入的点云数据。
    voxel_size (float): 体素的大小。

    返回:
    open3d.geometry.VoxelGrid: 生成的体素网格。
    """
    voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=voxel_size)
    return voxel_grid

def marching_cubes_from_voxel_grid(voxel_grid, level):
    """
    使用行进立方体算法从体素网格生成等值面。

    参数:
    voxel_grid (open3d.geometry.VoxelGrid): 输入的体素网格。
    level (float): 等值面的阈值。

    返回:
    open3d.geometry.TriangleMesh: 生成的三角网格。
    """
    # 将体素网格转换为3D numpy数组
    voxel_indices = np.asarray([v.grid_index for v in voxel_grid.get_voxels()])
    min_bound = voxel_grid.get_min_bound()
    max_bound = voxel_grid.get_max_bound()
    voxel_size = voxel_grid.voxel_size
    volume_shape = ((max_bound - min_bound) / voxel_size).astype(int) + 1
    volume = np.zeros(volume_shape, dtype=np.float32)

    for voxel in voxel_indices:
        grid_index = ((voxel * voxel_size - min_bound) / voxel_size).astype(int)
        if np.all(grid_index >= 0) and np.all(grid_index < volume_shape):
            volume[tuple(grid_index)] = 1

    # 使用scikit-image的marching_cubes方法生成等值面
    verts, faces, normals, _ = measure.marching_cubes(volume, level=level)

    # 将结果转换为Open3D的TriangleMesh
    mesh = o3d.geometry.TriangleMesh()
    mesh.vertices = o3d.utility.Vector3dVector(verts * voxel_size + min_bound)
    mesh.triangles = o3d.utility.Vector3iVector(faces)
    mesh.vertex_normals = o3d.utility.Vector3dVector(normals)
    return mesh

# 加载点云数据
pcd = o3d.io.read_point_cloud("hand.pcd")
o3d.visualization.draw_geometries([pcd], window_name="Cloud", width=800, height=600)
# 估计法向量
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))

# 生成体素网格
voxel_size = 0.05
voxel_grid = create_voxel_grid_from_point_cloud(pcd, voxel_size)

# 定义等值面阈值
isosurface_level = 0.5

# 使用行进立方体算法从体素网格生成等值面
mesh = marching_cubes_from_voxel_grid(voxel_grid, isosurface_level)

# 可视化生成的等值面
o3d.visualization.draw_geometries([mesh], window_name="Marching Cubes", width=800, height=600)

# 保存生成的等值面
# o3d.io.write_triangle_mesh("marching_cubes_mesh.ply", mesh)

三、实现效果

3.1原始点云

3.2重建后点云