FIX: Derive field from SyN displacements using EPI affine

sdcflows/interfaces/fmap.py

@@ -161,6 +161,7 @@ def _run_interface(self, runtime):
 
 class _DisplacementsField2FieldmapInputSpec(BaseInterfaceInputSpec):
     transform = File(exists=True, mandatory=True, desc="input displacements field")
+    epi = File(exists=True, mandatory=True, desc="source EPI image")
     ro_time = traits.Float(mandatory=True, desc="total readout time")
     pe_dir = traits.Enum(
         "j-", "j", "i", "i-", "k", "k-", mandatory=True, desc="phase encoding direction"
@@ -189,6 +190,7 @@ def _run_interface(self, runtime):
         )
         fmapnii = disp_to_fmap(
             nb.load(self.inputs.transform),
+            nb.load(self.inputs.epi),
             ro_time=self.inputs.ro_time,
             pe_dir=self.inputs.pe_dir,
             itk_format=self.inputs.itk_transform,

sdcflows/tests/test_transform.py

@@ -327,11 +327,11 @@ def test_conversions(tmpdir, testdata_dir, pe_dir):
             ro_time=0.2,
             pe_dir=pe_dir,
         ),
+        fmap_nii,
         ro_time=0.2,
         pe_dir=pe_dir,
     )
 
-    new_nii.to_filename("test.nii.gz")
     assert np.allclose(
         fmap_nii.get_fdata(dtype="float32"),
         new_nii.get_fdata(dtype="float32"),

sdcflows/transform.py

@@ -617,20 +617,22 @@ def fmap_to_disp(fmap_nii, ro_time, pe_dir, itk_format=True):
     return xyz_nii
 
 
-def disp_to_fmap(xyz_nii, ro_time, pe_dir, itk_format=True):
+def disp_to_fmap(xyz_nii, epi_nii, ro_time, pe_dir, itk_format=True):
     """
     Convert a displacements field into a fieldmap in Hz.
 
     This is the inverse operation to the previous function.
 
     Parameters
     ----------
-    xyz_nii : :obj:`os.pathlike`
-        Path to a displacements field in NIfTI format.
+    xyz_nii : :class:`nibabel.Nifti1Image`
+        Displacements field in NIfTI format.
+    epi_nii : :class:`nibabel.Nifti1Image`
+        Source EPI image.
     ro_time : :obj:`float`
-        The total readout time in seconds.
+        The total readout time of the EPI image in seconds.
     pe_dir : :obj:`str`
-        The ``PhaseEncodingDirection`` metadata value.
+        The ``PhaseEncodingDirection`` metadata value of the EPI image.
 
     Returns
     -------
@@ -644,12 +646,13 @@ def disp_to_fmap(xyz_nii, ro_time, pe_dir, itk_format=True):
         # ITK displacement vectors are in LPS orientation
         xyz_deltas[:, (0, 1)] *= -1
 
-    inv_aff = np.linalg.inv(xyz_nii.affine)
-    inv_aff[:3, 3] = 0  # Translations MUST NOT be applied.
+    # Use the EPI's affine to determine voxel spacing, axis ordering and flips
+    inv_aff = np.linalg.inv(epi_nii.affine)
+    inv_mat = inv_aff[:3, :3]
 
     # Convert displacements from mm to voxel units
     # Using the inverse affine accounts for reordering of axes, etc.
-    ijk_deltas = nb.affines.apply_affine(inv_aff, xyz_deltas).astype("float32")
+    ijk_deltas = (xyz_deltas @ inv_mat.T).astype("float32")
     pe_axis = "ijk".index(pe_dir[0])
     vsm = ijk_deltas[:, pe_axis].reshape(xyz_nii.shape[:3])
     scale_factor = -ro_time if pe_dir.endswith("-") else ro_time

sdcflows/workflows/fit/syn.py

@@ -247,7 +247,7 @@ def init_syn_sdc_wf(
 
     # Extract the corresponding fieldmap in Hz
     extract_field = pe.Node(
-        DisplacementsField2Fieldmap(demean=True), name="extract_field"
+        DisplacementsField2Fieldmap(), name="extract_field"
     )
 
     unwarp = pe.Node(ApplyCoeffsField(), name="unwarp")
@@ -267,14 +267,15 @@ def init_syn_sdc_wf(
     )
 
     # Regularize with B-Splines
-    bs_filter = pe.Node(BSplineApprox(), name="bs_filter")
+    bs_filter = pe.Node(
+        BSplineApprox(recenter=False, debug=debug, extrapolate=not debug),
+        name="bs_filter",
+    )
     bs_filter.interface._always_run = debug
     bs_filter.inputs.bs_spacing = (
         [DEFAULT_LF_ZOOMS_MM, DEFAULT_HF_ZOOMS_MM] if not sloppy else [DEFAULT_ZOOMS_MM]
     )
-    bs_filter.inputs.extrapolate = not debug
 
-    # fmt: off
     workflow.connect([
         (inputnode, readout_time, [(("epi_ref", _pop), "in_file"),
                                    (("epi_ref", _pull), "metadata")]),
@@ -314,6 +315,7 @@ def init_syn_sdc_wf(
         (epi_merge, syn, [("out", "moving_image")]),
         (moving_masks, syn, [("out", "moving_image_masks")]),
         (syn, extract_field, [(("forward_transforms", _pop), "transform")]),
+        (clip_epi, extract_field, [("out_file", "epi")]),
         (readout_time, extract_field, [("readout_time", "ro_time"),
                                        ("pe_direction", "pe_dir")]),
         (extract_field, zooms_field, [("out_file", "input_image")]),
@@ -329,8 +331,7 @@ def init_syn_sdc_wf(
         (bs_filter, outputnode, [("out_coeff", "fmap_coeff")]),
         (unwarp, outputnode, [("out_corrected", "fmap_ref"),
                               ("out_field", "fmap")]),
-    ])
-    # fmt: on
+    ])  # fmt:skip
 
     return workflow