Labeled volumes for MultiBlock with voxelize_volume and voxelize

[MosGeo]

Describe the feature you would like to be added.

It would be great to have the ability to voxelize_volume and voxelize to create labeled volume rather than a binary volume. I see this being used in conjection with MultiBlock. I want each block to have its own number.

If there is a workaround to do this right now, that would be very helpful.

Links to VTK Documentation, Examples, or Class Definitions.

No response

Pseudocode or Screenshots

I see the api as something like this:

mesh = pv.MultiBlock(...)
labeled_volume = pv.voxelize_volume(mesh, label=True)
binary_volume = pv.voxelize_volume(mesh, label=False)

Obviously, the default for label parameter is False.

[user27182]

Great idea! Here's a possible way to implement this for voxelize using the sample and split_values filters. See code/images below.

A few notes about this method:

• sample seems to only operate on point data, not sure how to make it work directly with cell data. Using ptc seems is a bit hacky because it means the data will be interpolated twice.
• the tolerance for sample must be sufficiently large, (maybe set it to inf? or just really large)
• split_values is a new filter (pyvista==0.44.0), currently only available from the main (dev) branch.

import pyvista as pv
from pyvista import examples
import numpy as np

# 1. Load mesh and add cell data labels
mesh = examples.download_bunny_coarse().clean()
labels = np.zeros((mesh.n_cells,))
labels[:700] = 1
mesh['labels'] = labels
mesh.plot()

# 2. Voxelize and sample
vox = pv.voxelize(mesh, density=0.01)
sampled = vox.sample(target=mesh, tolerance=1)
sampled.set_active_scalars('labels')
sampled.plot(show_edges=True)

# 3. We want cell data not point data
sampled = sampled.ptc()
sampled['labels'] = np.round(sampled['labels'])
sampled.plot(show_edges=True)

# 4. Split into multiblock
multiblock = sampled.split_values()

multiblock[0].plot()
# Also multiblock[1].plot(show_edges=True)

Output:

1. Inital labeled mesh:

2. Sampled voxelized mesh. The cell data is interpolated as point data (bad), but we want cell data

3. After fixing cell data

4. First block of the split voxelized labels

[MosGeo]

Thanks. I will check it out. I am fairly new to pyvista in terms of experinece so I don't know lots of capaiblities (e.g., sample function.

I also gave it a ago by modifying the voxelize_volume function. I was able to produce the desired result. I used a for loop so I don't to go over the blocks in the mesh so I don't know if it is the most efficient way of doing it. Here is an example of the result:

import collections
import numpy as np
import pyvista


def voxelize_volume_multiblock_label(mesh, density=None, check_surface=True):
    """Voxelize mesh to create a RectilinearGrid voxel volume.

    Creates a voxel volume that encloses the input mesh and discretizes the cells
    within the volume that intersect or are contained within the input mesh.
    ``InsideMesh``, an array in ``cell_data``, is ``1`` for cells inside and ``0`` outside.

    Parameters
    ----------
    mesh : pyvista.DataSet
        Mesh to voxelize.

    density : float | array_like[float]
        The uniform size of the voxels when single float passed.
        Nonuniform voxel size if a list of values are passed along x,y,z directions.
        Defaults to 1/100th of the mesh length.

    check_surface : bool, default: True
        Specify whether to check the surface for closure. If on, then the
        algorithm first checks to see if the surface is closed and
        manifold. If the surface is not closed and manifold, a runtime
        error is raised.

    Returns
    -------
    pyvista.RectilinearGrid
        RectilinearGrid as voxelized volume with discretized cells.

    See Also
    --------
    pyvista.voxelize
    pyvista.DataSetFilters.select_enclosed_points

    Examples
    --------
    Create an equal density voxel volume from input mesh.

    >>> import pyvista as pv
    >>> import numpy as np

    Load file from PyVista examples.

    >>> from pyvista import examples
    >>> mesh = examples.download_cow()

    Create an equal density voxel volume and plot the result.

    >>> vox = pv.voxelize_volume(mesh, density=0.15)
    >>> cpos = [(15, 3, 15), (0, 0, 0), (0, 0, 0)]
    >>> vox.plot(scalars='InsideMesh', show_edges=True, cpos=cpos)

    Slice the voxel volume to view ``InsideMesh``.

    >>> slices = vox.slice_orthogonal()
    >>> slices.plot(scalars='InsideMesh', show_edges=True)

    Create a voxel volume from unequal density dimensions and plot result.

    >>> vox = pv.voxelize_volume(mesh, density=[0.15, 0.15, 0.5])
    >>> vox.plot(scalars='InsideMesh', show_edges=True, cpos=cpos)

    Slice the unequal density voxel volume to view ``InsideMesh``.

    >>> slices = vox.slice_orthogonal()
    >>> slices.plot(scalars='InsideMesh', show_edges=True, cpos=cpos)

    """
    mesh = pyvista.wrap(mesh)
    if density is None:
        density = mesh.length / 100
    if isinstance(density, (int, float, np.number)):
        density_x, density_y, density_z = [density] * 3
    elif isinstance(density, (collections.abc.Sequence, np.ndarray)):
        density_x, density_y, density_z = density
    else:
        raise TypeError(f"Invalid density {density!r}, expected number or array-like.")

    x_min, x_max, y_min, y_max, z_min, z_max = mesh.bounds
    x = np.arange(x_min, x_max, density_x)
    y = np.arange(y_min, y_max, density_y)
    z = np.arange(z_min, z_max, density_z)

    # Initialize the voi
    voi = pyvista.RectilinearGrid(x, y, z)
    voi["InsideMesh"] = np.zeros(voi.n_cells)

    # check and pre-process input mesh
    sub_mesh: pyvista.PolyData
    for i, sub_mesh in enumerate(mesh):
        sub_mesh_number = i + 1
        surface: pyvista.PolyData = sub_mesh.extract_geometry()

        # filter preserves topology
        if not surface.faces.size:
            # we have a point cloud or an empty mesh
            raise ValueError("Input mesh must have faces for voxelization.")
        if not surface.is_all_triangles:
            # reduce chance for artifacts, see gh-1743
            surface.triangulate(inplace=True)

        # get part of the mesh within the mesh's bounding surface.
        selection = voi.select_enclosed_points(
            surface, tolerance=0.0, check_surface=check_surface
        )
        mask_vol = selection.point_data["SelectedPoints"].view(np.bool_)

        # Get voxels that fall within input mesh boundaries
        cell_ids = np.unique(
            voi.extract_points(np.argwhere(mask_vol))["vtkOriginalCellIds"]
        )

        # Create new element of grid where all cells _within_ mesh boundary are
        # given new name 'MeshCells' and a discrete value of 1
        voi["InsideMesh"][cell_ids] = sub_mesh_number

    return voi