[TEST] Test workflow prepare_phasediff_fmap (#900)

* init work

* add test

* some fixes to the test

* Update clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_workflows.py

Co-authored-by: Matthieu Joulot <85217698+MatthieuJoulot@users.noreply.github.com>

* Update clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_workflows.py

---------

Co-authored-by: Matthieu Joulot <85217698+MatthieuJoulot@users.noreply.github.com>

clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_utils.py

@@ -147,8 +147,28 @@ def print_end_pipeline(image_id, final_file):
     print_end_image(image_id)
 
 
-def rads2hz(in_file, delta_te, out_file=None):
-    """Convert input phase difference map to Hz."""
+def rads2hz(in_file: str, delta_te: float, out_file: str = None) -> str:
+    """Convert input phase difference map to Hz.
+
+    Parameters
+    ----------
+    in_file : str
+        Path to the phase difference map image.
+
+    delta_te : float
+        The DeltaEchoTime from BIDS specifications.
+
+    out_file : str, optional
+        Path to output file. If None, the output image
+        will be written in current folder. The file name
+        will have the same base name as in_file, but with
+        the "_radsec.nii.gz" suffix.
+
+    Returns
+    -------
+    out_file : str
+        Path to output file.
+    """
     import math
     import os
 
@@ -162,8 +182,9 @@ def rads2hz(in_file, delta_te, out_file=None):
         out_file = os.path.abspath(f"./{fname}_radsec.nii.gz")
 
     im = nb.load(in_file)
-    data = im.get_data().astype(np.float32) * (1.0 / (delta_te * 2 * math.pi))
+    data = im.get_data().astype(np.float32) * (1.0 / (float(delta_te) * 2 * math.pi))
     nb.Nifti1Image(data, im.affine, im.header).to_filename(out_file)
+
     return out_file
 
 

clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_workflows.py

@@ -3,25 +3,55 @@
 from nipype.pipeline.engine import Workflow
 
 
-def prepare_phasediff_fmap(name="prepare_phasediff_fmap"):
+def prepare_phasediff_fmap(
+    output_dir: Optional[str] = None,
+    name: Optional[str] = "prepare_phasediff_fmap",
+) -> Workflow:
     """This workflow adapts the fsl_prepare_fieldmap script from FSL for the FSL eddy command.
 
-    Please note that the step 3 converts the fieldmap into Hz instead of rad/s
-    in the initial script because FSL eddy --field expects the fieldmap to be in Hz.
+    Please note that the step 3 converts the field map into Hz instead of rad/s
+    in the initial script because FSL eddy --field expects the field map to be in Hz.
+
+    This workflow needs FSL.
 
-    Input node:
-        fmap_mask (str): Binary mask of the fieldmap.
-        fmap_phasediff (str): Phase difference fieldmap.
-        fmap_magnitude (str): Brain extracted magnitude fieldmap. Chose the fieldmap with the best contrast.
-        delta_echo_time (float): DeltaEchoTime from BIDS specifications.
+    The workflow performs the following steps:
+        - Step 1 - Convert the field map into radians
+        - Step 2 - Unwrap the field map with PRELUDE
+        - Step 3 - Convert the field map to Hz
+        - Step 4 - Call FUGUE to extrapolate from mask (fill holes, etc)
+        - Step 5 - Demean the field map to avoid gross shifting
+        - Step 6 - Clean up the edge voxels
 
-    Output node:
-        calibrated_fmap (str): Calibrated fieldmap for eddy --field command.
+    Parameters
+    ----------
+    output_dir: str, optional
+         Path to output directory.
+         If provided, the pipeline will write its output in this folder.
+         Default to None.
+
+    name : str, optional
+        Name of the workflow. Default="prepare_phasediff_fmap".
+
+    Returns
+    -------
+    Workflow :
+        The Nipype workflow.
+         This workflow has the following inputs:
+            - "fmap_mask" : str, the path to the binary mask of the field map.
+            - "fmap_phasediff" : str, the path to the phase difference field map.
+            - "fmap_magnitude" : str, the path to the brain extracted magnitude field map.
+              Chose the field map with the best contrast.
+            - "delta_echo_time" : float, the DeltaEchoTime from BIDS specifications.
+        And the following outputs:
+            - "calibrated_fmap" : str, the path to the calibrated field map for eddy --field command.
 
-    Warnings:
-        This workflow can not be used for PRELUDE. You would need to change the step 3 (conversion into rad/s instead of Hz).
+    Warnings
+    --------
+    This workflow can not be used for PRELUDE.
+    You would need to change the step 3 (conversion into rad/s instead of Hz).
     """
     import nipype.interfaces.fsl as fsl
+    import nipype.interfaces.io as nio
     import nipype.interfaces.utility as nutil
     import nipype.pipeline.engine as npe
     from niflow.nipype1.workflows.dmri.fsl.utils import (
@@ -34,79 +64,96 @@ def prepare_phasediff_fmap(name="prepare_phasediff_fmap"):
 
     input_node = npe.Node(
         nutil.IdentityInterface(
-            fields=["fmap_mask", "fmap_phasediff", "fmap_magnitude", "delta_echo_time"]
+            fields=[
+                "fmap_mask",
+                "fmap_phasediff",
+                "fmap_magnitude",
+                "delta_echo_time",
+            ]
         ),
         name="input_node",
     )
 
-    output_node = npe.Node(
-        nutil.IdentityInterface(fields=["calibrated_fmap"]), name="output_node"
-    )
-
-    # Step 1 - Convert the fmap into radians
-    pha2rads = npe.Node(
+    convert_fmap_to_rads = npe.Node(
         nutil.Function(
             input_names=["in_file"], output_names=["out_file"], function=siemens2rads
         ),
-        name="1-ConvertFMapToRads",
+        name="convert_fmap_to_rads",
     )
 
-    # Step 2 - Unwrap the fmap with PRELUDE
-    prelude = npe.Node(fsl.PRELUDE(process3d=True), name="2-PhaseUnwrap")
+    unwrap_fmap_with_prelude = npe.Node(
+        fsl.PRELUDE(process3d=True), name="unwrap_fmap_with_prelude"
+    )
 
-    # Step 3 - Convert the fmap to Hz
-    rads2hz = npe.Node(
+    convert_fmap_to_hz = npe.Node(
         nutil.Function(
             input_names=["in_file", "delta_te"],
             output_names=["out_file"],
             function=rads2hz,
         ),
-        name="3-ConvertFMapToHz",
+        name="convert_fmap_to_hz",
     )
 
-    # Step 4 - Call FUGUE to extrapolate from mask (fill holes, etc)
-    pre_fugue = npe.Node(fsl.FUGUE(save_fmap=True), name="4-FugueExtrapolationFromMask")
+    fugue_extrapolation_from_mask = npe.Node(
+        fsl.FUGUE(save_fmap=True), name="fugue_extrapolation_from_mask"
+    )
 
-    # Step 5 - Demean fmap to avoid gross shifting
-    demean = npe.Node(
+    demean_fmap = npe.Node(
         nutil.Function(
             input_names=["in_file", "in_mask"],
             output_names=["out_file"],
             function=demean_image,
         ),
-        name="5-DemeanFMap",
+        name="demean_fmap",
     )
 
-    # Step 6 - Clean up edge voxels
-    cleanup = cleanup_edge_pipeline(name="6-CleanUpEdgeVoxels")
+    cleanup_edge_voxels = cleanup_edge_pipeline(name="cleanup_edge_voxels")
+
+    output_node = npe.Node(
+        nutil.IdentityInterface(fields=["calibrated_fmap"]),
+        name="output_node",
+    )
+
+    if output_dir:
+        write_results = npe.Node(name="write_results", interface=nio.DataSink())
+        write_results.inputs.base_directory = output_dir
+        write_results.inputs.parameterization = False
 
     wf = npe.Workflow(name=name)
-    # fmt: off
-    wf.connect(
-        [
-            # Step 1 - Convert the fmap into radians
-            (input_node, pha2rads, [("fmap_phasediff", "in_file")]),
-            # Step 2 - Unwrap the fmap with PRELUDE
-            (pha2rads, prelude, [("out_file", "phase_file")]),
-            (input_node, prelude, [("fmap_magnitude", "magnitude_file")]),
-            (input_node, prelude, [("fmap_mask", "mask_file")]),
-            # Step 3 - Convert the fmap to Hz
-            (prelude, rads2hz, [("unwrapped_phase_file", "in_file")]),
-            (input_node, rads2hz, [("delta_echo_time", "delta_te")]),
-            # Step 4 - Call FUGUE to extrapolate from mask (fill holes, etc)
-            (rads2hz, pre_fugue, [("out_file", "fmap_in_file")]),
-            (input_node, pre_fugue, [("fmap_mask", "mask_file")]),
-            # Step 5 - Demean fmap to avoid gross shifting
-            (pre_fugue, demean, [("fmap_out_file", "in_file")]),
-            (input_node, demean, [("fmap_mask", "in_mask")]),
-            # Step 6 - Clean up edge voxels
-            (demean, cleanup, [("out_file", "inputnode.in_file")]),
-            (input_node, cleanup, [("fmap_mask", "inputnode.in_mask")]),
-            # Output node
-            (cleanup, output_node, [("outputnode.out_file", "calibrated_fmap")]),
+
+    connections = [
+        (input_node, convert_fmap_to_rads, [("fmap_phasediff", "in_file")]),
+        (convert_fmap_to_rads, unwrap_fmap_with_prelude, [("out_file", "phase_file")]),
+        (input_node, unwrap_fmap_with_prelude, [("fmap_magnitude", "magnitude_file")]),
+        (input_node, unwrap_fmap_with_prelude, [("fmap_mask", "mask_file")]),
+        (
+            unwrap_fmap_with_prelude,
+            convert_fmap_to_hz,
+            [("unwrapped_phase_file", "in_file")],
+        ),
+        (input_node, convert_fmap_to_hz, [("delta_echo_time", "delta_te")]),
+        (
+            convert_fmap_to_hz,
+            fugue_extrapolation_from_mask,
+            [("out_file", "fmap_in_file")],
+        ),
+        (input_node, fugue_extrapolation_from_mask, [("fmap_mask", "mask_file")]),
+        (fugue_extrapolation_from_mask, demean_fmap, [("fmap_out_file", "in_file")]),
+        (input_node, demean_fmap, [("fmap_mask", "in_mask")]),
+        (demean_fmap, cleanup_edge_voxels, [("out_file", "inputnode.in_file")]),
+        (input_node, cleanup_edge_voxels, [("fmap_mask", "inputnode.in_mask")]),
+        (
+            cleanup_edge_voxels,
+            output_node,
+            [("outputnode.out_file", "calibrated_fmap")],
+        ),
+    ]
+
+    if output_dir:
+        connections += [
+            (output_node, write_results, [("calibrated_fmap", "calibrated_fmap")]),
         ]
-    )
-    # fmt: on
+    wf.connect(connections)
 
     return wf
 

test/nonregression/pipelines/test_run_pipelines_dwi.py

@@ -129,6 +129,52 @@ def test_dwi_compute_reference_b0(cmdopt, tmp_path):
         assert similarity_measure(out_file, ref_file, 0.99)
 
 
+@pytest.mark.fast
+def test_prepare_phasediff_fmap(cmdopt, tmp_path):
+    """Test the pipeline prepare_phasediff_fmap which is part of
+     step 2 of the pipeline DWIPreprocessingUsingPhaseDiff.
+
+    This is a fast test which should run in less than a minute.
+    """
+    from clinica.pipelines.dwi_preprocessing_using_fmap.dwi_preprocessing_using_phasediff_fmap_workflows import (
+        prepare_phasediff_fmap,
+    )
+    from clinica.utils.filemanip import extract_metadata_from_json
+
+    base_dir = Path(cmdopt["input"])
+    input_dir, tmp_dir, ref_dir = configure_paths(
+        base_dir, tmp_path, "DWIPreparePhasediffFmap"
+    )
+    (tmp_path / "tmp").mkdir()
+
+    [echo_time1, echo_time2] = extract_metadata_from_json(
+        str(input_dir / "sub-01_ses-M000_phasediff.json"), ["EchoTime1", "EchoTime2"]
+    )
+    delta_echo_time = abs(echo_time2 - echo_time1)
+
+    wf = prepare_phasediff_fmap(output_dir=str(tmp_path / "tmp"))
+    wf.inputs.input_node.fmap_mask = str(
+        input_dir / "sub-01_ses-M000_magnitude1_corrected_brain_mask.nii.gz"
+    )
+    wf.inputs.input_node.fmap_phasediff = str(
+        input_dir / "sub-01_ses-M000_phasediff.nii.gz"
+    )
+    wf.inputs.input_node.fmap_magnitude = str(
+        input_dir / "sub-01_ses-M000_magnitude1_corrected_brain.nii.gz"
+    )
+    wf.inputs.input_node.delta_echo_time = str(delta_echo_time)
+
+    wf.run()
+
+    out_filename = (
+        "sub-01_ses-M000_phasediff_rads_unwrapped_radsec_fieldmap_demean_maths.nii.gz"
+    )
+    out_file = fspath(tmp_path / "tmp" / "calibrated_fmap" / out_filename)
+    ref_file = fspath(ref_dir / "calibrated_fmap" / out_filename)
+
+    assert similarity_measure(out_file, ref_file, 0.99)
+
+
 @pytest.mark.slow
 def test_dwi_preprocessing_using_phase_diff_field_map(cmdopt, tmp_path):
     base_dir = Path(cmdopt["input"])