ENH: Implement Jacobian weighting during unwarp [backport gh-391]

nipreps · Oct 4, 2023 · 6534591 · 6534591
1 parent 2b8400d
commit 6534591
Showing 1 changed file with 5 additions and 2 deletions.
diff --git a/sdcflows/transform.py b/sdcflows/transform.py
@@ -197,7 +197,10 @@ def apply(
         # fmap * ro_time is the voxel-shift map (VSM)
         # The VSM is just the displacements field given in index coordinates
         # voxcoords is the deformation field, i.e., the target position of each voxel
-        voxcoords[pe_axis, ...] += fmap.reshape(-1) * ro_time
+        vsm = fmap * ro_time
+        voxcoords[pe_axis, ...] += vsm.reshape(-1)
+
+        jacobian = 1 + np.gradient(vsm, axis=pe_axis)
 
         # Prepare data
         data = np.squeeze(np.asanyarray(spatialimage.dataobj))
@@ -212,7 +215,7 @@ def apply(
             mode=mode,
             cval=cval,
             prefilter=prefilter,
-        ).reshape(spatialimage.shape)
+        ).reshape(spatialimage.shape) * jacobian
 
         moved = spatialimage.__class__(
             resampled, spatialimage.affine, spatialimage.header