Labeled volumes for MultiBlock with voxelize_volume and voxelize about pyvista HOT 2 OPEN

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)

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