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)

Author: Julien Marabotto

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, ...]

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,
         )

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