diff --git a/README.md b/README.md
index 4ae2340f..2095f7a3 100644
--- a/README.md
+++ b/README.md
@@ -29,6 +29,11 @@ Compatibility between software tools and imaging formats is a necessary bridge
 to ensure the reproducibility of results and enable the optimization and evaluation
 of current image processing and analysis workflows.
 
+## BIDS' X5 format
+As of [the 25.0.0 release](https://github.com/nipy/nitransforms/releases/tag/25.0.0),
+*NiTransforms* experimentally supports writing X5 transform files, as drafted in the
+BIDS Extension Proposal 14 (BEP014).
+
 ## Integration with *NiBabel*
 *NiTransforms* started as a feature-repo spun off of *NiBabel*.
 Shortly after starting with [nipy/nibabel#656](https://github.com/nipy/nibabel/pull/656), it became apparent that it was going to build up in a humongous PR nobody would be able to review as thoroughly as it would require.
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index 6166386d..df237495 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -19,6 +19,7 @@
 from nibabel.loadsave import load as _nbload
 from nibabel.arrayproxy import get_obj_dtype
 from nibabel.spatialimages import SpatialImage
+from nibabel.funcs import squeeze_image
 from scipy import ndimage as ndi
 
 from nitransforms.base import (
@@ -233,6 +234,9 @@ def apply(
     if isinstance(spatialimage, (str, Path)):
         spatialimage = _nbload(str(spatialimage))
 
+    singleton_4d = spatialimage.ndim == 4 and spatialimage.shape[-1] == 1
+    spatialimage = squeeze_image(spatialimage)
+
     # Avoid opening the data array just yet
     input_dtype = cap_dtype(get_obj_dtype(spatialimage.dataobj), dtype_width)
 
@@ -255,22 +259,14 @@ def apply(
     targets = None
     ref_ndcoords = _ref.ndcoords
 
-    if hasattr(transform, "to_field") and callable(transform.to_field):
-        targets = ImageGrid(spatialimage).index(
-            _as_homogeneous(
-                transform.to_field(reference=reference).map(ref_ndcoords),
-                dim=_ref.ndim,
-            )
-        )
-    else:
-        # Targets' shape is (Nt, 3, Nv) with Nv = Num. voxels, Nt = Num. timepoints.
-        targets = (
-            ImageGrid(spatialimage).index(
-                _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim)
-            )
-            if targets is None
-            else targets
+    # Targets' shape is (Nt, 3, Nv) with Nv = Num. voxels, Nt = Num. timepoints.
+    targets = (
+        ImageGrid(spatialimage).index(
+            _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim)
         )
+        if targets is None
+        else targets
+    )
 
     if targets.ndim == 2:
         targets = targets.T[np.newaxis, ...]
@@ -367,7 +363,8 @@ def apply(
             with suppress(ValueError):
                 resampled = np.squeeze(resampled, axis=3)
 
-        moved = spatialimage.__class__(resampled, _ref.affine, hdr)
+        moved = spatialimage.__class__(
+            resampled[..., None] if singleton_4d else resampled, _ref.affine, hdr)
         return moved
 
     output_dtype = output_dtype or input_dtype
diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py
index d7538ea4..f9b88577 100644
--- a/nitransforms/tests/test_resampling.py
+++ b/nitransforms/tests/test_resampling.py
@@ -51,6 +51,20 @@
 }
 
 
+def test_apply_singleton_time_dimension():
+    """Resampling fails when the input image has a trailing singleton dimension (gh-270)."""
+
+    data = np.reshape(np.arange(27, dtype=np.uint8), (3, 3, 3, 1))
+    nii = nb.Nifti1Image(data, np.eye(4))
+    ref = nb.Nifti1Image(np.zeros((4, 4, 4)), np.eye(4))
+    xfm = nitl.Affine(np.eye(4), reference=ref)
+    movednii = apply(xfm, nii)
+    assert movednii.shape == ref.shape + (1, )
+
+    movednii = apply(xfm, nii, reference=ref)
+    assert movednii.shape == ref.shape + (1, )
+
+
 @pytest.mark.parametrize(
     "image_orientation",
     [