Repository: pyvista/pyvista

Issue:
Create a voxelize feature that takes (n_x, n_y, n_z) as parameters instead of just density
### Describe the feature you would like to be added.

Hi! It would be great if PyVista had a feature that allowed for the voxelization of a mesh according to a specified 3d-grid encompassing the object. The current voxelize feature allows for the adjustment of a density parameter to control the voxel size, but there is no way to directly create for example a 256 by 256 by 256 voxel grid centered on bounding box of a given mesh. For example, implementing something similar to the voxelization capabilities binvox or pyntcloud would be super useful in PyVista. It seems like perhaps this is already possible in vtk, but it's not entirely clear to me as I'm not too familiar with the inner workings of vtk. 

binvox: https://www.patrickmin.com/binvox/
pyntcloud: https://pyntcloud.readthedocs.io/en/latest/structures.html highlight=add_structure#pyntcloud.PyntCloud.add_structure

### Links to VTK Documentation, Examples, or Class Definitions.

https://vtk.org/doc/nightly/html/classvtkVoxelModeller.html



Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pyvista/pyvista/pull/6611?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) Report
Attention: Patch coverage is `90.00000%` with `3 lines` in your changes missing coverage. Please review.
> Project coverage is 96.99%. Comparing base [(`9e01845`)](https://app.codecov.io/gh/pyvista/pyvista/commit/9e0184566eda1b2644c8dc626a98e9d9171e52ec?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) to head [(`1605d46`)](https://app.codecov.io/gh/pyvista/pyvista/commit/1605d46c498d7e705c2e232d075156168f8ef566?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista).
> Report is 1 commits behind head on main.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #6611      +/-   ##
==========================================
- Coverage   97.35%   96.99%   -0.36%     
==========================================
  Files         143      143              
  Lines       27761    27809      +48     
==========================================
- Hits        27026    26973      -53     
- Misses        735      836     +101     
```

</details>
pyvista-bot: <!-- NETLIFY DEPLOY COMMENT GENERATED BY ACTIONS_NETLIFY - APP ID SHA256: d3574f413fe51079937885070831e027ca5b10f7facb9e609494e0a481a31f79 -->
🚀 Deployed on https://66d8b35f473d9a2b32ddf611--pyvista-dev.netlify.app
GuillaumeBroggi: I tried this new method to voxelize a mesh. As mentioned in #1743, it seems superior to the existing approaches as it does not produce stray voxels and is significantly faster.

Instead of raising an error when `reference_volume` is `None`, what about automatically generating a reference volume based on the bounds of the input mesh? I can write a commit if this seems like a good idea. 
tkoyama010: > I tried this new method to voxelize a mesh. As mentioned in #1743, it seems superior to the existing approaches as it does not produce stray voxels and is significantly faster.
> 
> Instead of raising an error when `reference_volume` is `None`, what about automatically generating a reference volume based on the bounds of the input mesh? I can write a commit if this seems like a good idea.

@GuillaumeBroggi Thanks for your thought. You can create Pull Request for this branch if you want :)
user27182: That's a good suggestion for improving the API. The current error is based on the Slicer implementation that was used as an initial reference. Still a few artifacts left over from porting that over. 

The hold up for this PR at the moment is with how images are transformed in PyVista , see #6616. Once resolved I'll work on improving this PR
user27182: @GuillaumeBroggi The filter has been updated to no longer require a reference volume, e.g.:

``` python
from pyvista import examples
import pyvista as pv
poly = examples.load_ant()
labelmap = poly.generate_labelmap()

pl = pv.Plotter()
pl.add_mesh(labelmap.points_to_cells().threshold(0.5), color='blue')
pl.add_mesh(poly, color='lime')
pl.show()
```
Voxelized is blue, input poly data is green:
<img width="624" alt="ant" src="https://github.com/user-attachments/assets/0df60501-38ae-4810-a89f-c77494f6ee24">

Feel free to review the code and suggest changes or improvements. (Also note that it's still a draft and I may also push changes)
user27182: With [d4eedba](https://github.com/pyvista/pyvista/pull/6611/commits/d4eedba1b83587e559cd517d295da341b6df3f8e), the filter now dynamically sets the spacing based on the input cell length. This differs from `voxelize_volume` which uses `mesh.length/100` by default.
 
E.g. Using the same code snippet above, but with `download_bunny_coarse`:
<img width="624" alt="image" src="https://github.com/user-attachments/assets/b5a40060-2f1d-4fd9-828c-87c1e25425f0">
And again but with full-res `download_bunny`:
<img width="624" alt="image" src="https://github.com/user-attachments/assets/d20de66f-d128-4400-8396-b6d279aad3c8">

And once again with the ant. This one is a good test because it has many large cells for the body, and tiny cells for the antennae:
<img width="624" alt="image" src="https://github.com/user-attachments/assets/0085570f-9ea9-4dc2-918a-30fe8ada06ef">



user27182: > @user27182 I reviewed the code, and it looks great, much better than I would have done.

Thanks for the feedback!

> In my implementation, I had chosen to use `**kwargs` in case the user wants to use a custom spacing or custom origin. I think your proposition brings more functionality.

Great, so the current option to specify `spacing` or `dimensions` separately is good? I originally considered adding `origin` explicitly too, but it gets too complicated since it also needs to be rounded for the grid to match the bounds correctly. If desired, users can specify `origin` implicitly using a reference volume instead.

> `spacing_bound` is a bit misleading IMO since the code rounds `dimensions`. See #6647 for a different implementation where `rounding_func` replaces `spacing_bound` and accepts any callable. This would provide more flexibility and is less verbose.

I agree with changing `spacing_bounds`.  #6647 looks good to me.

> Also, there is a floating point error when casting the results to a `RectilinearGrid`, but I am not sure it is possible to do anything about it

It appears the error stems from the casting to rectilinear grid, perhaps this function can be adjusted (in a separate PR) to remove the errror?
https://github.com/pyvista/pyvista/blob/e1b3f1559c1a95add0f9013d48005f79fb9d3905/pyvista/core/grid.py#L845-L849



GuillaumeBroggi: > Great, so the current option to specify `spacing` or `dimensions` separately is good?

Yes, it works nicely like that IMO.

> It appears the error stems from the casting to rectilinear grid, perhaps this function can be adjusted (in a separate PR) to remove the error?

Thanks for pointing it out, I'll try to correct the error.


user27182: I'd say this is ready for review. Docs _might_ still fail from minor RST formatting, but otherwise is feature complete.

Patch:
diff --git a/pyvista/core/_vtk_core.py b/pyvista/core/_vtk_core.py
index 32205387..f46b9ba5 100644
--- a/pyvista/core/_vtk_core.py
+++ b/pyvista/core/_vtk_core.py
@@ -363,6 +363,9 @@ from vtkmodules.vtkFiltersSources import vtkSphereSource
 from vtkmodules.vtkFiltersSources import vtkSuperquadricSource
 from vtkmodules.vtkFiltersSources import vtkTessellatedBoxSource
 from vtkmodules.vtkFiltersStatistics import vtkComputeQuartiles
+
+with contextlib.suppress(ImportError):
+    from vtkmodules.vtkFiltersStatistics import vtkLengthDistribution
 from vtkmodules.vtkFiltersTexture import vtkTextureMapToPlane
 from vtkmodules.vtkFiltersTexture import vtkTextureMapToSphere
 from vtkmodules.vtkFiltersVerdict import vtkCellQuality
@@ -384,6 +387,8 @@ from vtkmodules.vtkImagingGeneral import vtkImageGaussianSmooth
 from vtkmodules.vtkImagingGeneral import vtkImageMedian3D
 from vtkmodules.vtkImagingHybrid import vtkSampleFunction
 from vtkmodules.vtkImagingHybrid import vtkSurfaceReconstructionFilter
+from vtkmodules.vtkImagingStencil import vtkImageStencil
+from vtkmodules.vtkImagingStencil import vtkPolyDataToImageStencil
 from vtkmodules.vtkIOGeometry import vtkHoudiniPolyDataWriter
 from vtkmodules.vtkIOGeometry import vtkIVWriter
 from vtkmodules.vtkIOGeometry import vtkOBJWriter
diff --git a/pyvista/core/filters/image_data.py b/pyvista/core/filters/image_data.py
index 9753a066..c921ac1e 100644
--- a/pyvista/core/filters/image_data.py
+++ b/pyvista/core/filters/image_data.py
@@ -1270,7 +1270,7 @@ class ImageDataFilters(DataSetFilters):
             dims_mask,
         )
         new_image.spacing = self.spacing  # type: ignore[attr-defined]
-        new_image.SetDirectionMatrix(self.GetDirectionMatrix())  # type: ignore[attr-defined]
+        new_image.direction_matrix = self.direction_matrix  # type: ignore[attr-defined]
 
         # Copy field data
         new_image.field_data.update(self.field_data)  # type: ignore[attr-defined]
diff --git a/pyvista/core/filters/poly_data.py b/pyvista/core/filters/poly_data.py
index 67d6e7c6..958204be 100644
--- a/pyvista/core/filters/poly_data.py
+++ b/pyvista/core/filters/poly_data.py
@@ -2,12 +2,15 @@
 
 from __future__ import annotations
 
+from typing import TYPE_CHECKING
+from typing import Callable
 from typing import Sequence
 import warnings
 
 import numpy as np
 
 import pyvista
+from pyvista.core import _validation
 from pyvista.core import _vtk_core as _vtk
 from pyvista.core.errors import MissingDataError
 from pyvista.core.errors import NotAllTrianglesError
@@ -27,6 +30,9 @@ from pyvista.core.utilities.helpers import wrap
 from pyvista.core.utilities.misc import abstract_class
 from pyvista.core.utilities.misc import assert_empty_kwargs
 
+if TYPE_CHECKING:  # pragma: no cover
+    from pyvista.core._typing_core import VectorLike
+
 
 @abstract_class
 class PolyDataFilters(DataSetFilters):
@@ -3892,3 +3898,413 @@ class PolyDataFilters(DataSetFilters):
         alg.SetInputData(self)
         _update_alg(alg, progress_bar, "Generating Protein Ribbons")
         return _get_output(alg)
+
+    def voxelize_binary_mask(
+        self,
+        *,
+        background_value: int = 0,
+        foreground_value: int = 1,
+        reference_volume: pyvista.ImageData | None = None,
+        dimensions: VectorLike[int] | None = None,
+        spacing: float | VectorLike[float] | None = None,
+        rounding_func: Callable[[VectorLike[float]], VectorLike[int]] | None = None,
+        cell_length_percentile: float | None = None,
+        cell_length_sample_size: int | None = None,
+        mesh_length_fraction: float | None = None,
+        progress_bar: bool = False,
+    ):
+        """Voxelize :class:`~pyvista.PolyData` as a binary :class:`~pyvista.ImageData` mask.
+
+        Generate a voxelized representation of a surface. The voxelization can be
+        controlled in several ways.
+
+        #. Specify the output geometry using a ``reference_volume``.
+
+        #. Specify the ``spacing`` explicitly.
+
+        #. Specify the ``dimensions`` explicitly.
+
+        #. Specify the ``cell_length_percentile``. The spacing is estimated from the
+           surface's cells using the specified percentile.
+
+        #. Specify ``mesh_length_fraction``. The spacing is computed as a fraction of
+           the mesh's diagonal length.
+
+        Use ``reference_volume`` for full control of the output mask's geometry. For
+        all other options, the geometry is implicitly defined such that the generated
+        mask fits the bounds of the input surface.
+
+        If no inputs are provided, ``cell_length_percentile=0.1`` (10th percentile) is
+        used by default to estimate the spacing. On systems with VTK < 9.2, the default
+        spacing is set to ``mesh_length_fraction=1/100``.
+
+        .. note::
+            For best results, ensure the input surface is a closed surface. The
+            surface is considered closed if it has zero :attr:`~pyvista.PolyData.n_open_edges`.
+
+        .. note::
+            This filter returns voxels represented as point data, not :attr:`~pyvista.CellType.VOXEL` cells.
+            This differs from :func:`~pyvista.voxelize` and :func:`~pyvista.voxelize_volume`
+            which return meshes with voxel cells. See :ref:`image_representations_example`
+            for examples demonstrating the difference.
+
+        Parameters
+        ----------
+        background_value : int, default: 0
+            Background value of the generated mask.
+
+        foreground_value : int, default: 1
+            Foreground value of the generated mask.
+
+        reference_volume : pyvista.ImageData, optional
+            Volume to use as a reference. The output will have the same ``dimensions``,
+            ``origin``, ``spacing``, and ``direction_matrix`` as the reference.
+
+        dimensions : VectorLike[int], optional
+            Dimensions of the generated mask image. Set this value to control the
+            dimensions explicitly. If unset, the dimensions are defined implicitly
+            through other parameter. See summary and examples for details.
+
+        spacing : VectorLike[float], optional
+            Approximate spacing to use for the generated mask image. Set this value
+            to control the spacing explicitly. If unset, the spacing is defined
+            implicitly through other parameters. See summary and examples for details.
+
+        rounding_func : Callable[VectorLike[float], VectorLike[int]], optional
+            Control how the dimensions are rounded to integers based on the provided or
+            calculated ``spacing``. Should accept a length-3 vector containing the
+            dimension values along the three directions and return a length-3 vector.
+            :func:`numpy.round` is used by default.
+
+            Rounding the dimensions implies rounding the actual spacing.
+
+            Has no effect if ``reference_volume`` or ``dimensions`` are specified.
+
+        cell_length_percentile : float, optional
+            Cell length percentage ``p`` to use for computing the default ``spacing``.
+            Default is ``0.1`` (10th percentile) and must be between ``0`` and ``1``.
+            The ``p``-th percentile is computed from the cumulative distribution function
+            (CDF) of lengths which are representative of the cell length scales present
+            in the input. The CDF is computed by:
+
+            #. Triangulating the input cells.
+            #. Sampling a subset of up to ``cell_length_sample_size`` cells.
+            #. Computing the distance between two random points in each cell.
+            #. Inserting the distance into an ordered set to create the CDF.
+
+            Has no effect if ``dimension`` or ``reference_volume`` are specified.
+
+            .. note::
+                This option is only available for VTK 9.2 or greater.
+
+        cell_length_sample_size : int, optional
+            Number of samples to use for the cumulative distribution function (CDF)
+            when using the ``cell_length_percentile`` option. ``100 000`` samples are
+            used by default.
+
+        mesh_length_fraction : float, optional
+            Fraction of the surface mesh's length to use for computing the default
+            ``spacing``. Set this to any fractional value (e.g. ``1/100``) to enable
+            this option. This is used as an alternative to using ``cell_length_percentile``.
+
+        progress_bar : bool, default: False
+            Display a progress bar to indicate progress.
+
+        Returns
+        -------
+        pyvista.ImageData
+            Generated binary mask with a ``'mask'``  point data array. The data array
+            has dtype :class:`numpy.uint8` if the foreground and background values are
+            unsigned and less than 256.
+
+        See Also
+        --------
+        pyvista.voxelize
+            Similar function that returns a :class:`~pyvista.UnstructuredGrid` of
+            :attr:`~pyvista.CellType.VOXEL` cells.
+
+        pyvista.voxelize_volume
+            Similar function that returns a :class:`~pyvista.RectilinearGrid` with cell data.
+
+        pyvista.ImageDataFilters.contour_labeled
+            Filter that generates surface contours from labeled image data. Can be
+            loosely considered as an inverse of this filter.
+
+        pyvista.ImageDataFilters.points_to_cells
+            Convert voxels represented as points to :attr:`~pyvista.CellType.VOXEL`
+            cells.
+
+        pyvista.ImageData
+            Class used to build custom ``reference_volume``.
+
+        Examples
+        --------
+        Generate a binary mask from a coarse mesh.
+
+        >>> import numpy as np
+        >>> import pyvista as pv
+        >>> from pyvista import examples
+        >>> poly = examples.download_bunny_coarse()
+        >>> mask = poly.voxelize_binary_mask()
+
+        The mask is stored as :class:`~pyvista.ImageData` with point data scalars
+        (zeros for background, ones for foreground).
+
+        >>> mask
+        ImageData (...)
+          N Cells:      7056
+          N Points:     8228
+          X Bounds:     -1.245e-01, 1.731e-01
+          Y Bounds:     -1.135e-01, 1.807e-01
+          Z Bounds:     -1.359e-01, 9.140e-02
+          Dimensions:   22, 22, 17
+          Spacing:      1.417e-02, 1.401e-02, 1.421e-02
+          N Arrays:     1
+
+        >>> np.unique(mask.point_data['mask'])
+        pyvista_ndarray([0, 1], dtype=uint8)
+
+        To visualize it as voxel cells, use :meth:`~pyvista.ImageDataFilters.points_to_cells`,
+        then use :meth:`~pyvista.DataSetFilters.threshold` to extract the foreground.
+
+        We also plot the voxel cells in blue and the input poly data in green for
+        comparison.
+
+        >>> def mask_and_polydata_plotter(mask, poly):
+        ...     voxel_cells = mask.points_to_cells().threshold(0.5)
+        ...
+        ...     plot = pv.Plotter()
+        ...     _ = plot.add_mesh(voxel_cells, color='blue')
+        ...     _ = plot.add_mesh(poly, color='lime')
+        ...     plot.camera_position = 'xy'
+        ...     return plot
+        ...
+
+        >>> plot = mask_and_polydata_plotter(mask, poly)
+        >>> plot.show()
+
+        The spacing of the mask image is automatically adjusted to match the
+        density of the input.
+
+        Repeat the previous example with a finer mesh.
+
+        >>> poly = examples.download_bunny()
+        >>> mask = poly.voxelize_binary_mask()
+        >>> plot = mask_and_polydata_plotter(mask, poly)
+        >>> plot.show()
+
+        Control the spacing manually instead. Here, a very coarse spacing is used.
+
+        >>> mask = poly.voxelize_binary_mask(spacing=(0.01, 0.04, 0.02))
+        >>> plot = mask_and_polydata_plotter(mask, poly)
+        >>> plot.show()
+
+        Note that the spacing is only approximate. Check the mask's actual spacing.
+
+        >>> mask.spacing
+        (0.009731187485158443, 0.03858340159058571, 0.020112216472625732)
+
+        The actual values may be greater or less than the specified values. Use
+        ``rounding_func=np.floor`` to force all values to be greater.
+
+        >>> mask = poly.voxelize_binary_mask(
+        ...     spacing=(0.01, 0.04, 0.02), rounding_func=np.floor
+        ... )
+        >>> mask.spacing
+        (0.01037993331750234, 0.05144453545411428, 0.020112216472625732)
+
+        Set the dimensions instead of the spacing.
+
+        >>> mask = poly.voxelize_binary_mask(dimensions=(10, 20, 30))
+        >>> plot = mask_and_polydata_plotter(mask, poly)
+        >>> plot.show()
+
+        >>> mask.dimensions
+        (10, 20, 30)
+
+        Create a mask using a reference volume. First generate polydata from
+        an existing mask.
+
+        >>> volume = examples.load_frog_tissues()
+        >>> poly = volume.contour_labeled(smoothing=True)
+
+        Now create the mask from the polydata using the volume as a reference.
+
+        >>> mask = poly.voxelize_binary_mask(reference_volume=volume)
+        >>> plot = mask_and_polydata_plotter(mask, poly)
+        >>> plot.show()
+
+        """
+        _validation.check_greater_than(self.n_points, 1, name='n_points')  # type: ignore[attr-defined]
+        _validation.check_greater_than(self.n_cells, 1, name='n_cells')  # type: ignore[attr-defined]
+
+        def _preprocess_polydata(poly_in):
+            return poly_in.compute_normals().triangulate()
+
+        if reference_volume is not None:
+            if (
+                dimensions is not None
+                or spacing is not None
+                or rounding_func is not None
+                or cell_length_percentile is not None
+                or cell_length_sample_size is not None
+                or mesh_length_fraction is not None
+            ):
+                raise TypeError(
+                    'Cannot specify a reference volume with other geometry parameters. `reference_volume` must define the geometry exclusively.'
+                )
+            _validation.check_instance(reference_volume, pyvista.ImageData, name='reference volume')
+            # The image stencil filters do not support orientation, so we apply the
+            # inverse direction matrix to "remove" orientation from the polydata
+            poly_ijk = self.transform(reference_volume.direction_matrix.T, inplace=False)
+            poly_ijk = _preprocess_polydata(poly_ijk)
+        else:
+            # Compute reference volume geometry
+            if spacing is not None and dimensions is not None:
+                raise TypeError("Spacing and dimensions cannot both be set. Set one or the other.")
+
+            # Need to preprocess so that we have a triangle mesh for computing
+            # cell length percentile
+            poly_ijk = _preprocess_polydata(self)
+
+            if spacing is None:
+                if cell_length_percentile is not None and mesh_length_fraction is not None:
+                    raise TypeError(
+                        "Cell length percentile and mesh length fraction cannot both be set. Set one or the other."
+                    )
+
+                less_than_vtk92 = pyvista.vtk_version_info < (9, 2)
+                if (
+                    cell_length_percentile is not None or cell_length_sample_size is not None
+                ) and less_than_vtk92:
+                    raise TypeError(
+                        "Cell length percentile and sample size requires VTK 9.2 or greater."
+                    )
+
+                if mesh_length_fraction is not None or less_than_vtk92:
+                    # Compute spacing from mesh length
+                    mesh_length_fraction = (
+                        1 / 100
+                        if mesh_length_fraction is None
+                        else _validation.validate_number(
+                            mesh_length_fraction, must_have_dtype=float, must_be_in_range=[0.0, 1.0]
+                        )
+                    )
+                    spacing = self.length * mesh_length_fraction  # type: ignore[attr-defined]
+                else:
+                    # Estimate spacing from cell length percentile
+                    cell_length_percentile = (
+                        0.1 if cell_length_percentile is None else cell_length_percentile
+                    )
+                    cell_length_sample_size = (
+                        100_000 if cell_length_sample_size is None else cell_length_sample_size
+                    )
+                    spacing = _length_distribution_percentile(
+                        poly_ijk,
+                        cell_length_percentile,
+                        cell_length_sample_size,
+                        progress_bar=progress_bar,
+                    )
+            else:
+                # Spacing is specified directly. Make sure other params are not set.
+                if cell_length_percentile is not None:
+                    raise TypeError(
+                        "Spacing and cell length percentile cannot both be set. Set one or the other."
+                    )
+                if mesh_length_fraction is not None:
+                    raise TypeError(
+                        "Spacing and mesh length fraction cannot both be set. Set one or the other."
+                    )
+
+            # Get initial spacing (will be adjusted later)
+            initial_spacing = _validation.validate_array3(spacing, broadcast=True)
+
+            # Get size of poly data for computing dimensions
+            bnds = self.bounds  # type: ignore[attr-defined]
+            x_size = bnds.x_max - bnds.x_min
+            y_size = bnds.y_max - bnds.y_min
+            z_size = bnds.z_max - bnds.z_min
+            sizes = np.array((x_size, y_size, z_size))
+
+            if dimensions is None:
+                rounding_func = np.round if rounding_func is None else rounding_func
+                dimensions = np.array(rounding_func(sizes / initial_spacing), dtype=int)
+            elif rounding_func is not None:
+                raise TypeError(
+                    "Rounding func cannot be set when dimensions is specified. Set one or the other."
+                )
+
+            reference_volume = pyvista.ImageData()
+            reference_volume.dimensions = dimensions  # type: ignore[assignment]
+            # Dimensions are now fixed, now adjust spacing to match poly data bounds
+            # Since we are dealing with voxels as points, we want the bounds of the
+            # points to be 1/2 spacing width smaller than the polydata bounds
+            final_spacing = sizes / np.array(reference_volume.dimensions)
+            reference_volume.spacing = final_spacing
+            reference_volume.origin = np.array(self.bounds[::2]) + final_spacing / 2  # type: ignore[attr-defined]
+
+        # Init output structure. The image stencil filters do not support
+        # orientation, so we do not set the direction matrix
+        binary_mask = pyvista.ImageData()
+        binary_mask.dimensions = reference_volume.dimensions
+        binary_mask.spacing = reference_volume.spacing
+        binary_mask.origin = reference_volume.origin
+
+        # Init output scalars. Use uint8 dtype if possible.
+        scalars_shape = (binary_mask.n_points,)
+        scalars_dtype: type[np.uint8 | float | int]
+        if all(
+            isinstance(val, int) and val < 256 and val >= 0
+            for val in (background_value, foreground_value)
+        ):
+            scalars_dtype = np.uint8
+        elif all(
+            isinstance(val, (float, int)) and round(val) == val
+            for val in (background_value, foreground_value)
+        ):
+            scalars_dtype = int
+        else:
+            scalars_dtype = float
+        scalars = (  # Init with background value
+            np.zeros(scalars_shape, dtype=scalars_dtype)
+            if background_value == 0
+            else np.ones(scalars_shape, dtype=scalars_dtype) * background_value
+        )
+        binary_mask['mask'] = scalars  # type: ignore[assignment]
+
+        # Make sure that we have a clean triangle-strip polydata
+        # Note: Poly was partially pre-processed earlier
+        poly_ijk = poly_ijk.strip()
+
+        # Convert polydata to stencil
+        poly_to_stencil = _vtk.vtkPolyDataToImageStencil()
+        poly_to_stencil.SetInputData(poly_ijk)
+        poly_to_stencil.SetOutputSpacing(*reference_volume.spacing)
+        poly_to_stencil.SetOutputOrigin(*reference_volume.origin)
+        poly_to_stencil.SetOutputWholeExtent(*reference_volume.extent)
+        _update_alg(poly_to_stencil, progress_bar, "Converting polydata")
+
+        # Convert stencil to image
+        stencil = _vtk.vtkImageStencil()
+        stencil.SetInputData(binary_mask)
+        stencil.SetStencilConnection(poly_to_stencil.GetOutputPort())
+        stencil.ReverseStencilOn()
+        stencil.SetBackgroundValue(foreground_value)
+        _update_alg(stencil, progress_bar, "Generating binary mask")
+        output_volume = _get_output(stencil)
+
+        # Set the orientation of the output
+        output_volume.direction_matrix = reference_volume.direction_matrix
+
+        return pyvista.wrap(output_volume)
+
+
+def _length_distribution_percentile(poly, percentile, cell_length_sample_size, progress_bar):
+    percentile = _validation.validate_number(
+        percentile, must_be_in_range=[0.0, 1.0], name='percentile'
+    )
+    distribution = _vtk.vtkLengthDistribution()
+    distribution.SetInputData(poly)
+    distribution.SetSampleSize(cell_length_sample_size)
+    _update_alg(distribution, progress_bar, "Computing cell length distribution")
+    return distribution.GetLengthQuantile(percentile)
diff --git a/pyvista/core/utilities/features.py b/pyvista/core/utilities/features.py
index 80316807..0f9c2171 100644
--- a/pyvista/core/utilities/features.py
+++ b/pyvista/core/utilities/features.py
@@ -83,6 +83,14 @@ def voxelize(mesh, density=None, check_surface=True, enclosed=False, fit_bounds=
     Prior to version 0.39.0, this method improperly handled the order of
     structured coordinates.
 
+    See Also
+    --------
+    pyvista.voxelize_volume
+        Similar function that returns a :class:`pyvista.RectilinearGrid` with cell data.
+
+    pyvista.PolyDataFilters.voxelize_binary_mask
+        Similar function that returns a :class:`pyvista.ImageData` with point data.
+
     Examples
     --------
     Create an equal density voxelized mesh.
@@ -246,6 +254,12 @@ def voxelize_volume(mesh, density=None, check_surface=True, enclosed=False, fit_
     See Also
     --------
     pyvista.voxelize
+        Similar function that returns a :class:`pyvista.UnstructuredGrid` of
+        :attr:`~pyvista.CellType.VOXEL` cells.
+
+    pyvista.PolyDataFilters.voxelize_binary_mask
+        Similar function that returns a :class:`pyvista.ImageData` with point data.
+
     pyvista.DataSetFilters.select_enclosed_points
 
     Examples


Test patch:
diff --git a/tests/core/test_polydata_filters.py b/tests/core/test_polydata_filters.py
index e6e81045..5e486bcd 100644
--- a/tests/core/test_polydata_filters.py
+++ b/tests/core/test_polydata_filters.py
@@ -1,10 +1,12 @@
 from __future__ import annotations
 
+import numpy as np
 import pytest
 
 import pyvista as pv
 from pyvista import examples
 from pyvista.core.errors import MissingDataError
+from tests.conftest import flaky_test
 
 
 def test_contour_banded_raise(sphere):
@@ -74,3 +76,189 @@ def test_protein_ribbon():
     tgqp = examples.download_3gqp()
     ribbon = tgqp.protein_ribbon()
     assert ribbon.n_cells
+
+
+@pytest.fixture
+def frog_tissues_image():
+    return examples.load_frog_tissues()
+
+
+@pytest.fixture
+def frog_tissues_contour(frog_tissues_image):
+    return frog_tissues_image.contour_labeled(smoothing=False)
+
+
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_voxelize_binary_mask(frog_tissues_image, frog_tissues_contour):
+    mask = frog_tissues_contour.voxelize_binary_mask(
+        reference_volume=frog_tissues_image, progress_bar=True
+    )
+
+    expected_voxels = frog_tissues_image.points_to_cells().threshold(0.5)
+    actual_voxels = mask.points_to_cells().threshold(0.5)
+
+    assert expected_voxels.bounds == actual_voxels.bounds
+    assert expected_voxels.n_cells == actual_voxels.n_cells
+
+
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_voxelize_binary_mask_no_reference(frog_tissues_image, frog_tissues_contour):
+    mask = frog_tissues_contour.voxelize_binary_mask()
+    assert np.allclose(mask.points_to_cells().bounds, frog_tissues_contour.bounds)
+
+
+def test_voxelize_binary_mask_dimensions(sphere):
+    dims = (10, 11, 12)
+    mask = sphere.voxelize_binary_mask(dimensions=dims)
+    assert np.allclose(mask.points_to_cells().bounds, sphere.bounds)
+    assert mask.dimensions == dims
+
+
+def test_voxelize_binary_mask_auto_spacing(ant):
+    # Test default
+    mask_no_input = ant.voxelize_binary_mask()
+    if pv.vtk_version_info < (9, 2):
+        expected_mask = ant.voxelize_binary_mask(mesh_length_fraction=1 / 100)
+    else:
+        expected_mask = ant.voxelize_binary_mask(cell_length_percentile=0.1)
+    assert mask_no_input.spacing == expected_mask.spacing
+
+    # Test cell length
+    if pv.vtk_version_info < (9, 2):
+        match = "Cell length percentile and sample size requires VTK 9.2 or greater."
+        with pytest.raises(TypeError, match=match):
+            ant.voxelize_binary_mask(cell_length_percentile=0.2)
+    else:
+        mask_percentile_20 = ant.voxelize_binary_mask(cell_length_percentile=0.2)
+        mask_percentile_50 = ant.voxelize_binary_mask(cell_length_percentile=0.5)
+        assert np.all(np.array(mask_percentile_20.spacing) < mask_percentile_50.spacing)
+
+    # Test mesh length
+    mask_fraction_200 = ant.voxelize_binary_mask(mesh_length_fraction=1 / 200)
+    mask_fraction_500 = ant.voxelize_binary_mask(mesh_length_fraction=1 / 500)
+    assert np.all(np.array(mask_fraction_200.spacing) > mask_fraction_500.spacing)
+    # Check spacing matches mesh length. Use atol since spacing is approximate.
+    assert np.allclose(mask_fraction_500.spacing, ant.length / 500, atol=1e-3)
+
+
+# This test is flaky because of random sampling that cannot be controlled.
+# Sometimes the sampling produces the same output.
+@flaky_test
+def test_voxelize_binary_mask_cell_length_sample_size(ant):
+    if pv.vtk_version_info < (9, 2):
+        match = "Cell length percentile and sample size requires VTK 9.2 or greater."
+        with pytest.raises(TypeError, match=match):
+            ant.voxelize_binary_mask(cell_length_percentile=0.2)
+    else:
+        mask_samples_1 = ant.voxelize_binary_mask(cell_length_sample_size=100)
+        mask_samples_2 = ant.voxelize_binary_mask(cell_length_sample_size=200)
+        assert mask_samples_1.spacing != mask_samples_2.spacing
+
+        mask_samples_1 = ant.voxelize_binary_mask(cell_length_sample_size=ant.n_cells)
+        mask_samples_2 = ant.voxelize_binary_mask(cell_length_sample_size=ant.n_cells)
+        assert mask_samples_1.spacing == mask_samples_2.spacing
+
+
+@pytest.mark.parametrize(
+    'rounding_func',
+    [np.round, np.ceil, np.floor, lambda x: [np.round(x[0]), np.ceil(x[1]), np.floor(x[2])]],
+)
+def test_voxelize_binary_mask_rounding_func(sphere, rounding_func):
+    spacing = np.array((1.1, 1.2, 1.3))
+    mask = sphere.voxelize_binary_mask(spacing=spacing, rounding_func=rounding_func)
+    assert np.allclose(mask.points_to_cells().bounds, sphere.bounds)
+    if rounding_func == np.round:
+        assert np.any(mask.spacing > spacing)
+        assert np.any(mask.spacing < spacing)
+    elif rounding_func == np.ceil:
+        assert np.all(mask.spacing < spacing)
+    elif rounding_func == np.floor:
+        assert np.all(mask.spacing > spacing)
+    else:  # rounding_func == lambda x: [np.round(x[0]), np.ceil(x[1]), np.floor(x[2])]]
+        assert mask.spacing[1] < spacing[1]
+        assert mask.spacing[2] > spacing[2]
+
+
+@pytest.mark.parametrize('foreground', [1, 2.1])
+@pytest.mark.parametrize('background', [-1, 0])
+def test_voxelize_binary_mask_foreground_background(sphere, foreground, background):
+    mask = sphere.voxelize_binary_mask(foreground_value=foreground, background_value=background)
+    unique, counts = np.unique(mask['mask'], return_counts=True)
+    assert np.array_equal(unique, [background, foreground])
+    # Test we have more foreground than background (not always true, but is true for a sphere mesh)
+    assert counts[1] > counts[0]
+
+    # Test dtype
+    if (
+        isinstance(foreground, int)
+        and isinstance(background, int)
+        and foreground >= 0
+        and background >= 0
+    ):
+        assert mask['mask'].dtype == np.uint8
+    elif isinstance(foreground, int) and isinstance(background, int):
+        assert mask['mask'].dtype == int
+    else:
+        assert mask['mask'].dtype == float
+
+
+@pytest.fixture
+def oriented_image():
+    image = pv.ImageData()
+    image.spacing = (1.1, 1.2, 1.3)
+    image.dimensions = (10, 11, 12)
+    image.direction_matrix = pv.Transform().rotate_vector((4, 5, 6), 30).matrix[:3, :3]
+    image['scalars'] = np.ones((image.n_points,))
+    return image
+
+
+@pytest.fixture
+def oriented_polydata(oriented_image):
+    oriented_poly = oriented_image.pad_image().contour_labeled(smoothing=False)
+    assert np.allclose(oriented_poly.bounds, oriented_image.points_to_cells().bounds, atol=0.1)
+    return oriented_poly
+
+
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_voxelize_binary_mask_orientation(oriented_image, oriented_polydata):
+    mask = oriented_polydata.voxelize_binary_mask(reference_volume=oriented_image)
+    assert mask.bounds == oriented_image.bounds
+    mask_as_surface = mask.pad_image().contour_labeled(smoothing=False)
+    assert mask_as_surface.bounds == oriented_polydata.bounds
+
+
+def test_voxelize_binary_mask_raises(sphere):
+    match = 'Spacing and dimensions cannot both be set. Set one or the other.'
+    with pytest.raises(TypeError, match=match):
+        sphere.voxelize_binary_mask(dimensions=(1, 2, 3), spacing=(4, 5, 6))
+
+    match = 'Spacing and cell length percentile cannot both be set. Set one or the other.'
+    with pytest.raises(TypeError, match=match):
+        sphere.voxelize_binary_mask(spacing=(4, 5, 6), cell_length_percentile=0.2)
+
+    match = 'Spacing and mesh length fraction cannot both be set. Set one or the other.'
+    with pytest.raises(TypeError, match=match):
+        sphere.voxelize_binary_mask(spacing=(4, 5, 6), mesh_length_fraction=0.2)
+
+    match = (
+        'Cell length percentile and mesh length fraction cannot both be set. Set one or the other.'
+    )
+    with pytest.raises(TypeError, match=match):
+        sphere.voxelize_binary_mask(mesh_length_fraction=1 / 100, cell_length_percentile=0.2)
+
+    match = 'Rounding func cannot be set when dimensions is specified. Set one or the other.'
+    with pytest.raises(TypeError, match=match):
+        sphere.voxelize_binary_mask(dimensions=(1, 2, 3), rounding_func=np.round)
+
+    for parameter in [
+        'dimensions',
+        'spacing',
+        'rounding_func',
+        'cell_length_percentile',
+        'cell_length_sample_size',
+        'mesh_length_fraction',
+    ]:
+        kwargs = {parameter: 0}  # Give parameter any value for test
+        match = 'Cannot specify a reference volume with other geometry parameters. `reference_volume` must define the geometry exclusively.'
+        with pytest.raises(TypeError, match=match):
+            sphere.voxelize_binary_mask(reference_volume=pv.ImageData(), **kwargs)