diff --git a/.gitignore b/.gitignore index e4b5907d..fc63a3be 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ .cursor/.cursorenv +docs/superpowers/ # Test data qsiprep/tests/test_data/ diff --git a/docs/.gitignore b/docs/.gitignore index acbc6d5a..45d73c0a 100644 --- a/docs/.gitignore +++ b/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 diff --git a/docs/_static/denoising_example.png b/docs/_static/denoising_example.png new file mode 100644 index 00000000..c0c88239 Binary files /dev/null and b/docs/_static/denoising_example.png differ diff --git a/docs/_static/phasecorrection_example.png b/docs/_static/phasecorrection_example.png new file mode 100644 index 00000000..ccc4817f Binary files /dev/null and b/docs/_static/phasecorrection_example.png differ diff --git a/docs/conf.py b/docs/conf.py index 20fb6665..6615ff86 100644 --- a/docs/conf.py +++ b/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. diff --git a/docs/preprocessing.rst b/docs/preprocessing.rst index 107fdaf7..6282a233 100644 --- a/docs/preprocessing.rst +++ b/docs/preprocessing.rst @@ -125,66 +125,6 @@ You can enable regular expressions for more detailed filtering, for example:: will do a case-insensitive match of "mprage" within the "t1w" query. -.. _merge_denoise: - -Denoising and Merging Images -============================ - -The user can decide whether to do certain preprocessing steps and, if so, -whether they are performed *before* or *after* the DWI series are -concatenated. Specifically, image denoising (using ``dwidenoise`` or -``patch2self``) can be disabled with ``--denoise-method none``. Gibbs -unringing (using ``mrdegibbs``) is disabled by default but can be enabled -with ``--unringing-method mrdegibbs``. B1 bias field correction is applied by -default (using ``dwibiascorrect``) and can be disabled with the -``--dwi-no-biascorr`` option. The intensity of b=0 images is harmonized -across scans (i.e. scaled to an average value) by default, but this can be -turned off using ``--dwi-no-b0-harmonization``. - -.. tip:: - - If prescan normalization is enabled, - we recommend using ``--b1-biascorrect-stage none``. - This will skip B1 bias field correction, - which may introduce artifacts on normalized data. - -Together, denoising (MP-PCA or patch2self), Gibbs unringing B1 bias field -correction, and b=0 intensity normalization are referred to as *denoising* in -*QSIPrep*. Each of these image processing operations has assumptions about its -inputs and changes the distribution of noise in its outputs. Although the -inclusion of each operation can be decided by the user, the order in which -they are applied relative to one another is fixed. MP-PCA or patch2self are -applied directly to the BIDS inputs, which should be uninterpolated and as -"raw" as possible. Although Gibbs unringing should be performed on "raw" -data, it is recommended in the MRtrix3 documentation to apply MP-PCA before -Gibbs unringing. B1 bias field correction and b=0 intensity harmonization -do not have as specific requirements about their inputs so are run last. - -The last, and potentially very important decision, is whether the denoising -operations are applied to each input DWI series individually or whether the -denoising operations are applied to the concatenated input DWI files. At -present, there is little data to guide this choice. The more volumes -available, the more data MP-PCA/patch2self have to work with. However, if -there if the head is in a vastly different location in different scans, -denoising might be impacted in unpredictable ways. - -Consider MP-PCA. If a voxel contains CSF in one DWI series and the subject -repositions their head between scans so that the voxel contains corpus -callosum in the next DWI series, the non-noise signal will be very different -in the two series. Similarly, if the head is repositioned different areas -will be closer to the head coil and therefore be inconsistently affected by -B1 bias field. Similar problems can also occur *within* a DWI series due to -subject head motion, but these methods have been shown to work well even in -the presence of within-scan head movement. If the head position changes -across scans is of a similar magnitude to that of within-scan head motion, it -is likely fine to use the ``--denoise-after-combining`` option. To gauge how -much between-scan motion occurred, users can inspect the :ref:`qc_data` to see -whether Framewise Displacement is large where a new series begins. - -By default, the scans in the same warped space are individually denoised before -they are concatenated. When warped groups are concatenated an additional b=0 -image intensity normalization is performed. - Preprocessing HCP-style ======================= @@ -404,8 +344,9 @@ A single-line tsv file (``desc-image_qc.tsv``) is created for each output image. This file is particularly useful for comparing the relative quality across subjects before deciding who to include in a group analysis. The columns in this file come from DSI Studio's QC calculation and is described -in [Yeh2019]_. Columns prefixed by ``raw_`` reflect QC measurements from the -data before preprocessing. Columns prefixed by ``t1_`` or ``mni_`` contain QC +in :footcite:t:`yeh2019differential`. +Columns prefixed by ``raw_`` reflect QC measurements from the data before preprocessing. +Columns prefixed by ``t1_`` or ``mni_`` contain QC metrics calculated on the preprocessed data. Motion parameter summaries are also provided, such as the mean and max of framewise displacement (``mean_fd``, ``max_fd``). The max and mean absolute values for translation @@ -467,20 +408,10 @@ Processing the *Subject Anatomical Reference* T1w or T2w images from qsiprep.workflows.anatomical import init_anat_preproc_wf wf = init_anat_preproc_wf( - omp_nthreads=1, - reportlets_dir='.', - output_dir='.', - dwi_only=False, - infant_mode=False, - template='MNI152NLin2009cAsym', - output_spaces=['T1w'], - output_resolution=1.25, - skull_strip_template='OASIS', - force_spatial_normalization=True, - freesurfer=True, - debug=False, - hires=True, - num_t1w=1, + num_anat_images=1, + num_additional_t2ws=0, + has_rois=False, + anatomical_template='MNI152NLin2009cAsym', ) As of version 0.18 *QSIPrep* has been changed to be very flexible with anatomical @@ -520,7 +451,7 @@ Processing the *Anatomical Reference* images be bias corrected using N4 and aligned to one another. If ``--subject-anatomical-reference unbiased`` is specified they will be unbiasedly registered to each other using ANTs. Otherwise all the images are registered to the first (alphabetically) image (see - `Longitudinal T1w processing`_). + `Longitudinal anatomical processing`_). 3. Brain extraction is performed using ``SynthStrip``. .. figure:: _static/brainextraction_t1.svg @@ -555,7 +486,7 @@ Handling Lesions and abnormalities When processing images from patients with focal brain lesions (e.g., stroke, tumor resection), it is possible to provide a lesion mask to be used during spatial -normalization to MNI-space [Brett2001]_. +normalization to MNI-space :footcite:p:`brett2001`. ANTs will use this mask to minimize warping of healthy tissue into damaged areas (or vice-versa). Lesion masks should be binary NIfTI images (damaged areas = 1, everywhere else = 0) @@ -610,54 +541,268 @@ DWI preprocessing :func:`qsiprep.workflows.dwi.base.init_dwi_preproc_wf` -.. workflow:: - :graph2use: orig - :simple_form: yes +.. code-block:: python from qsiprep.workflows.dwi.base import init_dwi_preproc_wf wf = init_dwi_preproc_wf( - dwi_only=False, scan_groups={ - 'dwi_series': ['fake.nii'], + 'dwi_series': ['sub-1_dwi.nii.gz'], 'fieldmap_info': {'suffix': None}, 'dwi_series_pedir': 'j', }, + t2w_sdc=False, + output_prefix='sub-1', source_file='/data/sub-1/dwi/sub-1_dwi.nii.gz', - output_prefix='', - ignore=[], - b0_threshold=100, - motion_corr_to='iterative', - b0_to_t1w_transform='Rigid', - hmc_model='3dSHORE', - hmc_transform='Rigid', - shoreline_iters=2, - impute_slice_threshold=0, - eddy_config=None, - reportlets_dir='.', - output_spaces=['T1w'], - dwi_denoise_window=5, - denoise_method='dwidenoise', - unringing_method='mrdegibbs', - b1_biascorr_stage='final', - no_b0_harmonization=False, - denoise_before_combining=True, - template='MNI152NLin2009cAsym', - output_dir='.', - omp_nthreads=1, - fmap_bspline=False, - fmap_demean=True, - use_syn=True, - force_syn=False, - low_mem=False, - sloppy=True, - layout=None, + anatomical_template='MNI152NLin2009cAsym', ) Preprocessing of :abbr:`DWI (Diffusion Weighted Image)` files is split into multiple sub-workflows described below. +.. _merge_denoise: + +Denoising and Merging Images +---------------------------- + +*QSIPrep*'s "denoising" stage bundles together several image-processing operations +that reduce noise and remove acquisition artifacts before head motion correction. +Each can be enabled or disabled independently, +but the order in which they are applied is fixed (see :ref:`denoise_order`). +The sections below describe each operation -- its purpose, the relevant command-line +arguments and their effects, and when to use it: + +* :ref:`denoise_image` (``--denoise-method``, ``--dwi-denoise-window``) +* :ref:`complex_denoising` (``--dwi-phase-correction``) +* :ref:`denoise_gibbs` (``--unringing-method``) +* :ref:`denoise_biascorr` (``--b1-biascorrect-stage``) +* :ref:`denoise_b0harmonize` (``--no-b0-harmonization``) + +A separate choice -- whether these operations run on each input DWI series +individually or on the concatenated series -- is covered in :ref:`denoise_combine`. + + +.. _denoise_image: + +Image Denoising +^^^^^^^^^^^^^^^ + +Thermal noise corrupts diffusion-weighted images and biases downstream model +fits, an effect that is most severe at high *b*-values where the signal is weak. +*QSIPrep* reduces this noise with one of two algorithms, chosen with +``--denoise-method`` (default ``dwidenoise``; pass ``none`` to skip denoising). +The default, ``dwidenoise``, implements Marchenko-Pastur PCA (MP-PCA) +:footcite:p:`dwidenoise1,dwidenoise2` as implemented in MRtrix3 +:footcite:p:`mrtrix3`: it performs a local PCA in a sliding window and discards +the principal components whose variance is consistent with noise, exploiting the +redundancy across the volumes of a DWI series. The window size is controlled by +``--dwi-denoise-window``, whose default ``auto`` derives it from the number of +volumes (following the ``dwidenoise`` documentation); a manual value must be an +odd, positive integer, and a larger window gives the PCA more samples but tracks +spatially varying noise less well. The alternative, ``patch2self`` +:footcite:p:`patch2self` (implemented in DIPY :footcite:p:`dipy`), is a +self-supervised method that predicts each volume from the others via regression +and assumes no a-priori noise model; it ignores ``--dwi-denoise-window``. +``dwidenoise`` is a sound default for most acquisitions, whereas ``patch2self`` +can be advantageous when many volumes or directions are available. Because +MP-PCA assumes "raw", uninterpolated input, denoising is run as early as possible +in the pipeline (see :ref:`denoise_order`). + +When DWI **phase** data are available (a BIDS ``part-phase`` file accompanying +the ``part-mag`` file), ``dwidenoise`` denoises the **complex-valued** data +rather than the magnitude alone, which reduces the non-Gaussian (Rician) +noise-floor bias that affects magnitude data :footcite:p:`cordero2019complex`. +The magnitude and phase are combined into a single complex image, denoised +jointly, and -- unless :ref:`phase correction ` is requested +-- converted back to magnitude afterward. Phase is only used when it is available +and ``phase`` is not listed in ``--ignore``. Before being combined with the +magnitude, the phase is rescaled to radians in :math:`[-\pi, \pi)`, matching the +MRtrix3 complex-``dwidenoise`` convention. The denoising report always shows +magnitude images before and after denoising, even when complex denoising was +performed. Including phase also enables the optional +:ref:`phase-correction ` step, which can further suppress the +high-*b* noise floor. + +.. figure:: _static/denoising_example.png + :width: 100% + + The denoising reportlet, with both animation frames shown side by side: the + raw images (left) and the MP-PCA-denoised images (right), each for a low-*b* + and a high-*b* volume. + + Denoising should remove grainy, incoherent noise -- + most visibly at high *b* -- while preserving anatomical structure and edges. + The improvement is typically clearest in the high-*b* panels, where the raw + signal is weakest. + + The contour lines reflect the level of noise in the image. + These lines should not follow anatomical structures. + + +.. _complex_denoising: + +Phase Correction +^^^^^^^^^^^^^^^^ + +By default, after complex denoising the data are returned to magnitude. Phase +correction, enabled with ``--dwi-phase-correction`` (default ``none``), is an +optional alternative: instead of taking the magnitude, it estimates a smooth, +spatially varying *background phase* and rephases the complex signal so that the +diffusion signal is rotated onto the real axis. The **real channel** is then +kept for the rest of the pipeline, while the imaginary channel should contain +mostly noise. Retaining the real channel rather than the magnitude avoids the +rectified noise floor entirely, because Gaussian noise on the real axis stays +zero-mean -- the main benefit of having acquired phase data. Phase correction +requires phase data and ``--denoise-method dwidenoise``, and is ignored +otherwise. + +The three methods differ in how they estimate the background phase. ``tv`` +applies total-variation rephasing to the joint magnitude/phase image +:footcite:p:`eichner2015real`; this is the original formulation, but it operates +on the (wrapped) phase as a scalar field and is therefore sensitive to phase +wrapping. ``tvc`` is the paper-faithful, complex total-variation formulation +:footcite:p:`eichner2015real`, estimating the background phase from a +total-variation-denoised *complex* image, so it is wrap-free and +intensity-weighted by construction. ``dc`` uses decorrelated-phase rephasing +:footcite:p:`sprenger2017real`, smoothing the complex signal with a convolution +kernel and likewise wrap-free, with optional outlier rejection. The ``tvc`` and +``dc`` methods are invariant to a global phase offset, whereas ``tv`` is not. In +practice ``tvc`` is the most principled choice and ``dc`` a robust, faster +alternative, while the wrap-sensitive ``tv`` is a legacy option that is generally +not recommended. Whichever method you use, inspect the phase-correction report +(below) before trusting the result. + + +Interpreting the phase-correction report +"""""""""""""""""""""""""""""""""""""""" + +The phase-correction reportlet shows the retained **real channel** next to the +**imaginary residual** for a low-*b* and a high-*b* volume. The imaginary +residual is the diagnostic: it should look like featureless, noise-like +speckling with no recognizable brain anatomy (no sharp tissue boundaries, gyral +patterns, or contrast that mirrors the magnitude image). Coherent anatomy in the +imaginary residual indicates that real signal is leaking into the imaginary +channel -- i.e. the background-phase estimate is incomplete. + +It is normal, however, for the imaginary residual to show *higher* noise inside +the brain than outside. The methods remove only a smooth background phase, so +signal-dependent, non-smooth phase variation (for example from cardiac pulsation +or bulk motion during the diffusion-encoding gradients) legitimately leaves +additional, incoherent noise wherever there is tissue signal. A brain-shaped +envelope of stronger speckling is expected; readable anatomy is not. + +.. figure:: _static/phasecorrection_example.png + :width: 100% + + The phase-correction reportlet, with both animation frames shown side by + side: the retained real channel (left) and the imaginary residual (right), + each for a low-*b* and a high-*b* volume. The real channel should look like a + clean diffusion-weighted image. The imaginary residual is the diagnostic: it + should be noise-like speckling, possibly with a brain-shaped envelope of + stronger noise, but without recognizable anatomy. Coherent structure in the + imaginary residual indicates an incomplete background-phase estimate. + + +.. _denoise_gibbs: + +Gibbs Unringing +^^^^^^^^^^^^^^^ + +Gibbs ringing is the oscillating artifact that appears near sharp intensity +boundaries -- such as CSF/tissue interfaces -- because k-space is sampled only up +to a finite extent, truncating the Fourier series. Unringing is disabled by +default and is enabled with ``--unringing-method``, which offers two algorithms. +``mrdegibbs`` :footcite:p:`mrdegibbs` (implemented in MRtrix3 +:footcite:p:`mrtrix3`) uses the local subvoxel-shift method, resampling the image +at the points where the ringing crosses zero; because it assumes the image was +reconstructed from fully sampled k-space, it is **not** appropriate for +partial-Fourier acquisitions. For those, use ``rpg`` :footcite:p:`pfgibbs` +instead, which is designed to remove partial-Fourier-induced Gibbs ringing. +Gibbs unringing should be performed on "raw" data, but MRtrix3 recommends +applying MP-PCA *before* unringing, so *QSIPrep* enforces that order (see +:ref:`denoise_order`). + + +.. _denoise_biascorr: + +B1 Bias Field Correction +^^^^^^^^^^^^^^^^^^^^^^^^^ + +Spatial inhomogeneity of the receive coil's B1 field produces a smooth, +low-frequency intensity bias across the image -- some regions appear brighter +simply because they are closer to the coil. *QSIPrep* estimates and removes this +field using ANTs' N4 algorithm :footcite:p:`n4` (via ``dwibiascorrect``), and +does so by default. The ``--b1-biascorrect-stage`` option controls when, or +whether, this happens: ``final`` (the default) applies it after the data have +been resampled to their final space, ``none`` skips it, and ``legacy`` reproduces +the behavior of *QSIPrep* < 0.17. +(The older ``--dwi-no-biascorr`` flag is deprecated in favor of ``--b1-biascorrect-stage none``.) + +.. tip:: + + For prescan-normalized data, we recommend ``--b1-biascorrect-stage none``, + because bias correction may introduce artifacts on data that have already been + intensity-normalized. + + +.. _denoise_b0harmonize: + +b=0 Intensity Harmonization +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +When multiple DWI series are concatenated, differences in overall scaling +between series can introduce artificial intensity steps. *QSIPrep* harmonizes +the series by rescaling each so that its b=0 images share a common mean +intensity. This is done by default and can be turned off with +``--no-b0-harmonization``. + + +.. _denoise_order: + +Order of operations +^^^^^^^^^^^^^^^^^^^ + +Although each operation above can be toggled independently, the order in which +they are applied is fixed. Image denoising (MP-PCA or ``patch2self``) is applied +first, directly to the BIDS inputs, which should be uninterpolated and as "raw" +as possible. Gibbs unringing should also be performed on "raw" data, but the +MRtrix3 documentation recommends applying MP-PCA before unringing, so denoising +precedes unringing. B1 bias field correction and b=0 intensity harmonization +have no such strict input requirements and are run last. Each of these +operations has assumptions about its inputs and changes the distribution of +noise in its outputs, which is why the order is not user-configurable. + + +.. _denoise_combine: + +Denoising before or after merging +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +A consequential choice is whether the denoising operations are applied to each +input DWI series individually or to the concatenated DWI series. At present +there is little data to guide this choice. The more volumes available, the more +data MP-PCA/``patch2self`` have to work with. However, if the head is in a +vastly different location across scans, denoising may be affected in +unpredictable ways. + +Consider MP-PCA. If a voxel contains CSF in one DWI series and the subject +repositions their head between scans so that the voxel contains corpus callosum +in the next series, the non-noise signal will differ greatly between the two +series. Similarly, if the head is repositioned, different areas will be closer +to the head coil and therefore be inconsistently affected by the B1 bias field. +Similar problems can also occur *within* a DWI series due to head motion, but +these methods have been shown to work well even with within-scan movement. If +the across-scan head-position change is of a similar magnitude to within-scan +motion, it is likely fine to use ``--denoise-after-combining``, which runs +denoising on the concatenated series (before motion correction) rather than on +each series separately. By default, scans in the same warped space are denoised +individually before concatenation; when warped groups are concatenated, an +additional b=0 intensity normalization is performed. To gauge how much +between-scan motion occurred, inspect the :ref:`qc_data` to see whether +Framewise Displacement is large where a new series begins. + + .. _fsl_wf: Head-motion / Eddy Current / Distortion correction (FSL) @@ -702,13 +847,8 @@ and eddy current correction will be performed. The workflow looks like this: 'fieldmap_info': {'suffix': None}, 'dwi_series_pedir': 'j', }, - b0_threshold=100, - impute_slice_threshold=0., - fmap_demean=False, - fmap_bspline=False, - eddy_config=None, source_file='/path/to/dwi/sub-X_dwi.nii.gz', - omp_nthreads=1, + t2w_sdc=False, ) @@ -735,13 +875,8 @@ dedicated fieldmaps (in the ``fmap/`` directory) or DWI series }, 'concatenated_bids_name': 'sub-1', }, - b0_threshold=100, - impute_slice_threshold=0., - fmap_demean=False, - fmap_bspline=False, - eddy_config=None, source_file='/path/to/dwi/sub-X_dwi.nii.gz', - omp_nthreads=1, + t2w_sdc=False, ) @@ -752,9 +887,7 @@ If a GRE fieldmap or SyN-based fieldmapless distortion correction are detected, these will be performed on the outputs of ``eddy``. For details see :ref:`dwi_sdc`. -.. workflow:: - :graph2use: orig - :simple_form: yes +.. code-block:: python from qsiprep.workflows.dwi.fsl import init_fsl_hmc_wf @@ -772,13 +905,8 @@ For details see :ref:`dwi_sdc`. 'dwi_series_pedir': 'j', 'concatenated_bids_name': 'sub-1_dir-AP', }, - b0_threshold=100, - impute_slice_threshold=0., - fmap_demean=False, - fmap_bspline=False, - eddy_config=None, source_file='/path/to/dwi/sub-X_dwi.nii.gz', - omp_nthreads=1, + t2w_sdc=False, ) @@ -848,20 +976,8 @@ are saved for each slice for display in a carpet plot-like thing. 'dwi_series_pedir': 'j', }, source_file='/data/sub-1/dwi/sub-1_dwi.nii.gz', - b0_threshold=100, - hmc_transform='Affine', - hmc_model='3dSHORE', - hmc_align_to='iterative', - template='MNI152NLin2009cAsym', - shoreline_iters=1, - impute_slice_threshold=0, - omp_nthreads=1, - fmap_bspline=False, - fmap_demean=False, - use_syn=True, - force_syn=False, - name='qsiprep_hmcsdc_wf', - dwi_metadata={}, + t2w_sdc=False, + anatomical_template='MNI152NLin2009cAsym', ) @@ -878,7 +994,7 @@ The are three kinds of SDC available in *QSIPrep*: 1. :ref:`sdc_pepolar` (also called **blip-up/blip-down**): This is the implementation from sdcflows, using 3dQwarp to - correct a DWI series using a fieldmap in the fmaps directory [Jezzard1995]_. + correct a DWI series using a fieldmap in the fmaps directory :footcite:p:`jezzard1995`. The reverse phase encoding direction scan can come from the fieldmaps directory or the dwi directory. If using :ref:`fsl_wf`, then ``TOPUP`` is used for this correction. Also relevant is :ref:`best_b0`. @@ -953,13 +1069,11 @@ DWI reference image estimation :graph2use: orig :simple_form: yes - from qsiprep.workflows.anatomical import init_dwi_reference_wf + from qsiprep.workflows.dwi.util import init_dwi_reference_wf wf = init_dwi_reference_wf( - omp_nthreads=1, gen_report=True, source_file="sub-1_dwi.nii.gz", - register_t1=True, ) This workflow estimates a reference image for a DWI series. This @@ -983,13 +1097,10 @@ Pre-processed DWIs in a different space wf = init_dwi_trans_wf( source_file='sub-1_dwi.nii.gz', + mem_gb=3, template="ACPC", - output_resolution=1.2, use_compression=True, - to_mni=False, write_local_bvecs=True, - mem_gb=3, - omp_nthreads=1, ) A DWI series is resampled to an output space. The ``output_resolution`` is @@ -1019,8 +1130,6 @@ b0 to T1w registration from qsiprep.workflows.dwi.registration import init_b0_to_anat_registration_wf wf = init_b0_to_anat_registration_wf( - mem_gb=3, - omp_nthreads=1, transform_type="Rigid", write_report=False, ) diff --git a/docs/tools/svg_frames_to_png.py b/docs/tools/svg_frames_to_png.py new file mode 100644 index 00000000..8995a53a --- /dev/null +++ b/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 +```` holds two complete mosaics -- a ```` (the +"before" frame) and a ```` (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() diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index 7ffaed50..2ef9fe91 100644 --- a/qsiprep/cli/parser.py +++ b/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 diff --git a/qsiprep/config.py b/qsiprep/config.py index bac7ced4..f0c5a81f 100644 --- a/qsiprep/config.py +++ b/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 diff --git a/qsiprep/data/boilerplate.bib b/qsiprep/data/boilerplate.bib index 026d6f27..7bb1d5d8 100644 --- a/qsiprep/data/boilerplate.bib +++ b/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} +} diff --git a/qsiprep/data/reports-spec.yml b/qsiprep/data/reports-spec.yml index 6eaa8f2a..6074736a 100644 --- a/qsiprep/data/reports-spec.yml +++ b/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-b 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-b image." diff --git a/qsiprep/interfaces/dwi_merge.py b/qsiprep/interfaces/dwi_merge.py index 2fbfbfe3..c574bbae 100644 --- a/qsiprep/interfaces/dwi_merge.py +++ b/qsiprep/interfaces/dwi_merge.py @@ -846,10 +846,11 @@ 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. + -pi and pi, respectively, and scales the image to be between -pi and pi. - STATEMENT OF CHANGES: This class is derived from sources licensed under the Apache-2.0 terms, - and the code has not been changed. + STATEMENT OF CHANGES: This class is derived from sources licensed under the Apache-2.0 terms. + The rescaling target was changed from [0, 2*pi] to [-pi, pi] to match the phase convention + used by MRtrix's complex ``dwidenoise`` workflow (``mrcalc ... pi 4096 -div -mult -polar``). Notes ----- @@ -890,11 +891,11 @@ def _run_interface(self, runtime): 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())) + # Rescale to [-pi, pi] + data = (data - data.min()) * (2 * np.pi / (data.max() - data.min())) - np.pi # Round to float32 and clip - data = np.clip(np.float32(data), 0.0, 2 * np.pi) + data = np.clip(np.float32(data), -np.pi, np.pi) hdr.set_data_dtype(np.float32) hdr.set_xyzt_units('mm') diff --git a/qsiprep/interfaces/images.py b/qsiprep/interfaces/images.py index 16d408a8..dcf40633 100644 --- a/qsiprep/interfaces/images.py +++ b/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) diff --git a/qsiprep/interfaces/phase.py b/qsiprep/interfaces/phase.py new file mode 100644 index 00000000..fdf56a15 --- /dev/null +++ b/qsiprep/interfaces/phase.py @@ -0,0 +1,182 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Interfaces for complex-valued phase correction of DWI data. + +The array-based rephasing helpers driven by :class:`PhaseCorrect` live in +:mod:`qsiprep.utils.phase`. +""" + +import os + +import nibabel as nb +import numpy as np +from nilearn.image import threshold_img +from nipype import logging +from nipype.interfaces.base import ( + BaseInterfaceInputSpec, + File, + SimpleInterface, + TraitedSpec, + traits, +) +from nipype.utils.filemanip import fname_presuffix +from niworkflows.viz.utils import compose_view, cuts_from_bbox + +from ..utils.phase import rephase_dc, rephase_tv, rephase_tv_complex +from ..viz.utils import plot_denoise + +LOGGER = logging.getLogger('nipype.interface') + + +class _PhaseCorrectInputSpec(BaseInterfaceInputSpec): + complex_file = File(exists=True, mandatory=True, desc='complex-valued denoised DWI NIfTI') + method = traits.Enum('tv', 'tvc', 'dc', mandatory=True, desc='phase-correction method') + tv_weight = traits.Float(6.0, usedefault=True, desc='TV regularization weight') + dc_kernel = traits.Enum( + 'Opt5', + 'B3', + 'B5', + 'G3F1', + 'G5F2', + 'G3F1H', + 'G5F2H', + 'Opt3', + usedefault=True, + desc='DC convolution kernel', + ) + dc_outlier_detection = traits.Bool(False, usedefault=True, desc='enable DC outlier detection') + out_report = File( + 'phasecorrection_report.svg', + usedefault=True, + desc='filename for the visual report', + ) + + +class _PhaseCorrectOutputSpec(TraitedSpec): + out_file = File(exists=True, desc='real channel after phase correction') + out_report = File(desc='visual report') + + +class PhaseCorrect(SimpleInterface): + """Apply phase correction to complex DWI data and keep the real channel. + + Reads a complex-valued NIfTI (e.g. the output of complex ``dwidenoise``), + estimates a smooth background phase using the selected method, rephases the + data, and writes the real channel as a float32 NIfTI. The imaginary channel + should then contain mostly Gaussian noise. + + The background-phase estimators are implemented in + :mod:`qsiprep.utils.phase` (:func:`~qsiprep.utils.phase.rephase_tv`, + :func:`~qsiprep.utils.phase.rephase_tv_complex`, and + :func:`~qsiprep.utils.phase.rephase_dc`). + """ + + input_spec = _PhaseCorrectInputSpec + output_spec = _PhaseCorrectOutputSpec + + def _rephase_volume(self, cvol): + """Phase-correct a single 3D complex volume, returning (real, imag).""" + if self.inputs.method == 'tv': + real, imag, _ = rephase_tv(cvol, self.inputs.tv_weight) + elif self.inputs.method == 'tvc': + real, imag, _ = rephase_tv_complex(cvol, self.inputs.tv_weight) + else: + real, imag, _ = rephase_dc( + cvol, + self.inputs.dc_kernel, + self.inputs.dc_outlier_detection, + ) + return real, imag + + def _run_interface(self, runtime): + img = nb.load(self.inputs.complex_file) + # IMPORTANT: get_fdata() would coerce to float and drop the imaginary + # part; read the raw (complex) array instead. Keep it as complex64 (not + # complex128) and process one volume at a time: upcasting a whole 4D DWI + # series to complex128 would need tens of GB for a typical acquisition. + cdata = np.asanyarray(img.dataobj) + if not np.iscomplexobj(cdata): + cdata = cdata.astype(np.complex64) + elif cdata.dtype != np.complex64: + cdata = cdata.astype(np.complex64) + + hdr = img.header.copy() + hdr.set_data_dtype(np.float32) + out_file = fname_presuffix( + self.inputs.complex_file, + suffix='_real.nii.gz', + newpath=runtime.cwd, + use_ext=False, + ) + + if cdata.ndim == 3: + real, imag = self._rephase_volume(cdata) + out = real.astype(np.float32, copy=False) + report = (real, imag, real, imag) + elif cdata.ndim == 4: + nvol = cdata.shape[3] + out = np.empty(cdata.shape, dtype=np.float32) + vol_means = np.empty(nvol, dtype=np.float64) + for vol in range(nvol): + real_v, _imag_v = self._rephase_volume(cdata[..., vol]) + out[..., vol] = real_v + vol_means[vol] = real_v.mean() + # Recompute only the two report volumes to recover their imaginary + # residuals (cheaper than retaining the full imaginary series). + lowb_index = int(np.argmax(vol_means)) + highb_index = int(np.argmin(vol_means)) + real_low, imag_low = self._rephase_volume(cdata[..., lowb_index]) + real_high, imag_high = self._rephase_volume(cdata[..., highb_index]) + report = (real_low, imag_low, real_high, imag_high) + else: + raise ValueError('PhaseCorrect expects a 3D or 4D complex image') + + nb.Nifti1Image(out, img.affine, hdr).to_filename(out_file) + self._results['out_file'] = out_file + + try: + self._generate_report(report, img.affine) + self._results['out_report'] = self._out_report + except Exception as exc: # reports are non-critical + LOGGER.warning('Phase-correction reportlet failed: %s', exc) + + return runtime + + def _generate_report(self, report, affine): + """Build a real-channel vs imaginary-residual reportlet.""" + self._out_report = os.path.abspath(self.inputs.out_report) + + real_low, imag_low, real_high, imag_high = report + real_lowb = nb.Nifti1Image(real_low, affine) + real_highb = nb.Nifti1Image(real_high, affine) + imag_lowb = nb.Nifti1Image(imag_low, affine) + imag_highb = nb.Nifti1Image(imag_high, affine) + + mask_nii = threshold_img(real_lowb, 1e-3) + cuts = cuts_from_bbox(mask_nii, cuts=7) + + compose_view( + plot_denoise( + real_lowb, + real_highb, + 'moving-image', + estimate_brightness=True, + cuts=cuts, + label='Real (phase-corrected)', + lowb_contour=None, + highb_contour=None, + compress=False, + ), + plot_denoise( + imag_lowb, + imag_highb, + 'fixed-image', + estimate_brightness=True, + cuts=cuts, + label='Imaginary residual', + lowb_contour=None, + highb_contour=None, + compress=False, + ), + out_file=self._out_report, + ) diff --git a/qsiprep/tests/test_boilerplate_phase.py b/qsiprep/tests/test_boilerplate_phase.py new file mode 100644 index 00000000..79da6d6a --- /dev/null +++ b/qsiprep/tests/test_boilerplate_phase.py @@ -0,0 +1,9 @@ +"""Ensure phase-correction citations exist in the boilerplate bibliography.""" + +from qsiprep.data import load as load_data + + +def test_phase_citations_present(): + bib = load_data('boilerplate.bib').read_text() + assert '@article{eichner2015real' in bib + assert '@article{sprenger2017real' in bib diff --git a/qsiprep/tests/test_config_phase.py b/qsiprep/tests/test_config_phase.py new file mode 100644 index 00000000..53a69188 --- /dev/null +++ b/qsiprep/tests/test_config_phase.py @@ -0,0 +1,9 @@ +"""Tests for phase-correction config defaults.""" + +from qsiprep import config + + +def test_phase_correction_config_defaults(): + assert config.workflow.dwi_phase_correction is None + assert config.workflow.dwi_phase_tv_weight == 6.0 + assert config.workflow.dwi_phase_dc_kernel == 'Opt5' diff --git a/qsiprep/tests/test_interfaces_images.py b/qsiprep/tests/test_interfaces_images.py new file mode 100644 index 00000000..1758e534 --- /dev/null +++ b/qsiprep/tests/test_interfaces_images.py @@ -0,0 +1,40 @@ +"""Tests for qsiprep.interfaces.images.""" + +from qsiprep.interfaces.images import _find_gradient_file + + +def _touch(path): + path.write_text('0\n') + return str(path) + + +def test_find_gradient_file_part_specific(tmp_path): + """Prefers the part-specific gradient file when it exists.""" + dwi = str(tmp_path / 'sub-1_part-mag_dwi.nii.gz') + part_bval = _touch(tmp_path / 'sub-1_part-mag_dwi.bval') + _touch(tmp_path / 'sub-1_dwi.bval') # part-less also present + assert _find_gradient_file(dwi, '.bval') == part_bval + + +def test_find_gradient_file_strips_part_entity(tmp_path): + """Falls back to the part-stripped gradient file (complex DWI convention).""" + dwi = str(tmp_path / 'sub-1_dir-AP_part-mag_dwi.nii.gz') + # Only the part-less gradient files exist (shared by mag/phase). + bval = _touch(tmp_path / 'sub-1_dir-AP_dwi.bval') + bvec = _touch(tmp_path / 'sub-1_dir-AP_dwi.bvec') + assert _find_gradient_file(dwi, '.bval') == bval + assert _find_gradient_file(dwi, '.bvec') == bvec + + +def test_find_gradient_file_missing_returns_part_specific(tmp_path): + """When nothing exists, returns the part-specific candidate unchanged.""" + dwi = str(tmp_path / 'sub-1_part-mag_dwi.nii.gz') + expected = str(tmp_path / 'sub-1_part-mag_dwi.bval') + assert _find_gradient_file(dwi, '.bval') == expected + + +def test_find_gradient_file_no_part_entity(tmp_path): + """A DWI without a part entity resolves to the sibling gradient file.""" + dwi = str(tmp_path / 'sub-1_dwi.nii.gz') + bval = _touch(tmp_path / 'sub-1_dwi.bval') + assert _find_gradient_file(dwi, '.bval') == bval diff --git a/qsiprep/tests/test_interfaces_phase.py b/qsiprep/tests/test_interfaces_phase.py new file mode 100644 index 00000000..6f875695 --- /dev/null +++ b/qsiprep/tests/test_interfaces_phase.py @@ -0,0 +1,60 @@ +"""Tests for qsiprep.interfaces.phase.""" + +import numpy as np + + +def _synth_complex(shape=(32, 32, 8, 6), seed=0): + """Build synthetic complex data: positive true signal on the real axis, + rotated by a smooth background phase, plus Gaussian noise on both channels. + Returns (complex_data, true_real_signal).""" + rng = np.random.default_rng(seed) + nr, nc, ns, nb = shape + # True (real, positive) signal: a smooth blob per volume + yy, xx = np.mgrid[0:nr, 0:nc] + blob = np.exp(-(((xx - nc / 2) ** 2 + (yy - nr / 2) ** 2) / (2 * (nr / 5) ** 2))) + true = np.zeros(shape) + for v in range(nb): + for z in range(ns): + true[:, :, z, v] = (1.0 + 0.3 * v) * blob + # Smooth background phase (low spatial frequency), same across volumes + bg = 0.8 * np.sin(2 * np.pi * xx / nc) + 0.8 * np.cos(2 * np.pi * yy / nr) + bg = bg[:, :, None, None] * np.ones((1, 1, ns, nb)) + clean = true * np.exp(1j * bg) + # Moderate SNR: the regime where phase correction beats magnitude (the + # magnitude noise floor biases the signal). At very high SNR the magnitude + # already tracks the truth and the comparison is not meaningful. + noise = rng.normal(scale=0.1, size=shape) + 1j * rng.normal(scale=0.1, size=shape) + return clean + noise, true + + +def test_phasecorrect_interface_writes_real_channel(tmp_path): + import nibabel as nb + + from qsiprep.interfaces.phase import PhaseCorrect + + carr, true = _synth_complex() + affine = np.eye(4) + complex_file = tmp_path / 'complex.nii.gz' + nb.Nifti1Image(carr.astype(np.complex64), affine).to_filename(complex_file) + + interface = PhaseCorrect( + complex_file=str(complex_file), + method='tv', + tv_weight=6.0, + ) + results = interface.run(cwd=str(tmp_path)) + + import os + + out_file = results.outputs.out_file + assert os.path.isfile(out_file) + out_img = nb.load(out_file) + assert out_img.shape == carr.shape + # Output is real-valued float32 + assert np.issubdtype(out_img.dataobj.dtype, np.floating) + out_data = out_img.get_fdata() + err_real = np.mean((out_data - true) ** 2) + err_mag = np.mean((np.abs(carr) - true) ** 2) + assert err_real < err_mag + # A reportlet was produced + assert os.path.isfile(results.outputs.out_report) diff --git a/qsiprep/tests/test_parser_phase.py b/qsiprep/tests/test_parser_phase.py new file mode 100644 index 00000000..1399b577 --- /dev/null +++ b/qsiprep/tests/test_parser_phase.py @@ -0,0 +1,61 @@ +"""Tests for the --dwi-phase-correction CLI flag.""" + +import os +import tempfile + +import pytest + +from qsiprep.cli.parser import _build_parser + + +def test_phase_correction_flag_default(): + with tempfile.TemporaryDirectory() as tmpdir: + bids_dir = os.path.join(tmpdir, 'bids') + out_dir = os.path.join(tmpdir, 'out') + os.makedirs(bids_dir) + + parser = _build_parser() + opts = parser.parse_args([bids_dir, out_dir, 'participant', '--output-resolution', '1.2']) + assert opts.dwi_phase_correction == 'none' + + +def test_phase_correction_flag_set(): + with tempfile.TemporaryDirectory() as tmpdir: + bids_dir = os.path.join(tmpdir, 'bids') + out_dir = os.path.join(tmpdir, 'out') + os.makedirs(bids_dir) + + parser = _build_parser() + opts = parser.parse_args( + [ + bids_dir, + out_dir, + 'participant', + '--output-resolution', + '1.2', + '--dwi-phase-correction', + 'tv', + ] + ) + assert opts.dwi_phase_correction == 'tv' + + +def test_phase_correction_flag_rejects_bad_value(): + with tempfile.TemporaryDirectory() as tmpdir: + bids_dir = os.path.join(tmpdir, 'bids') + out_dir = os.path.join(tmpdir, 'out') + os.makedirs(bids_dir) + + parser = _build_parser() + with pytest.raises(SystemExit): + parser.parse_args( + [ + bids_dir, + out_dir, + 'participant', + '--output-resolution', + '1.2', + '--dwi-phase-correction', + 'bogus', + ] + ) diff --git a/qsiprep/tests/test_utils_phase.py b/qsiprep/tests/test_utils_phase.py new file mode 100644 index 00000000..2baa1bc5 --- /dev/null +++ b/qsiprep/tests/test_utils_phase.py @@ -0,0 +1,132 @@ +"""Tests for qsiprep.utils.phase.""" + +import numpy as np + + +def _synth_complex(shape=(32, 32, 8, 6), seed=0): + """Build synthetic complex data: positive true signal on the real axis, + rotated by a smooth background phase, plus Gaussian noise on both channels. + Returns (complex_data, true_real_signal).""" + rng = np.random.default_rng(seed) + nr, nc, ns, nb = shape + # True (real, positive) signal: a smooth blob per volume + yy, xx = np.mgrid[0:nr, 0:nc] + blob = np.exp(-(((xx - nc / 2) ** 2 + (yy - nr / 2) ** 2) / (2 * (nr / 5) ** 2))) + true = np.zeros(shape) + for v in range(nb): + for z in range(ns): + true[:, :, z, v] = (1.0 + 0.3 * v) * blob + # Smooth background phase (low spatial frequency), same across volumes + bg = 0.8 * np.sin(2 * np.pi * xx / nc) + 0.8 * np.cos(2 * np.pi * yy / nr) + bg = bg[:, :, None, None] * np.ones((1, 1, ns, nb)) + clean = true * np.exp(1j * bg) + # Moderate SNR: the regime where phase correction beats magnitude (the + # magnitude noise floor biases the signal). At very high SNR the magnitude + # already tracks the truth and the comparison is not meaningful. + noise = rng.normal(scale=0.1, size=shape) + 1j * rng.normal(scale=0.1, size=shape) + return clean + noise, true + + +def _anatomy_imprint_volume(): + """High-dynamic-range structured magnitude + smooth phase (single slice).""" + n = 64 + yy, xx = np.mgrid[0:n, 0:n] + mag = np.full((n, n), 400.0) + mag[16:48, 16:48] = 8000.0 + mag[26:38, 26:38] = 20000.0 + bg = 0.8 * np.sin(2 * np.pi * xx / n) + 0.8 * np.cos(2 * np.pi * yy / n) + return (mag * np.exp(1j * bg))[:, :, None], mag + + +def test_rephase_tv_recovers_real_signal(): + from qsiprep.utils.phase import rephase_tv + + carr, true = _synth_complex() + real, imag, phase_bg = rephase_tv(carr, weight=6.0) + + assert real.shape == carr.shape + assert imag.shape == carr.shape + # Real channel should track the true signal better than naive magnitude + mag = np.abs(carr) + err_real = np.mean((real - true) ** 2) + err_mag = np.mean((mag - true) ** 2) + assert err_real < err_mag + # Imaginary residual should be near zero-mean + assert abs(np.mean(imag)) < 0.05 + + +def test_rephase_tv_does_not_imprint_magnitude_anatomy(): + """High-dynamic-range magnitude + smooth phase must not leak anatomy. + + Regression guard for the background-phase estimator: with a structured, + high-dynamic-range magnitude and a purely smooth background phase (and no + noise), perfect correction leaves the imaginary channel ~zero. A vectorial + 2-channel TV (``channel_axis=-1``) couples the magnitude edges into the + phase estimate and leaks the magnitude structure into the imaginary channel; + the volumetric estimate used by ``rephase_tv`` must not. + """ + from qsiprep.utils.phase import rephase_tv + + carr, mag = _anatomy_imprint_volume() + _real, imag, _bg = rephase_tv(carr, weight=6.0) + + brain = mag > 4000.0 + leak = np.mean(np.abs(imag[:, :, 0][brain])) / np.mean(mag[brain]) + # channel_axis=-1 gives ~0.04 here; channel_axis=None gives <0.01. + assert leak < 0.02 + + +def test_rephase_tv_complex_recovers_real_signal(): + from qsiprep.utils.phase import rephase_tv_complex + + carr, true = _synth_complex() + real, imag, _phase_bg = rephase_tv_complex(carr, weight=6.0) + + assert real.shape == carr.shape + mag = np.abs(carr) + assert np.mean((real - true) ** 2) < np.mean((mag - true) ** 2) + assert abs(np.mean(imag)) < 0.05 + + +def test_rephase_tv_complex_does_not_imprint_magnitude_anatomy(): + """Paper-faithful complex TV must also not leak magnitude anatomy.""" + from qsiprep.utils.phase import rephase_tv_complex + + carr, mag = _anatomy_imprint_volume() + _real, imag, _bg = rephase_tv_complex(carr, weight=6.0) + brain = mag > 4000.0 + leak = np.mean(np.abs(imag[:, :, 0][brain])) / np.mean(mag[brain]) + assert leak < 0.02 + + +def test_dc_kernels_present_and_normalized(): + from qsiprep.utils.phase import DC_KERNELS + + expected = {'B3', 'B5', 'G3F1', 'G5F2', 'G3F1H', 'G5F2H', 'Opt3', 'Opt5'} + assert expected == set(DC_KERNELS) + # Boxcar kernels sum to 1 + assert abs(DC_KERNELS['B3'].sum() - 1.0) < 1e-9 + assert DC_KERNELS['B3'].shape == (3, 3) + assert DC_KERNELS['Opt5'].shape == (5, 5) + + +def test_rephase_dc_recovers_real_signal(): + from qsiprep.utils.phase import rephase_dc + + carr, true = _synth_complex() + real, imag, phase_bg = rephase_dc(carr, kernel_name='Opt5', outlier_detection=False) + + assert real.shape == carr.shape + mag = np.abs(carr) + err_real = np.mean((real - true) ** 2) + err_mag = np.mean((mag - true) ** 2) + assert err_real < err_mag + assert abs(np.mean(imag)) < 0.05 + + +def test_rephase_dc_outlier_detection_runs(): + from qsiprep.utils.phase import rephase_dc + + carr, _ = _synth_complex() + real, imag, phase_bg = rephase_dc(carr, kernel_name='Opt3', outlier_detection=True) + assert real.shape == carr.shape diff --git a/qsiprep/tests/test_workflows_phase.py b/qsiprep/tests/test_workflows_phase.py new file mode 100644 index 00000000..d85f85d2 --- /dev/null +++ b/qsiprep/tests/test_workflows_phase.py @@ -0,0 +1,83 @@ +"""Workflow-build tests for phase correction in the denoising workflow.""" + +import pytest + +from qsiprep import config +from qsiprep.workflows.dwi.merge import init_dwi_denoising_wf + + +@pytest.fixture +def _denoise_defaults(monkeypatch): + monkeypatch.setattr(config.workflow, 'denoise_method', 'dwidenoise', raising=False) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none', raising=False) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5, raising=False) + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True, raising=False) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1, raising=False) + + +def _node_names(wf): + # Top-level nodes in this workflow have flat names (no subworkflow prefix). + return set(wf.list_node_names()) + + +def _build(use_phase): + return init_dwi_denoising_wf( + source_file='/data/sub-01/dwi/sub-01_part-mag_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=10, + use_phase=use_phase, + do_biascorr=False, + name='test_denoise', + ) + + +@pytest.mark.usefixtures('_denoise_defaults') +def test_phase_correction_node_present_when_enabled(monkeypatch): + monkeypatch.setattr(config.workflow, 'dwi_phase_correction', 'tv', raising=False) + wf = _build(use_phase=True) + names = _node_names(wf) + assert 'phase_correct' in names + assert 'split_complex' not in names + + +@pytest.mark.usefixtures('_denoise_defaults') +def test_magnitude_path_unchanged_when_disabled(monkeypatch): + monkeypatch.setattr(config.workflow, 'dwi_phase_correction', 'none', raising=False) + wf = _build(use_phase=True) + names = _node_names(wf) + assert 'split_complex' in names + assert 'phase_correct' not in names + + +from nipype.interfaces.base import isdefined + +from qsiprep.workflows.dwi.derivatives import init_dwi_derivatives_wf + + +def _get_node(wf, name): + return wf.get_node(name) + + +def test_derivatives_tagged_part_real_when_enabled(monkeypatch): + monkeypatch.setattr(config.workflow, 'dwi_phase_correction', 'tv', raising=False) + monkeypatch.setattr(config.workflow, 'hmc_model', 'eddy', raising=False) + wf = init_dwi_derivatives_wf(source_file='/data/sub-01/dwi/sub-01_part-mag_dwi.nii.gz') + for node_name in ( + 'ds_dwi_t1', + 'ds_bvals_t1', + 'ds_bvecs_t1', + 'ds_gradient_table_t1', + 'ds_btable_t1', + ): + assert _get_node(wf, node_name).inputs.part == 'real' + # Derived references/maps are not tagged + assert not isdefined(_get_node(wf, 'ds_t1_b0_ref').inputs.part) + assert not isdefined(_get_node(wf, 'ds_dwi_mask_t1').inputs.part) + + +def test_derivatives_not_tagged_when_disabled(monkeypatch): + monkeypatch.setattr(config.workflow, 'dwi_phase_correction', 'none', raising=False) + monkeypatch.setattr(config.workflow, 'hmc_model', 'eddy', raising=False) + wf = init_dwi_derivatives_wf(source_file='/data/sub-01/dwi/sub-01_part-mag_dwi.nii.gz') + assert not isdefined(_get_node(wf, 'ds_dwi_t1').inputs.part) diff --git a/qsiprep/utils/phase.py b/qsiprep/utils/phase.py new file mode 100644 index 00000000..ed1d8540 --- /dev/null +++ b/qsiprep/utils/phase.py @@ -0,0 +1,498 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# STATEMENT OF CHANGES: The rephasing helpers in this file (``rephase_tv``, +# ``rephase_tv_complex``, ``rephase_dc``, ``_dc_outlier_real``, the per-slice +# convolution helper, and the ``DC_KERNELS`` coefficients) are derived from the +# BSD-2-Clause-licensed MRItools package by Francesco Grussu (University College +# London), specifically the scripts: +# - ``tools/imgrephaseTV.py`` (Eichner C et al, NeuroImage 2015, 122:373-384) +# - ``tools/imgrephaseDC.py`` (Sprenger T et al, Magnetic Resonance in +# Medicine 2017, 77:559-570) +# https://github.com/fragrussu/MRItools +# The original file-based command-line scripts were refactored into array-based +# functions operating on a single complex NIfTI, file I/O was removed, the +# per-slice loops were preserved, the ``rephase_tv_complex`` variant was added, +# and the scikit-image calls were updated from the deprecated ``multichannel`` +# argument to the modern ``channel_axis`` API. +# +# ORIGINAL WORK'S ATTRIBUTION NOTICE: +# +# Code released under BSD Two-Clause license. +# +# Copyright (c) 2020 University College London. +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +"""Utilities for complex-valued phase correction of DWI data. + +The rephasing helpers below are adapted from the BSD-2-Clause-licensed MRItools +package by Francesco Grussu (University College London): +``imgrephaseTV.py`` (Eichner C et al, NeuroImage 2015, 122:373-384) and +``imgrephaseDC.py`` (Sprenger T et al, Magnetic Resonance in Medicine 2017, +77:559-570). See the BSD 2-Clause attribution notice at the top of this file. + +These are pure-array utility functions; the Nipype interface that drives them +lives in :mod:`qsiprep.interfaces.phase`. +""" + +import numpy as np +from scipy import ndimage, signal +from skimage.restoration import denoise_tv_bregman + +DC_KERNELS = { + 'B3': np.full((3, 3), 1.0 / 9.0), + 'B5': np.full((5, 5), 1.0 / 25.0), + 'G3F1': np.array( + [ + [0.075113607954111, 0.123841403152974, 0.075113607954111], + [0.123841403152974, 0.204179955571658, 0.123841403152974], + [0.075113607954111, 0.123841403152974, 0.075113607954111], + ] + ), + 'G5F2': np.array( + [ + [ + 0.023246839878294, + 0.033823952439922, + 0.038327559383904, + 0.033823952439922, + 0.023246839878294, + ], + [ + 0.033823952439922, + 0.049213560408541, + 0.055766269846849, + 0.049213560408541, + 0.033823952439922, + ], + [ + 0.038327559383904, + 0.055766269846849, + 0.063191462410265, + 0.055766269846849, + 0.038327559383904, + ], + [ + 0.033823952439922, + 0.049213560408541, + 0.055766269846849, + 0.049213560408541, + 0.033823952439922, + ], + [ + 0.023246839878294, + 0.033823952439922, + 0.038327559383904, + 0.033823952439922, + 0.023246839878294, + ], + ] + ), + 'G3F1H': np.array( + [ + [0.075113607954111, 0.123841403152974, 0.075113607954111], + [0.123841403152974, 0.0, 0.123841403152974], + [0.075113607954111, 0.123841403152974, 0.075113607954111], + ] + ), + 'G5F2H': np.array( + [ + [ + 0.023246839878294, + 0.033823952439922, + 0.038327559383904, + 0.033823952439922, + 0.023246839878294, + ], + [ + 0.033823952439922, + 0.049213560408541, + 0.055766269846849, + 0.049213560408541, + 0.033823952439922, + ], + [0.038327559383904, 0.055766269846849, 0.0, 0.055766269846849, 0.038327559383904], + [ + 0.033823952439922, + 0.049213560408541, + 0.055766269846849, + 0.049213560408541, + 0.033823952439922, + ], + [ + 0.023246839878294, + 0.033823952439922, + 0.038327559383904, + 0.033823952439922, + 0.023246839878294, + ], + ] + ), + 'Opt3': np.array( + [ + [0.107235538162453, 0.142764461837547, 0.107235538162453], + [0.142764461837547, 0.0, 0.142764461837547], + [0.107235538162453, 0.142764461837547, 0.107235538162453], + ] + ), + 'Opt5': np.array( + [ + [ + 0.025441320175391, + 0.037016902431746, + 0.041945645727859, + 0.037016902431746, + 0.025441320175391, + ], + [ + 0.037016902431746, + 0.053859275233950, + 0.054719953999307, + 0.053859275233950, + 0.037016902431746, + ], + [0.041945645727859, 0.054719953999307, 0.0, 0.054719953999307, 0.041945645727859], + [ + 0.037016902431746, + 0.053859275233950, + 0.054719953999307, + 0.053859275233950, + 0.037016902431746, + ], + [ + 0.025441320175391, + 0.037016902431746, + 0.041945645727859, + 0.037016902431746, + 0.025441320175391, + ], + ] + ), +} + + +def rephase_tv(carr, weight): + """Rephase complex MRI data using total-variation phase smoothing. + + Parameters + ---------- + carr : :obj:`numpy.ndarray` of shape (X, Y, Z) or (X, Y, Z, V) + Complex-valued image, 3D (single volume) or 4D (a series of ``V`` + volumes). The array is internally cast to ``complex128``. For a 4D + input, each volume along the last axis is processed independently. + weight : :obj:`float` + Total-variation regularization weight (the ``lambda`` of + :footcite:t:`eichner2015real`). Larger values yield a smoother estimated + background phase. The weight passed to scikit-image's + :func:`~skimage.restoration.denoise_tv_bregman` is ``2.0 * weight``. + + Returns + ------- + real : :obj:`numpy.ndarray` + Real channel after rephasing, ``float64``, same shape as ``carr``. This + is the noise-floor-free channel retained for downstream processing. + imag : :obj:`numpy.ndarray` + Imaginary channel after rephasing, ``float64``, same shape as ``carr``. + After successful rephasing this should contain mostly Gaussian noise. + phase_bg : :obj:`numpy.ndarray` + Estimated smooth background phase in radians, ``float64``, same shape as + ``carr``. + + Notes + ----- + This implements the total-variation rephasing of + :footcite:t:`eichner2015real`. Diffusion-weighted images carry a smoothly + varying background phase (from eddy currents, B0 inhomogeneity, and motion) + on top of the desired real-valued signal. If that background phase is + removed, the true signal collapses onto the real channel and the noise + becomes Gaussian (rather than Rician), which avoids the magnitude noise + floor and enables true signal averaging. + + The background phase is estimated per in-plane slice. The magnitude and + wrapped phase are stacked and denoised jointly as a single image with + total-variation (Split-Bregman) regularization, and the denoised phase + channel is taken as the background phase ``phi_bg``. The data are then + rephased by ``carr * exp(-1j * phi_bg)``, rotating the smooth background + phase to zero. Because the phase is denoised directly, this estimator is + sensitive to phase wrapping; the wrap-free :func:`rephase_tv_complex` and + :func:`rephase_dc` variants avoid that limitation. + + References + ---------- + .. footbibliography:: + """ + carr = np.asarray(carr, dtype=np.complex128) + mag = np.abs(carr) + phase = np.angle(carr) + phase_bg = np.zeros(carr.shape, dtype=np.float64) + + # The magnitude/phase pair is denoised jointly as a single volumetric image + # (``channel_axis=None``), matching the reference ``imgrephaseTV.py`` (which + # relied on the old scikit-image ``multichannel=False`` default). Using + # ``channel_axis=-1`` instead performs *vectorial* 2-channel TV, which + # couples the (large dynamic range) magnitude edges into the phase estimate + # and imprints anatomy into the background phase, leaving coherent anatomical + # signal in the imaginary channel after rephasing. + if carr.ndim == 3: + for zz in range(carr.shape[2]): + buff = denoise_tv_bregman( + np.dstack((mag[:, :, zz], phase[:, :, zz])), + weight=2.0 * float(weight), + channel_axis=None, + ) + phase_bg[:, :, zz] = buff[:, :, 1] + elif carr.ndim == 4: + for vv in range(carr.shape[3]): + for zz in range(carr.shape[2]): + buff = denoise_tv_bregman( + np.dstack((mag[:, :, zz, vv], phase[:, :, zz, vv])), + weight=2.0 * float(weight), + channel_axis=None, + ) + phase_bg[:, :, zz, vv] = buff[:, :, 1] + else: + raise ValueError('rephase_tv expects a 3D or 4D array') + + rephased = carr * np.exp(-1j * phase_bg) + return np.real(rephased), np.imag(rephased), phase_bg + + +def rephase_tv_complex(carr, weight): + """Rephase complex MRI data using paper-faithful complex total variation. + + Parameters + ---------- + carr : :obj:`numpy.ndarray` of shape (X, Y, Z) or (X, Y, Z, V) + Complex-valued image, 3D (single volume) or 4D (a series of ``V`` + volumes). The array is internally cast to ``complex128``. For a 4D + input, each volume along the last axis is processed independently. + weight : :obj:`float` + Total-variation regularization weight (the ``lambda`` of + :footcite:t:`eichner2015real`). Larger values yield a smoother estimated + background phase. The weight passed to scikit-image's + :func:`~skimage.restoration.denoise_tv_bregman` is ``2.0 * weight``, + matching the convention used by :func:`rephase_tv`. + + Returns + ------- + real : :obj:`numpy.ndarray` + Real channel after rephasing, ``float64``, same shape as ``carr``. This + is the noise-floor-free channel retained for downstream processing. + imag : :obj:`numpy.ndarray` + Imaginary channel after rephasing, ``float64``, same shape as ``carr``. + After successful rephasing this should contain mostly Gaussian noise. + phase_bg : :obj:`numpy.ndarray` + Estimated smooth background phase in radians, ``float64``, same shape as + ``carr``. + + Notes + ----- + This is a paper-faithful variant of the total-variation rephasing of + :footcite:t:`eichner2015real`. Where :func:`rephase_tv` denoises a stacked + ``(magnitude, phase)`` image and reads off the phase channel, this function + estimates the background phase directly from a TV-denoised *complex* image, + ``phi_bg = angle(TV(complex signal))``. This corresponds to the complex + total-variation problem solved in the paper, + + .. math:: + + \\hat{x} = \\arg\\min_{x \\in \\mathbb{C}^n} + \\tfrac{1}{2}\\lVert y - x \\rVert_2^2 + + \\lambda \\lVert \\nabla x \\rVert_1, + + where ``y`` is the measured complex signal and ``angle(x_hat)`` is taken as + the background phase. Operating on the complex signal is wrap-free (no phase + unwrapping is required) and is intensity-weighted by construction, so it does + not imprint magnitude anatomy onto the phase estimate. + + The complex TV is realized per in-plane slice as a vectorial (2-channel) + total-variation denoising of the stacked real and imaginary parts. Before + denoising, each slice is normalized by a robust (95th-percentile) magnitude + so that the regularization ``weight`` is invariant to the large, + scanner-dependent magnitude scale; without this normalization scikit-image's + scale-dependent weight barely smooths DWI-scale data and ``phi_bg`` collapses + onto the noisy measured phase. The data are then rephased by + ``carr * exp(-1j * phi_bg)``. + + References + ---------- + .. footbibliography:: + """ + carr = np.asarray(carr, dtype=np.complex128) + phase_bg = np.zeros(carr.shape, dtype=np.float64) + sk_weight = 2.0 * float(weight) + + def _slice_phase_bg(cslice): + scale = np.percentile(np.abs(cslice), 95) + 1e-9 + cnorm = cslice / scale + denoised = denoise_tv_bregman( + np.dstack((cnorm.real, cnorm.imag)), + weight=sk_weight, + channel_axis=-1, # vectorial TV over (real, imag) == complex TV + ) + return np.angle(denoised[:, :, 0] + 1j * denoised[:, :, 1]) + + if carr.ndim == 3: + for zz in range(carr.shape[2]): + phase_bg[:, :, zz] = _slice_phase_bg(carr[:, :, zz]) + elif carr.ndim == 4: + for vv in range(carr.shape[3]): + for zz in range(carr.shape[2]): + phase_bg[:, :, zz, vv] = _slice_phase_bg(carr[:, :, zz, vv]) + else: + raise ValueError('rephase_tv_complex expects a 3D or 4D array') + + rephased = carr * np.exp(-1j * phase_bg) + return np.real(rephased), np.imag(rephased), phase_bg + + +def _convolve_per_slice(arr, kernel): + out = np.zeros_like(arr) + if arr.ndim == 3: + for zz in range(arr.shape[2]): + out[:, :, zz] = ndimage.convolve(arr[:, :, zz], kernel, mode='constant', cval=0.0) + elif arr.ndim == 4: + for vv in range(arr.shape[3]): + for zz in range(arr.shape[2]): + out[:, :, zz, vv] = ndimage.convolve( + arr[:, :, zz, vv], kernel, mode='constant', cval=0.0 + ) + else: + raise ValueError('expected a 3D or 4D array') + return out + + +def _dc_outlier_real(real, imag, ksize): + """Replace rephased-real outliers with the magnitude (Sprenger 2017). + + A voxel is an outlier when |magnitude - real| exceeds a local MAD-based + threshold computed on the imaginary channel within each in-plane slice. + """ + mag = np.sqrt(real * real + imag * imag) + delta = np.abs(mag - real) + thresh = np.zeros_like(real) + + def _slice_thresh(islice): + absdev = np.abs(islice - signal.medfilt(islice, ksize)) + medabsdev = signal.medfilt(absdev, ksize) + return 2.5 * 1.4826 * medabsdev + + if real.ndim == 3: + for zz in range(real.shape[2]): + thresh[:, :, zz] = _slice_thresh(imag[:, :, zz]) + elif real.ndim == 4: + for vv in range(real.shape[3]): + for zz in range(real.shape[2]): + thresh[:, :, zz, vv] = _slice_thresh(imag[:, :, zz, vv]) + else: + raise ValueError('expected a 3D or 4D array') + + real_thresh = mag.copy() # outliers fall back to magnitude + keep = delta < thresh + real_thresh[keep] = real[keep] + return real_thresh + + +def rephase_dc(carr, kernel_name, outlier_detection=False): + """Rephase complex MRI data using decorrelated phase filtering. + + Parameters + ---------- + carr : :obj:`numpy.ndarray` of shape (X, Y, Z) or (X, Y, Z, V) + Complex-valued image, 3D (single volume) or 4D (a series of ``V`` + volumes). The array is internally cast to ``complex128``. For a 4D + input, each volume along the last axis is processed independently. + kernel_name : :obj:`str` + Name of the 2D smoothing kernel used to estimate the background phase. + One of the keys of :data:`DC_KERNELS` (``'B3'``, ``'B5'``, ``'G3F1'``, + ``'G5F2'``, ``'G3F1H'``, ``'G5F2H'``, ``'Opt3'``, ``'Opt5'``). The + ``F``/``H`` and ``Opt`` kernels have a zeroed center, decoupling the + background-phase estimate at a voxel from that voxel's own value. + outlier_detection : :obj:`bool`, optional + If ``True``, voxels whose rephased real channel departs strongly from + the magnitude are flagged as outliers and replaced by the magnitude + (default: ``False``). See :func:`_dc_outlier_real`. + + Returns + ------- + real : :obj:`numpy.ndarray` + Real channel after rephasing, ``float64``, same shape as ``carr``. This + is the noise-floor-free channel retained for downstream processing. When + ``outlier_detection`` is ``True``, flagged voxels hold the magnitude + instead. + imag : :obj:`numpy.ndarray` + Imaginary channel after rephasing, ``float64``, same shape as ``carr``. + After successful rephasing this should contain mostly Gaussian noise. + phase_bg : :obj:`numpy.ndarray` + Estimated smooth background phase in radians, ``float64``, same shape as + ``carr``. + + Raises + ------ + ValueError + If ``kernel_name`` is not a key of :data:`DC_KERNELS`, or if ``carr`` is + neither 3D nor 4D. + + Notes + ----- + This implements the decorrelated-phase filtering of + :footcite:t:`sprenger2017real`. As with the total-variation methods, the goal + is to estimate and remove the smooth background phase so that the true signal + is recovered on the real channel with Gaussian (rather than Rician) noise, + avoiding the magnitude noise floor. + + Rather than total-variation denoising, the smooth background phase is + obtained by convolving the real and imaginary parts (separately, per + in-plane slice) with a small smoothing kernel and taking the phase of the + result, ``phi_bg = angle(conv(real) + 1j * conv(imag))``. Because real and + imaginary parts are smoothed independently, this is wrap-free. The + center-zeroed kernels (the ``H``- and ``Opt``-family) implement the + "decorrelation" of the paper: the phase estimate at each voxel is built only + from its neighbors, so the estimate is statistically decorrelated from the + voxel's own noise. The data are then rephased by ``carr * exp(-1j * + phi_bg)``. Optionally, an outlier step (:func:`_dc_outlier_real`) detects + voxels where rephasing failed and falls back to the magnitude there. + + References + ---------- + .. footbibliography:: + """ + carr = np.asarray(carr, dtype=np.complex128) + if kernel_name not in DC_KERNELS: + raise ValueError(f'Unsupported DC kernel: {kernel_name}') + kernel = DC_KERNELS[kernel_name] + + real_f = _convolve_per_slice(carr.real, kernel) + imag_f = _convolve_per_slice(carr.imag, kernel) + phase_bg = np.angle(real_f + 1j * imag_f) + + rephased = carr * np.exp(-1j * phase_bg) + real = np.real(rephased) + imag = np.imag(rephased) + + if outlier_detection: + real = _dc_outlier_real(real, imag, kernel.shape) + + return real, imag, phase_bg diff --git a/qsiprep/workflows/dwi/derivatives.py b/qsiprep/workflows/dwi/derivatives.py index 30e887de..d94869b2 100644 --- a/qsiprep/workflows/dwi/derivatives.py +++ b/qsiprep/workflows/dwi/derivatives.py @@ -168,6 +168,16 @@ def init_dwi_derivatives_wf(source_file) -> Workflow: mem_gb=DEFAULT_MEMORY_MIN_GB, ) + if config.workflow.dwi_phase_correction in ('tv', 'tvc', 'dc'): + for _node in ( + ds_dwi_t1, + ds_bvals_t1, + ds_bvecs_t1, + ds_gradient_table_t1, + ds_btable_t1, + ): + _node.inputs.part = 'real' + workflow.connect([ (inputnode, ds_dwi_t1, [('dwi_t1', 'in_file')]), (inputnode, ds_bvals_t1, [('bvals_t1', 'in_file')]), diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 3aae8e0c..5c860f87 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -29,6 +29,7 @@ PolarToComplex, ) from ...interfaces.nilearn import MaskEPI, Merge +from ...interfaces.phase import PhaseCorrect from ...interfaces.tortoise import Gibbs from ...utils.bids import IMPORTANT_DWI_FIELDS, update_metadata_from_nifti_header from .qc import init_modelfree_qc_wf @@ -282,6 +283,11 @@ def init_merge_and_denoise_wf( merge_confounds = pe.Node(niu.Merge(2), name='merge_confounds') hstack_confounds = pe.Node(StackConfounds(axis=1), name='hstack_confounds') n_volumes = dwi_df['NumVolumes'].sum() + if config.workflow.dwi_phase_correction in ('tv', 'tvc', 'dc'): + config.loggers.workflow.warning( + 'Phase correction is not applied to concatenated DWI series ' + '(denoising after combining); the magnitude will be used.' + ) denoising_wf = init_dwi_denoising_wf( partial_fourier=get_merged_parameter(dwi_df, 'PartialFourier', 'all'), phase_encoding_direction=get_merged_parameter(dwi_df, 'PhaseEncodingAxis', 'all'), @@ -462,8 +468,6 @@ def get_buffernode(): 'then denoised using the Marchenko-Pastur PCA method implemented in dwidenoise ' '[@mrtrix3; @dwidenoise1; @dwidenoise2; @cordero2019complex] ' f'with a{auto_str} window size of {dwi_denoise_window} voxels. ' - 'After denoising, the complex-valued data were split back into magnitude and ' - 'phase, and the denoised magnitude data were retained. ' ) last_step = 'After MP-PCA, ' @@ -501,16 +505,60 @@ def get_buffernode(): (denoiser, merge_confounds, [('nmse_text', f'in{step_num}')]), ]) # fmt:skip - split_complex = pe.Node( - ComplexToMagnitude(), - name='split_complex', - n_procs=omp_nthreads, - ) + phase_correction = config.workflow.dwi_phase_correction + if phase_correction in ('tv', 'tvc', 'dc'): + method_name = { + 'tv': 'total-variation rephasing [@eichner2015real]', + 'tvc': 'complex-valued total-variation rephasing [@eichner2015real]', + 'dc': 'decorrelated phase filtering [@sprenger2017real]', + }[phase_correction] + desc += ( + 'After denoising, the complex-valued data were phase-corrected ' + f'using {method_name}, and the real channel was retained for ' + 'subsequent processing. ' + ) - workflow.connect([ - (denoiser, split_complex, [('out_file', 'complex_file')]), - (split_complex, buffernodes[-1], [('out_file', 'dwi_file')]), - ]) # fmt:skip + phase_correct = pe.Node( + PhaseCorrect( + method=phase_correction, + tv_weight=config.workflow.dwi_phase_tv_weight, + dc_kernel=config.workflow.dwi_phase_dc_kernel, + ), + name='phase_correct', + n_procs=omp_nthreads, + ) + ds_report_phasecorr = pe.Node( + DerivativesDataSink( + datatype='figures', + desc='phasecorrection', + source_file=source_file, + ), + name=f'ds_report_{name}_phasecorrection', + run_without_submitting=True, + mem_gb=DEFAULT_MEMORY_MIN_GB, + ) + + workflow.connect([ + (denoiser, phase_correct, [('out_file', 'complex_file')]), + (phase_correct, buffernodes[-1], [('out_file', 'dwi_file')]), + (phase_correct, ds_report_phasecorr, [('out_report', 'in_file')]), + ]) # fmt:skip + else: + desc += ( + 'After denoising, the complex-valued data were split back into ' + 'magnitude and phase, and the denoised magnitude data were retained. ' + ) + + split_complex = pe.Node( + ComplexToMagnitude(), + name='split_complex', + n_procs=omp_nthreads, + ) + + workflow.connect([ + (denoiser, split_complex, [('out_file', 'complex_file')]), + (split_complex, buffernodes[-1], [('out_file', 'dwi_file')]), + ]) # fmt:skip elif denoise_method == 'dwidenoise': desc += (