Support complex-valued dwidenoising

qsiprep/cli/run.py

@@ -284,9 +284,13 @@ def get_parser():
         action="store",
         nargs="+",
         default=[],
-        choices=["fieldmaps"],
-        help="ignore selected aspects of the input dataset to disable "
-        "corresponding parts of the workflow (a space delimited list)",
+        choices=["fieldmaps", "phase"],
+        help=(
+            "ignore selected aspects of the input dataset to disable "
+            "corresponding parts of the workflow (a space delimited list). "
+            'Ignoring "phase" will disable complex-valued denoising, '
+            "when phase DWI data are available."
+        ),
     )
     g_conf.add_argument(
         "--longitudinal",

qsiprep/interfaces/dwi_merge.py

@@ -3,6 +3,7 @@
 import json
 import os.path as op
 
+import nibabel as nb
 import numpy as np
 import pandas as pd
 from nilearn.image import concat_imgs, index_img, iter_img, load_img, math_img
@@ -714,3 +715,77 @@ def create_provenance_dataframe(
     image_df = pd.concat(series_confounds, axis=0, ignore_index=True)
     image_df["original_file"] = bids_sources
     return image_df
+
+
+class _PhaseToRadInputSpec(BaseInterfaceInputSpec):
+    phase_file = File(exists=True, mandatory=True)
+
+
+class _PhaseToRadOutputSpec(TraitedSpec):
+    phase_file = File(exists=True)
+
+
+class PhaseToRad(SimpleInterface):
+    """Convert phase image from arbitrary units (au) to radians.
+
+    This method assumes that the phase image's minimum and maximum values correspond to
+    -pi and pi, respectively, and scales the image to be between 0 and 2*pi.
+
+    Notes
+    -----
+    The code is derived from
+    https://github.com/nipreps/sdcflows/blob/c6cd42944f4b6d638716ce020ffe51010e9eb58a/\
+    sdcflows/utils/phasemanip.py#L26.
+
+    License
+    -------
+    Copyright 2021 The NiPreps Developers <[email protected]>
+
+    Licensed under the Apache License, Version 2.0 (the "License");
+    you may not use this file except in compliance with the License.
+    You may obtain a copy of the License at
+
+        http://www.apache.org/licenses/LICENSE-2.0
+
+    Unless required by applicable law or agreed to in writing, software
+    distributed under the License is distributed on an "AS IS" BASIS,
+    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+    See the License for the specific language governing permissions and
+    limitations under the License.
+
+    We support and encourage derived works from this project, please read
+    about our expectations at
+
+        https://www.nipreps.org/community/licensing/
+
+    """
+
+    input_spec = _PhaseToRadInputSpec
+    output_spec = _PhaseToRadOutputSpec
+
+    def _run_interface(self, runtime):
+        im = nb.load(self.inputs.phase_file)
+        data = im.get_fdata(caching="unchanged")  # Read as float64 for safety
+        hdr = im.header.copy()
+
+        # Rescale to [0, 2*pi]
+        data = (data - data.min()) * (2 * np.pi / (data.max() - data.min()))
+
+        # Round to float32 and clip
+        data = np.clip(np.float32(data), 0.0, 2 * np.pi)
+
+        hdr.set_data_dtype(np.float32)
+        hdr.set_xyzt_units("mm")
+
+        # Set the output file name
+        self._results["phase_file"] = fname_presuffix(
+            self.inputs.phase_file,
+            suffix="_rad.nii.gz",
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+
+        # Save the output image
+        nb.Nifti1Image(data, None, hdr).to_filename(self._results["phase_file"])
+
+        return runtime

qsiprep/interfaces/mrtrix.py

@@ -1365,3 +1365,54 @@ class TransformHeader(CommandLine):
     input_spec = _TransformHeaderInputSpec
     output_spec = _TransformHeaderOutputSpec
     _cmd = "mrtransform -strides -1,-2,3"
+
+
+class _PolarToComplexInputSpec(CommandLineInputSpec):
+    mag_file = traits.File(exists=True, mandatory=True, position=0, argstr="%s")
+    phase_file = traits.File(exists=True, mandatory=True, position=1, argstr="%s")
+    out_file = traits.File(
+        exists=False,
+        name_source="mag_file",
+        name_template="%s_complex.nii.gz",
+        keep_extension=False,
+        position=-1,
+        argstr="-polar %s",
+    )
+
+
+class _PolarToComplexOutputSpec(TraitedSpec):
+    out_file = File(exists=True)
+
+
+class PolarToComplex(CommandLine):
+    """Convert a magnitude and phase image pair to a single complex image using mrcalc."""
+
+    input_spec = _PolarToComplexInputSpec
+    output_spec = _PolarToComplexOutputSpec
+
+    _cmd = "mrcalc"
+
+
+class _ComplexToMagnitudeInputSpec(CommandLineInputSpec):
+    complex_file = traits.File(exists=True, mandatory=True, position=0, argstr="%s")
+    out_file = traits.File(
+        exists=False,
+        name_source="complex_file",
+        name_template="%s_mag.nii.gz",
+        keep_extension=False,
+        position=-1,
+        argstr="-abs %s",
+    )
+
+
+class _ComplexToMagnitudeOutputSpec(TraitedSpec):
+    out_file = File(exists=True)
+
+
+class ComplexToMagnitude(CommandLine):
+    """Extract the magnitude portion of a complex image using mrcalc."""
+
+    input_spec = _ComplexToMagnitudeInputSpec
+    output_spec = _ComplexToMagnitudeOutputSpec
+
+    _cmd = "mrcalc"

qsiprep/utils/bids.py

@@ -157,10 +157,7 @@ def collect_participants(bids_dir, participant_label=None, strict=False, bids_va
 
 
 def collect_data(bids_dir, participant_label, filters=None, bids_validate=True):
-    """
-    Uses pybids to retrieve the input data for a given participant
-
-    """
+    """Use pybids to retrieve the input data for a given participant."""
     if isinstance(bids_dir, BIDSLayout):
         layout = bids_dir
     else:
@@ -173,7 +170,7 @@ def collect_data(bids_dir, participant_label, filters=None, bids_validate=True):
         "t2w": {"datatype": "anat", "suffix": "T2w"},
         "t1w": {"datatype": "anat", "suffix": "T1w"},
         "roi": {"datatype": "anat", "suffix": "roi"},
-        "dwi": {"datatype": "dwi", "suffix": "dwi"},
+        "dwi": {"datatype": "dwi", "part": ["mag", None], "suffix": "dwi"},
     }
     bids_filters = filters or {}
     for acq, entities in bids_filters.items():

qsiprep/utils/grouping.py

@@ -51,7 +51,6 @@ def group_dwi_scans(
         A dict where the keys are the BIDS derivatives name of the output file after
         concatenation. The values are lists of dwi files in that group.
     """
-
     # Handle the grouping of multiple dwi files within a session
     dwi_session_groups = get_session_groups(bids_layout, subject_data, combine_scans)
 

qsiprep/workflows/base.py

@@ -158,8 +158,7 @@ def init_qsiprep_wf(
     anatomical_contrst : str
         Which contrast to use for the anatomical reference
     ignore : list
-        Preprocessing steps to skip (may include "slicetiming",
-        "fieldmaps")
+        Preprocessing steps to skip (may include "slicetiming", "fieldmaps", "phase").
     low_mem : bool
         Write uncompressed .nii files in some cases to reduce memory usage
     anat_only : bool
@@ -401,7 +400,7 @@ def init_single_subject_wf(
     name : str
         Name of workflow
     ignore : list
-        Preprocessing steps to skip (may include "sbref", "fieldmaps")
+        Preprocessing steps to skip (may include "sbref", "fieldmaps", "phase")
     debug : bool
         Do inaccurate but fast normalization
     low_mem : bool
@@ -496,7 +495,10 @@ def init_single_subject_wf(
         LOGGER.warning("Building a test workflow")
     else:
         subject_data, layout = collect_data(
-            bids_dir, subject_id, filters=bids_filters, bids_validate=False
+            bids_dir,
+            subject_id,
+            filters=bids_filters,
+            bids_validate=False,
         )
 
     # Warn about --dwi-only and non-none --anat-modality

qsiprep/workflows/dwi/base.py

@@ -373,6 +373,8 @@ def init_dwi_preproc_wf(
         source_file=source_file,
         low_mem=low_mem,
         denoise_before_combining=denoise_before_combining,
+        layout=layout,
+        ignore=ignore,
         omp_nthreads=omp_nthreads,
     )
     test_pre_hmc_connect = pe.Node(TestInput(), name="test_pre_hmc_connect")