`resampling.apply` chokes when transforming niftis that have fourth dimension equal to 1

[psadil]

import numpy as np
import nibabel as nb
import nitransforms as nt

dataobj=np.reshape(np.arange(27, dtype=np.uint8), (3,3,3))
nii = nb.nifti1.Nifti1Image(dataobj=dataobj, affine=np.eye(4))
transform=nt.linear.Affine(nii.affine, nii)
nt.resampling.apply(transform, spatialimage=nii)
print("success!")
nii2 = nb.nifti1.Nifti1Image(dataobj=dataobj[..., np.newaxis], affine=np.eye(4))
transform=nt.linear.Affine(nii.affine, nii2)
# this will fail
nt.resampling.apply(transform, spatialimage=nii2)

gives

success!
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[1], line 12
     10 nii2 = nb.nifti1.Nifti1Image(dataobj=dataobj[..., np.newaxis], affine=np.eye(4))
     11 transform=nt.linear.Affine(nii.affine, nii2)
---> 12 nt.resampling.apply(transform, spatialimage=nii2)

File ~/git/neuro/qcapp/django-qcapp-ratings/.venv/lib/python3.12/site-packages/nitransforms/resampling.py:333, in apply(transform, spatialimage, reference, order, mode, cval, prefilter, output_dtype, dtype_width, serialize_nvols, max_concurrent)
    330         targets = np.vstack((time_row, ijk_targets))
    331         data = np.rollaxis(data, data.ndim - 1, 0)
--> 333     resampled = ndi.map_coordinates(
    334         data,
    335         targets,
    336         order=order,
    337         mode=mode,
    338         cval=cval,
    339         prefilter=prefilter,
    340     )
    342 if isinstance(_ref, ImageGrid):  # If reference is grid, reshape
    343     hdr = (
    344         _ref.header.copy()
    345         if _ref.header is not None
    346         else spatialimage.header.__class__()
    347     )

File ~/git/neuro/qcapp/django-qcapp-ratings/.venv/lib/python3.12/site-packages/scipy/ndimage/_interpolation.py:456, in map_coordinates(input, coordinates, output, order, mode, cval, prefilter)
    454     raise RuntimeError('input and output rank must be > 0')
    455 if coordinates.shape[0] != input.ndim:
--> 456     raise RuntimeError('invalid shape for coordinate array')
    457 complex_output = np.iscomplexobj(input)
    458 output = _ni_support._get_output(output, input, shape=output_shape,
    459                                  complex_output=complex_output)

RuntimeError: invalid shape for coordinate array

A workaround is to remove the final dimension, e.g. with nb.funcs.squeeze_image, but it would be nice if that were handled automatically. Not sure whether the result should be 4d or 3d.

[oesteban]

Thanks for reporting! I'll add tests for this

[effigies]

Just going to quickly reopen this to resolve this discussion. I'd rather not continue it on a closed PR...

In #278 (comment), @oesteban wrote:

Okay tests do think it's okay to squeeze: there are these potential cases:

Input	Transform	Output	Example
3D (including 4D singleton)	3D	3D	A T1w to MNI T1w transform
4D	4D	4D (the number of timepoints in input and transform must match)	BOLD head-motion correction
3D (including 4D singleton)	4D	4D (the 3D input resampled through each of the transformations)	A T1w to each of BOLD timepoints, including head-motion
4D	3D	4D (the same transform applied to each timepoint)	Moving BOLD timepoints into T1w space, without considering head-motion

There could be an argument as to whether the output should be 4D singleton in the particular case this is fixing, but I would honestly place that burden on the user to match the input image.

I'm confused why the singleton case doesn't fall into the (input=4D, transform=3D) case. That is exactly what I would expect to get, and it seems that we don't actually support it.

[oesteban]

It's nice to see some paying attention! The only reason I closed this prematurely is laziness (or lack of bandwidth to write three more LoCs).