Open3D Voxelization 體素化

Voxelization 體素化

點雲和三角形網格是非常靈活但不規則的幾何類型。體素網格（voxel grid）是在 3D 網格上定義的 3D 中的另一種幾何類型，體素可以被視為 2D 中像素（pixel）的 3D 對應物。Open3D 具有可用于處理體素網格的幾何類型VoxelGrid。

From triangle mesh 從三角形網格

Open3D 提供了從三角形網格建立體素網格的方法create_from_triangle_mesh。它傳回一個體素網格，其中所有與三角形相交的體素都設定為1 ，所有其他體素都設定為 0。參數voxel_size定義體素網格的分辨率。

print('input')
bunny = o3d.data.BunnyMesh()
mesh = o3d.io.read_triangle_mesh(bunny.path)

# fit to unit cube
mesh.scale(1 / np.max(mesh.get_max_bound() - mesh.get_min_bound()),
           center=mesh.get_center())
o3d.visualization.draw_geometries([mesh])

print('voxelization')
voxel_grid = o3d.geometry.VoxelGrid.create_from_triangle_mesh(mesh,
                                                              voxel_size=0.05)
o3d.visualization.draw_geometries([voxel_grid])      

From point cloud 從點雲

體素網格也可以使用create_from_point_cloud方法從點雲建立。如果點雲的至少一個點在體素内，則體素被占用。體素的顔色是體素内所有點的平均值。參數voxel_size定義體素網格的分辨率。

print('input')
armadillo = o3d.data.ArmadilloMesh()
mesh = o3d.io.read_triangle_mesh(armadillo.path)

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])

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

Inclusion test 包含測試

體素網格還可用于測試點是否在占用的體素内。check_if_included方法将(n,3)數組作為輸入并輸出bool數組。

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

Voxel carving 體素雕刻

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,
                 use_depth=True,
                 surface_method='pointcloud'):
   mesh.compute_vertex_normals()
   camera_sphere = o3d.geometry.TriangleMesh.create_sphere()

   # 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])

   # 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
   pcd_agg = o3d.geometry.PointCloud()
   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)
       pcd_agg += o3d.geometry.PointCloud.create_from_depth_image(
           o3d.geometry.Image(depth),
           param.intrinsic,
           param.extrinsic,
           depth_scale=1)

       # depth map carving method
       if use_depth:
           voxel_carving.carve_depth_map(o3d.geometry.Image(depth), param)
       else:
           voxel_carving.carve_silhouette(o3d.geometry.Image(depth), param)
       print("Carve view %03d/%03d" % (cid + 1, len(camera_sphere.vertices)))
   vis.destroy_window()

   # add voxel grid survace
   print('Surface voxel grid from %s' % surface_method)
   if surface_method == 'pointcloud':
       voxel_surface = o3d.geometry.VoxelGrid.create_from_point_cloud_within_bounds(
           pcd_agg,
           voxel_size=cubic_size / voxel_resolution,
           min_bound=(-cubic_size / 2, -cubic_size / 2, -cubic_size / 2),
           max_bound=(cubic_size / 2, cubic_size / 2, cubic_size / 2))
   elif surface_method == 'mesh':
       voxel_surface = o3d.geometry.VoxelGrid.create_from_triangle_mesh_within_bounds(
           mesh,
           voxel_size=cubic_size / voxel_resolution,
           min_bound=(-cubic_size / 2, -cubic_size / 2, -cubic_size / 2),
           max_bound=(cubic_size / 2, cubic_size / 2, cubic_size / 2))
   else:
       raise Exception('invalid surface method')
   voxel_carving_surface = voxel_surface + voxel_carving

   return voxel_carving_surface, voxel_carving, voxel_surface

armadillo_path = r'../data/ArmadilloMesh.ply'
mesh = o3d.io.read_triangle_mesh(armadillo_path)

visualization = True
cubic_size = 2.0
voxel_resolution = 128.0

voxel_grid, voxel_carving, voxel_surface = voxel_carving(
   mesh, cubic_size, voxel_resolution)

print("surface voxels")
print(voxel_surface)
o3d.visualization.draw_geometries([voxel_surface])

print("carved voxels")
print(voxel_carving)
o3d.visualization.draw_geometries([voxel_carving])

print("combined voxels (carved + surface)")
print(voxel_grid)
o3d.visualization.draw_geometries([voxel_grid])