voxelization to 3D numpy.ndarray

[loisvazquezind]

Hello everybody.

I am new with pyvista. I would like to export the voxelization that pyvista does of a geometry to a numpy 3D array. I was looking for related posts and found the discussions #141 and #4339.

In my case, I modified the code as follows:

from pyvista import examples
import pyvista
import numpy as np

mesh = examples.download_foot_bones()
density = mesh.length / 100

voxels = pyvista.voxelize(
    mesh.triangulate(), density=density, check_surface=False
).connectivity()

x_min, x_max, y_min, y_max, z_min, z_max = mesh.bounds
x = np.arange(x_min, x_max, density)
y = np.arange(y_min, y_max, density)
z = np.arange(z_min, z_max, density)
x, y, z = np.meshgrid(x, y, z)

# # Create unstructured grid from the structured grid
# grid = pyvista.StructuredGrid(x, y, z)
# ugrid = pyvista.UnstructuredGrid(grid)

# get part of the mesh within the mesh's bounding surface.
selection = voxels.select_enclosed_points(mesh.extract_surface(),
                                          tolerance=0.0,
                                          check_surface=False)
mask = selection.point_arrays['SelectedPoints'].view(bool)
mask = mask.reshape(x.shape)
mask = pyvista.wrap(mask).threshold(1)

# pyvista.plot(grid.points, scalars=mask, notebook=False)
p = pyvista.Plotter(notebook=False)
p.add_mesh(mask)
p.show()

However, I always get the same error on the mask = mask.reshape(x.shape):

ValueError: cannot reshape array of size 18973 into shape (30,90,34)

I've been trying different options but I don't think I quite understand this library (sorry). I am using python 3.8 and pyvista v0.43.3.

[dveni]

I've been trying to do this the whole morning and finally got to a working solution. I'm also new to pyvista so maybe the menta model I have is not accurate and I might be wrong, but this is what I think is going on:

The first thing is that you cannot put the array of mask into the grid since they have different dimensions as the error says. The problem is that voxels is just an array of points defining the coordinates of the voxels extracted from the mesh (I think that each point can have additional attributes) . So when you select the points of voxels that are inside the mesh:

selection = voxels.select_enclosed_points(mesh.extract_surface(),
                                          tolerance=0.0,
                                          check_surface=False)

the result is that all the points are inside (as could be expected, since all the points in voxels are extracted from mesh) except for the actual surface (interesting method to get the skin of the volume enclosed by the mesh):

Actually, if you change the tolerance slightly to 0.1, most of the points in selection become true.

Now grid is a 3D space volume (defined from a 3D np meshgrid) within the mesh bounds that you are sampling every density steps. What you want is to set the points of grid =True if they corresponde to the mask, and False otherwise:

# Create unstructured grid from the structured grid
grid = pyvista.StructuredGrid(x, y, z)
ugrid = pyvista.UnstructuredGrid(grid)

# get part of the mesh within the mesh's bounding surface.
selection = ugrid.select_enclosed_points(mesh.extract_surface(),
                                         tolerance=0.0,
                                         check_surface=False)
mask = selection.point_data['SelectedPoints'].view(bool)
pyvista.plot(ugrid.points, scalars=mask.astype(int))

Now you can plot the grid with the scalar values from selection and should look somehting like this since you are plotting the sampled 3D volume:

So what you want is to get the array mask (an array of True/False corresponding to the points in grid) into the shape of the 3D grid and then extract the np.ndarray. I have been very confused with this, because if you do this:

points_C_order = mask.reshape(grid.dimensions)

your result will look something like this (as reported in pyvista/pyvista-support#141 and #5778):

The problem is the reshape order as hinted here: #4504:

points_F_order = mask.reshape(grid.dimensions, order='F')

will result in:

where points_F_order is a pyvista_ndarray that you can get into a np.ndarray with my_np_array = np.array(points_F_order).

The complete code with some additionall plotting:

from pyvista import examples
import pyvista
import numpy as np
import matplotlib.pyplot as plt

def get_array(mesh: pyvista.PolyData, density:float):
    x_min, x_max, y_min, y_max, z_min, z_max = mesh.bounds
    x = np.arange(x_min, x_max, density)
    y = np.arange(y_min, y_max, density)
    z = np.arange(z_min, z_max, density)
    x, y, z = np.meshgrid(x, y, z)

    # Create unstructured grid from the structured grid
    grid = pyvista.StructuredGrid(x, y, z)
    ugrid = pyvista.UnstructuredGrid(grid)

    # get part of the mesh within the mesh's bounding surface.
    selection = ugrid.select_enclosed_points(mesh.extract_surface(),
                                            tolerance=0.0,
                                            check_surface=False)
    mask = selection['SelectedPoints'].view(bool)
    mask = mask.reshape(x.shape, order='F')
    mask = np.array(mask)
    return mask


mesh = examples.download_foot_bones()
mask = get_array(mesh, mesh.length/100)
print(mask.shape, type(mask))


fig, axes = plt.subplots(1, 3)

axes[0].imshow(mask.sum(0))
axes[0].set_title('Sum along X-axis')
axes[1].imshow(mask.sum(1))
axes[1].set_title('Sum along Y-axis')
axes[2].imshow(mask.sum(2))
axes[2].set_title('Sum along Z-axis')

plt.show()

# # Plot mask in 3D with matplotlib
mask = np.transpose(mask, (1,2,0))
xs, ys, zs = np.where(mask)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.voxels(mask)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_box_aspect((np.ptp(xs), np.ptp(ys), np.ptp(zs)))  # aspect ratio is 1:1:1 in data space

plt.show()

[loisvazquezind]

This solves the problem! Thank you very much for your time!