diff --git a/clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_utils.py b/clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_utils.py index 755c509d5..d99195c64 100644 --- a/clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_utils.py +++ b/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 diff --git a/clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_workflows.py b/clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_workflows.py index 2558d5d48..b2e4c699f 100644 --- a/clinica/pipelines/dwi_preprocessing_using_fmap/dwi_preprocessing_using_phasediff_fmap_workflows.py +++ b/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 diff --git a/test/nonregression/pipelines/test_run_pipelines_dwi.py b/test/nonregression/pipelines/test_run_pipelines_dwi.py index 14b9c2123..d7944c53e 100644 --- a/test/nonregression/pipelines/test_run_pipelines_dwi.py +++ b/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"])