FIX: Derive field from SyN displacements using EPI affine (#421)

Theoretically, a fieldmap in Hz is the number of voxels/second that signal is shifted from the true location. SyN-SDC calculates the shift in mm from the true location. In order to find the equivalent fieldmap, we need the size of the voxels, readout time and phase-encoding direction of the EPI image that is to be corrected.

Previously, we were deriving the fieldmap from the readout time and phase-encoding direction of the EPI image, but the voxel size of the displacement field. This causes an overestimate of the number of voxels to be shifted, and thus exaggerated correction.

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