Merge pull request #239 from nipreps/fix/syn-revision

FIX: SDC-SyN ("fieldmap-less") overhaul

Author: oesteban

sdcflows/data/affine.json

@@ -1,22 +1,24 @@
 {
     "collapse_output_transforms": true,
-    "convergence_threshold": [ 1e-06, 1e-06 ],
-    "convergence_window_size": [ 10, 10 ],
+    "convergence_threshold": [ 1e-06, 1e-07 ],
+    "convergence_window_size": [ 10, 5 ],
     "dimension": 3,
     "initial_moving_transform_com": 1,
     "interpolation": "Linear",
     "metric": [ "Mattes", "Mattes" ],
     "metric_weight": [ 1.0, 1.0 ],
-    "number_of_iterations": [ [ 1000, 500, 250, 100 ], [ 1000, 500, 250, 100 ] ],
+    "number_of_iterations": [ [ 1000, 0 ], [ 250, 100 ] ],
     "radius_or_number_of_bins": [ 32, 32 ],
-    "sampling_percentage": [ 0.25, 0.25 ],
+    "sampling_percentage": [ 1.0, 0.4 ],
     "sampling_strategy": [ "Regular", "Random" ],
-    "shrink_factors": [[ 8, 4, 2, 1] , [ 8, 4, 2, 1 ]],
-    "sigma_units": [ "vox", "vox" ],
-    "smoothing_sigmas": [ [ 4, 2, 1, 0 ], [ 4, 2, 1, 0 ] ],
-    "transform_parameters": [ [ 1.0 ], [ 0.6 ] ],
+    "shrink_factors": [[ 4, 1] , [ 2, 1 ]],
+    "sigma_units": [ "mm", "vox" ],
+    "smoothing_sigmas": [ [ 6, 0 ], [ 2, 1 ] ],
+    "transform_parameters": [ [ 2.0 ], [ 1.0 ] ],
     "transforms": [ "Rigid", "Affine" ],
     "use_estimate_learning_rate_once": [ false, false ],
     "use_histogram_matching": [ true, true ],
+    "winsorize_lower_quantile": 0.001,
+    "winsorize_upper_quantile": 0.998,
     "write_composite_transform": false
 }

sdcflows/data/sd_syn.json

@@ -1,25 +1,24 @@
 {
     "collapse_output_transforms": true,
-    "convergence_threshold": [ 1e-06, 1e-06, 1e-08 ],
-    "convergence_window_size": [ 5, 3, 2 ],
+    "convergence_threshold": [ 1e-06, 1e-08 ],
+    "convergence_window_size": [ 5, 2 ],
     "dimension": 3,
     "interpolation": "Linear",
-    "metric": [ "Mattes", "CC", "CC" ],
-    "metric_weight": [ 1, 1, 1 ],
-    "number_of_iterations": [ [ 100, 0 ], [ 20, 10, 0 ], [ 10 ] ],
+    "metric": [ ["Mattes", "Mattes"], ["Mattes", "Mattes"] ],
+    "metric_weight": [ [0.5, 0.5], [0.5, 0.5] ],
+    "number_of_iterations": [ [ 200, 100 ], [ 10 ] ],
     "output_transform_prefix": "fmap_syn",
-    "radius_or_number_of_bins": [ 56, 8, 5 ],
-    "sampling_percentage": [ 0.80, 0.80, 1.0 ],
-    "sampling_strategy": [ "Random", "Random", "Random" ],
-    "shrink_factors": [ [ 1, 1 ], [ 1, 1, 1 ], [ 1 ] ],
-    "sigma_units": [ "vox", "vox", "vox" ],
-    "smoothing_sigmas": [ [ 2, 0 ], [ 2, 1, 0 ], [ 0 ] ],
-    "transform_parameters": [ [ 0.2, 2.8, 2.8 ], [ 0.2, 2.8, 1.8 ], [ 0.8, 2.0, 1.0 ] ],
-    "transforms": [ "SyN", "SyN", "SyN" ],
-    "use_estimate_learning_rate_once": [ false, true, true ],
-    "use_histogram_matching": [ true, false, false ],
-    "verbose": true,
-    "winsorize_lower_quantile": 0.02,
-    "winsorize_upper_quantile": 1.0,
+    "radius_or_number_of_bins": [ [48, 48], [48, 48] ],
+    "sampling_percentage": [ [0.8, 0.8], [1.0, 1.0] ],
+    "sampling_strategy": [["Random", "Random"], [null, null] ],
+    "shrink_factors": [ [ 1, 1 ], [ 1 ] ],
+    "sigma_units": [ "vox", "vox" ],
+    "smoothing_sigmas": [ [ 2, 0 ], [ 0 ] ],
+    "transform_parameters": [ [ 0.8, 6.0, 3.0 ], [ 0.8, 2.0, 1.0 ] ],
+    "transforms": [ "SyN", "SyN" ],
+    "use_estimate_learning_rate_once": [ true, true ],
+    "use_histogram_matching": [ true, true ],
+    "winsorize_lower_quantile": 0.001,
+    "winsorize_upper_quantile": 0.998,
     "write_composite_transform": false
 }

sdcflows/data/sd_syn_sloppy.json

@@ -1,25 +1,25 @@
 {
     "collapse_output_transforms": true,
-    "convergence_threshold": [ 1e-06, 1e-06, 1e-08 ],
-    "convergence_window_size": [ 5, 3, 2 ],
+    "convergence_threshold": [ 1e-06, 1e-08 ],
+    "convergence_window_size": [ 5, 2 ],
     "dimension": 3,
     "interpolation": "Linear",
-    "metric": [ "Mattes", "CC", "CC" ],
-    "metric_weight": [ 1, 1, 1 ],
-    "number_of_iterations": [ [ 10, 0 ], [ 10, 5, 0 ], [ 2 ] ],
+    "metric": [ ["Mattes", "Mattes"], ["Mattes", "Mattes"] ],
+    "metric_weight": [ [0.5, 0.5], [0.5, 0.5] ],
+    "number_of_iterations": [ [ 20, 10 ], [ 2 ] ],
     "output_transform_prefix": "fmap_syn",
-    "radius_or_number_of_bins": [ 56, 8, 5 ],
-    "sampling_percentage": [ 0.80, 0.80, 1.0 ],
-    "sampling_strategy": [ "Random", "Random", "Random" ],
-    "shrink_factors": [ [ 1, 1 ], [ 1, 1, 1 ], [ 1 ] ],
-    "sigma_units": [ "vox", "vox", "vox" ],
-    "smoothing_sigmas": [ [ 2, 0 ], [ 2, 1, 0 ], [ 0 ] ],
-    "transform_parameters": [ [ 0.2, 2.8, 2.8 ], [ 0.2, 2.8, 1.8 ], [ 0.8, 2.0, 1.0 ] ],
-    "transforms": [ "SyN", "SyN", "SyN" ],
-    "use_estimate_learning_rate_once": [ false, true, true ],
-    "use_histogram_matching": [ true, false, false ],
+    "radius_or_number_of_bins": [ [48, 48], [48, 48] ],
+    "sampling_percentage": [ [0.8, 0.8], [1.0, 1.0] ],
+    "sampling_strategy": [["Random", "Random"], [null, null] ],
+    "shrink_factors": [ [ 1, 1 ], [ 1 ] ],
+    "sigma_units": [ "vox", "vox" ],
+    "smoothing_sigmas": [ [ 2, 0 ], [ 0 ] ],
+    "transform_parameters": [ [ 0.8, 6.0, 3.0 ], [ 0.8, 2.0, 1.0 ] ],
+    "transforms": [ "SyN", "SyN" ],
+    "use_estimate_learning_rate_once": [ true, true ],
+    "use_histogram_matching": [ true, true ],
     "verbose": true,
-    "winsorize_lower_quantile": 0.02,
-    "winsorize_upper_quantile": 1.0,
+    "winsorize_lower_quantile": 0.001,
+    "winsorize_upper_quantile": 0.998,
     "write_composite_transform": false
 }

sdcflows/interfaces/bspline.py

@@ -52,7 +52,7 @@
 
 class _BSplineApproxInputSpec(BaseInterfaceInputSpec):
     in_data = File(exists=True, mandatory=True, desc="path to a fieldmap")
-    in_mask = File(exists=True, mandatory=True, desc="path to a brain mask")
+    in_mask = File(exists=True, desc="path to a brain mask")
     bs_spacing = InputMultiObject(
         [DEFAULT_ZOOMS_MM],
         traits.Tuple(traits.Float, traits.Float, traits.Float),
@@ -127,7 +127,11 @@ def _run_interface(self, runtime):
         # Load in the fieldmap
         fmapnii = nb.load(self.inputs.in_data)
         data = fmapnii.get_fdata(dtype="float32")
-        mask = nb.load(self.inputs.in_mask).get_fdata() > 0
+        mask = (
+            nb.load(self.inputs.in_mask).get_fdata() > 0
+            if isdefined(self.inputs.in_mask)
+            else np.ones_like(data, dtype=bool)
+        )
         bs_spacing = [np.array(sp, dtype="float32") for sp in self.inputs.bs_spacing]
 
         # Recenter the fieldmap

sdcflows/interfaces/fmap.py

@@ -152,3 +152,52 @@ def _run_interface(self, runtime):
             self._results["out_file"]
         )
         return runtime
+
+
+class _DisplacementsField2FieldmapInputSpec(BaseInterfaceInputSpec):
+    transform = File(exists=True, mandatory=True, desc="input displacements field")
+    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"
+    )
+    demean = traits.Bool(False, usedefault=True, desc="regress field to the mean")
+    itk_transform = traits.Bool(
+        True, usedefault=True, desc="whether this is an ITK/ANTs transform"
+    )
+
+
+class _DisplacementsField2FieldmapOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="output fieldmap in Hz")
+
+
+class DisplacementsField2Fieldmap(SimpleInterface):
+    """Convert from a transform to a B0 fieldmap in Hz."""
+
+    input_spec = _DisplacementsField2FieldmapInputSpec
+    output_spec = _DisplacementsField2FieldmapOutputSpec
+
+    def _run_interface(self, runtime):
+        from sdcflows.transform import disp_to_fmap
+
+        self._results["out_file"] = fname_presuffix(
+            self.inputs.transform, suffix="_Hz", newpath=runtime.cwd
+        )
+        fmapnii = disp_to_fmap(
+            nb.load(self.inputs.transform),
+            ro_time=self.inputs.ro_time,
+            pe_dir=self.inputs.pe_dir,
+            itk_format=self.inputs.itk_transform,
+        )
+
+        if self.inputs.demean:
+            data = np.asanyarray(fmapnii.dataobj)
+            data -= np.median(data)
+
+            fmapnii = fmapnii.__class__(
+                data.astype("float32"),
+                fmapnii.affine,
+                fmapnii.header,
+            )
+
+        fmapnii.to_filename(self._results["out_file"])
+        return runtime

sdcflows/tests/test_transform.py

@@ -134,3 +134,26 @@ def test_displacements_field(tmpdir, testdata_dir, outdir, pe_dir, rotation, fli
                 f"_y-{rotation[1] or 0}_z-{rotation[2] or 0}.svg"
             ),
         ).run()
+
+
+@pytest.mark.parametrize("pe_dir", ["j", "j-", "i", "i-", "k", "k-"])
+def test_conversions(tmpdir, testdata_dir, pe_dir):
+    """Check inverse functions."""
+    tmpdir.chdir()
+
+    fmap_nii = nb.load(testdata_dir / "topup-field.nii.gz")
+    new_nii = tf.disp_to_fmap(
+        tf.fmap_to_disp(
+            fmap_nii,
+            ro_time=0.2,
+            pe_dir=pe_dir,
+        ),
+        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

@@ -89,12 +89,14 @@ def fit(self, spatialimage):
 
         # Interpolate the VSM (voxel-shift map)
         vsm = np.zeros(spatialimage.shape[:3], dtype="float32")
-        vsm = (np.squeeze(np.vstack(coeffs).T) @ sparse_vstack(weights)).reshape(
+        vsm = (np.squeeze(np.hstack(coeffs).T) @ sparse_vstack(weights)).reshape(
             vsm.shape
         )
 
         # Cache
         self.shifts = nb.Nifti1Image(vsm, spatialimage.affine, None)
+        self.shifts.header.set_intent("estimate", name="Voxel shift")
+        self.shifts.header.set_xyzt_units("mm")
         return True
 
     def apply(
@@ -215,36 +217,99 @@ def to_displacements(self, ro_time, pe_dir, itk_format=True):
             A NIfTI 1.0 object containing the distortion.
 
         """
-        # Set polarity & scale VSM (voxel-shift-map) by readout time
-        vsm = self.shifts.get_fdata().copy()
-        pe_axis = "ijk".index(pe_dir[0])
-        vsm *= -1.0 if pe_dir.endswith("-") else 1.0
-        vsm *= ro_time
-
-        # Shape of displacements field
-        # Note that ITK NIfTI fields are 5D (have an empty 4th dimension)
-        fieldshape = tuple(list(vsm.shape[:3]) + [1, 3])
-
-        # Convert VSM to voxel displacements
-        ijk_deltas = np.zeros((vsm.size, 3), dtype="float32")
-        ijk_deltas[:, pe_axis] = vsm.reshape(-1)
-
-        # To convert from VSM to RAS field we just apply the affine
-        aff = self.shifts.affine.copy()
-        aff[:3, 3] = 0  # Translations MUST NOT be applied, though.
-        xyz_deltas = nb.affines.apply_affine(aff, ijk_deltas)
-        if itk_format:
-            # ITK displacement vectors are in LPS orientation
-            xyz_deltas[..., (0, 1)] *= -1.0
-
-        xyz_nii = nb.Nifti1Image(
-            xyz_deltas.reshape(fieldshape),
-            self.shifts.affine,
-            None,
-        )
-        xyz_nii.header.set_intent("vector", (), "")
-        xyz_nii.header.set_xyzt_units("mm")
-        return xyz_nii
+        return fmap_to_disp(self.shifts, ro_time, pe_dir, itk_format=itk_format)
+
+
+def fmap_to_disp(fmap_nii, ro_time, pe_dir, itk_format=True):
+    """
+    Convert a fieldmap in Hz into an ITK/ANTs-compatible displacements field.
+
+    The displacements field can be calculated following
+    `Eq. (2) in the fieldmap fitting section
+    <sdcflows.workflows.fit.fieldmap.html#mjx-eqn-eq%3Afieldmap-2>`__.
+
+    Parameters
+    ----------
+    fmap_nii : :obj:`os.pathlike`
+        Path to a voxel-shift-map (VSM) in NIfTI format
+    ro_time : :obj:`float`
+        The total readout time in seconds
+    pe_dir : :obj:`str`
+        The ``PhaseEncodingDirection`` metadata value
+
+    Returns
+    -------
+    spatialimage : :obj:`nibabel.nifti.Nifti1Image`
+        A NIfTI 1.0 object containing the distortion.
+
+    """
+    # Set polarity & scale VSM (voxel-shift-map) by readout time
+    vsm = fmap_nii.get_fdata().copy() * (-ro_time if pe_dir.endswith("-") else ro_time)
+
+    # Shape of displacements field
+    # Note that ITK NIfTI fields are 5D (have an empty 4th dimension)
+    fieldshape = vsm.shape[:3] + (1, 3)
+
+    # Convert VSM to voxel displacements
+    pe_axis = "ijk".index(pe_dir[0])
+    ijk_deltas = np.zeros((vsm.size, 3), dtype="float32")
+    ijk_deltas[:, pe_axis] = vsm.reshape(-1)
+
+    # To convert from VSM to RAS field we just apply the affine
+    aff = fmap_nii.affine.copy()
+    aff[:3, 3] = 0  # Translations MUST NOT be applied, though.
+    xyz_deltas = nb.affines.apply_affine(aff, ijk_deltas)
+    if itk_format:
+        # ITK displacement vectors are in LPS orientation
+        xyz_deltas[..., (0, 1)] *= -1.0
+
+    xyz_nii = nb.Nifti1Image(xyz_deltas.reshape(fieldshape), fmap_nii.affine)
+    xyz_nii.header.set_intent("vector", name="SDC")
+    xyz_nii.header.set_xyzt_units("mm")
+    return xyz_nii
+
+
+def disp_to_fmap(xyz_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.
+    ro_time : :obj:`float`
+        The total readout time in seconds.
+    pe_dir : :obj:`str`
+        The ``PhaseEncodingDirection`` metadata value.
+
+    Returns
+    -------
+    spatialimage : :obj:`nibabel.nifti.Nifti1Image`
+        A NIfTI 1.0 object containing the field in Hz.
+
+    """
+    xyz_deltas = np.squeeze(xyz_nii.get_fdata(dtype="float32")).reshape((-1, 3))
+
+    if itk_format:
+        # 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.
+
+    # 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")
+    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
+
+    fmap_nii = nb.Nifti1Image(vsm / scale_factor, xyz_nii.affine)
+    fmap_nii.header.set_intent("estimate", name="Delta_B0 [Hz]")
+    fmap_nii.header.set_xyzt_units("mm")
+    return fmap_nii
 
 
 def _cubic_bspline(d):