From e38c16a461d44dff932a732072cd2548693676ce Mon Sep 17 00:00:00 2001
From: Julien Marabotto <166002186+jmarabotto@users.noreply.github.com>
Date: Fri, 5 Apr 2024 10:00:33 +0200
Subject: [PATCH 01/15] enh: outsource the apply function

---
 nitransforms/base.py       |  95 -------------------------------
 nitransforms/resampling.py | 114 +++++++++++++++++++++++++++++++++++++
 2 files changed, 114 insertions(+), 95 deletions(-)
 create mode 100644 nitransforms/resampling.py

diff --git a/nitransforms/base.py b/nitransforms/base.py
index 96f00ed..b970434 100644
--- a/nitransforms/base.py
+++ b/nitransforms/base.py
@@ -222,101 +222,6 @@ def ndim(self):
         """Access the dimensions of the reference space."""
         raise TypeError("TransformBase has no dimensions")
 
-    def apply(
-        self,
-        spatialimage,
-        reference=None,
-        order=3,
-        mode="constant",
-        cval=0.0,
-        prefilter=True,
-        output_dtype=None,
-    ):
-        """
-        Apply a transformation to an image, resampling on the reference spatial object.
-
-        Parameters
-        ----------
-        spatialimage : `spatialimage`
-            The image object containing the data to be resampled in reference
-            space
-        reference : spatial object, optional
-            The image, surface, or combination thereof containing the coordinates
-            of samples that will be sampled.
-        order : int, optional
-            The order of the spline interpolation, default is 3.
-            The order has to be in the range 0-5.
-        mode : {'constant', 'reflect', 'nearest', 'mirror', 'wrap'}, optional
-            Determines how the input image is extended when the resamplings overflows
-            a border. Default is 'constant'.
-        cval : float, optional
-            Constant value for ``mode='constant'``. Default is 0.0.
-        prefilter: bool, optional
-            Determines if the image's data array is prefiltered with
-            a spline filter before interpolation. The default is ``True``,
-            which will create a temporary *float64* array of filtered values
-            if *order > 1*. If setting this to ``False``, the output will be
-            slightly blurred if *order > 1*, unless the input is prefiltered,
-            i.e. it is the result of calling the spline filter on the original
-            input.
-        output_dtype: dtype specifier, optional
-            The dtype of the returned array or image, if specified.
-            If ``None``, the default behavior is to use the effective dtype of
-            the input image. If slope and/or intercept are defined, the effective
-            dtype is float64, otherwise it is equivalent to the input image's
-            ``get_data_dtype()`` (on-disk type).
-            If ``reference`` is defined, then the return value is an image, with
-            a data array of the effective dtype but with the on-disk dtype set to
-            the input image's on-disk dtype.
-
-        Returns
-        -------
-        resampled : `spatialimage` or ndarray
-            The data imaged after resampling to reference space.
-
-        """
-        if reference is not None and isinstance(reference, (str, Path)):
-            reference = _nbload(str(reference))
-
-        _ref = (
-            self.reference if reference is None else SpatialReference.factory(reference)
-        )
-
-        if _ref is None:
-            raise TransformError("Cannot apply transform without reference")
-
-        if isinstance(spatialimage, (str, Path)):
-            spatialimage = _nbload(str(spatialimage))
-
-        data = np.asanyarray(spatialimage.dataobj)
-        targets = ImageGrid(spatialimage).index(  # data should be an image
-            _as_homogeneous(self.map(_ref.ndcoords.T), dim=_ref.ndim)
-        )
-
-        resampled = ndi.map_coordinates(
-            data,
-            targets.T,
-            output=output_dtype,
-            order=order,
-            mode=mode,
-            cval=cval,
-            prefilter=prefilter,
-        )
-
-        if isinstance(_ref, ImageGrid):  # If reference is grid, reshape
-            hdr = None
-            if _ref.header is not None:
-                hdr = _ref.header.copy()
-                hdr.set_data_dtype(output_dtype or spatialimage.get_data_dtype())
-            moved = spatialimage.__class__(
-                resampled.reshape(_ref.shape),
-                _ref.affine,
-                hdr,
-            )
-            return moved
-
-        return resampled
-
     def map(self, x, inverse=False):
         r"""
         Apply :math:`y = f(x)`.
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
new file mode 100644
index 0000000..7d2765a
--- /dev/null
+++ b/nitransforms/resampling.py
@@ -0,0 +1,114 @@
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+"""Resampling utilities."""
+from pathlib import Path
+import numpy as np
+import h5py
+import warnings
+from nibabel.loadsave import load as _nbload
+from nibabel import funcs as _nbfuncs
+from nibabel.nifti1 import intent_codes as INTENT_CODES
+from nibabel.cifti2 import Cifti2Image
+from scipy import ndimage as ndi
+
+
+def apply(
+    transform,
+    spatialimage,
+    reference=None,
+    order=3,
+    mode="constant",
+    cval=0.0,
+    prefilter=True,
+    output_dtype=None,
+):
+    """
+    Apply a transformation to an image, resampling on the reference spatial object.
+
+    Parameters
+    ----------
+    spatialimage : `spatialimage`
+        The image object containing the data to be resampled in reference
+        space
+    reference : spatial object, optional
+        The image, surface, or combination thereof containing the coordinates
+        of samples that will be sampled.
+    order : int, optional
+        The order of the spline interpolation, default is 3.
+        The order has to be in the range 0-5.
+    mode : {'constant', 'reflect', 'nearest', 'mirror', 'wrap'}, optional
+        Determines how the input image is extended when the resamplings overflows
+        a border. Default is 'constant'.
+    cval : float, optional
+        Constant value for ``mode='constant'``. Default is 0.0.
+    prefilter: bool, optional
+        Determines if the image's data array is prefiltered with
+        a spline filter before interpolation. The default is ``True``,
+        which will create a temporary *float64* array of filtered values
+        if *order > 1*. If setting this to ``False``, the output will be
+        slightly blurred if *order > 1*, unless the input is prefiltered,
+        i.e. it is the result of calling the spline filter on the original
+        input.
+    output_dtype: dtype specifier, optional
+        The dtype of the returned array or image, if specified.
+        If ``None``, the default behavior is to use the effective dtype of
+        the input image. If slope and/or intercept are defined, the effective
+        dtype is float64, otherwise it is equivalent to the input image's
+        ``get_data_dtype()`` (on-disk type).
+        If ``reference`` is defined, then the return value is an image, with
+        a data array of the effective dtype but with the on-disk dtype set to
+        the input image's on-disk dtype.
+
+    Returns
+    -------
+    resampled : `spatialimage` or ndarray
+        The data imaged after resampling to reference space.
+
+    """
+    if reference is not None and isinstance(reference, (str, Path)):
+        reference = _nbload(str(reference))
+
+    _ref = (
+        transform.reference if reference is None else SpatialReference.factory(reference)
+    )
+
+    if _ref is None:
+        raise TransformError("Cannot apply transform without reference")
+
+    if isinstance(spatialimage, (str, Path)):
+        spatialimage = _nbload(str(spatialimage))
+
+    data = np.asanyarray(spatialimage.dataobj)
+    targets = ImageGrid(spatialimage).index(  # data should be an image
+        _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
+    )
+
+    resampled = ndi.map_coordinates(
+        data,
+        targets.T,
+        output=output_dtype,
+        order=order,
+        mode=mode,
+        cval=cval,
+        prefilter=prefilter,
+    )
+
+    if isinstance(_ref, ImageGrid):  # If reference is grid, reshape
+        hdr = None
+        if _ref.header is not None:
+            hdr = _ref.header.copy()
+            hdr.set_data_dtype(output_dtype or spatialimage.get_data_dtype())
+        moved = spatialimage.__class__(
+            resampled.reshape(_ref.shape),
+            _ref.affine,
+            hdr,
+        )
+        return moved
+
+    return resampled

From f396f94daa6bd542104d49dadb6e744f7f22b322 Mon Sep 17 00:00:00 2001
From: Julien Marabotto <166002186+jmarabotto@users.noreply.github.com>
Date: Fri, 5 Apr 2024 10:06:05 +0200
Subject: [PATCH 02/15] sty: pacify flake8

---
 nitransforms/resampling.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index 7d2765a..a876bb1 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -9,14 +9,14 @@
 """Resampling utilities."""
 from pathlib import Path
 import numpy as np
-import h5py
-import warnings
 from nibabel.loadsave import load as _nbload
-from nibabel import funcs as _nbfuncs
-from nibabel.nifti1 import intent_codes as INTENT_CODES
-from nibabel.cifti2 import Cifti2Image
-from scipy import ndimage as ndi
 
+from nitransforms.base import (
+    ImageGrid,
+    TransformError,
+    SpatialReference,
+    _as_homogeneous,
+)
 
 def apply(
     transform,

From a150da488b48e038c898d0392d5d8bc33633b342 Mon Sep 17 00:00:00 2001
From: Julien Marabotto <166002186+jmarabotto@users.noreply.github.com>
Date: Fri, 5 Apr 2024 10:09:33 +0200
Subject: [PATCH 03/15] sty: fix imports

---
 nitransforms/base.py       | 1 -
 nitransforms/resampling.py | 2 ++
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/nitransforms/base.py b/nitransforms/base.py
index b970434..e7f67e6 100644
--- a/nitransforms/base.py
+++ b/nitransforms/base.py
@@ -15,7 +15,6 @@
 from nibabel import funcs as _nbfuncs
 from nibabel.nifti1 import intent_codes as INTENT_CODES
 from nibabel.cifti2 import Cifti2Image
-from scipy import ndimage as ndi
 
 EQUALITY_TOL = 1e-5
 
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index a876bb1..e89b081 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -10,6 +10,7 @@
 from pathlib import Path
 import numpy as np
 from nibabel.loadsave import load as _nbload
+from scipy import ndimage as ndi
 
 from nitransforms.base import (
     ImageGrid,
@@ -18,6 +19,7 @@
     _as_homogeneous,
 )
 
+
 def apply(
     transform,
     spatialimage,

From 3ff7407c8603f64c6a364115ad94aa40a2097fbc Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Fri, 5 Apr 2024 10:35:24 +0200
Subject: [PATCH 04/15] fix: update many apply() calls

---
 nitransforms/cli.py               |  4 ++-
 nitransforms/tests/test_base.py   |  9 ++++---
 nitransforms/tests/test_io.py     | 41 ++++++++++++++++++++-----------
 nitransforms/tests/test_linear.py | 16 ++++++------
 4 files changed, 45 insertions(+), 25 deletions(-)

diff --git a/nitransforms/cli.py b/nitransforms/cli.py
index 63b8bed..8f8f5ce 100644
--- a/nitransforms/cli.py
+++ b/nitransforms/cli.py
@@ -5,6 +5,7 @@
 
 from .linear import load as linload
 from .nonlinear import load as nlinload
+from .resampling import apply
 
 
 def cli_apply(pargs):
@@ -38,7 +39,8 @@ def cli_apply(pargs):
     # ensure a reference is set
     xfm.reference = pargs.ref or pargs.moving
 
-    moved = xfm.apply(
+    moved = apply(
+        xfm,
         pargs.moving,
         order=pargs.order,
         mode=pargs.mode,
diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py
index 07a7e4e..a1402ba 100644
--- a/nitransforms/tests/test_base.py
+++ b/nitransforms/tests/test_base.py
@@ -6,6 +6,7 @@
 
 from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase
 from .. import linear as nitl
+from ..resampling import apply
 
 
 def test_SpatialReference(testdata_path):
@@ -94,9 +95,10 @@ def _to_hdf5(klass, x5_root):
     # Test identity transform
     xfm = TransformBase()
     xfm.reference = fname
+
     with pytest.raises(TypeError):
         _ = xfm.ndim
-    moved = xfm.apply(fname, order=0)
+    moved = apply(xfm, fname, order=0)
     assert np.all(
         imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())
     )
@@ -104,9 +106,10 @@ def _to_hdf5(klass, x5_root):
     # Test identity transform - setting reference
     xfm = TransformBase()
     xfm.reference = fname
+
     with pytest.raises(TypeError):
         _ = xfm.ndim
-    moved = xfm.apply(str(fname), reference=fname, order=0)
+    moved = apply(xfm, str(fname), reference=fname, order=0)
     assert np.all(
         imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())
     )
@@ -126,7 +129,7 @@ def _to_hdf5(klass, x5_root):
             )
         ]
     )
-    giimoved = xfm.apply(fname, reference=gii, order=0)
+    giimoved = apply(xfm, fname, reference=gii, order=0)
     assert np.allclose(giimoved.reshape(xfm.reference.shape), moved.get_fdata())
 
     # Test to_filename
diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py
index bcee919..0cc79d1 100644
--- a/nitransforms/tests/test_io.py
+++ b/nitransforms/tests/test_io.py
@@ -28,6 +28,8 @@
 )
 from nitransforms.io.base import LinearParameters, TransformIOError, TransformFileError
 from nitransforms.conftest import _datadir, _testdir
+from nitransforms.resampling import apply
+
 
 LPS = np.diag([-1, -1, 1, 1])
 ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS))
@@ -497,10 +499,13 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
     assert np.allclose(card_aff, nb.load("deob_3drefit.nii.gz").affine)
 
     # Check that nitransforms can emulate 3drefit -deoblique
-    nt3drefit = Affine(
-        afni._cardinal_rotation(img.affine, False),
-        reference="deob_3drefit.nii.gz",
-    ).apply("orig.nii.gz")
+    nt3drefit = apply(
+        Affine(
+            afni._cardinal_rotation(img.affine, False),
+            reference="deob_3drefit.nii.gz",
+        ),
+        "orig.nii.gz",
+    )
 
     diff = (
         np.asanyarray(img.dataobj, dtype="uint8")
@@ -509,10 +514,13 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
     assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1
 
     # Check that nitransforms can revert 3drefit -deoblique
-    nt_undo3drefit = Affine(
-        afni._cardinal_rotation(img.affine, True),
-        reference="orig.nii.gz",
-    ).apply("deob_3drefit.nii.gz")
+    nt_undo3drefit = apply(
+        Affine(
+            afni._cardinal_rotation(img.affine, True),
+            reference="orig.nii.gz",
+        ),
+        "deob_3drefit.nii.gz",
+    )
 
     diff = (
         np.asanyarray(img.dataobj, dtype="uint8")
@@ -531,16 +539,21 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
     assert np.allclose(deobaff, deobnii.affine)
 
     # Check resampling in deobliqued grid
-    ntdeobnii = Affine(np.eye(4), reference=deobnii.__class__(
-        np.zeros(deobshape, dtype="uint8"),
-        deobaff,
-        deobnii.header
-    )).apply(img, order=0)
+    ntdeobnii = apply(
+        Affine(np.eye(4), reference=deobnii.__class__(
+            np.zeros(deobshape, dtype="uint8"),
+            deobaff,
+            deobnii.header
+        )),
+        img,
+        order=0,
+    )
 
     # Generate an internal box to exclude border effects
     box = np.zeros(img.shape, dtype="uint8")
     box[10:-10, 10:-10, 10:-10] = 1
-    ntdeobmask = Affine(np.eye(4), reference=ntdeobnii).apply(
+    ntdeobmask = apply(
+        Affine(np.eye(4), reference=ntdeobnii),
         nb.Nifti1Image(box, img.affine, img.header),
         order=0,
     )
diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py
index 2957f59..9a06fe3 100644
--- a/nitransforms/tests/test_linear.py
+++ b/nitransforms/tests/test_linear.py
@@ -13,6 +13,7 @@
 from nibabel.affines import from_matvec
 from nitransforms import linear as nitl
 from nitransforms import io
+from nitransforms.resampling import apply
 from .utils import assert_affines_by_filename
 
 RMSE_TOL = 0.1
@@ -285,7 +286,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient
     assert exit_code == 0
     sw_moved_mask = nb.load("resampled_brainmask.nii.gz")
 
-    nt_moved_mask = xfm.apply(msk, order=0)
+    nt_moved_mask = apply(xfm, msk, order=0)
     nt_moved_mask.set_data_dtype(msk.get_data_dtype())
     nt_moved_mask.to_filename("ntmask.nii.gz")
     diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj)
@@ -305,7 +306,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient
     sw_moved = nb.load("resampled.nii.gz")
     sw_moved.set_data_dtype(img.get_data_dtype())
 
-    nt_moved = xfm.apply(img, order=0)
+    nt_moved = apply(xfm, img, order=0)
     diff = (
         np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype())
         - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
@@ -314,7 +315,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient
     # A certain tolerance is necessary because of resampling at borders
     assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL
 
-    nt_moved = xfm.apply("img.nii.gz", order=0)
+    nt_moved = apply(xfm, "img.nii.gz", order=0)
     diff = (
         np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype())
         - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
@@ -343,8 +344,8 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path):
     assert isinstance(hmc, nitl.LinearTransformsMapping)
 
     # Test-case: realign functional data on to sbref
-    nii = hmc.apply(
-        testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz"
+    nii = apply(
+        hmc, testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz"
     )
     assert nii.dataobj.shape[-1] == len(hmc)
 
@@ -352,13 +353,14 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path):
     hmcinv = nitl.LinearTransformsMapping(
         np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz"
     )
-    nii = hmcinv.apply(testdata_path / "fmap.nii.gz", order=1)
+    nii = apply(hmcinv, testdata_path / "fmap.nii.gz", order=1)
     assert nii.dataobj.shape[-1] == len(hmc)
 
     # Ensure a ValueError is issued when trying to do weird stuff
     hmc = nitl.LinearTransformsMapping(hmc.matrix[:1, ...])
     with pytest.raises(ValueError):
-        hmc.apply(
+        apply(
+            hmc,
             testdata_path / "func.nii.gz",
             order=1,
             reference=testdata_path / "sbref.nii.gz",

From cc3b21e205f027e0ab7d7f360de397befbab4298 Mon Sep 17 00:00:00 2001
From: Julien Marabotto <julienmarabotto@MacBook-Pro-de-Henri.local>
Date: Wed, 17 Apr 2024 11:18:02 +0200
Subject: [PATCH 05/15] Updated outsource Apply

outsourced Apply(); fixed resampled.py (line 101); implemented np.tensordor to _apply_affine() in main.py (line 287). Left to do: fix test_linear RunTime error (line 358)
---
 nitransforms/base.py              |  8 ++++++--
 nitransforms/resampling.py        | 15 ++++++++++++---
 nitransforms/tests/test_linear.py |  6 +++++-
 3 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/nitransforms/base.py b/nitransforms/base.py
index e7f67e6..99a0ee9 100644
--- a/nitransforms/base.py
+++ b/nitransforms/base.py
@@ -283,7 +283,11 @@ def _as_homogeneous(xyz, dtype="float32", dim=3):
 
     return np.hstack((xyz, np.ones((xyz.shape[0], 1), dtype=dtype)))
 
-
+#import pdb; pdb.set_trace()
 def _apply_affine(x, affine, dim):
     """Get the image array's indexes corresponding to coordinates."""
-    return affine.dot(_as_homogeneous(x, dim=dim).T)[:dim, ...].T
+    return np.tensordot(
+        affine,
+        _as_homogeneous(x, dim=dim).T,
+        axes=1,
+    )[:dim, ...]
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index e89b081..f621f7a 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -87,27 +87,36 @@ def apply(
         spatialimage = _nbload(str(spatialimage))
 
     data = np.asanyarray(spatialimage.dataobj)
+    
     targets = ImageGrid(spatialimage).index(  # data should be an image
         _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
     )
 
+    if data.ndim < targets.shape[-1]:
+        data = data[..., np.newaxis]
+
+    import pdb; pdb.set_trace()
+     
+    #import pdb; pdb.set_trace()
     resampled = ndi.map_coordinates(
         data,
-        targets.T,
+        #targets.T,
+        #Reshape targets (516096, 3, 8) --> (4, 4128768) : 
+        _as_homogeneous(targets.reshape(-2, targets.shape[0])).T,
         output=output_dtype,
         order=order,
         mode=mode,
         cval=cval,
         prefilter=prefilter,
     )
-
+     
     if isinstance(_ref, ImageGrid):  # If reference is grid, reshape
         hdr = None
         if _ref.header is not None:
             hdr = _ref.header.copy()
             hdr.set_data_dtype(output_dtype or spatialimage.get_data_dtype())
         moved = spatialimage.__class__(
-            resampled.reshape(_ref.shape),
+            resampled.reshape(_ref.shape if data.ndim < 4 else _ref.shape + (-1, )),
             _ref.affine,
             hdr,
         )
diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py
index 9a06fe3..b7f7384 100644
--- a/nitransforms/tests/test_linear.py
+++ b/nitransforms/tests/test_linear.py
@@ -353,7 +353,11 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path):
     hmcinv = nitl.LinearTransformsMapping(
         np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz"
     )
-    nii = apply(hmcinv, testdata_path / "fmap.nii.gz", order=1)
+
+    import pdb; pdb.set_trace()
+    nii = apply(
+        hmcinv, testdata_path / "fmap.nii.gz", order=1
+    )
     assert nii.dataobj.shape[-1] == len(hmc)
 
     # Ensure a ValueError is issued when trying to do weird stuff

From f16b737dfc7c3b63c722e1af772deacdbf59abff Mon Sep 17 00:00:00 2001
From: Julien Marabotto <julienmarabotto@MacBook-Pro-de-Henri.local>
Date: Wed, 17 Apr 2024 11:26:00 +0200
Subject: [PATCH 06/15] Updated changes - see previous commit for details

---
 nitransforms/resampling.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index f621f7a..88f7e6b 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -92,11 +92,9 @@ def apply(
         _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
     )
 
-    if data.ndim < targets.shape[-1]:
-        data = data[..., np.newaxis]
-
-    import pdb; pdb.set_trace()
-     
+    #if data.ndim < targets.shape[-1]:
+    #    data = data[..., np.newaxis]
+ 
     #import pdb; pdb.set_trace()
     resampled = ndi.map_coordinates(
         data,

From be7e9a9449944b7e3f9b9ad50e7ccb084a070767 Mon Sep 17 00:00:00 2001
From: Julien Marabotto
 <julienmarabotto@macbook-pro-de-henri.vpn.intranet.chuv>
Date: Mon, 22 Apr 2024 10:44:26 +0200
Subject: [PATCH 07/15] FIX: Outsource Apply

Outsourced apply, test_linear.py successful
---
 nitransforms/base.py              |   1 -
 nitransforms/linear.py            | 121 +-----------------------------
 nitransforms/resampling.py        |  10 +--
 nitransforms/tests/test_linear.py |   1 -
 4 files changed, 9 insertions(+), 124 deletions(-)

diff --git a/nitransforms/base.py b/nitransforms/base.py
index 99a0ee9..6a6ae7e 100644
--- a/nitransforms/base.py
+++ b/nitransforms/base.py
@@ -283,7 +283,6 @@ def _as_homogeneous(xyz, dtype="float32", dim=3):
 
     return np.hstack((xyz, np.ones((xyz.shape[0], 1), dtype=dtype)))
 
-#import pdb; pdb.set_trace()
 def _apply_affine(x, affine, dim):
     """Get the image array's indexes corresponding to coordinates."""
     return np.tensordot(
diff --git a/nitransforms/linear.py b/nitransforms/linear.py
index af14f39..cf16104 100644
--- a/nitransforms/linear.py
+++ b/nitransforms/linear.py
@@ -112,6 +112,10 @@ def __invert__(self):
 
         """
         return self.__class__(self._inverse)
+    
+    def __len__(self):
+        """Enable using len()."""
+        return 1 if self._matrix.ndim == 2 else len(self._matrix)
 
     def __matmul__(self, b):
         """
@@ -330,10 +334,6 @@ def __getitem__(self, i):
         """Enable indexed access to the series of matrices."""
         return Affine(self.matrix[i, ...], reference=self._reference)
 
-    def __len__(self):
-        """Enable using len()."""
-        return len(self._matrix)
-
     def map(self, x, inverse=False):
         r"""
         Apply :math:`y = f(x)`.
@@ -402,119 +402,6 @@ def to_filename(self, filename, fmt="X5", moving=None):
             ).to_filename(filename)
         return filename
 
-    def apply(
-        self,
-        spatialimage,
-        reference=None,
-        order=3,
-        mode="constant",
-        cval=0.0,
-        prefilter=True,
-        output_dtype=None,
-    ):
-        """
-        Apply a transformation to an image, resampling on the reference spatial object.
-
-        Parameters
-        ----------
-        spatialimage : `spatialimage`
-            The image object containing the data to be resampled in reference
-            space
-        reference : spatial object, optional
-            The image, surface, or combination thereof containing the coordinates
-            of samples that will be sampled.
-        order : int, optional
-            The order of the spline interpolation, default is 3.
-            The order has to be in the range 0-5.
-        mode : {"constant", "reflect", "nearest", "mirror", "wrap"}, optional
-            Determines how the input image is extended when the resamplings overflows
-            a border. Default is "constant".
-        cval : float, optional
-            Constant value for ``mode="constant"``. Default is 0.0.
-        prefilter: bool, optional
-            Determines if the image's data array is prefiltered with
-            a spline filter before interpolation. The default is ``True``,
-            which will create a temporary *float64* array of filtered values
-            if *order > 1*. If setting this to ``False``, the output will be
-            slightly blurred if *order > 1*, unless the input is prefiltered,
-            i.e. it is the result of calling the spline filter on the original
-            input.
-
-        Returns
-        -------
-        resampled : `spatialimage` or ndarray
-            The data imaged after resampling to reference space.
-
-        """
-
-        if reference is not None and isinstance(reference, (str, Path)):
-            reference = _nbload(str(reference))
-
-        _ref = (
-            self.reference if reference is None else SpatialReference.factory(reference)
-        )
-
-        if isinstance(spatialimage, (str, Path)):
-            spatialimage = _nbload(str(spatialimage))
-
-        # Avoid opening the data array just yet
-        input_dtype = get_obj_dtype(spatialimage.dataobj)
-        output_dtype = output_dtype or input_dtype
-
-        # Prepare physical coordinates of input (grid, points)
-        xcoords = _ref.ndcoords.astype("f4").T
-
-        # Invert target's (moving) affine once
-        ras2vox = ~Affine(spatialimage.affine)
-
-        if spatialimage.ndim == 4 and (len(self) != spatialimage.shape[-1]):
-            raise ValueError(
-                "Attempting to apply %d transforms on a file with "
-                "%d timepoints" % (len(self), spatialimage.shape[-1])
-            )
-
-        # Order F ensures individual volumes are contiguous in memory
-        # Also matches NIfTI, making final save more efficient
-        resampled = np.zeros(
-            (xcoords.shape[0], len(self)), dtype=output_dtype, order="F"
-        )
-
-        dataobj = (
-            np.asanyarray(spatialimage.dataobj, dtype=input_dtype)
-            if spatialimage.ndim in (2, 3)
-            else None
-        )
-
-        for t, xfm_t in enumerate(self):
-            # Map the input coordinates on to timepoint t of the target (moving)
-            ycoords = xfm_t.map(xcoords)[..., : _ref.ndim]
-
-            # Calculate corresponding voxel coordinates
-            yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim]
-
-            # Interpolate
-            resampled[..., t] = ndi.map_coordinates(
-                (
-                    dataobj
-                    if dataobj is not None
-                    else spatialimage.dataobj[..., t].astype(input_dtype, copy=False)
-                ),
-                yvoxels.T,
-                output=output_dtype,
-                order=order,
-                mode=mode,
-                cval=cval,
-                prefilter=prefilter,
-            )
-
-        if isinstance(_ref, ImageGrid):  # If reference is grid, reshape
-            newdata = resampled.reshape(_ref.shape + (len(self),))
-            moved = spatialimage.__class__(newdata, _ref.affine, spatialimage.header)
-            moved.header.set_data_dtype(output_dtype)
-            return moved
-
-        return resampled
-
 
 def load(filename, fmt=None, reference=None, moving=None):
     """
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index 88f7e6b..1a1e223 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -92,14 +92,14 @@ def apply(
         _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
     )
 
-    #if data.ndim < targets.shape[-1]:
-    #    data = data[..., np.newaxis]
+    if data.ndim == 4 and data.shape[-1] != len(transform):
+        raise ValueError("The fourth dimension of the data does not match the tranform's shape.")
+
+    if data.ndim < transform.ndim:
+        data = data[..., np.newaxis]
  
-    #import pdb; pdb.set_trace()
     resampled = ndi.map_coordinates(
         data,
-        #targets.T,
-        #Reshape targets (516096, 3, 8) --> (4, 4128768) : 
         _as_homogeneous(targets.reshape(-2, targets.shape[0])).T,
         output=output_dtype,
         order=order,
diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py
index b7f7384..50cc537 100644
--- a/nitransforms/tests/test_linear.py
+++ b/nitransforms/tests/test_linear.py
@@ -354,7 +354,6 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path):
         np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz"
     )
 
-    import pdb; pdb.set_trace()
     nii = apply(
         hmcinv, testdata_path / "fmap.nii.gz", order=1
     )

From c5b86e1d5bae736d75b6f0cf486268d161801b50 Mon Sep 17 00:00:00 2001
From: Julien Marabotto
 <julienmarabotto@macbook-pro-de-henri.vpn.intranet.chuv>
Date: Mon, 22 Apr 2024 14:50:33 +0200
Subject: [PATCH 08/15] ENH: update outsoucre apply

---
 nitransforms/resampling.py      |  1 +
 nitransforms/tests/test_base.py | 17 ++++++++++-------
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index 1a1e223..fa02c2e 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -98,6 +98,7 @@ def apply(
     if data.ndim < transform.ndim:
         data = data[..., np.newaxis]
  
+    import pdb; pdb.set_trace()
     resampled = ndi.map_coordinates(
         data,
         _as_homogeneous(targets.reshape(-2, targets.shape[0])).T,
diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py
index a1402ba..fe7855a 100644
--- a/nitransforms/tests/test_base.py
+++ b/nitransforms/tests/test_base.py
@@ -4,7 +4,7 @@
 import pytest
 import h5py
 
-from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase
+from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase, _as_homogeneous
 from .. import linear as nitl
 from ..resampling import apply
 
@@ -42,11 +42,13 @@ def test_ImageGrid(get_testdata, image_orientation):
     # Test ras2vox and vox2ras conversions
     ijk = [[10, 10, 10], [40, 4, 20], [0, 0, 0], [s - 1 for s in im.shape[:3]]]
     xyz = [img._affine.dot(idx + [1])[:-1] for idx in ijk]
+    # xyz = np.array([np.tensordot(img._affine, idx + [1], axes=1)[:-1] for idx in ijk])
 
-    assert np.allclose(img.ras(ijk[0]), xyz[0])
+    # import pdb; pdb.set_trace()
+    assert np.allclose(np.squeeze(img.ras(ijk[0])), xyz[0])
     assert np.allclose(np.round(img.index(xyz[0])), ijk[0])
-    assert np.allclose(img.ras(ijk), xyz)
-    assert np.allclose(np.round(img.index(xyz)), ijk)
+    assert np.allclose(img.ras(ijk).T, xyz)
+    assert np.allclose(np.round(img.index(xyz)).T, ijk)
 
     # nd index / coords
     idxs = img.ndindex
@@ -92,12 +94,13 @@ def _to_hdf5(klass, x5_root):
     img = nb.load(fname)
     imgdata = np.asanyarray(img.dataobj, dtype=img.get_data_dtype())
 
-    # Test identity transform
     xfm = TransformBase()
-    xfm.reference = fname
-
     with pytest.raises(TypeError):
         _ = xfm.ndim
+
+    # Test identity transform
+    xfm = nitl.Affine()
+    xfm.reference = fname
     moved = apply(xfm, fname, order=0)
     assert np.all(
         imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())

From 95a215edf5e3e93a8a02669900c0167c34034f37 Mon Sep 17 00:00:00 2001
From: Julien Marabotto <julienmarabotto@MacBook-Pro-de-Henri.local>
Date: Thu, 25 Apr 2024 15:33:41 +0200
Subject: [PATCH 09/15] Updated: offsource apply

---
 nitransforms/nonlinear.py            | 13 +++++++++----
 nitransforms/resampling.py           |  6 ++++--
 nitransforms/tests/test_base.py      | 17 +++++------------
 nitransforms/tests/test_nonlinear.py |  5 +++--
 4 files changed, 21 insertions(+), 20 deletions(-)

diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 79c3aa4..52d854e 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -14,6 +14,7 @@
 from nitransforms import io
 from nitransforms.io.base import _ensure_image
 from nitransforms.interp.bspline import grid_bspline_weights, _cubic_bspline
+from nitransforms.resampling import apply
 from nitransforms.base import (
     TransformBase,
     TransformError,
@@ -257,7 +258,7 @@ def __init__(self, coefficients, reference=None, order=3):
         if reference is not None:
             self.reference = reference
 
-            if coefficients.shape[-1] != self.ndim:
+            if coefficients.shape[-1] != self.reference.ndim:
                 raise TransformError(
                     'Number of components of the coefficients does '
                     'not match the number of dimensions')
@@ -310,19 +311,23 @@ def apply(
         spatialimage = _ensure_image(spatialimage)
 
         # If locations to be interpolated are not on a grid, run map()
+        #import pdb; pdb.set_trace()
         if not isinstance(_ref, ImageGrid):
-            return super().apply(
+            return apply(
+                super(),
                 spatialimage,
                 reference=_ref,
+                output_dtype=output_dtype,
                 order=order,
                 mode=mode,
                 cval=cval,
                 prefilter=prefilter,
-                output_dtype=output_dtype,
+                
             )
 
         # If locations to be interpolated are on a grid, generate a displacements field
-        return self.to_field(reference=reference).apply(
+        return apply(
+            self.to_field(reference=reference),
             spatialimage,
             reference=reference,
             order=order,
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index fa02c2e..7cbdd9b 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -98,10 +98,12 @@ def apply(
     if data.ndim < transform.ndim:
         data = data[..., np.newaxis]
  
-    import pdb; pdb.set_trace()
+    if transform.ndim == 4:
+        targets = _as_homogeneous(targets.reshape(-2, targets.shape[0])).T
+
     resampled = ndi.map_coordinates(
         data,
-        _as_homogeneous(targets.reshape(-2, targets.shape[0])).T,
+        targets,
         output=output_dtype,
         order=order,
         mode=mode,
diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py
index fe7855a..fe08d2e 100644
--- a/nitransforms/tests/test_base.py
+++ b/nitransforms/tests/test_base.py
@@ -94,10 +94,14 @@ def _to_hdf5(klass, x5_root):
     img = nb.load(fname)
     imgdata = np.asanyarray(img.dataobj, dtype=img.get_data_dtype())
 
+    # Test identity transform - setting reference
     xfm = TransformBase()
     with pytest.raises(TypeError):
         _ = xfm.ndim
 
+    # Test to_filename
+    xfm.to_filename("data.x5")
+
     # Test identity transform
     xfm = nitl.Affine()
     xfm.reference = fname
@@ -106,17 +110,6 @@ def _to_hdf5(klass, x5_root):
         imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())
     )
 
-    # Test identity transform - setting reference
-    xfm = TransformBase()
-    xfm.reference = fname
-
-    with pytest.raises(TypeError):
-        _ = xfm.ndim
-    moved = apply(xfm, str(fname), reference=fname, order=0)
-    assert np.all(
-        imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())
-    )
-
     # Test ndim returned by affine
     assert nitl.Affine().ndim == 3
     assert nitl.LinearTransformsMapping(
@@ -136,7 +129,7 @@ def _to_hdf5(klass, x5_root):
     assert np.allclose(giimoved.reshape(xfm.reference.shape), moved.get_fdata())
 
     # Test to_filename
-    xfm.to_filename("data.x5")
+    xfm.to_filename("data.xfm", fmt='itk')
 
 
 def test_SampledSpatialData(testdata_path):
diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py
index 93d3fd4..dd4cbf9 100644
--- a/nitransforms/tests/test_nonlinear.py
+++ b/nitransforms/tests/test_nonlinear.py
@@ -8,6 +8,7 @@
 
 import numpy as np
 import nibabel as nb
+from nitransforms.resampling import apply
 from nitransforms.base import TransformError
 from nitransforms.io.base import TransformFileError
 from nitransforms.nonlinear import (
@@ -247,8 +248,8 @@ def test_bspline(tmp_path, testdata_path):
     bsplxfm = BSplineFieldTransform(bs_name, reference=img_name)
     dispxfm = DenseFieldTransform(disp_name)
 
-    out_disp = dispxfm.apply(img_name)
-    out_bspl = bsplxfm.apply(img_name)
+    out_disp = apply(dispxfm,img_name)
+    out_bspl = apply(bsplxfm,img_name)
 
     out_disp.to_filename("resampled_field.nii.gz")
     out_bspl.to_filename("resampled_bsplines.nii.gz")

From 9f93a67177504ad25fd96ef69c7c0af2133dffee Mon Sep 17 00:00:00 2001
From: Julien Marabotto <julienmarabotto@MacBook-Pro-de-Henri.local>
Date: Fri, 26 Apr 2024 10:35:37 +0200
Subject: [PATCH 10/15] enh: removed straneous comments, update nonlinear

---
 nitransforms/nonlinear.py       | 2 +-
 nitransforms/tests/test_base.py | 2 --
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 52d854e..488d01c 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -168,7 +168,7 @@ def map(self, x, inverse=False):
         indexes = np.round(ijk).astype("int")
 
         if np.all(np.abs(ijk - indexes) < 1e-3):
-            indexes = tuple(tuple(i) for i in indexes.T)
+            indexes = tuple(tuple(i) for i in indexes)
             return self._field[indexes]
 
         return np.vstack(tuple(
diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py
index fe08d2e..4a34526 100644
--- a/nitransforms/tests/test_base.py
+++ b/nitransforms/tests/test_base.py
@@ -42,9 +42,7 @@ def test_ImageGrid(get_testdata, image_orientation):
     # Test ras2vox and vox2ras conversions
     ijk = [[10, 10, 10], [40, 4, 20], [0, 0, 0], [s - 1 for s in im.shape[:3]]]
     xyz = [img._affine.dot(idx + [1])[:-1] for idx in ijk]
-    # xyz = np.array([np.tensordot(img._affine, idx + [1], axes=1)[:-1] for idx in ijk])
 
-    # import pdb; pdb.set_trace()
     assert np.allclose(np.squeeze(img.ras(ijk[0])), xyz[0])
     assert np.allclose(np.round(img.index(xyz[0])), ijk[0])
     assert np.allclose(img.ras(ijk).T, xyz)

From f59720d2f2dae6e166c9e8416b3080aa6864ff28 Mon Sep 17 00:00:00 2001
From: Julien Marabotto <julienmarabotto@MacBook-Pro-de-Henri.local>
Date: Thu, 2 May 2024 14:42:23 +0200
Subject: [PATCH 11/15] FIX: Offsource Apply

Apply function offsourced. Tests:  139 passed, 163 Skipped, 15 Warnings
---
 nitransforms/nonlinear.py            | 67 +++++++---------------------
 nitransforms/resampling.py           | 13 ++++--
 nitransforms/tests/test_nonlinear.py | 13 +++---
 3 files changed, 31 insertions(+), 62 deletions(-)

diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 488d01c..9b67b81 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -30,6 +30,11 @@ class DenseFieldTransform(TransformBase):
 
     __slots__ = ("_field", "_deltas")
 
+    @property
+    def ndim(self):
+        """Access the dimensions of this Desne Field Transform."""
+        return self._field.ndim - 1
+
     def __init__(self, field=None, is_deltas=True, reference=None):
         """
         Create a dense field transform.
@@ -82,11 +87,10 @@ def __init__(self, field=None, is_deltas=True, reference=None):
                 "Reference is not a spatial image"
             )
 
-        ndim = self._field.ndim - 1
-        if self._field.shape[-1] != ndim:
+        if self._field.shape[-1] != self.ndim:
             raise TransformError(
                 "The number of components of the field (%d) does not match "
-                "the number of dimensions (%d)" % (self._field.shape[-1], ndim)
+                "the number of dimensions (%d)" % (self._field.shape[-1], self.ndim)
             )
 
         if is_deltas:
@@ -245,6 +249,12 @@ class BSplineFieldTransform(TransformBase):
 
     __slots__ = ['_coeffs', '_knots', '_weights', '_order', '_moving']
 
+    @property
+    def ndim(self):
+        """Access the dimensions of this BSpline."""
+        #return ndim = self._coeffs.shape[-1]
+        return self._coeffs.ndim - 1
+
     def __init__(self, coefficients, reference=None, order=3):
         """Create a smooth deformation field using B-Spline basis."""
         super().__init__()
@@ -277,14 +287,12 @@ def to_field(self, reference=None, dtype="float32"):
         if _ref is None:
             raise TransformError("A reference must be defined")
 
-        ndim = self._coeffs.shape[-1]
-
         if self._weights is None:
             self._weights = grid_bspline_weights(_ref, self._knots)
 
-        field = np.zeros((_ref.npoints, ndim))
+        field = np.zeros((_ref.npoints, self.ndim))
 
-        for d in range(ndim):
+        for d in range(self.ndim):
             #  1 x Nvox :                          (1 x K) @ (K x Nvox)
             field[:, d] = self._coeffs[..., d].reshape(-1) @ self._weights
 
@@ -292,51 +300,6 @@ def to_field(self, reference=None, dtype="float32"):
             field.astype(dtype).reshape(*_ref.shape, -1), reference=_ref
         )
 
-    def apply(
-        self,
-        spatialimage,
-        reference=None,
-        order=3,
-        mode="constant",
-        cval=0.0,
-        prefilter=True,
-        output_dtype=None,
-    ):
-        """Apply a B-Spline transform on input data."""
-
-        _ref = (
-            self.reference if reference is None else
-            SpatialReference.factory(_ensure_image(reference))
-        )
-        spatialimage = _ensure_image(spatialimage)
-
-        # If locations to be interpolated are not on a grid, run map()
-        #import pdb; pdb.set_trace()
-        if not isinstance(_ref, ImageGrid):
-            return apply(
-                super(),
-                spatialimage,
-                reference=_ref,
-                output_dtype=output_dtype,
-                order=order,
-                mode=mode,
-                cval=cval,
-                prefilter=prefilter,
-                
-            )
-
-        # If locations to be interpolated are on a grid, generate a displacements field
-        return apply(
-            self.to_field(reference=reference),
-            spatialimage,
-            reference=reference,
-            order=order,
-            mode=mode,
-            cval=cval,
-            prefilter=prefilter,
-            output_dtype=output_dtype,
-        )
-
     def map(self, x, inverse=False):
         r"""
         Apply the transformation to a list of physical coordinate points.
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index 7cbdd9b..942ab07 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -87,16 +87,21 @@ def apply(
         spatialimage = _nbload(str(spatialimage))
 
     data = np.asanyarray(spatialimage.dataobj)
-    
-    targets = ImageGrid(spatialimage).index(  # data should be an image
-        _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
-    )
 
     if data.ndim == 4 and data.shape[-1] != len(transform):
         raise ValueError("The fourth dimension of the data does not match the tranform's shape.")
 
     if data.ndim < transform.ndim:
         data = data[..., np.newaxis]
+    
+    if hasattr(transform, 'to_field') and callable(transform.to_field):
+        targets = ImageGrid(spatialimage).index(
+            _as_homogeneous(transform.to_field(reference=reference).map(_ref.ndcoords.T), dim=_ref.ndim)
+        )
+    else:
+        targets = ImageGrid(spatialimage).index(  # data should be an image
+            _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
+        )
  
     if transform.ndim == 4:
         targets = _as_homogeneous(targets.reshape(-2, targets.shape[0])).T
diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py
index dd4cbf9..4a802b5 100644
--- a/nitransforms/tests/test_nonlinear.py
+++ b/nitransforms/tests/test_nonlinear.py
@@ -97,13 +97,14 @@ def test_bsplines_references(testdata_path):
         ).to_field()
 
     with pytest.raises(TransformError):
-        BSplineFieldTransform(
-            testdata_path / "someones_bspline_coefficients.nii.gz"
-        ).apply(testdata_path / "someones_anatomy.nii.gz")
+        apply(
+            BSplineFieldTransform(testdata_path / "someones_bspline_coefficients.nii.gz"),
+            testdata_path / "someones_anatomy.nii.gz",
+        )
 
-    BSplineFieldTransform(
-        testdata_path / "someones_bspline_coefficients.nii.gz"
-    ).apply(
+    apply(
+        BSplineFieldTransform(
+        testdata_path / "someones_bspline_coefficients.nii.gz"),
         testdata_path / "someones_anatomy.nii.gz",
         reference=testdata_path / "someones_anatomy.nii.gz"
     )

From 0837e912b710e107f7a7b27970b9ec5b142a1c2d Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Thu, 16 May 2024 08:05:19 -0400
Subject: [PATCH 12/15] sty: pacify flake8

---
 nitransforms/base.py                 |  6 ++-
 nitransforms/linear.py               |  6 +--
 nitransforms/nonlinear.py            | 57 ++++++++++++++--------------
 nitransforms/resampling.py           | 23 +++++++----
 nitransforms/tests/test_base.py      | 17 +++++----
 nitransforms/tests/test_nonlinear.py | 49 +++++++++++++-----------
 6 files changed, 85 insertions(+), 73 deletions(-)

diff --git a/nitransforms/base.py b/nitransforms/base.py
index 6a6ae7e..26c0d47 100644
--- a/nitransforms/base.py
+++ b/nitransforms/base.py
@@ -177,7 +177,10 @@ def __ne__(self, other):
 class TransformBase:
     """Abstract image class to represent transforms."""
 
-    __slots__ = ("_reference", "_ndim",)
+    __slots__ = (
+        "_reference",
+        "_ndim",
+    )
 
     def __init__(self, reference=None):
         """Instantiate a transform."""
@@ -283,6 +286,7 @@ def _as_homogeneous(xyz, dtype="float32", dim=3):
 
     return np.hstack((xyz, np.ones((xyz.shape[0], 1), dtype=dtype)))
 
+
 def _apply_affine(x, affine, dim):
     """Get the image array's indexes corresponding to coordinates."""
     return np.tensordot(
diff --git a/nitransforms/linear.py b/nitransforms/linear.py
index cf16104..71df6a1 100644
--- a/nitransforms/linear.py
+++ b/nitransforms/linear.py
@@ -10,16 +10,12 @@
 import warnings
 import numpy as np
 from pathlib import Path
-from scipy import ndimage as ndi
 
-from nibabel.loadsave import load as _nbload
 from nibabel.affines import from_matvec
-from nibabel.arrayproxy import get_obj_dtype
 
 from nitransforms.base import (
     ImageGrid,
     TransformBase,
-    SpatialReference,
     _as_homogeneous,
     EQUALITY_TOL,
 )
@@ -112,7 +108,7 @@ def __invert__(self):
 
         """
         return self.__class__(self._inverse)
-    
+
     def __len__(self):
         """Enable using len()."""
         return 1 if self._matrix.ndim == 2 else len(self._matrix)
diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 9b67b81..93f891f 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -14,12 +14,10 @@
 from nitransforms import io
 from nitransforms.io.base import _ensure_image
 from nitransforms.interp.bspline import grid_bspline_weights, _cubic_bspline
-from nitransforms.resampling import apply
 from nitransforms.base import (
     TransformBase,
     TransformError,
     ImageGrid,
-    SpatialReference,
     _as_homogeneous,
 )
 from scipy.ndimage import map_coordinates
@@ -77,14 +75,12 @@ def __init__(self, field=None, is_deltas=True, reference=None):
             is_deltas = True
 
         try:
-            self.reference = ImageGrid(
-                reference if reference is not None else field
-            )
+            self.reference = ImageGrid(reference if reference is not None else field)
         except AttributeError:
             raise TransformError(
                 "Field must be a spatial image if reference is not provided"
-                if reference is None else
-                "Reference is not a spatial image"
+                if reference is None
+                else "Reference is not a spatial image"
             )
 
         if self._field.shape[-1] != self.ndim:
@@ -175,16 +171,19 @@ def map(self, x, inverse=False):
             indexes = tuple(tuple(i) for i in indexes)
             return self._field[indexes]
 
-        return np.vstack(tuple(
-            map_coordinates(
-                self._field[..., i],
-                ijk.T,
-                order=3,
-                mode="constant",
-                cval=0,
-                prefilter=True,
-            ) for i in range(self.reference.ndim)
-        )).T
+        return np.vstack(
+            tuple(
+                map_coordinates(
+                    self._field[..., i],
+                    ijk.T,
+                    order=3,
+                    mode="constant",
+                    cval=0,
+                    prefilter=True,
+                )
+                for i in range(self.reference.ndim)
+            )
+        ).T
 
     def __matmul__(self, b):
         """
@@ -206,9 +205,9 @@ def __matmul__(self, b):
         True
 
         """
-        retval = b.map(
-            self._field.reshape((-1, self._field.shape[-1]))
-        ).reshape(self._field.shape)
+        retval = b.map(self._field.reshape((-1, self._field.shape[-1]))).reshape(
+            self._field.shape
+        )
         return DenseFieldTransform(retval, is_deltas=False, reference=self.reference)
 
     def __eq__(self, other):
@@ -247,12 +246,12 @@ def from_filename(cls, filename, fmt="X5"):
 class BSplineFieldTransform(TransformBase):
     """Represent a nonlinear transform parameterized by BSpline basis."""
 
-    __slots__ = ['_coeffs', '_knots', '_weights', '_order', '_moving']
+    __slots__ = ["_coeffs", "_knots", "_weights", "_order", "_moving"]
 
     @property
     def ndim(self):
         """Access the dimensions of this BSpline."""
-        #return ndim = self._coeffs.shape[-1]
+        # return ndim = self._coeffs.shape[-1]
         return self._coeffs.ndim - 1
 
     def __init__(self, coefficients, reference=None, order=3):
@@ -270,8 +269,9 @@ def __init__(self, coefficients, reference=None, order=3):
 
             if coefficients.shape[-1] != self.reference.ndim:
                 raise TransformError(
-                    'Number of components of the coefficients does '
-                    'not match the number of dimensions')
+                    "Number of components of the coefficients does "
+                    "not match the number of dimensions"
+                )
 
     @property
     def ndim(self):
@@ -281,8 +281,7 @@ def ndim(self):
     def to_field(self, reference=None, dtype="float32"):
         """Generate a displacements deformation field from this B-Spline field."""
         _ref = (
-            self.reference if reference is None else
-            ImageGrid(_ensure_image(reference))
+            self.reference if reference is None else ImageGrid(_ensure_image(reference))
         )
         if _ref is None:
             raise TransformError("A reference must be defined")
@@ -350,9 +349,9 @@ def _map_xyz(x, reference, knots, coeffs):
     # Probably this will change if the order of the B-Spline is different
     w_start, w_end = np.ceil(ijk - 2).astype(int), np.floor(ijk + 2).astype(int)
     # Generate a grid of indexes corresponding to the window
-    nonzero_knots = tuple([
-        np.arange(start, end + 1) for start, end in zip(w_start, w_end)
-    ])
+    nonzero_knots = tuple(
+        [np.arange(start, end + 1) for start, end in zip(w_start, w_end)]
+    )
     nonzero_knots = tuple(np.meshgrid(*nonzero_knots, indexing="ij"))
     window = np.array(nonzero_knots).reshape((ndim, -1))
 
diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index 942ab07..e1ac154 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -77,7 +77,9 @@ def apply(
         reference = _nbload(str(reference))
 
     _ref = (
-        transform.reference if reference is None else SpatialReference.factory(reference)
+        transform.reference
+        if reference is None
+        else SpatialReference.factory(reference)
     )
 
     if _ref is None:
@@ -89,20 +91,25 @@ def apply(
     data = np.asanyarray(spatialimage.dataobj)
 
     if data.ndim == 4 and data.shape[-1] != len(transform):
-        raise ValueError("The fourth dimension of the data does not match the tranform's shape.")
+        raise ValueError(
+            "The fourth dimension of the data does not match the tranform's shape."
+        )
 
     if data.ndim < transform.ndim:
         data = data[..., np.newaxis]
-    
-    if hasattr(transform, 'to_field') and callable(transform.to_field):
+
+    if hasattr(transform, "to_field") and callable(transform.to_field):
         targets = ImageGrid(spatialimage).index(
-            _as_homogeneous(transform.to_field(reference=reference).map(_ref.ndcoords.T), dim=_ref.ndim)
+            _as_homogeneous(
+                transform.to_field(reference=reference).map(_ref.ndcoords.T),
+                dim=_ref.ndim,
+            )
         )
     else:
         targets = ImageGrid(spatialimage).index(  # data should be an image
             _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim)
         )
- 
+
     if transform.ndim == 4:
         targets = _as_homogeneous(targets.reshape(-2, targets.shape[0])).T
 
@@ -115,14 +122,14 @@ def apply(
         cval=cval,
         prefilter=prefilter,
     )
-     
+
     if isinstance(_ref, ImageGrid):  # If reference is grid, reshape
         hdr = None
         if _ref.header is not None:
             hdr = _ref.header.copy()
             hdr.set_data_dtype(output_dtype or spatialimage.get_data_dtype())
         moved = spatialimage.__class__(
-            resampled.reshape(_ref.shape if data.ndim < 4 else _ref.shape + (-1, )),
+            resampled.reshape(_ref.shape if data.ndim < 4 else _ref.shape + (-1,)),
             _ref.affine,
             hdr,
         )
diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py
index 4a34526..d32ce7f 100644
--- a/nitransforms/tests/test_base.py
+++ b/nitransforms/tests/test_base.py
@@ -4,7 +4,12 @@
 import pytest
 import h5py
 
-from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase, _as_homogeneous
+from ..base import (
+    SpatialReference,
+    SampledSpatialData,
+    ImageGrid,
+    TransformBase,
+)
 from .. import linear as nitl
 from ..resampling import apply
 
@@ -104,15 +109,11 @@ def _to_hdf5(klass, x5_root):
     xfm = nitl.Affine()
     xfm.reference = fname
     moved = apply(xfm, fname, order=0)
-    assert np.all(
-        imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())
-    )
+    assert np.all(imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()))
 
     # Test ndim returned by affine
     assert nitl.Affine().ndim == 3
-    assert nitl.LinearTransformsMapping(
-        [nitl.Affine(), nitl.Affine()]
-    ).ndim == 4
+    assert nitl.LinearTransformsMapping([nitl.Affine(), nitl.Affine()]).ndim == 4
 
     # Test applying to Gifti
     gii = nb.gifti.GiftiImage(
@@ -127,7 +128,7 @@ def _to_hdf5(klass, x5_root):
     assert np.allclose(giimoved.reshape(xfm.reference.shape), moved.get_fdata())
 
     # Test to_filename
-    xfm.to_filename("data.xfm", fmt='itk')
+    xfm.to_filename("data.xfm", fmt="itk")
 
 
 def test_SampledSpatialData(testdata_path):
diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py
index 4a802b5..cfaa12c 100644
--- a/nitransforms/tests/test_nonlinear.py
+++ b/nitransforms/tests/test_nonlinear.py
@@ -29,7 +29,7 @@
 3dNwarpApply -nwarp {transform} -source {moving} \
 -master {reference} -interp NN -prefix {output} {extra}\
 """.format,
-    'fsl': """\
+    "fsl": """\
 applywarp -i {moving} -r {reference} -o {output} {extra}\
 -w {transform} --interp=nn""".format,
 }
@@ -39,7 +39,9 @@
 def test_itk_disp_load(size):
     """Checks field sizes."""
     with pytest.raises(TransformFileError):
-        ITKDisplacementsField.from_image(nb.Nifti1Image(np.zeros(size), np.eye(4), None))
+        ITKDisplacementsField.from_image(
+            nb.Nifti1Image(np.zeros(size), np.eye(4), None)
+        )
 
 
 @pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 2, 3), (20, 20, 20, 1, 4)])
@@ -98,15 +100,16 @@ def test_bsplines_references(testdata_path):
 
     with pytest.raises(TransformError):
         apply(
-            BSplineFieldTransform(testdata_path / "someones_bspline_coefficients.nii.gz"),
+            BSplineFieldTransform(
+                testdata_path / "someones_bspline_coefficients.nii.gz"
+            ),
             testdata_path / "someones_anatomy.nii.gz",
         )
 
     apply(
-        BSplineFieldTransform(
-        testdata_path / "someones_bspline_coefficients.nii.gz"),
+        BSplineFieldTransform(testdata_path / "someones_bspline_coefficients.nii.gz"),
         testdata_path / "someones_anatomy.nii.gz",
-        reference=testdata_path / "someones_anatomy.nii.gz"
+        reference=testdata_path / "someones_anatomy.nii.gz",
     )
 
 
@@ -170,7 +173,7 @@ def test_displacements_field1(
     nt_moved_mask.set_data_dtype(msk.get_data_dtype())
     diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj)
 
-    assert np.sqrt((diff ** 2).mean()) < RMSE_TOL
+    assert np.sqrt((diff**2).mean()) < RMSE_TOL
     brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool)
 
     # Then apply the transform and cross-check with software
@@ -179,7 +182,7 @@ def test_displacements_field1(
         reference=tmp_path / "reference.nii.gz",
         moving=tmp_path / "reference.nii.gz",
         output=tmp_path / "resampled.nii.gz",
-        extra="--output-data-type uchar" if sw_tool == "itk" else ""
+        extra="--output-data-type uchar" if sw_tool == "itk" else "",
     )
 
     exit_code = check_call([cmd], shell=True)
@@ -190,10 +193,9 @@ def test_displacements_field1(
     nt_moved.set_data_dtype(nii.get_data_dtype())
     nt_moved.to_filename("nt_resampled.nii.gz")
     sw_moved.set_data_dtype(nt_moved.get_data_dtype())
-    diff = (
-        np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype())
-        - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
-    )
+    diff = np.asanyarray(
+        sw_moved.dataobj, dtype=sw_moved.get_data_dtype()
+    ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
     # A certain tolerance is necessary because of resampling at borders
     assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL
 
@@ -230,12 +232,11 @@ def test_displacements_field2(tmp_path, testdata_path, sw_tool):
     nt_moved = xfm.apply(img_fname, order=0)
     nt_moved.to_filename("nt_resampled.nii.gz")
     sw_moved.set_data_dtype(nt_moved.get_data_dtype())
-    diff = (
-        np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype())
-        - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
-    )
+    diff = np.asanyarray(
+        sw_moved.dataobj, dtype=sw_moved.get_data_dtype()
+    ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
     # A certain tolerance is necessary because of resampling at borders
-    assert np.sqrt((diff ** 2).mean()) < RMSE_TOL
+    assert np.sqrt((diff**2).mean()) < RMSE_TOL
 
 
 def test_bspline(tmp_path, testdata_path):
@@ -249,12 +250,16 @@ def test_bspline(tmp_path, testdata_path):
     bsplxfm = BSplineFieldTransform(bs_name, reference=img_name)
     dispxfm = DenseFieldTransform(disp_name)
 
-    out_disp = apply(dispxfm,img_name)
-    out_bspl = apply(bsplxfm,img_name)
+    out_disp = apply(dispxfm, img_name)
+    out_bspl = apply(bsplxfm, img_name)
 
     out_disp.to_filename("resampled_field.nii.gz")
     out_bspl.to_filename("resampled_bsplines.nii.gz")
 
-    assert np.sqrt(
-        (out_disp.get_fdata(dtype="float32") - out_bspl.get_fdata(dtype="float32")) ** 2
-    ).mean() < 0.2
+    assert (
+        np.sqrt(
+            (out_disp.get_fdata(dtype="float32") - out_bspl.get_fdata(dtype="float32"))
+            ** 2
+        ).mean()
+        < 0.2
+    )

From c7a958cf50f11110ca463288b631d1df9d200d77 Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Thu, 16 May 2024 08:18:54 -0400
Subject: [PATCH 13/15] fix: remove double definition

---
 nitransforms/nonlinear.py | 11 -----------
 1 file changed, 11 deletions(-)

diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 93f891f..db29633 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -28,11 +28,6 @@ class DenseFieldTransform(TransformBase):
 
     __slots__ = ("_field", "_deltas")
 
-    @property
-    def ndim(self):
-        """Access the dimensions of this Desne Field Transform."""
-        return self._field.ndim - 1
-
     def __init__(self, field=None, is_deltas=True, reference=None):
         """
         Create a dense field transform.
@@ -248,12 +243,6 @@ class BSplineFieldTransform(TransformBase):
 
     __slots__ = ["_coeffs", "_knots", "_weights", "_order", "_moving"]
 
-    @property
-    def ndim(self):
-        """Access the dimensions of this BSpline."""
-        # return ndim = self._coeffs.shape[-1]
-        return self._coeffs.ndim - 1
-
     def __init__(self, coefficients, reference=None, order=3):
         """Create a smooth deformation field using B-Spline basis."""
         super().__init__()

From bfe592d59b666bea9d9ae9564bae81cb5a5e4fa5 Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <effigies@gmail.com>
Date: Thu, 16 May 2024 14:38:53 -0400
Subject: [PATCH 14/15] Fix bad merge

---
 nitransforms/nonlinear.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py
index 0aa6d36..f4b9514 100644
--- a/nitransforms/nonlinear.py
+++ b/nitransforms/nonlinear.py
@@ -166,7 +166,7 @@ def map(self, x, inverse=False):
             indexes = tuple(tuple(i) for i in indexes)
             return self._field[indexes]
 
-        return np.vstack(
+        new_map = np.vstack(
             tuple(
                 map_coordinates(
                     self._field[..., i],

From 5b1736bf2bbda7c737e8c7a7ec8806dd7510a1ad Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Fri, 17 May 2024 11:45:54 -0400
Subject: [PATCH 15/15] Update nitransforms/resampling.py

---
 nitransforms/resampling.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py
index e1ac154..9de0d2d 100644
--- a/nitransforms/resampling.py
+++ b/nitransforms/resampling.py
@@ -98,6 +98,7 @@ def apply(
     if data.ndim < transform.ndim:
         data = data[..., np.newaxis]
 
+    # For model-based nonlinear transforms, generate the corresponding dense field
     if hasattr(transform, "to_field") and callable(transform.to_field):
         targets = ImageGrid(spatialimage).index(
             _as_homogeneous(