Affine transformation from 2D to 3D using SciPy leads to empty arrays

Question

Using scipy.ndimage.affine_transform, I am trying to apply an affine transformation on a 3D array with one degenerate dimension, e.g. with shape (10, 1, 10), and get a non-degenerate 3D output shape, e.g. (10, 10, 10). The purpose of this operation is to transform 2D medical images to the coordinate system of a 3D template.

Let's generate an example slice:

import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import affine_transform

rng = np.random.default_rng(42)
random_slice = rng.random((10, 1, 10))

plt.imshow(random_slice[:, 0])
plt.axis("off")

Now, applying a simple offset using affine_transform and outputting to shape (10, 10, 10) works as intended:

affine = np.eye(3)
transformed_slice = affine_transform(
    random_slice, affine, offset=-1.0, output_shape=(10, 10, 10)
)

plt.imshow(transformed_slice[:, 1, :])
plt.axis("off")

However, as soon as scaling components become non-unitary, the resulting array is empty, e.g.

affine = np.eye(3) * 0.999
transformed_slice = affine_transform(
    random_slice, affine, offset=-1.0, output_shape=(10, 10, 10)
)

>>> transformed_slice.sum()
0.0

What's causing this?

What I've tested:

• Using mode="grid-constant" solves the problem by allowing interpolation outside of the original array bounds. However, it creates a lot of ringing artifacts, making this approach unusable. The problem might come from SciPy's computation of the original array bounds though, but I'm not sure how to test this.
• Setting the scaling component along the degenerate dimension to 1 solves the problem. This is the best solution for now, but I would like to understand what's causing affine_transform to fail even for scaling components very close to 1.