Error in creating voxelgrid from triangulated mesh

[fnoi]

Checklist

• [X] I have searched for similar issues.
• [x] For Python issues, I have tested with the latest development wheel.
• [X] I have checked the release documentation and the latest documentation (for master branch).

Describe the issue

I am facing problems with an occupancy grid from a triangulated mesh using the o3d.geometry.VoxelGrid.create_from_triangle_mesh function to create a voxelgrid of occupied voxels. The resulting grid does not include all faces of the mesh. It looks like the voxels are shifted or something, but I couldn't figure out the reason/ logic of what is happening.

Steps to reproduce the bug

import open3d as o3d
import meshio
import plotly.graph_objects as go
import numpy as np

mesh_monki = o3d.io.read_triangle_mesh('untitled.obj')
trimodel_mesh = meshio.read('untitled.obj')
vertices_trimodel = trimodel_mesh.points
triangles_trimodel = trimodel_mesh.cells_dict['triangle']

x=vertices_trimodel[:, 0]
y=vertices_trimodel[:, 1]
z=vertices_trimodel[:, 2]
i=triangles_trimodel[:, 0]
j=triangles_trimodel[:, 1]
k=triangles_trimodel[:, 2]

edges = []
for face in triangles_trimodel:
    for ind in range(len(face)):
        if ind < len(face)-1:
            edges.append(
                (face[ind], face[ind+1])
            )
        else:
            edges.append(
                (face[ind], face[0])
            )

xe=[]
ye=[]
ze=[]

for edge in edges:
    xe.extend([x[edge[0]], x[edge[1]], None])
    ye.extend([y[edge[0]], y[edge[1]], None])
    ze.extend([z[edge[0]], z[edge[1]], None])

trace_monki_mesh = go.Mesh3d(x=x, y=y, z=z, color='white', alphahull=3, #flatshading=True,
                               i=i,
                               j=j,
                               k=k,)

trace_monki_edges = go.Scatter3d(x=xe, y=ye, z=ze, mode='lines', line=dict(color='pink', width=1))

vs = 0.25
vg_monki = o3d.geometry.VoxelGrid.create_from_triangle_mesh(input=mesh_monki, voxel_size=vs)
monki_vox = vg_monki.get_voxels()
monki_vox_idx = [_.grid_index for _ in monki_vox]
monki_vox_ptx = np.asarray([vg_monki.get_voxel_bounding_points(_) for _ in monki_vox_idx])
edge_x = []
edge_y = []
edge_z = []

for box in monki_vox_ptx:
    for pair in [(0,1), (0,2), (0,4), (1,3), (1,5), (2,6), (3,7), (4,5), (4,6), (5,7), (6, 7), (2,3)]:
        edge_x.extend([box[pair[0]][0], box[pair[1]][0], None])
        edge_y.extend([box[pair[0]][1], box[pair[1]][1], None])
        edge_z.extend([box[pair[0]][2], box[pair[1]][2], None])

trace_edges = go.Scatter3d(x=edge_x, y=edge_y, z=edge_z, mode='lines')

data = [trace_monki_mesh, trace_monki_edges, trace_edges]

layout = go.Layout(width=1000,
                   height=1000,
                   showlegend=False,
                   hovermode='closest')

fig = go.Figure(data=data, layout=layout)

fig.update_scenes(camera_projection_type='orthographic')

fig.show()

Error message

No response

Expected behavior

o3d.geometry.VoxelGrid.create_from_triangle_mesh should return a grid of occupied voxels, and all faces of the mesh should lie within the occupied voxels; this is not true in my tests (see attached screenshots of the visualization).

Open3D, Python and System information

- Operating system: macOS 12.5.1
- Python version: Python 3.9.9
- Open3D version: 0.14.1
- System architecture: apple-silicon
- Is this a remote workstation?: no
- How did you install Open3D?: pip
- Compiler version (if built from source): NA

Additional information

[FlorianBertonBrightClue]

I have a similar problem in open3d with the function o3d.geometry.VoxelGrid.create_from_triangle_mesh. The voxel grid I obtained from a mesh is either shifted or it's missing lot's of voxel. Furthermore, if I sampled points on the mesh, I found that most of them are not present in the voxel grid (using the function check_if_included). I also noticed that I have a better results and voxel grid if I sampled points on the mesh surface and then use these points to create a voxel grid.

You can check it by using the example I attached in this comment.

Steps to reproduce the bug

# load mesh
path_file = f"base.obj"
mesh = o3d.io.read_triangle_mesh(path_file)

for nb_sampling in [500,1000]:
    print(f"Number of sampling point {nb_sampling}")
    # get surface point cloud
    pcd = mesh.sample_points_poisson_disk(number_of_points=int(nb_sampling))

    # estimate voxel size
    knn = NearestNeighbors(n_neighbors=2, n_jobs=-1)
    knn.fit(np.asarray(pcd.points))
    dists, _ = knn.kneighbors(np.asarray(pcd.points), n_neighbors=2, return_distance=True)
    voxel_size = 1.5*np.mean(dists[:,1])

    voxel_grid_from_mesh = o3d.geometry.VoxelGrid.create_from_triangle_mesh(mesh,voxel_size=voxel_size)
    voxel_grid_from_point_cloud = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd,voxel_size=voxel_size)

    for text,voxel_grid in zip(["mesh","point cloud"],[voxel_grid_from_mesh,voxel_grid_from_point_cloud]):
        print(f"    Test points on {text} surface in voxel grid from point cloud")
        is_points_in_grid = voxel_grid.check_if_included(o3d.utility.Vector3dVector(pcd.points))
        nb_points, nb_total = (np.sum(is_points_in_grid),len(is_points_in_grid))
        print(f"        number points in voxel grid: {nb_points}/{nb_total} ({np.round(100*nb_points/nb_total,1)}%)")

Ouput

Number of sampling point 500
    Test points on mesh surface in voxel grid from point cloud
        number points in voxel grid: 0/500 (0.0%)
    Test points on point cloud surface in voxel grid from point cloud
        number points in voxel grid: 500/500 (100.0%)
Number of sampling point 1000
    Test points on mesh surface in voxel grid from point cloud
        number points in voxel grid: 35/1000 (3.5%)
    Test points on point cloud surface in voxel grid from point cloud
        number points in voxel grid: 1000/1000 (100.0%)

Graphic display

• Using 500 points to sample

• Using 1000 points to sample

base.zip

[Fjanks]

I have the same problem and tried to find a minimal example. It creates a mesh made of one triangle and uses this mesh to create the voxelgrid.

import numpy as np
import open3d as o3d

z = 0

def create_voxelgrid(z, voxel_size):
    p1 = np.array([0,0,0])
    p2 = np.array([100,0,0])
    p3 = np.array([100,100,z])
    vertices = o3d.utility.Vector3dVector(np.array([p1,p2,p3]))
    triangles = o3d.utility.Vector3iVector([(0,1,2)])
    mesh = o3d.geometry.TriangleMesh(vertices, triangles)
    voxelgrid = o3d.geometry.VoxelGrid.create_from_triangle_mesh(mesh, voxel_size)
    return voxelgrid

vg1 = create_voxelgrid(1, 3)
print("vg1 has {} voxel.".format(len(vg1.get_voxels())))

vg2 = create_voxelgrid(0, 3)
print("vg2 has {} voxel.".format(len(vg2.get_voxels())))

vg3 = create_voxelgrid(1, 2)
print("vg3 has {} voxel.".format(len(vg3.get_voxels())))

Output

vg1 has 68 voxel.
vg2 has 628 voxel.
vg3 has 1478 voxel.

In vg1 the most voxel are missing. Moving the third vertex to z=0 (vg2) or using a smaller voxel size (vg3) avoids the problem.