Add phase correction options

qsiprep/cli/parser.py

@@ -395,6 +395,19 @@ def _bids_filter(value, parser):
         help='Image-based denoising method. Either "dwidenoise" (MRtrix), '
         '"patch2self" (DIPY) or "none". (default: dwidenoise)',
     )
+    g_conf.add_argument(
+        '--dwi-phase-correction',
+        action='store',
+        choices=['none', 'tv', 'tvc', 'dc'],
+        default='none',
+        help='Phase-correction method for complex-valued DWI data. '
+        'When "tv" (Eichner 2015, TV on the magnitude/phase stack), "tvc" '
+        '(Eichner 2015, paper-faithful TV on the complex signal) or "dc" '
+        '(Sprenger 2017), the real channel after rephasing is used for '
+        'downstream processing instead of the magnitude. Requires phase data '
+        '(a BIDS part-phase file) and --denoise-method dwidenoise; otherwise it '
+        'is ignored. (default: none)',
+    )
     g_conf.add_argument(
         '--unringing-method',
         action='store',
@@ -755,6 +768,18 @@ def parse_args(args=None, namespace=None):
                 'The --dwi-denoise-window option is not used when --denoise-method=none'
             )
 
+    if config.workflow.dwi_phase_correction != 'none':
+        if config.workflow.denoise_method != 'dwidenoise':
+            config.loggers.cli.warning(
+                'The --dwi-phase-correction option requires --denoise-method '
+                'dwidenoise; it will be ignored.'
+            )
+        if 'phase' in config.workflow.ignore:
+            config.loggers.cli.warning(
+                'The --dwi-phase-correction option has no effect when phase '
+                'data are ignored (--ignore phase); it will be ignored.'
+            )
+
     bids_dir = config.execution.bids_dir
     output_dir = config.execution.output_dir
     work_dir = config.execution.work_dir

qsiprep/config.py

@@ -569,6 +569,18 @@ class workflow(_Config):
     denoise_method = None
     """Image-based denoising method. Either "dwidenoise" (MRtrix), "patch2self" (DIPY)
     or "none"."""
+    dwi_phase_correction = None
+    """Phase-correction method applied to complex denoised DWI data.
+    One of "none", "tv" (Eichner 2015) or "dc" (Sprenger 2017). When "tv" or
+    "dc", the real channel after rephasing is passed downstream in place of the
+    magnitude. Only takes effect when phase data are present and
+    ``denoise_method`` is "dwidenoise"."""
+    dwi_phase_tv_weight = 6.0
+    """Total-variation regularization weight for the "tv" phase-correction
+    method (Eichner 2015)."""
+    dwi_phase_dc_kernel = 'Opt5'
+    """Convolution kernel name for the "dc" phase-correction method
+    (Sprenger 2017). One of B3, B5, G3F1, G5F2, G3F1H, G5F2H, Opt3, Opt5."""
     distortion_group_merge = None
     """How to combine images across distortion groups (concatenate, average or none)."""
     dwi_denoise_window = None

qsiprep/data/boilerplate.bib

@@ -525,3 +525,26 @@ @article{tortoisev4
 	journal      = {Imaging Neuroscience},
 	doi          = {10.1162/IMAG.a.948},
 }
+
+@article{eichner2015real,
+	title        = {Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast},
+	author       = {Eichner, Cornelius and Cauley, Stephen F. and Cohen-Adad, Julien and Möller, Harald E. and Turner, Robert and Setsompop, Kawin and Wald, Lawrence L.},
+	year         = {2015},
+	journal      = {NeuroImage},
+	volume       = {122},
+	pages        = {373--384},
+	doi          = {10.1016/j.neuroimage.2015.07.074},
+	pmid         = {26241680},
+}
+
+@article{sprenger2017real,
+	title        = {Real valued diffusion-weighted imaging using decorrelated phase filtering},
+	author       = {Sprenger, Tim and Sperl, Jonathan I. and Fernandez, Brice and Haase, Axel and Menzel, Marion I.},
+	year         = {2017},
+	journal      = {Magnetic Resonance in Medicine},
+	volume       = {77},
+	number       = {2},
+	pages        = {559--570},
+	doi          = {10.1002/mrm.26138},
+	pmid         = {26869192},
+}

qsiprep/data/reports-spec.yml

@@ -186,6 +186,15 @@ sections:
     static: false
     subtitle: DWI Denoising
 
+  - bids: {datatype: figures, desc: phasecorrection, suffix: dwi}
+    caption: |
+      "Phase correction of the complex-valued data. The real channel (retained for
+      downstream processing) is shown alongside the imaginary residual for a low and
+      high-<em>b</em> image. After successful correction the imaginary residual should
+      contain only noise, with no coherent anatomical structure."
+    static: false
+    subtitle: Complex Phase Correction
+
   - bids: {datatype: figures, desc: unringing, suffix: dwi}
     caption: |
       "Effect of removing Gibbs ringing on a low and high-<em>b</em> image."

qsiprep/interfaces/images.py

@@ -9,6 +9,7 @@
 
 import glob
 import os
+import re
 from subprocess import PIPE, Popen
 from textwrap import indent
 
@@ -396,6 +397,30 @@ class ConformDwiOutputSpec(TraitedSpec):
     out_report = File(exists=True, desc='HTML segment containing warning')
 
 
+def _find_gradient_file(dwi_file, ext):
+    """Locate the bval/bvec file (``ext``) associated with a DWI image.
+
+    The gradient file is normally named like the DWI image but with a ``.bval``
+    or ``.bvec`` extension. For complex-valued acquisitions the magnitude and
+    phase images share a single pair of gradient files that omit the ``part-``
+    entity (e.g. ``..._dwi.bval`` for ``..._part-mag_dwi.nii.gz``), so fall back
+    to a ``part``-stripped name when the part-specific file is absent. If
+    neither exists, the part-specific candidate is returned unchanged so the
+    caller's existing ``os.path.exists`` handling applies.
+    """
+    candidate = fname_presuffix(dwi_file, suffix=ext, use_ext=False)
+    if os.path.exists(candidate):
+        return candidate
+
+    stripped = re.sub(r'_part-[A-Za-z0-9]+', '', dwi_file)
+    if stripped != dwi_file:
+        stripped_candidate = fname_presuffix(stripped, suffix=ext, use_ext=False)
+        if os.path.exists(stripped_candidate):
+            return stripped_candidate
+
+    return candidate
+
+
 class ConformDwi(SimpleInterface):
     """Conform a series of dwi images to enable merging.
     Performs three basic functions:
@@ -419,12 +444,12 @@ def _run_interface(self, runtime):
         if isdefined(self.inputs.bval_file):
             bval_fname = self.inputs.bval_file
         else:
-            bval_fname = fname_presuffix(fname, suffix='.bval', use_ext=False)
+            bval_fname = _find_gradient_file(fname, '.bval')
 
         if isdefined(self.inputs.bvec_file):
             bvec_fname = self.inputs.bvec_file
         else:
-            bvec_fname = fname_presuffix(fname, suffix='.bvec', use_ext=False)
+            bvec_fname = _find_gradient_file(fname, '.bvec')
 
         out_bvec_fname = fname_presuffix(bvec_fname, suffix=suffix, newpath=runtime.cwd)
         validator = ValidateImage(in_file=fname)