Support complex-valued dwidenoising

qsiprep/cli/run.py

@@ -254,9 +254,14 @@ 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',
         action='store_true',

qsiprep/interfaces/bids.py

@@ -115,6 +115,7 @@ class BIDSDataGrabberOutputSpec(TraitedSpec):
     t2w = OutputMultiPath(desc='output T2w images')
     flair = OutputMultiPath(desc='output FLAIR images')
     dwi = OutputMultiPath(desc='output DWI images')
+    dwi_phase = OutputMultiPath(desc='output DWI images')
 
 
 class BIDSDataGrabber(SimpleInterface):

qsiprep/interfaces/dwi_merge.py

@@ -1,6 +1,8 @@
 """Handle merging and spliting of DSI files."""
 import os.path as op
 import json
+
+import nibabel as nb
 import numpy as np
 import pandas as pd
 from nilearn.image import concat_imgs, load_img, index_img, math_img, iter_img
@@ -9,6 +11,7 @@
 from nipype.utils.filemanip import fname_presuffix
 from nipype.interfaces import ants
 from nipype import logging
+
 from .fmap import get_distortion_grouping
 from ..workflows.dwi.util import _get_concatenated_bids_name
 LOGGER = logging.getLogger('nipype.workflow')
@@ -659,3 +662,36 @@ def create_provenance_dataframe(bids_sources, harmonized_niis, b0_means,
     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):
+    input_spec = _PhaseToRadInputSpec
+    output_spec = _PhaseToRadOutputSpec
+
+    def _run_interface(self, runtime):
+        phase_img = nb.load(self.inputs.phase_file)
+
+        phase_data = phase_img.get_fdata()
+        imax = phase_data.max()
+        imin = phase_data.min()
+        scaled = (phase_data - imin) / (imax - imin)
+        rad_data = 2 * np.pi * scaled
+        out_img = nb.Nifti1Image(rad_data, phase_img.affine, phase_img.header)
+
+        self._results['phase_file'] = fname_presuffix(
+            self.inputs.phase_file,
+            suffix='_rad.nii.gz',
+            newpath=runtime.cwd,
+            use_ext=False,
+        )
+        out_img.to_filename(self._results['phase_file'])
+
+        return runtime

qsiprep/interfaces/mrtrix.py

@@ -1356,3 +1356,61 @@ 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,
+        mandatory=True,
+        position=-1,
+        argstr="-polar %s"
+    )
+
+
+class _PolarToComplexOutputSpec(TraitedSpec):
+    out_file = File(exists=True)
+
+
+class PolarToComplex(CommandLine):
+    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,
+        mandatory=True,
+        position=-1,
+        argstr="-abs %s"
+    )
+
+
+class _ComplexToMagnitudeOutputSpec(TraitedSpec):
+    out_file = File(exists=True)
+
+
+class ComplexToMagnitude(CommandLine):
+    input_spec = _ComplexToMagnitudeInputSpec
+    output_spec = _ComplexToMagnitudeOutputSpec
+
+    _cmd = "mrcalc"

qsiprep/utils/bids.py

@@ -135,10 +135,7 @@ def collect_participants(bids_dir, participant_label=None, strict=False,
 
 
 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:
@@ -151,7 +148,8 @@ 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'},
+        'dwi_phase': {'datatype': 'dwi', 'part': 'phase', 'suffix': 'dwi'},
     }
     bids_filters = filters or {}
     for acq, entities in bids_filters.items():

qsiprep/utils/grouping.py

@@ -41,7 +41,6 @@ def group_dwi_scans(bids_layout, subject_data, using_fsl=False, combine_scans=Tr
             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

@@ -132,8 +132,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
@@ -369,7 +368,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
@@ -479,12 +478,17 @@ 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,
         )
 
+    if 'phase' in ignore:
+        subject_data['dwi_phase'] = []
+
+    phase_available = bool(subject_data['dwi_phase'])
+
     # Warn about --dwi-only and non-none --anat-modality
     if dwi_only and anatomical_contrast != "none":
         LOGGER.warn(
@@ -741,6 +745,7 @@ def init_single_subject_wf(
             fmap_bspline=fmap_bspline,
             fmap_demean=fmap_demean,
             t2w_sdc=bool(subject_data.get('t2w')),
+            phase_available=phase_available,
             sloppy=debug,
             source_file=source_file
         )

qsiprep/workflows/dwi/base.py

@@ -64,6 +64,7 @@ def init_dwi_preproc_wf(dwi_only,
                         low_mem,
                         sloppy,
                         source_file,
+                        phase_available,
                         layout=None):
     """
     This workflow controls the dwi preprocessing stages of qsiprep.
@@ -104,6 +105,7 @@ def init_dwi_preproc_wf(dwi_only,
                                   low_mem=False,
                                   sloppy=True,
                                   source_file='/data/bids/sub-1/dwi/sub-1_dwi.nii.gz',
+                                  phase_available=False,
                                   layout=None)
 
     **Parameters**
@@ -172,6 +174,10 @@ def init_dwi_preproc_wf(dwi_only,
             (default is 1)
         sloppy : bool
             Use low-quality settings for motion correction
+        phase_available : bool
+            True if phase data are available for the DWI scan.
+            If True, and ``denoise_method`` is ``dwidenoise``, then ``dwidenoise``
+            will be run on the complex-valued data.
         source_file : str
             The file name template used for derivatives
 
@@ -324,6 +330,7 @@ def init_dwi_preproc_wf(dwi_only,
                                      source_file=source_file,
                                      low_mem=low_mem,
                                      denoise_before_combining=denoise_before_combining,
+                                     phase_available=phase_available,
                                      omp_nthreads=omp_nthreads)
     test_pre_hmc_connect = pe.Node(TestInput(), name='test_pre_hmc_connect')
     if hmc_model in ('none', '3dSHORE','tensor'):