Add phase correction options

.gitignore

@@ -1,4 +1,5 @@
 .cursor/.cursorenv
+docs/superpowers/
 
 # Test data
 qsiprep/tests/test_data/

docs/.gitignore

@@ -2,3 +2,7 @@
 *.dot
 etc/
 args.txt
+
+# Example figures referenced from the documentation (kept despite the *.png rule)
+!_static/denoising_example.png
+!_static/phasecorrection_example.png

docs/conf.py

@@ -26,6 +26,22 @@
 
 import qsiprep
 
+# Seed the global qsiprep config with sensible defaults so that the example
+# workflows rendered by the ``.. workflow::`` directives can be built. Most
+# workflow functions read parameters (b0_threshold, output_resolution, etc.)
+# from ``qsiprep.config`` rather than from their arguments, so without this the
+# example workflows fail to instantiate. ``init=False`` only populates the
+# settings; it does not validate or create the (CI-specific) paths in the file.
+from qsiprep import config as _qsiprep_config  # noqa: E402
+
+try:
+    _qsiprep_config.load(
+        os.path.join(os.path.dirname(qsiprep.__file__), 'data', 'tests', 'config.toml'),
+        init=False,
+    )
+except Exception as _exc:  # pragma: no cover - docs build best-effort
+    print(f'Could not seed qsiprep config for workflow examples: {_exc}')
+
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.

docs/tools/svg_frames_to_png.py

@@ -0,0 +1,192 @@
+"""Flatten a niworkflows "flicker" report SVG into a static side-by-side PNG.
+
+*QSIPrep*'s denoising and phase-correction reportlets are animated SVGs: a single
+``<svg>`` holds two complete mosaics -- a ``<g class="background-svg">`` (the
+"before" frame) and a ``<g class="foreground-svg">`` (the "after" frame) -- that
+a CSS opacity animation flickers between. That animation cannot be shown in the
+static HTML/PDF documentation, so this script renders each frame on its own and
+stacks them side by side into a single PNG (background on the left, foreground on
+the right), letting readers compare both frames at a glance.
+
+Rasterization is delegated to ``rsvg-convert`` (librsvg), which is found on
+``PATH`` or supplied with ``--rsvg``. The two frames are recombined with Pillow.
+
+Examples
+--------
+Regenerate the two figures bundled in the docs (writes to ``docs/_static``)::
+
+    python docs/tools/svg_frames_to_png.py
+
+Convert an arbitrary report SVG::
+
+    python docs/tools/svg_frames_to_png.py --pair INPUT.svg OUTPUT.png
+"""
+
+import argparse
+import copy
+import os
+import shutil
+import subprocess
+import sys
+import tempfile
+
+from lxml import etree
+from PIL import Image
+
+SVG_NS = 'http://www.w3.org/2000/svg'
+NSMAP = {'s': SVG_NS}
+
+# Default conversions: the example figures referenced from docs/preprocessing.rst.
+# Paths are relative to the repository's ``docs`` directory.
+DEFAULT_PAIRS = [
+    (
+        '/mnt/c/Users/tsalo/Documents/datasets/qsiprep-development/dwidenoise-out/'
+        'sub-20188/ses-1/figures/'
+        'sub-20188_ses-1_dir-AP_run-01_part-mag_desc-denoising_dwi.svg',
+        '_static/denoising_example.png',
+    ),
+    (
+        '/mnt/c/Users/tsalo/Documents/datasets/qsiprep-development/tvc-out-pi/'
+        'sub-20188/ses-1/figures/'
+        'sub-20188_ses-1_dir-AP_run-01_part-mag_desc-phasecorrection_dwi.svg',
+        '_static/phasecorrection_example.png',
+    ),
+]
+
+
+def _find_rsvg(explicit=None):
+    """Locate the ``rsvg-convert`` executable."""
+    if explicit:
+        if not os.path.isfile(explicit):
+            raise FileNotFoundError(f'rsvg-convert not found at {explicit}')
+        return explicit
+    found = shutil.which('rsvg-convert')
+    if found:
+        return found
+    # Fall back to common micromamba/conda environments that ship librsvg.
+    for env in ('fsl', 'mapenv'):
+        candidate = os.path.expanduser(f'~/micromamba/envs/{env}/bin/rsvg-convert')
+        if os.path.isfile(candidate):
+            return candidate
+    raise FileNotFoundError(
+        'Could not find "rsvg-convert". Install librsvg or pass --rsvg /path/to/rsvg-convert.'
+    )
+
+
+def _frame_svg(root, keep_class):
+    """Return a deep copy of ``root`` containing only the requested frame group.
+
+    The sibling frame is removed, and the kept group is forced to full opacity
+    with its animation class stripped, so ``rsvg-convert`` renders it as an
+    opaque still frame.
+    """
+    clone = copy.deepcopy(root)
+    drop_class = 'foreground-svg' if keep_class == 'background-svg' else 'background-svg'
+
+    for el in clone.xpath(f".//s:g[@class='{drop_class}']", namespaces=NSMAP):
+        el.getparent().remove(el)
+
+    for el in clone.xpath(f".//s:g[@class='{keep_class}']", namespaces=NSMAP):
+        # Remove the class (it triggers the CSS flicker animation) and pin opacity.
+        el.attrib.pop('class', None)
+        el.set('style', 'opacity:1')
+
+    return clone
+
+
+def _rasterize(svg_root, rsvg, zoom):
+    """Rasterize an SVG element tree to a Pillow image via ``rsvg-convert``."""
+    with tempfile.TemporaryDirectory() as tmp:
+        svg_path = os.path.join(tmp, 'frame.svg')
+        png_path = os.path.join(tmp, 'frame.png')
+        etree.ElementTree(svg_root).write(svg_path, xml_declaration=True, encoding='utf-8')
+        subprocess.run(
+            [
+                rsvg,
+                '--zoom',
+                str(zoom),
+                '--background-color',
+                'white',
+                '--output',
+                png_path,
+                svg_path,
+            ],
+            check=True,
+        )
+        return Image.open(png_path).convert('RGB')
+
+
+def _stack_side_by_side(left, right, gap, bg):
+    """Place two images side by side on a single canvas, vertically centered."""
+    height = max(left.height, right.height)
+    width = left.width + gap + right.width
+    canvas = Image.new('RGB', (width, height), bg)
+    canvas.paste(left, (0, (height - left.height) // 2))
+    canvas.paste(right, (left.width + gap, (height - right.height) // 2))
+    return canvas
+
+
+def convert(in_svg, out_png, rsvg, zoom=1.0, gap=16, bg=(255, 255, 255)):
+    """Flatten one flicker SVG into a side-by-side PNG."""
+    root = etree.parse(in_svg).getroot()
+
+    found = {
+        c
+        for c in ('background-svg', 'foreground-svg')
+        if root.xpath(f".//s:g[@class='{c}']", namespaces=NSMAP)
+    }
+    if found != {'background-svg', 'foreground-svg'}:
+        raise ValueError(
+            f'{in_svg}: expected both background-svg and foreground-svg frames, found {found}'
+        )
+
+    left = _rasterize(_frame_svg(root, 'background-svg'), rsvg, zoom)
+    right = _rasterize(_frame_svg(root, 'foreground-svg'), rsvg, zoom)
+
+    combined = _stack_side_by_side(left, right, gap, bg)
+    os.makedirs(os.path.dirname(os.path.abspath(out_png)) or '.', exist_ok=True)
+    combined.save(out_png, optimize=True)
+    print(f'{in_svg}\n  -> {out_png} ({combined.width}x{combined.height})')
+
+
+def _parse_args(argv):
+    parser = argparse.ArgumentParser(
+        description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
+    )
+    parser.add_argument(
+        '--pair',
+        nargs=2,
+        metavar=('IN_SVG', 'OUT_PNG'),
+        action='append',
+        help='An input flicker SVG and the output PNG path. May be repeated. '
+        'If omitted, the example denoising and phase-correction figures are regenerated.',
+    )
+    parser.add_argument(
+        '--rsvg', help='Path to the rsvg-convert executable (default: search PATH).'
+    )
+    parser.add_argument(
+        '--zoom', type=float, default=1.0, help='Rasterization zoom factor (default: 1.0).'
+    )
+    parser.add_argument(
+        '--gap', type=int, default=16, help='Pixel gap between the two frames (default: 16).'
+    )
+    return parser.parse_args(argv)
+
+
+def main(argv=None):
+    args = _parse_args(argv if argv is not None else sys.argv[1:])
+    rsvg = _find_rsvg(args.rsvg)
+
+    if args.pair:
+        pairs = args.pair
+    else:
+        # Resolve default output paths rel. to docs directory (this file's parent's parent)
+        docs_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+        pairs = [(src, os.path.join(docs_dir, dst)) for src, dst in DEFAULT_PAIRS]
+
+    for in_svg, out_png in pairs:
+        convert(in_svg, out_png, rsvg, zoom=args.zoom, gap=args.gap)
+
+
+if __name__ == '__main__':
+    main()

qsiprep/cli/parser.py

@@ -392,8 +392,26 @@ def _bids_filter(value, parser):
         action='store',
         choices=['dwidenoise', 'patch2self', 'none'],
         default='dwidenoise',
-        help='Image-based denoising method. Either "dwidenoise" (MRtrix), '
-        '"patch2self" (DIPY) or "none". (default: dwidenoise)',
+        help=(
+            'Image-based denoising method. '
+            'Either "dwidenoise" (MRtrix), "patch2self" (DIPY) or "none".'
+        ),
+    )
+    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. '
+            'For any method other than "none", the real channel after rephasing is '
+            'used for downstream processing instead of the magnitude. The available '
+            'methods are "tv" (Eichner 2015, total variation on the magnitude/phase '
+            'stack), "tvc" (Eichner 2015, paper-faithful total variation on the '
+            'complex signal) and "dc" (Sprenger 2017, decorrelated-phase filtering). '
+            'Requires phase data (a BIDS part-phase file) and --denoise-method '
+            'dwidenoise; otherwise it is ignored. (default: none)'
+        ),
     )
     g_conf.add_argument(
         '--unringing-method',
@@ -755,6 +773,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), "tvc" (Eichner 2015) or "dc"
+    (Sprenger 2017). For any method other than "none", 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,60 @@ @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},
+}
+
+@article{jezzard1995,
+	title        = {Correction for geometric distortion in echo planar images from {B0} field variations},
+	author       = {Jezzard, Peter and Balaban, Robert S.},
+	year         = {1995},
+	journal      = {Magnetic Resonance in Medicine},
+	volume       = {34},
+	number       = {1},
+	pages        = {65--73},
+	doi          = {10.1002/mrm.1910340111},
+	pmid         = {7674900},
+}
+
+@article{brett2001,
+	title        = {Spatial normalization of brain images with focal lesions using cost function masking},
+	author       = {Brett, Matthew and Leff, Alexander P. and Rorden, Chris and Ashburner, John},
+	year         = {2001},
+	journal      = {NeuroImage},
+	volume       = {14},
+	number       = {2},
+	pages        = {486--500},
+	doi          = {10.1006/nimg.2001.0845},
+	pmid         = {11467921},
+}
+
+@article{yeh2019differential,
+	title={Differential tractography as a track-based biomarker for neuronal injury},
+	author={Yeh, Fang-Cheng and Zaydan, Islam M and Suski, Valerie R and Lacomis, David and Richardson, R Mark and Maroon, Joseph C and Barrios-Martinez, Jessica},
+	journal={Neuroimage},
+	volume={202},
+	pages={116131},
+	year={2019},
+	publisher={Elsevier}
+}

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."