Repository: pyvista/pyvista

Issue:
Add vtkImageResample  filter for upscaling
## Description

How would we go about properly upscaling a mesh much like the image from [this SO post](https://stackoverflow.com/questions/40988864/creating-a-grid-using-vtk-activiz-c-sharp?rq=1):

![tuVWd](https://user-images.githubusercontent.com/22067021/61654710-44193100-ac7a-11e9-8798-68373728fedd.jpg)


## Example Data

I started this using the `sample` filter, but this isn't exactly what I'm looking for as the sample filter can only sample onto the point data of the output grid. I'd like a way to properly upscale a uniform/ImageData grid on the cell data. Also, the coarse grid that I created does not perfectly align with the fine grid which we should fix.

```py
import pyvista as pv
from pyvista import examples

fine = examples.load_channels()
# Perfrom an upscaling of sorts
coarse = pv.create_grid(fine, dimensions=(50,50,45)).sample(fine)

# Show it
dp = dict(cmap='Reds', show_edges=True, lighting=False)

p = pv.Plotter(shape=(1,2), notebook=False)
p.add_mesh(fine, **dp)

p.subplot(0,1)
p.add_mesh(coarse, **dp)

p.link_views()
p.show()
```

<img width="1574" alt="Screen Shot 2019-07-22 at 12 10 47 PM" src="https://user-images.githubusercontent.com/22067021/61654618-0e744800-ac7a-11e9-91c3-15425b71f220.png">

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pyvista/pyvista/pull/7068?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 `62.06897%` with `44 lines` in your changes missing coverage. Please review.
> Project coverage is 91.35%. Comparing base [(`2921752`)](https://app.codecov.io/gh/pyvista/pyvista/commit/29217527ea90dd26fba3dd30859380624319e1f4?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) to head [(`f026ebf`)](https://app.codecov.io/gh/pyvista/pyvista/commit/f026ebffec5e954031e5b407c801ef0ee73810d3?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista).

> :exclamation:  There is a different number of reports uploaded between BASE (2921752) and HEAD (f026ebf). Click for more details.
> 
> <details><summary>HEAD has 5 uploads less than BASE</summary>
>
>| Flag | BASE (2921752) | HEAD (f026ebf) |
>|------|------|------|
>||6|1|
></details>

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


```diff
@@            Coverage Diff             @@
##             main    #7068      +/-   ##
==========================================
- Coverage   97.53%   91.35%   -6.18%     
==========================================
  Files         144      144              
  Lines       28848    28962     +114     
==========================================
- Hits        28137    26459    -1678     
- Misses        711     2503    +1792     
```

</details>
pyvista-bot: <!-- NETLIFY DEPLOY COMMENT GENERATED BY ACTIONS_NETLIFY - APP ID SHA256: d3574f413fe51079937885070831e027ca5b10f7facb9e609494e0a481a31f79 -->
🚀 Deployed on https://677e219c08aa6d1250ab1734--pyvista-dev.netlify.app

Patch:
diff --git a/pyvista/core/_vtk_core.py b/pyvista/core/_vtk_core.py
index e99d65c8..bfb4633f 100644
--- a/pyvista/core/_vtk_core.py
+++ b/pyvista/core/_vtk_core.py
@@ -451,6 +451,7 @@ from vtkmodules.vtkImagingCore import vtkImageDifference as vtkImageDifference
 from vtkmodules.vtkImagingCore import vtkImageExtractComponents as vtkImageExtractComponents
 from vtkmodules.vtkImagingCore import vtkImageFlip as vtkImageFlip
 from vtkmodules.vtkImagingCore import vtkImageMirrorPad as vtkImageMirrorPad
+from vtkmodules.vtkImagingCore import vtkImageResample as vtkImageResample
 from vtkmodules.vtkImagingCore import vtkImageThreshold as vtkImageThreshold
 from vtkmodules.vtkImagingCore import vtkImageWrapPad as vtkImageWrapPad
 from vtkmodules.vtkImagingCore import vtkRTAnalyticSource as vtkRTAnalyticSource
diff --git a/pyvista/core/filters/data_set.py b/pyvista/core/filters/data_set.py
index 7f174c82..256c905b 100644
--- a/pyvista/core/filters/data_set.py
+++ b/pyvista/core/filters/data_set.py
@@ -3966,6 +3966,10 @@ class DataSetFilters:
         See Also
         --------
         pyvista.DataSetFilters.interpolate
+            Interpolate values from one mesh onto another.
+
+        :meth:`pyvista.ImageDataFilters.resample`
+            Resample image data to modify its dimensions and spacing.
 
         Examples
         --------
@@ -4110,6 +4114,10 @@ class DataSetFilters:
         See Also
         --------
         pyvista.DataSetFilters.sample
+            Resample array data from one mesh onto another.
+
+        :meth:`pyvista.ImageDataFilters.resample`
+            Resample image data to modify its dimensions and spacing.
 
         Examples
         --------
diff --git a/pyvista/core/filters/image_data.py b/pyvista/core/filters/image_data.py
index 84bf9d86..e6c7a5e6 100644
--- a/pyvista/core/filters/image_data.py
+++ b/pyvista/core/filters/image_data.py
@@ -2927,3 +2927,517 @@ class ImageDataFilters(DataSetFilters):
             )
 
         return dimensions_mask, dimensions_result
+
+    def resample(  # type: ignore[misc]
+        self: ImageData,
+        sample_rate: float | VectorLike[float] | None = None,
+        interpolation: Literal['nearest', 'linear', 'cubic'] = 'nearest',
+        *,
+        reference_image: ImageData | None = None,
+        dimensions: VectorLike[int] | None = None,
+        extend_border: bool | None = None,
+        scalars: str | None = None,
+        preference: Literal['point', 'cell'] = 'point',
+        inplace: bool = False,
+        progress_bar: bool = False,
+    ):
+        """Resample the image to modify its dimensions and spacing.
+
+        The resampling can be controlled in several ways:
+
+        #. Specify the output geometry using a ``reference_image``.
+
+        #. Specify the ``dimensions`` explicitly.
+
+        #. Specify the ``sample_rate`` explicitly.
+
+        Use ``reference_image`` for full control of the resampled geometry. For
+        all other options, the geometry is implicitly defined such that the resampled
+        image fits the bounds of the input.
+
+        This filter may be used to resample either point or cell data. For point data,
+        this filter assumes the data is from discrete samples in space which represent
+        pixels or voxels; the resampled bounds are therefore extended by 1/2 voxel
+        spacing by default though this may be disabled.
+
+        .. note::
+
+            Singleton dimensions are not resampled by this filter, e.g. 2D images
+            will remain 2D.
+
+        .. versionadded:: 0.45
+
+        Parameters
+        ----------
+        sample_rate : float | VectorLike[float], optional
+            Sampling rate(s) to use. Can be a single value or vector of three values
+            for each axis. Values greater than ``1.0`` will up-sample the axis and
+            values less than ``1.0`` will down-sample it. Values must be greater than ``0``.
+
+        interpolation : 'nearest' | 'linear' | 'cubic', default: 'nearest'
+            Interpolation mode to use. By default, ``'nearest'`` is used which
+            duplicates (if upsampling) or removes (if downsampling) values but
+            does not modify them. The ``'linear'`` and ``'cubic'`` options use linear
+            and cubic interpolation, respectively, and may modify the values.
+
+            .. note:
+
+                Integer scalars are converted to floats if ``'linear'`` or
+                ``'cubic'`` interpolation is used.
+
+        reference_image : ImageData, optional
+            Reference image to use. If specified, the input is resampled
+            to match the geometry of the reference. The :attr:`~pyvista.ImageData.dimensions`,
+            :attr:`~pyvista.ImageData.spacing`, :attr:`~pyvista.ImageData.origin`,
+            :attr:`~pyvista.ImageData.offset`, and :attr:`~pyvista.ImageData.direction_matrix`
+            of the resampled image will all match the reference image.
+
+        dimensions : VectorLike[int]
+            Set the output :attr:`~pyvista.ImageData.dimensions` of the resampled image.
+
+            .. note::
+
+                Dimensions is the number of `points` along each axis. If resampling
+                `cell` data, each dimension should be one more than the number of
+                desired output cells (since there are ``N`` cells and ``N+1`` points
+                along each axis). See examples.
+
+        extend_border : bool, optional
+            Extend the apparent input border by approximately half the
+            :attr:`~pyvista.ImageData.spacing`. If enabled, the bounds of the
+            resampled points will be larger than the input image bounds.
+            Enabling this option also has the effect that the re-sampled spacing
+            will directly correlate with the resampled dimensions, e.g. if
+            the dimensions are doubled the spacing will be halved. See examples.
+
+            This option is enabled by default when resampling point data. Has no effect
+            when resampling cell data or when a ``reference_image`` is provided.
+
+        scalars : str, optional
+            Name of scalars to resample. Defaults to currently active scalars.
+
+        preference : str, default: 'point'
+            When scalars is specified, this is the preferred array type to search
+            for in the dataset.  Must be either ``'point'`` or ``'cell'``.
+
+        inplace : bool, default: False
+            If ``True``, resample the image inplace. By default, a new
+            :class:`~pyvista.ImageData` instance is returned.
+
+        progress_bar : bool, default: False
+            Display a progress bar to indicate progress.
+
+        Returns
+        -------
+        ImageData
+            Resampled image.
+
+        See Also
+        --------
+        :meth:`~pyvista.DataSetFilters.sample`
+            Resample array data from one mesh onto another.
+
+        :meth:`~pyvista.DataSetFilters.interpolate`
+            Interpolate values from one mesh onto another.
+
+        :ref:`image_representations_example`
+            Compare images represented as points vs. cells.
+
+        Examples
+        --------
+        Create a small 2D grayscale image with dimensions ``3 x 2`` for demonstration.
+
+        >>> import pyvista as pv
+        >>> import numpy as np
+        >>> from pyvista import examples
+        >>> image = pv.ImageData(dimensions=(3, 2, 1))
+        >>> image.point_data['data'] = np.linspace(0, 255, 6, dtype=np.uint8)
+
+        Define a custom plotter to show the image. Although the image data is defined
+        as point data, we use uses :meth:`points_to_cells` to display the image as
+        :attr:`~pyvista.CellType.PIXEL` (or :attr:`~pyvista.CellType.VOXEL`) cells
+        instead. Grayscale coloring is used and the camera is adjusted to fit the image.
+
+        >>> def image_plotter(image: pv.ImageData) -> pv.Plotter:
+        ...     pl = pv.Plotter()
+        ...     image = image.points_to_cells()
+        ...     pl.add_mesh(
+        ...         image,
+        ...         lighting=False,
+        ...         show_edges=True,
+        ...         cmap='grey',
+        ...         show_scalar_bar=False,
+        ...     )
+        ...     pl.view_xy()
+        ...     pl.camera.tight()
+        ...     return pl
+
+        Show the image.
+
+        >>> plot = image_plotter(image)
+        >>> plot.show()
+
+        Use ``sample_rate`` to up-sample the image. ``'nearest'`` interpolation is
+        used by default.
+
+        >>> upsampled = image.resample(sample_rate=2)
+        >>> plot = image_plotter(upsampled)
+        >>> plot.show()
+
+        Use ``'linear'`` interpolation instead.
+
+        >>> upsampled = image.resample(sample_rate=2, interpolation='linear')
+        >>> plot = image_plotter(upsampled)
+        >>> plot.show()
+
+        Use ``'cubic'`` interpolation instead. Here we also specify the output
+        ``dimensions`` explicitly instead of using ``sample_rate``.
+
+        >>> upsampled = image.resample(dimensions=(6, 4, 1), interpolation='cubic')
+        >>> plot = image_plotter(upsampled)
+        >>> plot.show()
+
+        Compare the relative physical size of the image before and after resampling.
+
+        >>> image
+        ImageData (...)
+          N Cells:      2
+          N Points:     6
+          X Bounds:     0.000e+00, 2.000e+00
+          Y Bounds:     0.000e+00, 1.000e+00
+          Z Bounds:     0.000e+00, 0.000e+00
+          Dimensions:   3, 2, 1
+          Spacing:      1.000e+00, 1.000e+00, 1.000e+00
+          N Arrays:     1
+
+        >>> upsampled
+        ImageData (...)
+          N Cells:      15
+          N Points:     24
+          X Bounds:     -2.500e-01, 2.250e+00
+          Y Bounds:     -2.500e-01, 1.250e+00
+          Z Bounds:     0.000e+00, 0.000e+00
+          Dimensions:   6, 4, 1
+          Spacing:      5.000e-01, 5.000e-01, 1.000e+00
+          N Arrays:     1
+
+        Note that the upsampled :attr:`~pyvista.ImageData.dimensions` are doubled and
+        the :attr:`~pyvista.ImageData.spacing` is halved (as expected). Also note,
+        however, that the physical bounds of the input differ from the output.
+        The upsampled :attr:`~pyvista.ImageData.origin` also differs:
+
+        >>> image.origin
+        (0.0, 0.0, 0.0)
+        >>> upsampled.origin
+        (-0.25, -0.25, 0.0)
+
+        This is because the resampling is done with ``extend_border`` enabled by default
+        which adds a half cell-width border to the image and adjusts the origin and
+        spacing such that the bounds match when the image is represented as cells.
+
+        Apply :meth:`points_to_cells` to the input and resampled images and show that
+        the bounds match.
+
+        >>> image_as_cells = image.points_to_cells()
+        >>> image_as_cells.bounds
+        BoundsTuple(x_min=-0.5, x_max=2.5, y_min=-0.5, y_max=1.5, z_min=0.0, z_max=0.0)
+
+        >>> upsampled_as_cells = upsampled.points_to_cells()
+        >>> upsampled_as_cells.bounds
+        BoundsTuple(x_min=-0.5, x_max=2.5, y_min=-0.5, y_max=1.5, z_min=0.0, z_max=0.0)
+
+        Plot the two images together as wireframe to visualize them. The original is in
+        red, and the resampled image is in black.
+
+        >>> plt = pv.Plotter()
+        >>> _ = plt.add_mesh(
+        ...     image_as_cells, style='wireframe', color='red', line_width=10
+        ... )
+        >>> _ = plt.add_mesh(
+        ...     upsampled_as_cells, style='wireframe', color='black', line_width=2
+        ... )
+        >>> plt.view_xy()
+        >>> plt.camera.tight()
+        >>> plt.show()
+
+        Disable ``extend_border`` to force the input and output bounds of the points
+        to be the same instead.
+
+        >>> upsampled = image.resample(sample_rate=2, extend_border=False)
+
+        Compare the two images again.
+
+        >>> image
+        ImageData (...)
+          N Cells:      2
+          N Points:     6
+          X Bounds:     0.000e+00, 2.000e+00
+          Y Bounds:     0.000e+00, 1.000e+00
+          Z Bounds:     0.000e+00, 0.000e+00
+          Dimensions:   3, 2, 1
+          Spacing:      1.000e+00, 1.000e+00, 1.000e+00
+          N Arrays:     1
+
+        >>> upsampled
+        ImageData (...)
+          N Cells:      15
+          N Points:     24
+          X Bounds:     0.000e+00, 2.000e+00
+          Y Bounds:     0.000e+00, 1.000e+00
+          Z Bounds:     0.000e+00, 0.000e+00
+          Dimensions:   6, 4, 1
+          Spacing:      4.000e-01, 3.333e-01, 1.000e+00
+          N Arrays:     1
+
+        This time the input and output bounds match without any further processing.
+        Like before, the dimensions have doubled; unlike before, however, the spacing is
+        not halved, but is instead smaller than half which is necessaru to ensure the
+        bounds remain the same. Also unlike before, the origin is unaffected:
+
+        >>> image.origin
+        (0.0, 0.0, 0.0)
+        >>> upsampled.origin
+        (0.0, 0.0, 0.0)
+
+        All the above examples are with 2D images with point data. However, the filter
+        also works with 3D volumes and will also work with cell data.
+
+        Convert the 2D image with point data into a 3D volume with cell data and plot
+        it for context.
+
+        >>> volume = image.points_to_cells(dimensionality='3D')
+        >>> volume.plot(show_edges=True, cmap='grey')
+
+        Up-sample the volume. Set the sampling rate for each axis separately.
+
+        >>> resampled = volume.resample(sample_rate=(3.0, 2.0, 1.0))
+        >>> resampled.plot(show_edges=True, cmap='grey')
+
+        Alternatively, we could have set the dimensions explicitly. Since we want
+        ``9 x 4 x 1`` cells along the x-y-z axes (respectively), we set the dimensions
+        to ``(10, 5, 2)``, i.e. one more than the desired number of cells.
+
+        >>> resampled = volume.resample(dimensions=(10, 5, 2))
+        >>> resampled.plot(show_edges=True, cmap='grey')
+
+        Compare the bounds before and after resampling. Unlike with point data, the
+        bounds are not (and cannot be) extended.
+
+        >>> volume.bounds
+        BoundsTuple(x_min=-0.5, x_max=2.5, y_min=-0.5, y_max=1.5, z_min=-0.5, z_max=0.5)
+        >>> resampled.bounds
+        BoundsTuple(x_min=-0.5, x_max=2.5, y_min=-0.5, y_max=1.5, z_min=-0.5, z_max=0.5)
+
+        Use a reference image to control the resampling instead. Here we load two
+        images with different dimensions:
+        :func:`~pyvista.examples.downloads.download_puppy` and
+        :func:`~pyvista.examples.downloads.download_gourds`.
+
+        >>> puppy = examples.download_puppy()
+        >>> puppy.dimensions
+        (1600, 1200, 1)
+
+        >>> gourds = examples.download_gourds()
+        >>> gourds.dimensions
+        (640, 480, 1)
+
+        Use ``reference_image`` to resample the puppy to match the gourds geometry or
+        vice-versa.
+
+        >>> puppy_resampled = puppy.resample(reference_image=gourds)
+        >>> puppy_resampled.dimensions
+        (640, 480, 1)
+
+        >>> gourds_resampled = gourds.resample(reference_image=puppy)
+        >>> gourds_resampled.dimensions
+        (1600, 1200, 1)
+
+        Downsample the puppy image to 1/10th its original resolution.
+
+        >>> downsampled = puppy.resample(0.1)
+        >>> downsampled.dimensions
+        (160, 120, 1)
+
+        Compare the downsampled image to the original and zoom in to show detail.
+
+        >>> plt = pv.Plotter(shape=(1, 2))
+        >>> _ = plt.add_mesh(puppy, rgba=True, show_edges=False, lighting=False)
+        >>> plt.subplot(0, 1)
+        >>> _ = plt.add_mesh(downsampled, rgba=True, show_edges=False, lighting=False)
+        >>> plt.link_views()
+        >>> plt.view_xy()
+        >>> plt.camera.zoom(3.0)
+        >>> plt.show()
+
+        """
+        # Process scalars
+        if scalars is None:
+            field, name = set_default_active_scalars(self)
+        else:
+            name = scalars
+            field = self.get_array_association(scalars, preference=preference)
+
+        active_scalars = self.get_array(name, preference=field.name.lower())  # type: ignore[arg-type]
+
+        # Validate interpolation and modify scalars as needed
+        input_dtype = active_scalars.dtype
+        has_int_scalars = np.issubdtype(input_dtype, np.integer)
+        _validation.check_contains(
+            ['linear', 'nearest', 'cubic'], must_contain=interpolation, name='interpolation'
+        )
+        if has_int_scalars:
+            # Need floats for interpolation but also to avoid crashing in some cases
+            input_image = self.copy(deep=False)
+            input_image[name] = active_scalars.astype(float)
+        else:
+            input_image = self
+
+        # Make sure we have point scalars
+        processing_cell_scalars = field == FieldAssociation.CELL
+        if processing_cell_scalars:
+            if extend_border:
+                raise ValueError('`extend_border` cannot be set when resampling cell data.')
+            dimensionality = input_image.dimensionality
+            input_image = input_image.cells_to_points(scalars=scalars, copy=False)
+
+        # Set default extend_border value
+        if extend_border is None:
+            # Only extend border with point data
+            extend_border = not processing_cell_scalars
+        elif extend_border and reference_image is not None:
+            raise ValueError('`extend_border` cannot be set when a `image_reference` is provided.')
+
+        # Setup reference image
+        if reference_image is None:
+            # Use the input as a reference
+            reference_image = pyvista.ImageData()
+            reference_image.copy_structure(input_image)
+            reference_image_provided = False
+        else:
+            if dimensions is not None or sample_rate is not None:
+                raise ValueError(
+                    'Cannot specify a reference image along with `dimensions` or `sample_rate` parameters.\n'
+                    '`reference_image` must define the geometry exclusively.'
+                )
+            _validation.check_instance(reference_image, pyvista.ImageData, name='reference_image')
+            reference_image_provided = True
+
+        # Ideally vtkImageResample would directly support setting output dimensions
+        # (e.g. via SetOutputExtent) but this doesn't work, so instead we are
+        # stuck using SetMagnificationFactors to indirectly set the dimensions.
+        # To compute the magnification factors we first define input (old) and output
+        # (new) dimensions.
+        old_dimensions = np.array(input_image.dimensions)
+        if sample_rate is not None:
+            if reference_image_provided or dimensions is not None:
+                raise ValueError(
+                    'Cannot specify a sample rate along with `reference_image` or `sample_rate` parameters.\n'
+                    '`sample_rate` must define the sampling geometry exclusively.'
+                )
+            # Set reference dimensions from the sample rate
+            sample_rate_ = _validation.validate_array3(
+                sample_rate,
+                broadcast=True,
+                must_be_in_range=[0, np.inf],
+                strict_lower_bound=True,
+                name='sample_rate',
+            )
+            new_dimensions = old_dimensions * sample_rate_
+        else:
+            if dimensions is not None:
+                dimensions_ = np.array(dimensions)
+                dimensions_ = dimensions_ - 1 if processing_cell_scalars else dimensions_
+                reference_image.dimensions = dimensions_  # type: ignore[assignment]
+            new_dimensions = np.array(reference_image.dimensions)
+
+        # Compute the magnification factors to use with vtkImageResample
+        # Note that SetMagnificationFactors will multiply the factors by the extent
+        # but we want to multiply the dimensions. These values are off by one.
+        singleton_dims = old_dimensions == 1
+        with np.errstate(divide='ignore', invalid='ignore'):
+            # Ignore division by zero, this is fixed with singleton_dims on the next line
+            magnification_factors = (new_dimensions - 1) / (old_dimensions - 1)
+        magnification_factors[singleton_dims] = 1
+
+        resample = _vtk.vtkImageResample()
+        resample.SetInputData(input_image)
+        resample.SetMagnificationFactors(*magnification_factors)
+
+        # Set interpolation mode
+        if interpolation == 'linear':
+            resample.SetInterpolationModeToLinear()
+        elif interpolation == 'cubic':
+            resample.SetInterpolationModeToCubic()
+        else:  # 'nearest'
+            resample.SetInterpolationModeToNearestNeighbor()
+
+        # Get output
+        _update_alg(resample, progress_bar=progress_bar)
+        output_image = _get_output(resample).copy(deep=False)
+
+        # Set geometry from the reference
+        output_image.direction_matrix = reference_image.direction_matrix
+        output_image.origin = reference_image.origin
+        output_image.offset = reference_image.offset
+
+        if reference_image_provided:
+            output_image.spacing = reference_image.spacing
+        else:
+            # Need to fixup the spacing
+            old_spacing = np.array(input_image.spacing)
+            output_dimensions = np.array(output_image.dimensions)
+
+            if extend_border and not processing_cell_scalars:
+                # Compute spacing to have the same effective sample rate as the dimensions
+                actual_sample_rate = output_dimensions / old_dimensions
+                new_spacing = old_spacing / actual_sample_rate
+
+                # This will enlarge the image, so we need to shift the origin accordingly
+                # Shift the origin by 1/2 of the old and new spacing, but keep the spacing
+                # unchanged for singleton dimensions.
+                shift_old = old_spacing[~singleton_dims] / 2
+                shift_new = new_spacing[~singleton_dims] / 2
+                new_origin = np.array(input_image.origin)
+                new_origin[~singleton_dims] += shift_new - shift_old
+
+                output_image.origin = new_origin
+            else:
+                # Compute spacing to match bounds of input and dimensions of output
+                bnds = input_image.bounds
+                size = np.array(
+                    (bnds.x_max - bnds.x_min, bnds.y_max - bnds.y_min, bnds.z_max - bnds.z_min)
+                )
+                if processing_cell_scalars:
+                    new_spacing = (size + input_image.spacing) / (output_dimensions)
+                else:
+                    with np.errstate(divide='ignore', invalid='ignore'):
+                        # Ignore division by zero, this is fixed with singleton_dims on the next line
+                        new_spacing = size / (output_dimensions - 1)
+
+            # For singleton dimensions, keep the original spacing value
+            new_spacing[singleton_dims] = old_spacing[singleton_dims]
+            output_image.spacing = new_spacing
+
+        if output_image.active_scalars_name == 'ImageScalars':
+            output_image.rename_array('ImageScalars', name)
+
+        if has_int_scalars and interpolation == 'nearest':
+            # Can safely cast to int to match input
+            output_image.point_data[name] = output_image.point_data[name].astype(input_dtype)
+
+        if processing_cell_scalars:
+            # Convert back to cells. This modifies origin so we need to reset it.
+            output_image = output_image.points_to_cells(
+                scalars=name, copy=False, dimensionality=dimensionality
+            )
+            output_image.origin = (
+                reference_image.origin if reference_image_provided else self.origin
+            )
+            output_image.point_data.clear()
+        else:
+            output_image.cell_data.clear()
+
+        if inplace:
+            self.copy_from(output_image)
+            return self
+        return output_image


Test patch:
diff --git a/tests/core/test_imagedata_filters.py b/tests/core/test_imagedata_filters.py
index 4a0d4425..3a9b23d2 100644
--- a/tests/core/test_imagedata_filters.py
+++ b/tests/core/test_imagedata_filters.py
@@ -1091,3 +1091,161 @@ def test_validate_dim_operation_invalid_parameters(
         image._validate_dimensional_operation(
             operation_mask=operation_mask, operator=operator, operation_size=(1, 3, 5)
         )
+
+
+@pytest.mark.parametrize('spacing', [None, [0.3, 0.4, 0.5]])
+@pytest.mark.parametrize('direction_matrix', [None, np.diag((-1, 1, 1))])
+@pytest.mark.parametrize('origin', [None, (1.1, 2.2, 3.3)])
+@pytest.mark.parametrize('dimensions', [None, (5, 6, 7)])
+@pytest.mark.parametrize('offset', [None, (-100, -101, -102)])
+def test_resample_reference_image(uniform, spacing, direction_matrix, origin, dimensions, offset):
+    uniform = pv.ImageData(dimensions=(4, 4, 4))
+    uniform['data'] = np.arange(uniform.n_points, dtype=float)
+    reference = uniform.copy()
+    if spacing is not None:
+        reference.spacing = spacing
+    if direction_matrix is not None:
+        reference.direction_matrix = direction_matrix
+    if origin is not None:
+        reference.origin = origin
+    if dimensions is not None:
+        reference.dimensions = dimensions
+    if offset is not None:
+        reference.offset = offset
+    reference['data'] = range(reference.n_points)
+
+    resampled = uniform.resample(reference_image=reference, progress_bar=True)
+    assert isinstance(resampled, pv.ImageData)
+    assert resampled is not uniform
+    assert resampled is not reference
+    assert np.array_equal(resampled.dimensions, reference.dimensions)
+    assert np.array_equal(resampled.offset, reference.offset)
+    assert np.allclose(resampled.index_to_physical_matrix, reference.index_to_physical_matrix)
+    assert len(resampled.active_scalars) == resampled.n_points
+    assert resampled.active_scalars_name == uniform.active_scalars_name
+    assert np.allclose(resampled.bounds, reference.bounds)
+
+
+@pytest.mark.parametrize(
+    ('name', 'value'),
+    [
+        ('sample_rate', 1.5),
+        ('sample_rate', 2.0),
+        ('sample_rate', 0.5),
+        ('dimensions', (15, 15, 15)),
+    ],
+)
+@pytest.mark.parametrize('extend_border', [True, False])
+def test_resample_extend_border(uniform, extend_border, name, value):
+    kwarg = {name: value}
+    resampled = uniform.resample(**kwarg, extend_border=extend_border)
+    expected_dimensions = uniform.dimensions * np.array(value) if name == 'sample_rate' else value
+
+    assert np.array_equal(resampled.dimensions, expected_dimensions)
+    assert len(resampled.active_scalars) == resampled.n_points
+    assert resampled.active_scalars_name == uniform.active_scalars_name
+
+    if extend_border:
+        sample_rate = np.array(resampled.dimensions) / uniform.dimensions
+        expected_spacing = uniform.spacing / sample_rate
+        assert np.allclose(resampled.spacing, expected_spacing)
+        assert not np.allclose(resampled.bounds, uniform.bounds)
+
+        # Test bounds are the same when represented as cells
+        expected_cell_bounds = uniform.points_to_cells().bounds
+        actual_cell_bounds = resampled.points_to_cells().bounds
+        assert np.allclose(actual_cell_bounds, expected_cell_bounds)
+    else:
+        assert np.allclose(resampled.bounds, uniform.bounds)
+
+
+@pytest.mark.parametrize('dtype', ['uint8', 'int', 'float'])
+@pytest.mark.parametrize('interpolation', ['linear', 'nearest', 'cubic'])
+def test_resample_interpolation(uniform, interpolation, dtype):
+    array = uniform.active_scalars
+    uniform[uniform.active_scalars_name] = array.astype(dtype)
+    resampled = uniform.resample(interpolation=interpolation)
+
+    expected_dtype = float if interpolation in ['linear', 'cubic'] else dtype
+    actual_dtype = resampled.active_scalars.dtype
+    assert actual_dtype == expected_dtype
+
+
+@pytest.mark.parametrize(
+    ('name', 'value'),
+    [
+        ('scalars', 'Spatial Point Data'),
+        ('scalars', 'Spatial Cell Data'),
+        ('preference', 'point'),
+        ('preference', 'cell'),
+    ],
+)
+def test_resample_scalars(uniform, name, value):
+    kwargs = {name: value}
+    if name == 'preference':
+        # Make array names ambiguous
+        scalars = 'data'
+        uniform.rename_array('Spatial Point Data', scalars)
+        uniform.rename_array('Spatial Cell Data', scalars)
+        kwargs['scalars'] = scalars
+
+    resampled = uniform.resample(**kwargs)
+
+    if 'cell' in value.lower():
+        assert len(resampled.cell_data) == 1
+        assert len(resampled.point_data) == 0
+    else:
+        assert len(resampled.cell_data) == 0
+        assert len(resampled.point_data) == 1
+
+
+def test_resample_cell_data(uniform):
+    uniform.point_data.clear()
+    input_cell_dimensions = np.array(uniform.dimensions) - 1
+    sample_rate = np.array((2, 3, 4))
+
+    # Test sample rate
+    resampled = uniform.resample(sample_rate)
+    resample_cell_dimensions = np.array(resampled.dimensions) - 1
+    expected_cell_dimensions = input_cell_dimensions * sample_rate
+    assert np.array_equal(resample_cell_dimensions, expected_cell_dimensions)
+    assert np.allclose(uniform.bounds, resampled.bounds)
+
+    # Test dimensions
+    output_dimensions = (12, 13, 14)
+    resampled = uniform.resample(dimensions=output_dimensions)
+    assert np.array_equal(resampled.dimensions, output_dimensions)
+    assert np.allclose(uniform.bounds, resampled.bounds)
+
+
+def test_resample_inplace(uniform):
+    resampled = uniform.resample()
+    assert resampled is not uniform
+    resampled = uniform.resample(inplace=True)
+    assert resampled is uniform
+
+
+def test_resample_raises(uniform):
+    match = (
+        'Cannot specify a reference image along with `dimensions` or `sample_rate` parameters.\n'
+        '`reference_image` must define the geometry exclusively.'
+    )
+    with pytest.raises(ValueError, match=re.escape(match)):
+        uniform.resample(sample_rate=2, reference_image=uniform)
+    with pytest.raises(ValueError, match=re.escape(match)):
+        uniform.resample(dimensions=(2, 2, 2), reference_image=uniform)
+
+    match = (
+        'Cannot specify a sample rate along with `reference_image` or `sample_rate` parameters.\n'
+        '`sample_rate` must define the sampling geometry exclusively.'
+    )
+    with pytest.raises(ValueError, match=re.escape(match)):
+        uniform.resample(sample_rate=2, dimensions=(2, 2, 2))
+
+    match = '`extend_border` cannot be set when resampling cell data.'
+    with pytest.raises(ValueError, match=re.escape(match)):
+        uniform.resample(scalars='Spatial Cell Data', extend_border=True)
+
+    match = '`extend_border` cannot be set when a `image_reference` is provided.'
+    with pytest.raises(ValueError, match=re.escape(match)):
+        uniform.resample(reference_image=uniform, extend_border=True)