ENH: Allow --ignore fmap-jacobian to disable Jacobian determinant modulation during fieldmap correction

fmriprep/cli/parser.py

@@ -285,7 +285,7 @@ def _slice_time_ref(value, parser):
         action="store",
         nargs="+",
         default=[],
-        choices=["fieldmaps", "slicetiming", "sbref", "t2w", "flair"],
+        choices=["fieldmaps", "slicetiming", "sbref", "t2w", "flair", "fmap-jacobian"],
         help="Ignore selected aspects of the input dataset to disable corresponding "
         "parts of the workflow (a space delimited list)",
     )

fmriprep/interfaces/resampling.py

@@ -49,6 +49,7 @@
         "k-",
         desc="the phase-encoding direction corresponding to in_data",
     )
+    jacobian = traits.Bool(mandatory=True, desc="Whether to apply Jacobian correction")
     num_threads = traits.Int(1, usedefault=True, desc="Number of threads to use for resampling")
     output_data_type = traits.Str("float32", usedefault=True, desc="Data type of output image")
     order = traits.Int(3, usedefault=True, desc="Order of interpolation (0=nearest, 3=cubic)")
@@ -105,6 +106,7 @@
             transforms=transforms,
             fieldmap=fieldmap,
             pe_info=pe_info,
+            jacobian=self.inputs.jacobian,
             nthreads=self.inputs.num_threads,
             output_dtype=self.inputs.output_data_type,
             order=self.inputs.order,
@@ -217,6 +219,7 @@
     data: np.ndarray,
     coordinates: np.ndarray,
     pe_info: tuple[int, float],
+    jacobian: bool,
     hmc_xfm: np.ndarray | None,
     fmap_hz: np.ndarray,
     output: np.dtype | np.ndarray | None = None,
@@ -282,8 +285,6 @@
     vsm = fmap_hz * pe_info[1]
     coordinates[pe_info[0], ...] += vsm
 
-    jacobian = 1 + np.gradient(vsm, axis=pe_info[0])
-
     result = ndi.map_coordinates(
         data,
         coordinates,
@@ -293,14 +294,18 @@
         cval=cval,
         prefilter=prefilter,
     )
-    result *= jacobian
+
+    if jacobian:
+        result *= 1 + np.gradient(vsm, axis=pe_info[0])
+
     return result
 
 
 async def resample_series_async(
     data: np.ndarray,
     coordinates: np.ndarray,
     pe_info: list[tuple[int, float]],
+    jacobian: bool,
     hmc_xfms: list[np.ndarray] | None,
     fmap_hz: np.ndarray,
     output_dtype: np.dtype | None = None,
@@ -361,6 +366,7 @@
             data,
             coordinates,
             pe_info[0],
+            jacobian,
             hmc_xfms[0] if hmc_xfms else None,
             fmap_hz,
             output_dtype,
@@ -384,6 +390,7 @@
                     data=volume,
                     coordinates=coordinates,
                     pe_info=pe_info[volid],
+                    jacobian=jacobian,
                     hmc_xfm=hmc_xfms[volid] if hmc_xfms else None,
                     fmap_hz=fmap_hz,
                     output=out_array[..., volid],
@@ -407,6 +414,7 @@
     data: np.ndarray,
     coordinates: np.ndarray,
     pe_info: list[tuple[int, float]],
+    jacobian: bool,
     hmc_xfms: list[np.ndarray] | None,
     fmap_hz: np.ndarray,
     output_dtype: np.dtype | None = None,
@@ -467,6 +475,7 @@
             data=data,
             coordinates=coordinates,
             pe_info=pe_info,
+            jacobian=jacobian,
             hmc_xfms=hmc_xfms,
             fmap_hz=fmap_hz,
             output_dtype=output_dtype,
@@ -485,6 +494,7 @@
     transforms: nt.TransformChain,
     fieldmap: nb.Nifti1Image | None,
     pe_info: list[tuple[int, float]] | None,
+    jacobian: bool = True,
     nthreads: int = 1,
     output_dtype: np.dtype | str | None = 'f4',
     order: int = 3,
@@ -566,6 +576,7 @@
         data=source.get_fdata(dtype='f4'),
         coordinates=mapped_coordinates.T.reshape((3, *target.shape[:3])),
         pe_info=pe_info,
+        jacobian=jacobian,
         hmc_xfms=hmc_xfms,
         fmap_hz=fieldmap.get_fdata(dtype='f4'),
         output_dtype=output_dtype,

fmriprep/workflows/bold/apply.py

@@ -16,6 +16,7 @@ def init_bold_volumetric_resample_wf(
     *,
     metadata: dict,
     mem_gb: dict[str, float],
+    jacobian: bool,
     fieldmap_id: str | None = None,
     omp_nthreads: int = 1,
     name: str = 'bold_volumetric_resample_wf',
@@ -123,7 +124,7 @@ def init_bold_volumetric_resample_wf(
     boldref2target = pe.Node(niu.Merge(2), name='boldref2target', run_without_submitting=True)
     bold2target = pe.Node(niu.Merge(2), name='bold2target', run_without_submitting=True)
     resample = pe.Node(
-        ResampleSeries(),
+        ResampleSeries(jacobian=jacobian),
         name="resample",
         n_procs=omp_nthreads,
         mem_gb=mem_gb['resampled'],

fmriprep/workflows/bold/base.py

@@ -380,6 +380,7 @@ def init_bold_wf(
         fieldmap_id=fieldmap_id if not multiecho else None,
         omp_nthreads=omp_nthreads,
         mem_gb=mem_gb,
+        jacobian='fmap-jacobian' not in config.workflow.ignore,
         name='bold_anat_wf',
     )
     bold_anat_wf.inputs.inputnode.resolution = "native"
@@ -437,6 +438,7 @@ def init_bold_wf(
             fieldmap_id=fieldmap_id if not multiecho else None,
             omp_nthreads=omp_nthreads,
             mem_gb=mem_gb,
+            jacobian='fmap-jacobian' not in config.workflow.ignore,
             name='bold_std_wf',
         )
         ds_bold_std_wf = init_ds_volumes_wf(
@@ -521,6 +523,7 @@ def init_bold_wf(
             fieldmap_id=fieldmap_id if not multiecho else None,
             omp_nthreads=omp_nthreads,
             mem_gb=mem_gb,
+            jacobian='fmap-jacobian' not in config.workflow.ignore,
             name='bold_MNI6_wf',
         )
 

fmriprep/workflows/bold/fit.py

@@ -838,7 +838,7 @@ def init_bold_native_wf(
 
     # Resample to boldref
     boldref_bold = pe.Node(
-        ResampleSeries(),
+        ResampleSeries(jacobian="fmap-jacobian" not in config.workflow.ignore),
         name="boldref_bold",
         n_procs=omp_nthreads,
         mem_gb=mem_gb["resampled"],