From 99f8f28bf44d72b03f6faf8c1ac27a7ca9202340 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 13 Jul 2026 16:09:27 -0400 Subject: [PATCH 01/14] Draft DWIDenoise2 interface. --- Dockerfile | 44 ++++++++ qsiprep/interfaces/mrtrix.py | 191 +++++++++++++++++++++++++++++++++-- 2 files changed, 227 insertions(+), 8 deletions(-) diff --git a/Dockerfile b/Dockerfile index ca1cc36f8..b4a398879 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,35 @@ ARG BASE_IMAGE=pennlinc/qsiprep-base:20260415 +ARG DWIDENOISE2_COMMIT=892d5e8dd8f453ce1a561878f7dcb3998ae258ba +ARG MRTRIX3_DWIDENOISE2_COMMIT=fa6ee952913fbc1df79aeca745600852155f533f + +FROM buildpack-deps:bookworm AS dwidenoise2-build +ARG DWIDENOISE2_COMMIT +ARG MRTRIX3_DWIDENOISE2_COMMIT +RUN apt-get update && \ + apt-get install -y --no-install-recommends \ + cmake \ + libfftw3-dev \ + ninja-build \ + pkg-config \ + zlib1g-dev && \ + apt-get clean && rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* + +WORKDIR /src/dwidenoise2 +RUN git init . && \ + git remote add origin https://github.com/tsalo/dwidenoise2.git && \ + git fetch --depth 1 origin ${DWIDENOISE2_COMMIT} && \ + git checkout --detach FETCH_HEAD + +WORKDIR /src/mrtrix3 +RUN git clone --filter=blob:none --no-checkout https://github.com/MRtrix3/mrtrix3.git . && \ + git checkout --detach ${MRTRIX3_DWIDENOISE2_COMMIT} +RUN cp /src/dwidenoise2/cpp/cmd/dwidenoise2.cpp cpp/cmd/dwidenoise2.cpp && \ + cp -r /src/dwidenoise2/cpp/core/denoise cpp/core/denoise && \ + cmake -B build -GNinja \ + -DMRTRIX_BUILD_GUI=OFF \ + -DCMAKE_COMPILE_WARNING_AS_ERROR=ON \ + --preset=release && \ + cmake --build build --target dwidenoise2 FROM ghcr.io/prefix-dev/pixi:0.58.0 AS build RUN apt-get update && \ @@ -37,6 +68,19 @@ FROM ${BASE_IMAGE} AS base WORKDIR /home/qsiprep ENV HOME="/home/qsiprep" +COPY --from=dwidenoise2-build \ + /src/mrtrix3/build/bin/dwidenoise2 \ + /opt/dwidenoise2/bin/dwidenoise2 +COPY --from=dwidenoise2-build \ + /src/mrtrix3/build/cpp/core/libmrtrix-core.so \ + /opt/dwidenoise2/lib/libmrtrix-core.so +COPY --from=dwidenoise2-build \ + /src/dwidenoise2/LICENSE \ + /opt/dwidenoise2/LICENSE +ENV PATH="/opt/dwidenoise2/bin:$PATH" \ + LD_LIBRARY_PATH="/opt/dwidenoise2/lib:$LD_LIBRARY_PATH" +RUN dwidenoise2 -version + RUN chmod -R go=u $HOME WORKDIR /tmp diff --git a/qsiprep/interfaces/mrtrix.py b/qsiprep/interfaces/mrtrix.py index 1affb1f0a..9dc5ece02 100644 --- a/qsiprep/interfaces/mrtrix.py +++ b/qsiprep/interfaces/mrtrix.py @@ -115,19 +115,169 @@ def _run_interface(self, runtime): class DWIDenoiseInputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): in_file = File(exists=True, argstr='%s', position=-2, mandatory=True, desc='input DWI image') - mask = File(exists=True, argstr='-mask %s', position=1, desc='mask image') - extent = traits.Tuple( - (traits.Int, traits.Int, traits.Int), - argstr='-extent %d,%d,%d', - desc='set the window size of the denoising filter. (default = 5,5,5)', + mask = File(exists=True, desc='mask image used only to define the visual report contour') + onepass = traits.Bool(argstr='-onepass', desc='estimate noise and denoise in one pass') + datatype = traits.Enum( + 'float32', + 'float64', + argstr='-datatype %s', + desc='eigendecomposition datatype', + ) + decomposition = traits.Enum( + 'bdcsvd', + 'selfadjoint', + argstr='-decomposition %s', + desc='patch decomposition method', + ) + estimator = traits.Enum( + 'Exp1', + 'Exp2', + 'Med', + 'MRM2023', + argstr='-estimator %s', + desc='noise level estimator', + ) + noise_in = traits.Either( + traits.Float, + File(exists=True), + argstr='-noise_in %s', + xor=('fixed_rank',), + desc='scalar noise level or pre-estimated noise map', + ) + fixed_rank = traits.Int( + argstr='-fixed_rank %d', xor=('noise_in',), desc='fixed input signal rank' + ) + shape = traits.Enum( + 'cuboid', + 'sphere', + argstr='-shape %s', + desc='sliding spatial window shape', + ) + radius = traits.Float(argstr='-radius %g', desc='absolute spherical kernel radius in mm') + aspect_ratio = traits.Float( + argstr='-aspect_ratio %g', + desc='ratio of kernel voxels to input volumes', + ) + minvoxels = traits.Int(argstr='-minvoxels %d', desc='minimum voxels in a spherical kernel') + extent = traits.Either( + traits.Int, + traits.Tuple(traits.Int, traits.Int, traits.Int), + argstr='-extent %s', + desc='cuboid window size as one integer or a triplet', + ) + subsample = traits.Either( + traits.Int, + traits.Tuple(traits.Int, traits.Int, traits.Int), + argstr='-subsample %s', + desc='PCA kernel subsampling factor as one integer or a triplet', + ) + demodulate = traits.Enum( + 'none', + 'linear', + 'nonlinear', + argstr='-demodulate %s', + desc='phase demodulation mode', + ) + demod_axes = traits.Str( + argstr='-demod_axes %s', + desc='comma-separated FFT axes for phase demodulation', + ) + demean = traits.Enum( + 'none', + 'volume_groups', + 'shells', + 'all', + argstr='-demean %s', + desc='demeaning method before PCA', + ) + vst = File( + exists=True, + argstr='-vst %s', + desc='noise map for variance-stabilising transformation', + ) + preconditioned_input = File( + argstr='-preconditioned_input %s', + desc='export preconditioned PCA input', + ) + preconditioned_output = File( + argstr='-preconditioned_output %s', + desc='export output before reversing preconditioning', + ) + filter_method = traits.Enum( + 'optshrink', + 'optthresh', + 'truncate', + argstr='-filter %s', + desc='eigenvalue filtering method', + ) + aggregator = traits.Enum( + 'exclusive', + 'gaussian', + 'invl0', + 'rank', + 'uniform', + argstr='-aggregator %s', + desc='overlapping-patch aggregation method', ) noise_image = File( - argstr='-noise %s', + argstr='-noise_out %s', name_template='%s_noise.nii.gz', name_source=['in_file'], keep_extension=False, desc='the output noise map', ) + lamplus = File(argstr='-lamplus %s', desc='estimated upper noise eigenspectrum bound') + rank_pcanonzero = File(argstr='-rank_pcanonzero %s', desc='non-zero PCA rank before denoising') + rank_input = File(argstr='-rank_input %s', desc='estimated input rank per denoising patch') + rank_output = File( + argstr='-rank_output %s', + desc='estimated output rank after patch aggregation', + ) + variance_removed = File( + argstr='-variance_removed %s', + desc='fraction of variance removed by PCA', + ) + eigenspectra = File( + argstr='-eigenspectra %s', + desc='matrix of eigenvalue spectra across patches', + ) + residual_statistics = traits.Tuple( + File(), + File(), + File(), + argstr='-residual_statistics %s %s %s', + desc='residual mean, variance, and maximum-absolute-value images', + ) + max_dist = File(argstr='-max_dist %s', desc='maximum within-patch voxel distance') + voxelcount = File(argstr='-voxelcount %s', desc='voxels contributing to each PCA') + patchcount = File(argstr='-patchcount %s', desc='unique patches containing each voxel') + sum_aggregation = File( + argstr='-sum_aggregation %s', + desc='sum of aggregation weights per voxel', + ) + sum_optshrink = File( + argstr='-sum_optshrink %s', + desc='sum of optimal-shrinkage weights per patch', + ) + grad_file = File( + exists=True, + argstr='-grad %s', + xor=('bvec_file', 'bval_file'), + desc='MRtrix-format diffusion gradient scheme', + ) + bvec_file = File( + exists=True, + argstr='-fslgrad %s %s', + requires=('bval_file',), + xor=('grad_file',), + desc='FSL-format diffusion gradient b-vector file', + ) + bval_file = File( + exists=True, + requires=('bvec_file',), + xor=('grad_file',), + desc='FSL-format diffusion gradient b-value file', + ) out_file = File( name_template='%s_denoised.nii.gz', name_source=['in_file'], @@ -137,13 +287,31 @@ class DWIDenoiseInputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): desc='the output denoised DWI image', ) out_report = File( - 'dwidenoise_report.svg', usedefault=True, desc='filename for the visual report' + 'dwidenoise_report.svg', + usedefault=True, + desc='filename for the visual report', ) class DWIDenoiseOutputSpec(SeriesPreprocReportOutputSpec): noise_image = File(desc='the output noise map', exists=True) out_file = File(desc='the output denoised DWI image', exists=True) + preconditioned_input = File(exists=True, desc='preconditioned PCA input') + preconditioned_output = File(exists=True, desc='output before reversal of preconditioning') + lamplus = File(exists=True, desc='estimated upper noise eigenspectrum bound') + rank_pcanonzero = File(exists=True, desc='non-zero PCA rank before denoising') + rank_input = File(exists=True, desc='estimated input rank per denoising patch') + rank_output = File(exists=True, desc='estimated output rank after patch aggregation') + variance_removed = File(exists=True, desc='fraction of variance removed by PCA') + eigenspectra = File(exists=True, desc='matrix of eigenvalue spectra across patches') + residual_statistics = traits.Tuple( + File(exists=True), File(exists=True), File(exists=True), desc='residual statistic images' + ) + max_dist = File(exists=True, desc='maximum within-patch voxel distance') + voxelcount = File(exists=True, desc='voxels contributing to each PCA') + patchcount = File(exists=True, desc='unique patches containing each voxel') + sum_aggregation = File(exists=True, desc='sum of aggregation weights per voxel') + sum_optshrink = File(exists=True, desc='sum of optimal-shrinkage weights per patch') class DWIDenoise(SeriesPreprocReport, MRTrix3Base): @@ -167,10 +335,17 @@ class DWIDenoise(SeriesPreprocReport, MRTrix3Base): """ - _cmd = 'dwidenoise' + _cmd = 'dwidenoise2' input_spec = DWIDenoiseInputSpec output_spec = DWIDenoiseOutputSpec + def _format_arg(self, name, spec, value): + if name in ('extent', 'subsample') and not isinstance(value, int): + value = ','.join(str(item) for item in value) + elif name == 'bvec_file': + value = (value, self.inputs.bval_file) + return super()._format_arg(name, spec, value) + def _get_plotting_images(self): input_dwi = load_img(self.inputs.in_file) outputs = self._list_outputs() From fb4ae67227b2ed0075bccb17bd8d2fcbdf03cdd6 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 14 Jul 2026 10:07:10 -0400 Subject: [PATCH 02/14] Work on dwidenoise parser. --- qsiprep/cli/parser.py | 22 +++-- qsiprep/config.py | 3 +- qsiprep/interfaces/mrtrix.py | 17 +++- qsiprep/tests/test_interfaces_mrtrix3.py | 61 +++++++++++++ qsiprep/tests/test_utils_misc.py | 70 ++++++++++++++- qsiprep/utils/misc.py | 109 +++++++++++++++++++++++ qsiprep/workflows/dwi/merge.py | 25 ++++-- 7 files changed, 289 insertions(+), 18 deletions(-) diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index acd640850..354061b81 100644 --- a/qsiprep/cli/parser.py +++ b/qsiprep/cli/parser.py @@ -27,6 +27,7 @@ import sys from .. import config +from ..utils.misc import parse_denoise_method def _build_parser(**kwargs): @@ -121,6 +122,13 @@ def _int_or_auto(value, parser): return value + def _denoise_method(value, parser): + try: + parse_denoise_method(value) + except ValueError as exc: + parser.error(f'Invalid --denoise-method specification: {exc}') + return value + def _to_gb(value): scale = {'G': 1, 'T': 10**3, 'M': 1e-3, 'K': 1e-6, 'B': 1e-9} digits = ''.join([c for c in value if c.isdigit()]) @@ -177,6 +185,7 @@ def _bids_filter(value, parser): IsFile = partial(_is_file, parser=parser) PositiveInt = partial(_min_one, parser=parser) IntOrAuto = partial(_int_or_auto, parser=parser) + DenoiseMethod = partial(_denoise_method, parser=parser) BIDSFilter = partial(_bids_filter, parser=parser) # Arguments as specified by BIDS-Apps @@ -394,10 +403,12 @@ def _bids_filter(value, parser): g_conf.add_argument( '--denoise-method', action='store', - choices=['dwidenoise', 'patch2self', 'none'], + type=DenoiseMethod, default='dwidenoise', - help='Image-based denoising method. Either "dwidenoise" (MRtrix), ' - '"patch2self" (DIPY) or "none". (default: dwidenoise)', + help='Image-based denoising method: "dwidenoise" (MRtrix), "patch2self" (DIPY), ' + 'or "none". DWIDenoise parameters may follow the method as semicolon-delimited ' + 'name:value pairs, for example ' + '"dwidenoise;demodulate:nonlinear;decomposition:bdcsvd".', ) g_conf.add_argument( '--unringing-method', @@ -749,12 +760,13 @@ def parse_args(args=None, namespace=None): ) # Validate the tricky options here + denoise_method, _ = parse_denoise_method(config.workflow.denoise_method) if config.workflow.dwi_denoise_window != 'auto': - if config.workflow.denoise_method == 'patch2self': + if denoise_method == 'patch2self': config.loggers.cli.error( 'The --dwi-denoise-window option is not used when --denoise-method=patch2self' ) - elif config.workflow.denoise_method == 'none': + elif denoise_method == 'none': config.loggers.cli.warning( 'The --dwi-denoise-window option is not used when --denoise-method=none' ) diff --git a/qsiprep/config.py b/qsiprep/config.py index bac7ced4c..0d1338719 100644 --- a/qsiprep/config.py +++ b/qsiprep/config.py @@ -568,7 +568,8 @@ class workflow(_Config): """Run ``dwidenoise`` after combining dwis, but before motion correction.""" denoise_method = None """Image-based denoising method. Either "dwidenoise" (MRtrix), "patch2self" (DIPY) - or "none".""" + or "none". DWIDenoise parameters may be appended as semicolon-delimited name:value + pairs.""" distortion_group_merge = None """How to combine images across distortion groups (concatenate, average or none).""" dwi_denoise_window = None diff --git a/qsiprep/interfaces/mrtrix.py b/qsiprep/interfaces/mrtrix.py index 9dc5ece02..c0ff69fbe 100644 --- a/qsiprep/interfaces/mrtrix.py +++ b/qsiprep/interfaces/mrtrix.py @@ -148,12 +148,16 @@ class DWIDenoiseInputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): argstr='-fixed_rank %d', xor=('noise_in',), desc='fixed input signal rank' ) shape = traits.Enum( - 'cuboid', 'sphere', + 'cuboid', argstr='-shape %s', desc='sliding spatial window shape', ) - radius = traits.Float(argstr='-radius %g', desc='absolute spherical kernel radius in mm') + radius = traits.Float( + argstr='-radius %g', + xor=('extent',), + desc='absolute spherical kernel radius in mm', + ) aspect_ratio = traits.Float( argstr='-aspect_ratio %g', desc='ratio of kernel voxels to input volumes', @@ -163,6 +167,7 @@ class DWIDenoiseInputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): traits.Int, traits.Tuple(traits.Int, traits.Int, traits.Int), argstr='-extent %s', + xor=('radius',), desc='cuboid window size as one integer or a triplet', ) subsample = traits.Either( @@ -346,6 +351,14 @@ def _format_arg(self, name, spec, value): value = (value, self.inputs.bval_file) return super()._format_arg(name, spec, value) + def _parse_inputs(self, skip=None): + shape = self.inputs.shape if isdefined(self.inputs.shape) else 'sphere' + if shape == 'sphere' and isdefined(self.inputs.extent): + raise ValueError("'extent' cannot be used when 'shape' is 'sphere'") + if shape == 'cuboid' and isdefined(self.inputs.radius): + raise ValueError("'radius' cannot be used when 'shape' is 'cuboid'") + return super()._parse_inputs(skip=skip) + def _get_plotting_images(self): input_dwi = load_img(self.inputs.in_file) outputs = self._list_outputs() diff --git a/qsiprep/tests/test_interfaces_mrtrix3.py b/qsiprep/tests/test_interfaces_mrtrix3.py index 84e4de55c..695a7882b 100644 --- a/qsiprep/tests/test_interfaces_mrtrix3.py +++ b/qsiprep/tests/test_interfaces_mrtrix3.py @@ -3,8 +3,11 @@ import os import nibabel as nb +import pytest +from qsiprep import config from qsiprep.interfaces import mrtrix +from qsiprep.workflows.dwi.merge import init_dwi_denoising_wf def test_dwidenoise(datasets, tmp_path_factory): @@ -16,6 +19,7 @@ def test_dwidenoise(datasets, tmp_path_factory): in_img = nb.load(in_file) interface = mrtrix.DWIDenoise( + shape='cuboid', extent=(5, 5, 5), in_file=in_file, nthreads=1, @@ -33,3 +37,60 @@ def test_dwidenoise(datasets, tmp_path_factory): assert os.path.isfile(results.outputs.out_report) assert os.path.isfile(results.outputs.nmse_text) + + +@pytest.mark.parametrize( + ('shape', 'kernel_option', 'error'), + [ + ('sphere', {'extent': (5, 5, 5)}, "'extent' cannot be used"), + ('cuboid', {'radius': 2.5}, "'radius' cannot be used"), + ], +) +def test_dwidenoise_kernel_shape_validation(tmp_path, shape, kernel_option, error): + """Reject kernel options that do not apply to the selected shape.""" + in_file = tmp_path / 'dwi.nii.gz' + in_file.touch() + interface = mrtrix.DWIDenoise(in_file=in_file, shape=shape, **kernel_option) + + with pytest.raises(ValueError, match=error): + _ = interface.cmdline + + +def test_dwidenoise_kernel_options_are_mutually_exclusive(tmp_path): + """Reject simultaneous spherical and cuboid kernel size options.""" + in_file = tmp_path / 'dwi.nii.gz' + in_file.touch() + + with pytest.raises(OSError, match='mutually exclusive'): + mrtrix.DWIDenoise( + in_file=in_file, + shape='sphere', + radius=2.5, + extent=(5, 5, 5), + ) + + +def test_dwidenoise_cli_parameters_reach_workflow(monkeypatch): + """Forward parsed DWIDenoise parameters to the workflow node.""" + monkeypatch.setattr( + config.workflow, + 'denoise_method', + 'dwidenoise;demodulate:nonlinear;decomposition:bdcsvd', + ) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=False, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + assert denoiser.inputs.demodulate == 'nonlinear' + assert denoiser.inputs.decomposition == 'bdcsvd' diff --git a/qsiprep/tests/test_utils_misc.py b/qsiprep/tests/test_utils_misc.py index c22c25479..a0a03f745 100644 --- a/qsiprep/tests/test_utils_misc.py +++ b/qsiprep/tests/test_utils_misc.py @@ -3,8 +3,10 @@ import logging import numpy as np +import pytest -from qsiprep.utils.misc import safe_unit_vector +from qsiprep.cli.parser import _build_parser +from qsiprep.utils.misc import parse_denoise_method, safe_unit_vector def test_safe_unit_vector_zero_magnitude_substitutes_x_axis(): @@ -45,3 +47,69 @@ def test_angle_between_finite_for_zero_vector(): angle = angle_between(np.array([0.0, 0.0, 0.0]), np.array([1.0, 0.0, 0.0])) assert np.isfinite(angle) + + +def test_parse_denoise_method_parameters(): + method, parameters = parse_denoise_method( + 'dwidenoise;demodulate:nonlinear;decomposition:bdcsvd;' + 'onepass:true;radius:2.5;subsample:2,2,2' + ) + + assert method == 'dwidenoise' + assert parameters == { + 'demodulate': 'nonlinear', + 'decomposition': 'bdcsvd', + 'onepass': True, + 'radius': 2.5, + 'subsample': (2, 2, 2), + } + + +@pytest.mark.parametrize( + 'spec', + [ + 'unknown', + 'patch2self;decomposition:bdcsvd', + 'dwidenoise;decomposition', + 'dwidenoise;unknown:value', + 'dwidenoise;decomposition:bdcsvd;decomposition:selfadjoint', + 'dwidenoise;decomposition:invalid', + 'dwidenoise;onepass:maybe', + 'dwidenoise;extent:1,2', + ], +) +def test_parse_denoise_method_rejects_invalid_specs(spec): + with pytest.raises(ValueError, match='.'): + parse_denoise_method(spec) + + +def test_denoise_method_cli_parameter(tmp_path): + spec = 'dwidenoise;demodulate:nonlinear;decomposition:bdcsvd' + opts = _build_parser().parse_args( + [ + str(tmp_path), + str(tmp_path / 'out'), + 'participant', + '--output-resolution', + '2', + '--denoise-method', + spec, + ] + ) + + assert opts.denoise_method == spec + + +def test_denoise_method_cli_rejects_invalid_parameter(tmp_path): + with pytest.raises(SystemExit): + _build_parser().parse_args( + [ + str(tmp_path), + str(tmp_path / 'out'), + 'participant', + '--output-resolution', + '2', + '--denoise-method', + 'dwidenoise;decomposition:invalid', + ] + ) diff --git a/qsiprep/utils/misc.py b/qsiprep/utils/misc.py index b64d80633..485536ca8 100644 --- a/qsiprep/utils/misc.py +++ b/qsiprep/utils/misc.py @@ -8,6 +8,115 @@ LOGGER = logging.getLogger('nipype.interface') +_DWIDENOISE_ENUM_PARAMETERS = { + 'datatype': ('float32', 'float64'), + 'decomposition': ('bdcsvd', 'selfadjoint'), + 'estimator': ('Exp1', 'Exp2', 'Med', 'MRM2023'), + 'shape': ('sphere', 'cuboid'), + 'demodulate': ('none', 'linear', 'nonlinear'), + 'demean': ('none', 'volume_groups', 'shells', 'all'), + 'filter_method': ('optshrink', 'optthresh', 'truncate'), + 'aggregator': ('exclusive', 'gaussian', 'invl0', 'rank', 'uniform'), +} +_DWIDENOISE_STRING_PARAMETERS = { + 'demod_axes', + 'vst', + 'preconditioned_input', + 'preconditioned_output', + 'noise_image', + 'lamplus', + 'rank_pcanonzero', + 'rank_input', + 'rank_output', + 'variance_removed', + 'eigenspectra', + 'max_dist', + 'voxelcount', + 'patchcount', + 'sum_aggregation', + 'sum_optshrink', + 'grad_file', + 'bvec_file', + 'bval_file', +} +_DWIDENOISE_PARAMETERS = ( + set(_DWIDENOISE_ENUM_PARAMETERS) + | _DWIDENOISE_STRING_PARAMETERS + | { + 'onepass', + 'noise_in', + 'fixed_rank', + 'radius', + 'aspect_ratio', + 'minvoxels', + 'extent', + 'subsample', + 'residual_statistics', + } +) + + +def parse_denoise_method(spec): + """Parse a denoising method and semicolon-delimited parameters. + + Parameters use ``name:value`` syntax, for example + ``dwidenoise;demodulate:nonlinear;decomposition:bdcsvd``. + """ + elements = spec.split(';') + method = elements[0].strip() + if method not in ('dwidenoise', 'patch2self', 'none'): + raise ValueError(f'Unknown denoising method: {method!r}') + if len(elements) > 1 and method != 'dwidenoise': + raise ValueError(f'{method!r} does not accept DWIDenoise parameters') + + parameters = {} + for element in elements[1:]: + name, separator, value = element.partition(':') + name = name.strip() + value = value.strip() + if not separator or not name or not value: + raise ValueError(f'Invalid DWIDenoise parameter: {element!r}') + if name not in _DWIDENOISE_PARAMETERS: + raise ValueError(f'Unknown DWIDenoise parameter: {name!r}') + if name in parameters: + raise ValueError(f'Duplicate DWIDenoise parameter: {name!r}') + + if name in _DWIDENOISE_ENUM_PARAMETERS: + choices = _DWIDENOISE_ENUM_PARAMETERS[name] + if value not in choices: + raise ValueError(f'Invalid value for {name!r}: {value!r}; choose from {choices}') + parsed_value = value + elif name == 'onepass': + bool_values = {'true': True, 'false': False, '1': True, '0': False} + try: + parsed_value = bool_values[value.lower()] + except KeyError as exc: + raise ValueError(f'Invalid boolean value for {name!r}: {value!r}') from exc + elif name in ('fixed_rank', 'minvoxels'): + parsed_value = int(value) + elif name in ('radius', 'aspect_ratio'): + parsed_value = float(value) + elif name == 'noise_in': + try: + parsed_value = float(value) + except ValueError: + parsed_value = value + elif name in ('extent', 'subsample'): + values = tuple(int(item.strip()) for item in value.split(',')) + if len(values) not in (1, 3): + raise ValueError(f'{name!r} must contain one or three integers') + parsed_value = values[0] if len(values) == 1 else values + elif name == 'residual_statistics': + parsed_value = tuple(item.strip() for item in value.split(',')) + if len(parsed_value) != 3 or not all(parsed_value): + raise ValueError(f'{name!r} must contain three file names') + else: + parsed_value = value + + parameters[name] = parsed_value + + return method, parameters + def safe_unit_vector(vector): """Return the unit vector of ``vector``. diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 60c5384a0..f26322c78 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -31,6 +31,7 @@ from ...interfaces.nilearn import MaskEPI, Merge from ...interfaces.tortoise import Gibbs from ...utils.bids import IMPORTANT_DWI_FIELDS, update_metadata_from_nifti_header +from ...utils.misc import parse_denoise_method from .qc import init_modelfree_qc_wf from .util import _get_wf_name @@ -415,7 +416,7 @@ def get_buffernode(): ]) # fmt:skip # Which steps to apply? - denoise_method = config.workflow.denoise_method + denoise_method, dwidenoise_params = parse_denoise_method(config.workflow.denoise_method) unringing_method = config.workflow.unringing_method do_denoise = denoise_method in ('patch2self', 'dwidenoise') do_unringing = config.workflow.unringing_method in ('mrdegibbs', 'rpg') @@ -487,11 +488,10 @@ def get_buffernode(): (phase_to_radians, combine_complex, [('phase_file', 'phase_file')]), ]) # fmt:skip + dwidenoise_inputs = {'shape': 'sphere', 'nthreads': omp_nthreads} + dwidenoise_inputs.update(dwidenoise_params) denoiser = pe.Node( - DWIDenoise( - extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), - nthreads=omp_nthreads, - ), + DWIDenoise(**dwidenoise_inputs), name='denoiser', n_procs=omp_nthreads, ) @@ -522,11 +522,16 @@ def get_buffernode(): ) last_step = 'After MP-PCA, ' + dwidenoise_inputs = { + 'shape': 'cuboid', + 'extent': (dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), + 'nthreads': omp_nthreads, + } + if dwidenoise_params.get('shape') == 'sphere' and 'extent' not in dwidenoise_params: + dwidenoise_inputs.pop('extent') + dwidenoise_inputs.update(dwidenoise_params) denoiser = pe.Node( - DWIDenoise( - extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), - nthreads=omp_nthreads, - ), + DWIDenoise(**dwidenoise_inputs), name='denoiser', n_procs=omp_nthreads, ) @@ -541,6 +546,8 @@ def get_buffernode(): name='denoiser', n_procs=omp_nthreads, ) + + if denoise_method in ('dwidenoise', 'patch2self'): workflow.connect([(inputnode, denoiser, [('bval_file', 'bval_file')])]) if (denoise_method in ('dwidenoise', 'patch2self')) and not use_phase: From 7bd8580a503f4ac95b05f2824f6bea0ef228e4fb Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 14 Jul 2026 10:54:35 -0400 Subject: [PATCH 03/14] Try fixing. --- qsiprep/cli/parser.py | 14 ++++++++++---- qsiprep/tests/test_cli.py | 4 +++- qsiprep/tests/test_interfaces_mrtrix3.py | 4 ++++ qsiprep/workflows/dwi/merge.py | 10 ++++++++-- 4 files changed, 25 insertions(+), 7 deletions(-) diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index 354061b81..c8da2e801 100644 --- a/qsiprep/cli/parser.py +++ b/qsiprep/cli/parser.py @@ -405,10 +405,16 @@ def _bids_filter(value, parser): action='store', type=DenoiseMethod, default='dwidenoise', - help='Image-based denoising method: "dwidenoise" (MRtrix), "patch2self" (DIPY), ' - 'or "none". DWIDenoise parameters may follow the method as semicolon-delimited ' - 'name:value pairs, for example ' - '"dwidenoise;demodulate:nonlinear;decomposition:bdcsvd".', + help=( + 'Image-based denoising method: "dwidenoise" (MRtrix), "patch2self" (DIPY), ' + 'or "none".\n' + 'dwidenoise parameters may follow the method as semicolon-delimited ' + 'name:value pairs, for example ' + '"dwidenoise;demodulate:nonlinear;decomposition:bdcsvd".\n' + 'To approximate legacy "dwidenoise", use ' + '"dwidenoise;shape:cuboid;subsample:1;demodulate:none;demean:none;' + 'filter_method:truncate;aggregator:exclusive".' + ), ) g_conf.add_argument( '--unringing-method', diff --git a/qsiprep/tests/test_cli.py b/qsiprep/tests/test_cli.py index d68386916..9fa7503f3 100644 --- a/qsiprep/tests/test_cli.py +++ b/qsiprep/tests/test_cli.py @@ -60,6 +60,8 @@ def test_dsdti_fmap(data_dir, output_dir, working_dir): '--write-graph', '--mem-mb=4096', '--output-resolution=5', + '--denoise-method=dwidenoise;shape:cuboid;subsample:1;demodulate:none;demean:none;' + 'filter_method:truncate;aggregator:exclusive', ] _run_and_generate(TEST_NAME, parameters, test_main=False) @@ -566,7 +568,7 @@ def test_maternal_brain_project(data_dir, output_dir, working_dir): @pytest.mark.integration @pytest.mark.forrest_gump def test_forrest_gump(data_dir, output_dir, working_dir): - """Run QSIPrep on Forrest Gump data with dwidenoise denoising. + """Run QSIPrep on Forrest Gump data without dwidenoise denoising. The dataset was built from the Forrest Gump dataset: https://openneuro.org/datasets/ds000113/versions/1.3.0 diff --git a/qsiprep/tests/test_interfaces_mrtrix3.py b/qsiprep/tests/test_interfaces_mrtrix3.py index 695a7882b..90bf69fd3 100644 --- a/qsiprep/tests/test_interfaces_mrtrix3.py +++ b/qsiprep/tests/test_interfaces_mrtrix3.py @@ -21,6 +21,8 @@ def test_dwidenoise(datasets, tmp_path_factory): interface = mrtrix.DWIDenoise( shape='cuboid', extent=(5, 5, 5), + onepass=True, + subsample=1, in_file=in_file, nthreads=1, ) @@ -94,3 +96,5 @@ def test_dwidenoise_cli_parameters_reach_workflow(monkeypatch): assert denoiser.inputs.demodulate == 'nonlinear' assert denoiser.inputs.decomposition == 'bdcsvd' + assert denoiser.inputs.onepass is True + assert denoiser.inputs.subsample == 1 diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index f26322c78..fa01b62f8 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -525,10 +525,16 @@ def get_buffernode(): dwidenoise_inputs = { 'shape': 'cuboid', 'extent': (dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), + # Fixed odd extents are incompatible with the changing subsampling parity of + # iterative mode, so reproduce legacy fixed-window behavior in one pass. + 'onepass': True, + 'subsample': 1, 'nthreads': omp_nthreads, } - if dwidenoise_params.get('shape') == 'sphere' and 'extent' not in dwidenoise_params: - dwidenoise_inputs.pop('extent') + if dwidenoise_params.get('shape') == 'sphere': + for parameter in ('extent', 'onepass', 'subsample'): + if parameter not in dwidenoise_params: + dwidenoise_inputs.pop(parameter) dwidenoise_inputs.update(dwidenoise_params) denoiser = pe.Node( DWIDenoise(**dwidenoise_inputs), From 62af2d3abc4cc59ff767afe659a1d1ea7b5c2c80 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 14 Jul 2026 11:15:14 -0400 Subject: [PATCH 04/14] Add regular dwidenoise back. --- qsiprep/cli/parser.py | 14 ++-- qsiprep/interfaces/mrtrix.py | 76 ++++++++++++++++++++-- qsiprep/tests/test_cli.py | 6 +- qsiprep/tests/test_interfaces_mrtrix3.py | 42 ++++++++++-- qsiprep/utils/misc.py | 16 ++--- qsiprep/workflows/dwi/merge.py | 83 ++++++++++++++---------- 6 files changed, 175 insertions(+), 62 deletions(-) diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index c8da2e801..d53272b5b 100644 --- a/qsiprep/cli/parser.py +++ b/qsiprep/cli/parser.py @@ -406,13 +406,13 @@ def _bids_filter(value, parser): type=DenoiseMethod, default='dwidenoise', help=( - 'Image-based denoising method: "dwidenoise" (MRtrix), "patch2self" (DIPY), ' - 'or "none".\n' - 'dwidenoise parameters may follow the method as semicolon-delimited ' + 'Image-based denoising method: "dwidenoise" (MRtrix), "dwidenoise2", ' + '"patch2self" (DIPY), or "none".\n' + 'dwidenoise2 parameters may follow the method as semicolon-delimited ' 'name:value pairs, for example ' - '"dwidenoise;demodulate:nonlinear;decomposition:bdcsvd".\n' - 'To approximate legacy "dwidenoise", use ' - '"dwidenoise;shape:cuboid;subsample:1;demodulate:none;demean:none;' + '"dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd".\n' + 'To approximate legacy "dwidenoise" with "dwidenoise2", use ' + '"dwidenoise2;shape:cuboid;subsample:1;demodulate:none;demean:none;' 'filter_method:truncate;aggregator:exclusive".' ), ) @@ -452,7 +452,7 @@ def _bids_filter(value, parser): g_conf.add_argument( '--denoise-after-combining', action='store_true', - help='run ``dwidenoise`` after combining dwis, but before motion correction', + help='run denoising after combining dwis, but before motion correction', ) g_conf.add_argument( '--separate-all-dwis', diff --git a/qsiprep/interfaces/mrtrix.py b/qsiprep/interfaces/mrtrix.py index c0ff69fbe..26927d1ec 100644 --- a/qsiprep/interfaces/mrtrix.py +++ b/qsiprep/interfaces/mrtrix.py @@ -114,6 +114,74 @@ def _run_interface(self, runtime): class DWIDenoiseInputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): + in_file = File(exists=True, argstr='%s', position=-2, mandatory=True, desc='input DWI image') + mask = File(exists=True, argstr='-mask %s', position=1, desc='mask image') + extent = traits.Tuple( + (traits.Int, traits.Int, traits.Int), + argstr='-extent %d,%d,%d', + desc='set the window size of the denoising filter. (default = 5,5,5)', + ) + noise_image = File( + argstr='-noise %s', + name_template='%s_noise.nii.gz', + name_source=['in_file'], + keep_extension=False, + desc='the output noise map', + ) + out_file = File( + name_template='%s_denoised.nii.gz', + name_source=['in_file'], + keep_extension=False, + argstr='%s', + position=-1, + desc='the output denoised DWI image', + ) + out_report = File( + 'dwidenoise_report.svg', usedefault=True, desc='filename for the visual report' + ) + + +class DWIDenoiseOutputSpec(SeriesPreprocReportOutputSpec): + noise_image = File(desc='the output noise map', exists=True) + out_file = File(desc='the output denoised DWI image', exists=True) + + +class DWIDenoise(SeriesPreprocReport, MRTrix3Base): + """ + Denoise DWI data and estimate the noise level based on the optimal + threshold for PCA. + + DWI data denoising and noise map estimation by exploiting data redundancy + in the PCA domain using the prior knowledge that the eigenspectrum of + random covariance matrices is described by the universal Marchenko Pastur + distribution. + + Important note: image denoising must be performed as the first step of the + image processing pipeline. The routine will fail if interpolation or + smoothing has been applied to the data prior to denoising. + + Note that this function does not correct for non-Gaussian noise biases. + + For more information, see + + + """ + + _cmd = 'dwidenoise' + input_spec = DWIDenoiseInputSpec + output_spec = DWIDenoiseOutputSpec + + def _get_plotting_images(self): + input_dwi = load_img(self.inputs.in_file) + outputs = self._list_outputs() + ref_name = outputs.get('out_file') + denoised_nii = load_img(ref_name) + noise_name = outputs['noise_image'] + noisenii = load_img(noise_name) + return input_dwi, denoised_nii, noisenii + + +class DWIDenoise2InputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): in_file = File(exists=True, argstr='%s', position=-2, mandatory=True, desc='input DWI image') mask = File(exists=True, desc='mask image used only to define the visual report contour') onepass = traits.Bool(argstr='-onepass', desc='estimate noise and denoise in one pass') @@ -298,7 +366,7 @@ class DWIDenoiseInputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): ) -class DWIDenoiseOutputSpec(SeriesPreprocReportOutputSpec): +class DWIDenoise2OutputSpec(SeriesPreprocReportOutputSpec): noise_image = File(desc='the output noise map', exists=True) out_file = File(desc='the output denoised DWI image', exists=True) preconditioned_input = File(exists=True, desc='preconditioned PCA input') @@ -319,7 +387,7 @@ class DWIDenoiseOutputSpec(SeriesPreprocReportOutputSpec): sum_optshrink = File(exists=True, desc='sum of optimal-shrinkage weights per patch') -class DWIDenoise(SeriesPreprocReport, MRTrix3Base): +class DWIDenoise2(SeriesPreprocReport, MRTrix3Base): """ Denoise DWI data and estimate the noise level based on the optimal threshold for PCA. @@ -341,8 +409,8 @@ class DWIDenoise(SeriesPreprocReport, MRTrix3Base): """ _cmd = 'dwidenoise2' - input_spec = DWIDenoiseInputSpec - output_spec = DWIDenoiseOutputSpec + input_spec = DWIDenoise2InputSpec + output_spec = DWIDenoise2OutputSpec def _format_arg(self, name, spec, value): if name in ('extent', 'subsample') and not isinstance(value, int): diff --git a/qsiprep/tests/test_cli.py b/qsiprep/tests/test_cli.py index 9fa7503f3..9a209b530 100644 --- a/qsiprep/tests/test_cli.py +++ b/qsiprep/tests/test_cli.py @@ -38,7 +38,7 @@ def test_dsdti_fmap(data_dir, output_dir, working_dir): This tests the following features: - Blip-up + Blip-down DWI series for TOPUP/Eddy - Eddy is run on a CPU - - dwidenoise is enabled implicitly. + - dwidenoise is enabled explicitly. Inputs ------ @@ -60,7 +60,7 @@ def test_dsdti_fmap(data_dir, output_dir, working_dir): '--write-graph', '--mem-mb=4096', '--output-resolution=5', - '--denoise-method=dwidenoise;shape:cuboid;subsample:1;demodulate:none;demean:none;' + '--denoise-method=dwidenoise2;shape:cuboid;subsample:1;demodulate:none;demean:none;' 'filter_method:truncate;aggregator:exclusive', ] @@ -100,7 +100,7 @@ def test_dscsdsi_fmap(data_dir, output_dir, working_dir): f'-w={work_dir}', '--boilerplate', '--sloppy', - '--denoise-method=dwidenoise', + '--denoise-method=dwidenoise2', '--b0-motion-corr-to=first', '--write-graph', '--mem-mb=4096', diff --git a/qsiprep/tests/test_interfaces_mrtrix3.py b/qsiprep/tests/test_interfaces_mrtrix3.py index 90bf69fd3..0bced4046 100644 --- a/qsiprep/tests/test_interfaces_mrtrix3.py +++ b/qsiprep/tests/test_interfaces_mrtrix3.py @@ -19,6 +19,34 @@ def test_dwidenoise(datasets, tmp_path_factory): in_img = nb.load(in_file) interface = mrtrix.DWIDenoise( + extent=(5, 5, 5), + in_file=in_file, + nthreads=1, + ) + results = interface.run(cwd=tmpdir) + + assert os.path.isfile(results.outputs.out_file) + denoised_img = nb.load(results.outputs.out_file) + assert denoised_img.shape == in_img.shape + + assert os.path.isfile(results.outputs.noise_image) + noise_img = nb.load(results.outputs.noise_image) + assert noise_img.shape == in_img.shape[:3] + assert noise_img.ndim == 3 + + assert os.path.isfile(results.outputs.out_report) + assert os.path.isfile(results.outputs.nmse_text) + + +def test_dwidenoise2(datasets, tmp_path_factory): + """Test qsiprep.interfaces.mrtrix.DWIDenoise2.""" + tmpdir = tmp_path_factory.mktemp('test_dwidenoise') + + in_dir = datasets['forrest_gump'] + in_file = os.path.join(in_dir, 'sub-01/ses-forrestgump/dwi/sub-01_ses-forrestgump_dwi.nii.gz') + in_img = nb.load(in_file) + + interface = mrtrix.DWIDenoise2( shape='cuboid', extent=(5, 5, 5), onepass=True, @@ -48,23 +76,23 @@ def test_dwidenoise(datasets, tmp_path_factory): ('cuboid', {'radius': 2.5}, "'radius' cannot be used"), ], ) -def test_dwidenoise_kernel_shape_validation(tmp_path, shape, kernel_option, error): +def test_dwidenoise2_kernel_shape_validation(tmp_path, shape, kernel_option, error): """Reject kernel options that do not apply to the selected shape.""" in_file = tmp_path / 'dwi.nii.gz' in_file.touch() - interface = mrtrix.DWIDenoise(in_file=in_file, shape=shape, **kernel_option) + interface = mrtrix.DWIDenoise2(in_file=in_file, shape=shape, **kernel_option) with pytest.raises(ValueError, match=error): _ = interface.cmdline -def test_dwidenoise_kernel_options_are_mutually_exclusive(tmp_path): +def test_dwidenoise2_kernel_options_are_mutually_exclusive(tmp_path): """Reject simultaneous spherical and cuboid kernel size options.""" in_file = tmp_path / 'dwi.nii.gz' in_file.touch() with pytest.raises(OSError, match='mutually exclusive'): - mrtrix.DWIDenoise( + mrtrix.DWIDenoise2( in_file=in_file, shape='sphere', radius=2.5, @@ -72,12 +100,12 @@ def test_dwidenoise_kernel_options_are_mutually_exclusive(tmp_path): ) -def test_dwidenoise_cli_parameters_reach_workflow(monkeypatch): - """Forward parsed DWIDenoise parameters to the workflow node.""" +def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): + """Forward parsed DWIDenoise2 parameters to the workflow node.""" monkeypatch.setattr( config.workflow, 'denoise_method', - 'dwidenoise;demodulate:nonlinear;decomposition:bdcsvd', + 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd', ) monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) monkeypatch.setattr(config.workflow, 'unringing_method', 'none') diff --git a/qsiprep/utils/misc.py b/qsiprep/utils/misc.py index 485536ca8..5c34ca1bb 100644 --- a/qsiprep/utils/misc.py +++ b/qsiprep/utils/misc.py @@ -59,15 +59,15 @@ def parse_denoise_method(spec): """Parse a denoising method and semicolon-delimited parameters. - Parameters use ``name:value`` syntax, for example - ``dwidenoise;demodulate:nonlinear;decomposition:bdcsvd``. + Parameters for dwidenoise2 use ``name:value`` syntax, for example + ``dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd``. """ elements = spec.split(';') method = elements[0].strip() - if method not in ('dwidenoise', 'patch2self', 'none'): + if method not in ('dwidenoise', 'dwidenoise2', 'patch2self', 'none'): raise ValueError(f'Unknown denoising method: {method!r}') - if len(elements) > 1 and method != 'dwidenoise': - raise ValueError(f'{method!r} does not accept DWIDenoise parameters') + if len(elements) > 1 and method != 'dwidenoise2': + raise ValueError(f'{method!r} does not accept DWIDenoise2 parameters') parameters = {} for element in elements[1:]: @@ -75,11 +75,11 @@ def parse_denoise_method(spec): name = name.strip() value = value.strip() if not separator or not name or not value: - raise ValueError(f'Invalid DWIDenoise parameter: {element!r}') + raise ValueError(f'Invalid DWIDenoise2 parameter: {element!r}') if name not in _DWIDENOISE_PARAMETERS: - raise ValueError(f'Unknown DWIDenoise parameter: {name!r}') + raise ValueError(f'Unknown DWIDenoise2 parameter: {name!r}') if name in parameters: - raise ValueError(f'Duplicate DWIDenoise parameter: {name!r}') + raise ValueError(f'Duplicate DWIDenoise2 parameter: {name!r}') if name in _DWIDENOISE_ENUM_PARAMETERS: choices = _DWIDENOISE_ENUM_PARAMETERS[name] diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index fa01b62f8..517fafc7d 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -25,6 +25,7 @@ ComplexToMagnitude, DWIBiasCorrect, DWIDenoise, + DWIDenoise2, MRDeGibbs, PolarToComplex, ) @@ -446,7 +447,7 @@ def get_buffernode(): dwi_denoise_window = config.workflow.dwi_denoise_window auto_str = '' - if denoise_method == 'dwidenoise' and dwi_denoise_window == 'auto': + if (denoise_method == 'dwidenoise') and dwi_denoise_window == 'auto': # Configure the denoising window import numpy as np @@ -458,7 +459,7 @@ def get_buffernode(): ) auto_str = 'n automatically-determined' - if (denoise_method == 'dwidenoise') and use_phase: + if denoise_method.startswith('dwidenoise') and use_phase: desc += ( 'Magnitude and phase DWI data were combined into a complex-valued file, ' 'then denoised using the Marchenko-Pastur PCA method implemented in dwidenoise ' @@ -488,13 +489,20 @@ def get_buffernode(): (phase_to_radians, combine_complex, [('phase_file', 'phase_file')]), ]) # fmt:skip - dwidenoise_inputs = {'shape': 'sphere', 'nthreads': omp_nthreads} - dwidenoise_inputs.update(dwidenoise_params) - denoiser = pe.Node( - DWIDenoise(**dwidenoise_inputs), - name='denoiser', - n_procs=omp_nthreads, - ) + if denoise_method == 'dwidenoise2': + dwidenoise_inputs = {'shape': 'sphere', 'nthreads': omp_nthreads} + dwidenoise_inputs.update(dwidenoise_params) + denoiser = pe.Node( + DWIDenoise2(**dwidenoise_inputs), + name='denoiser', + n_procs=omp_nthreads, + ) + else: + denoiser = pe.Node( + DWIDenoise(extent=dwi_denoise_window), + name='denoiser', + n_procs=omp_nthreads, + ) workflow.connect([ (combine_complex, denoiser, [('out_file', 'in_file')]), @@ -507,40 +515,46 @@ def get_buffernode(): 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': + elif denoise_method.startswith('dwidenoise2'): desc += ( 'DWI data were ' - 'denoised using the Marchenko-Pastur PCA method implemented in dwidenoise ' + 'denoised using the Marchenko-Pastur PCA method implemented in dwidenoise2 ' '[@mrtrix3; @dwidenoise1; @dwidenoise2; @cordero2019complex] ' f'with a{auto_str} window size of {dwi_denoise_window} voxels. ' ) last_step = 'After MP-PCA, ' - dwidenoise_inputs = { - 'shape': 'cuboid', - 'extent': (dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), - # Fixed odd extents are incompatible with the changing subsampling parity of - # iterative mode, so reproduce legacy fixed-window behavior in one pass. - 'onepass': True, - 'subsample': 1, - 'nthreads': omp_nthreads, - } - if dwidenoise_params.get('shape') == 'sphere': - for parameter in ('extent', 'onepass', 'subsample'): - if parameter not in dwidenoise_params: - dwidenoise_inputs.pop(parameter) - dwidenoise_inputs.update(dwidenoise_params) - denoiser = pe.Node( - DWIDenoise(**dwidenoise_inputs), - name='denoiser', - n_procs=omp_nthreads, - ) + if denoise_method == 'dwidenoise2': + dwidenoise_inputs = { + 'shape': 'cuboid', + 'extent': (dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), + # Fixed odd extents are incompatible with the changing subsampling parity of + # iterative mode, so reproduce legacy fixed-window behavior in one pass. + 'onepass': True, + 'subsample': 1, + 'nthreads': omp_nthreads, + } + if dwidenoise_params.get('shape') == 'sphere': + for parameter in ('extent', 'onepass', 'subsample'): + if parameter not in dwidenoise_params: + dwidenoise_inputs.pop(parameter) + dwidenoise_inputs.update(dwidenoise_params) + denoiser = pe.Node( + DWIDenoise2(**dwidenoise_inputs), + name='denoiser', + n_procs=omp_nthreads, + ) + else: + denoiser = pe.Node( + DWIDenoise(extent=dwi_denoise_window), + name='denoiser', + n_procs=omp_nthreads, + ) else: desc += ( "DWI data were denoised using DiPy's Patch2Self algorithm [@dipy; @patch2self] " @@ -553,10 +567,13 @@ def get_buffernode(): n_procs=omp_nthreads, ) - if denoise_method in ('dwidenoise', 'patch2self'): + if denoise_method in ('dwidenoise2', 'patch2self'): workflow.connect([(inputnode, denoiser, [('bval_file', 'bval_file')])]) - if (denoise_method in ('dwidenoise', 'patch2self')) and not use_phase: + if denoise_method == 'dwidenoise2': + workflow.connect([(inputnode, denoiser, [('bvec_file', 'bvec_file')])]) + + if (denoise_method != 'none') and not use_phase: workflow.connect([ (buffernodes[-2], denoiser, [('dwi_file', 'in_file')]), (denoiser, ds_report_denoising, [('out_report', 'in_file')]), From 70ccdabd4f14010ec3e0e84c996c72d95b311fd0 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 14 Jul 2026 11:50:50 -0400 Subject: [PATCH 05/14] Fix some things. --- qsiprep/tests/test_interfaces_mrtrix3.py | 2 +- qsiprep/tests/test_utils_misc.py | 17 +++++++++-------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/qsiprep/tests/test_interfaces_mrtrix3.py b/qsiprep/tests/test_interfaces_mrtrix3.py index 0bced4046..d5d30e935 100644 --- a/qsiprep/tests/test_interfaces_mrtrix3.py +++ b/qsiprep/tests/test_interfaces_mrtrix3.py @@ -40,7 +40,7 @@ def test_dwidenoise(datasets, tmp_path_factory): def test_dwidenoise2(datasets, tmp_path_factory): """Test qsiprep.interfaces.mrtrix.DWIDenoise2.""" - tmpdir = tmp_path_factory.mktemp('test_dwidenoise') + tmpdir = tmp_path_factory.mktemp('test_dwidenoise2') in_dir = datasets['forrest_gump'] in_file = os.path.join(in_dir, 'sub-01/ses-forrestgump/dwi/sub-01_ses-forrestgump_dwi.nii.gz') diff --git a/qsiprep/tests/test_utils_misc.py b/qsiprep/tests/test_utils_misc.py index a0a03f745..5ab2e4c07 100644 --- a/qsiprep/tests/test_utils_misc.py +++ b/qsiprep/tests/test_utils_misc.py @@ -55,7 +55,7 @@ def test_parse_denoise_method_parameters(): 'onepass:true;radius:2.5;subsample:2,2,2' ) - assert method == 'dwidenoise' + assert method == 'dwidenoise2' assert parameters == { 'demodulate': 'nonlinear', 'decomposition': 'bdcsvd', @@ -70,12 +70,13 @@ def test_parse_denoise_method_parameters(): [ 'unknown', 'patch2self;decomposition:bdcsvd', - 'dwidenoise;decomposition', - 'dwidenoise;unknown:value', - 'dwidenoise;decomposition:bdcsvd;decomposition:selfadjoint', - 'dwidenoise;decomposition:invalid', - 'dwidenoise;onepass:maybe', - 'dwidenoise;extent:1,2', + 'dwidenoise;decomposition:bdcsvd', + 'dwidenoise2;decomposition', + 'dwidenoise2;unknown:value', + 'dwidenoise2;decomposition:bdcsvd;decomposition:selfadjoint', + 'dwidenoise2;decomposition:invalid', + 'dwidenoise2;onepass:maybe', + 'dwidenoise2;extent:1,2', ], ) def test_parse_denoise_method_rejects_invalid_specs(spec): @@ -84,7 +85,7 @@ def test_parse_denoise_method_rejects_invalid_specs(spec): def test_denoise_method_cli_parameter(tmp_path): - spec = 'dwidenoise;demodulate:nonlinear;decomposition:bdcsvd' + spec = 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd' opts = _build_parser().parse_args( [ str(tmp_path), From 8a02ff6e01d8c3787230000fbe1009df3a216a76 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 14 Jul 2026 12:11:06 -0400 Subject: [PATCH 06/14] Update. --- qsiprep/tests/test_utils_misc.py | 2 +- qsiprep/workflows/dwi/merge.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/qsiprep/tests/test_utils_misc.py b/qsiprep/tests/test_utils_misc.py index 5ab2e4c07..043f0c300 100644 --- a/qsiprep/tests/test_utils_misc.py +++ b/qsiprep/tests/test_utils_misc.py @@ -51,7 +51,7 @@ def test_angle_between_finite_for_zero_vector(): def test_parse_denoise_method_parameters(): method, parameters = parse_denoise_method( - 'dwidenoise;demodulate:nonlinear;decomposition:bdcsvd;' + 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd;' 'onepass:true;radius:2.5;subsample:2,2,2' ) diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 517fafc7d..5c7c8c84f 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -419,7 +419,7 @@ def get_buffernode(): # Which steps to apply? denoise_method, dwidenoise_params = parse_denoise_method(config.workflow.denoise_method) unringing_method = config.workflow.unringing_method - do_denoise = denoise_method in ('patch2self', 'dwidenoise') + do_denoise = denoise_method in ('patch2self', 'dwidenoise', 'dwidenoise2') do_unringing = config.workflow.unringing_method in ('mrdegibbs', 'rpg') harmonize_b0s = not config.workflow.no_b0_harmonization From d59d2b80a86e5b60400e1498dfa264854ff92fa2 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 4 Aug 2026 10:20:22 -0400 Subject: [PATCH 07/14] Try fixing failing tests. --- ...s_mrtrix3.py => test_interfaces_mrtrix.py} | 53 +++++++++++++++++-- qsiprep/workflows/dwi/merge.py | 39 +++++++++----- 2 files changed, 75 insertions(+), 17 deletions(-) rename qsiprep/tests/{test_interfaces_mrtrix3.py => test_interfaces_mrtrix.py} (66%) diff --git a/qsiprep/tests/test_interfaces_mrtrix3.py b/qsiprep/tests/test_interfaces_mrtrix.py similarity index 66% rename from qsiprep/tests/test_interfaces_mrtrix3.py rename to qsiprep/tests/test_interfaces_mrtrix.py index d5d30e935..b5777802c 100644 --- a/qsiprep/tests/test_interfaces_mrtrix3.py +++ b/qsiprep/tests/test_interfaces_mrtrix.py @@ -1,4 +1,4 @@ -"""Tests for the qsiprep.interfaces.dipy module.""" +"""Tests for the qsiprep.interfaces.mrtrix module.""" import os @@ -47,8 +47,8 @@ def test_dwidenoise2(datasets, tmp_path_factory): in_img = nb.load(in_file) interface = mrtrix.DWIDenoise2( - shape='cuboid', - extent=(5, 5, 5), + shape='sphere', + radius=3, onepass=True, subsample=1, in_file=in_file, @@ -100,6 +100,53 @@ def test_dwidenoise2_kernel_options_are_mutually_exclusive(tmp_path): ) +@pytest.mark.parametrize('use_phase', [False, True]) +def test_dwidenoise_workflow_uses_dwidenoise(monkeypatch, use_phase): + """Build a DWIDenoise node, not Patch2Self, when ``dwidenoise`` is requested.""" + monkeypatch.setattr(config.workflow, 'denoise_method', 'dwidenoise') + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=use_phase, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + assert isinstance(denoiser.interface, mrtrix.DWIDenoise) + assert denoiser.inputs.extent == (5, 5, 5) + assert denoiser.inputs.nthreads == 1 + + +@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2']) +def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method): + """Resolve the default ``auto`` window size for every dwidenoise variant.""" + monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 'auto') + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=False, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + # cbrt(30) rounded up to the closest odd integer + assert denoiser.inputs.extent == (5, 5, 5) + + def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): """Forward parsed DWIDenoise2 parameters to the workflow node.""" monkeypatch.setattr( diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 5c7c8c84f..eab5f89ca 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -447,7 +447,7 @@ def get_buffernode(): dwi_denoise_window = config.workflow.dwi_denoise_window auto_str = '' - if (denoise_method == 'dwidenoise') and dwi_denoise_window == 'auto': + if denoise_method.startswith('dwidenoise') and dwi_denoise_window == 'auto': # Configure the denoising window import numpy as np @@ -462,7 +462,8 @@ def get_buffernode(): if denoise_method.startswith('dwidenoise') and use_phase: desc += ( 'Magnitude and phase DWI data were combined into a complex-valued file, ' - 'then denoised using the Marchenko-Pastur PCA method implemented in dwidenoise ' + 'then denoised using the Marchenko-Pastur PCA method implemented in ' + f'{denoise_method} ' '[@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 ' @@ -499,7 +500,10 @@ def get_buffernode(): ) else: denoiser = pe.Node( - DWIDenoise(extent=dwi_denoise_window), + DWIDenoise( + extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), + nthreads=omp_nthreads, + ), name='denoiser', n_procs=omp_nthreads, ) @@ -520,10 +524,10 @@ def get_buffernode(): (split_complex, buffernodes[-1], [('out_file', 'dwi_file')]), ]) # fmt:skip - elif denoise_method.startswith('dwidenoise2'): + elif denoise_method.startswith('dwidenoise'): desc += ( 'DWI data were ' - 'denoised using the Marchenko-Pastur PCA method implemented in dwidenoise2 ' + f'denoised using the Marchenko-Pastur PCA method implemented in {denoise_method} ' '[@mrtrix3; @dwidenoise1; @dwidenoise2; @cordero2019complex] ' f'with a{auto_str} window size of {dwi_denoise_window} voxels. ' ) @@ -531,18 +535,22 @@ def get_buffernode(): if denoise_method == 'dwidenoise2': dwidenoise_inputs = { - 'shape': 'cuboid', - 'extent': (dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), - # Fixed odd extents are incompatible with the changing subsampling parity of - # iterative mode, so reproduce legacy fixed-window behavior in one pass. - 'onepass': True, - 'subsample': 1, + 'shape': 'sphere', + 'radius': dwi_denoise_window, 'nthreads': omp_nthreads, } - if dwidenoise_params.get('shape') == 'sphere': - for parameter in ('extent', 'onepass', 'subsample'): + if dwidenoise_params.get('shape') == 'cuboid': + for parameter in ('radius',): if parameter not in dwidenoise_params: dwidenoise_inputs.pop(parameter) + + # cuboid uses extent instead of radius + dwidenoise_inputs['extent'] = ( + dwi_denoise_window, + dwi_denoise_window, + dwi_denoise_window, + ) + dwidenoise_inputs.update(dwidenoise_params) denoiser = pe.Node( DWIDenoise2(**dwidenoise_inputs), @@ -551,7 +559,10 @@ def get_buffernode(): ) else: denoiser = pe.Node( - DWIDenoise(extent=dwi_denoise_window), + DWIDenoise( + extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), + nthreads=omp_nthreads, + ), name='denoiser', n_procs=omp_nthreads, ) From 9332fabd2365365ed9e049e6e7c81a6ec30611b3 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 4 Aug 2026 10:21:44 -0400 Subject: [PATCH 08/14] Update lint.yml --- .github/workflows/lint.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml index 9d4e6832e..81f276d67 100644 --- a/.github/workflows/lint.yml +++ b/.github/workflows/lint.yml @@ -24,8 +24,8 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v7 - - run: pipx run ruff check . - - run: pipx run ruff format --diff . + - run: pipx run --spec ruff==0.14.11 ruff check . + - run: pipx run --spec ruff==0.14.11 ruff format --diff . codespell: name: Check for spelling errors From bae2d1cb6548dccc9da9c77eec5fdd9e117c16a4 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Tue, 4 Aug 2026 11:26:02 -0400 Subject: [PATCH 09/14] Keep working. --- Dockerfile | 2 +- qsiprep/tests/test_interfaces_mrtrix.py | 77 ----------------------- qsiprep/tests/test_workflows_merge.py | 83 +++++++++++++++++++++++++ 3 files changed, 84 insertions(+), 78 deletions(-) create mode 100644 qsiprep/tests/test_workflows_merge.py diff --git a/Dockerfile b/Dockerfile index b4a398879..7116adf32 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,4 +1,4 @@ -ARG BASE_IMAGE=pennlinc/qsiprep-base:20260415 +ARG BASE_IMAGE=pennlinc/qsiprep-base:20260804 ARG DWIDENOISE2_COMMIT=892d5e8dd8f453ce1a561878f7dcb3998ae258ba ARG MRTRIX3_DWIDENOISE2_COMMIT=fa6ee952913fbc1df79aeca745600852155f533f diff --git a/qsiprep/tests/test_interfaces_mrtrix.py b/qsiprep/tests/test_interfaces_mrtrix.py index b5777802c..e3ff331c3 100644 --- a/qsiprep/tests/test_interfaces_mrtrix.py +++ b/qsiprep/tests/test_interfaces_mrtrix.py @@ -5,9 +5,7 @@ import nibabel as nb import pytest -from qsiprep import config from qsiprep.interfaces import mrtrix -from qsiprep.workflows.dwi.merge import init_dwi_denoising_wf def test_dwidenoise(datasets, tmp_path_factory): @@ -98,78 +96,3 @@ def test_dwidenoise2_kernel_options_are_mutually_exclusive(tmp_path): radius=2.5, extent=(5, 5, 5), ) - - -@pytest.mark.parametrize('use_phase', [False, True]) -def test_dwidenoise_workflow_uses_dwidenoise(monkeypatch, use_phase): - """Build a DWIDenoise node, not Patch2Self, when ``dwidenoise`` is requested.""" - monkeypatch.setattr(config.workflow, 'denoise_method', 'dwidenoise') - monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) - monkeypatch.setattr(config.workflow, 'unringing_method', 'none') - monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) - monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) - - workflow = init_dwi_denoising_wf( - source_file='sub-01_dwi.nii.gz', - partial_fourier=1.0, - phase_encoding_direction='j', - n_volumes=30, - use_phase=use_phase, - do_biascorr=False, - ) - denoiser = workflow.get_node('denoiser') - - assert isinstance(denoiser.interface, mrtrix.DWIDenoise) - assert denoiser.inputs.extent == (5, 5, 5) - assert denoiser.inputs.nthreads == 1 - - -@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2']) -def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method): - """Resolve the default ``auto`` window size for every dwidenoise variant.""" - monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) - monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 'auto') - monkeypatch.setattr(config.workflow, 'unringing_method', 'none') - monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) - monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) - - workflow = init_dwi_denoising_wf( - source_file='sub-01_dwi.nii.gz', - partial_fourier=1.0, - phase_encoding_direction='j', - n_volumes=30, - use_phase=False, - do_biascorr=False, - ) - denoiser = workflow.get_node('denoiser') - - # cbrt(30) rounded up to the closest odd integer - assert denoiser.inputs.extent == (5, 5, 5) - - -def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): - """Forward parsed DWIDenoise2 parameters to the workflow node.""" - monkeypatch.setattr( - config.workflow, - 'denoise_method', - 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd', - ) - monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) - monkeypatch.setattr(config.workflow, 'unringing_method', 'none') - monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) - monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) - - workflow = init_dwi_denoising_wf( - source_file='sub-01_dwi.nii.gz', - partial_fourier=1.0, - phase_encoding_direction='j', - n_volumes=30, - use_phase=False, - do_biascorr=False, - ) - denoiser = workflow.get_node('denoiser') - - assert denoiser.inputs.demodulate == 'nonlinear' - assert denoiser.inputs.decomposition == 'bdcsvd' - assert denoiser.inputs.onepass is True - assert denoiser.inputs.subsample == 1 diff --git a/qsiprep/tests/test_workflows_merge.py b/qsiprep/tests/test_workflows_merge.py new file mode 100644 index 000000000..4616be74d --- /dev/null +++ b/qsiprep/tests/test_workflows_merge.py @@ -0,0 +1,83 @@ +"""Tests for the qsiprep.workflows.dwi.merge module.""" + +import pytest + +from qsiprep import config + +from qsiprep.interfaces import mrtrix +from qsiprep.workflows.dwi.merge import init_dwi_denoising_wf + + +@pytest.mark.parametrize('use_phase', [False, True]) +def test_dwidenoise_workflow_uses_dwidenoise(monkeypatch, use_phase): + """Build a DWIDenoise node, not Patch2Self, when ``dwidenoise`` is requested.""" + monkeypatch.setattr(config.workflow, 'denoise_method', 'dwidenoise') + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=use_phase, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + assert isinstance(denoiser.interface, mrtrix.DWIDenoise) + assert denoiser.inputs.extent == (5, 5, 5) + assert denoiser.inputs.nthreads == 1 + + +@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2']) +def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method): + """Resolve the default ``auto`` window size for every dwidenoise variant.""" + monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 'auto') + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=False, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + # cbrt(30) rounded up to the closest odd integer + assert denoiser.inputs.extent == (5, 5, 5) + + +def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): + """Forward parsed DWIDenoise2 parameters to the workflow node.""" + monkeypatch.setattr( + config.workflow, + 'denoise_method', + 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd', + ) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=False, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + assert denoiser.inputs.demodulate == 'nonlinear' + assert denoiser.inputs.decomposition == 'bdcsvd' + assert denoiser.inputs.onepass is True + assert denoiser.inputs.subsample == 1 From 19ba74bb8923c80e7054e0565b010ede169ed387 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 5 Aug 2026 09:25:44 -0400 Subject: [PATCH 10/14] Pass mask into denoiser interfaces. This won't actually have any effect bc dwidenoise2 doesn't use the mask yet. --- qsiprep/interfaces/mrtrix.py | 6 +- qsiprep/tests/conftest.py | 31 ++ qsiprep/tests/test_interfaces_mrtrix.py | 17 ++ qsiprep/tests/test_utils_misc.py | 10 + qsiprep/tests/test_workflows_merge.py | 366 +++++++++++++++++++++++- qsiprep/utils/misc.py | 4 + qsiprep/workflows/dwi/merge.py | 62 +++- 7 files changed, 476 insertions(+), 20 deletions(-) diff --git a/qsiprep/interfaces/mrtrix.py b/qsiprep/interfaces/mrtrix.py index 26927d1ec..9c3763794 100644 --- a/qsiprep/interfaces/mrtrix.py +++ b/qsiprep/interfaces/mrtrix.py @@ -183,7 +183,7 @@ def _get_plotting_images(self): class DWIDenoise2InputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): in_file = File(exists=True, argstr='%s', position=-2, mandatory=True, desc='input DWI image') - mask = File(exists=True, desc='mask image used only to define the visual report contour') + mask = File(exists=True, desc='mask image') onepass = traits.Bool(argstr='-onepass', desc='estimate noise and denoise in one pass') datatype = traits.Enum( 'float32', @@ -416,7 +416,9 @@ def _format_arg(self, name, spec, value): if name in ('extent', 'subsample') and not isinstance(value, int): value = ','.join(str(item) for item in value) elif name == 'bvec_file': - value = (value, self.inputs.bval_file) + # -fslgrad takes both files, so format them here rather than passing a tuple + # to a File trait, which nipype would try to shell-quote as a single value. + return spec.argstr % (value, self.inputs.bval_file) return super()._format_arg(name, spec, value) def _parse_inputs(self, skip=None): diff --git a/qsiprep/tests/conftest.py b/qsiprep/tests/conftest.py index 00ecb08ff..fcf41d88d 100644 --- a/qsiprep/tests/conftest.py +++ b/qsiprep/tests/conftest.py @@ -1,6 +1,7 @@ """Fixtures for the CircleCI tests.""" import os +from pathlib import Path import pytest @@ -40,4 +41,34 @@ def datasets(data_dir): """Locate downloaded datasets.""" dsets = {} dsets['forrest_gump'] = os.path.join(data_dir, 'forrest_gump') + dsets['nibs-ci'] = os.path.join(data_dir, 'nibs-ci') return dsets + + +@pytest.fixture(scope='session') +def nibs_dwi(data_dir): + """Locate the nibs-ci DWI series used to test the denoising workflow. + + The series is small (48x48x29x76) and has both magnitude and phase data, so it can + exercise the complex-valued denoising paths without a long runtime. + + Tests using this fixture are skipped when the dataset is unavailable, which keeps + them runnable outside of the container. + """ + if not data_dir: + pytest.skip('--data_dir was not provided') + + dwi_dir = Path(data_dir) / 'nibs-ci' / 'sub-22449' / 'ses-01' / 'dwi' + stem = 'sub-22449_ses-01_acq-HBCD75_rec-norm_dir-AP_run-01' + files = { + 'dwi_file': dwi_dir / f'{stem}_part-mag_dwi.nii.gz', + 'phase_file': dwi_dir / f'{stem}_part-phase_dwi.nii.gz', + 'bval_file': dwi_dir / f'{stem}_part-mag_dwi.bval', + 'bvec_file': dwi_dir / f'{stem}_part-mag_dwi.bvec', + 'json_file': dwi_dir / f'{stem}_part-mag_dwi.json', + } + missing = sorted(str(f) for f in files.values() if not f.is_file()) + if missing: + pytest.skip(f'nibs-ci dataset is unavailable; missing {missing}') + + return {key: str(value) for key, value in files.items()} diff --git a/qsiprep/tests/test_interfaces_mrtrix.py b/qsiprep/tests/test_interfaces_mrtrix.py index e3ff331c3..1f3f67bc4 100644 --- a/qsiprep/tests/test_interfaces_mrtrix.py +++ b/qsiprep/tests/test_interfaces_mrtrix.py @@ -96,3 +96,20 @@ def test_dwidenoise2_kernel_options_are_mutually_exclusive(tmp_path): radius=2.5, extent=(5, 5, 5), ) + + +def test_dwidenoise2_formats_fslgrad(tmp_path): + """Pass the bvec and bval files to dwidenoise2 as a single -fslgrad option.""" + in_file = tmp_path / 'dwi.nii.gz' + bvec_file = tmp_path / 'dwi.bvec' + bval_file = tmp_path / 'dwi.bval' + for path in (in_file, bvec_file, bval_file): + path.touch() + + interface = mrtrix.DWIDenoise2( + in_file=in_file, + bvec_file=bvec_file, + bval_file=bval_file, + ) + + assert f'-fslgrad {bvec_file} {bval_file}' in interface.cmdline diff --git a/qsiprep/tests/test_utils_misc.py b/qsiprep/tests/test_utils_misc.py index 043f0c300..73feb2bf8 100644 --- a/qsiprep/tests/test_utils_misc.py +++ b/qsiprep/tests/test_utils_misc.py @@ -84,6 +84,16 @@ def test_parse_denoise_method_rejects_invalid_specs(spec): parse_denoise_method(spec) +def test_parse_denoise_method_rejects_cuboid_shape(): + """Reject cuboid kernels, which need an even extent that QSIPrep never produces.""" + with pytest.raises(ValueError, match='not supported yet'): + parse_denoise_method('dwidenoise2;shape:cuboid') + + method, parameters = parse_denoise_method('dwidenoise2;shape:sphere') + assert method == 'dwidenoise2' + assert parameters == {'shape': 'sphere'} + + def test_denoise_method_cli_parameter(tmp_path): spec = 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd' opts = _build_parser().parse_args( diff --git a/qsiprep/tests/test_workflows_merge.py b/qsiprep/tests/test_workflows_merge.py index 4616be74d..854411f0f 100644 --- a/qsiprep/tests/test_workflows_merge.py +++ b/qsiprep/tests/test_workflows_merge.py @@ -1,10 +1,20 @@ """Tests for the qsiprep.workflows.dwi.merge module.""" +import json +import os +from pathlib import Path + +import nibabel as nb +import numpy as np +import pandas as pd import pytest +from nipype.interfaces import io as nio +from nipype.interfaces.base import isdefined +from nipype.pipeline import engine as pe from qsiprep import config - from qsiprep.interfaces import mrtrix +from qsiprep.interfaces.dipy import Patch2Self from qsiprep.workflows.dwi.merge import init_dwi_denoising_wf @@ -15,6 +25,7 @@ def test_dwidenoise_workflow_uses_dwidenoise(monkeypatch, use_phase): monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) monkeypatch.setattr(config.workflow, 'unringing_method', 'none') monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) workflow = init_dwi_denoising_wf( @@ -32,13 +43,22 @@ def test_dwidenoise_workflow_uses_dwidenoise(monkeypatch, use_phase): assert denoiser.inputs.nthreads == 1 -@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2']) -def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method): +@pytest.mark.parametrize( + ('denoise_method', 'kernel_input'), + [ + # dwidenoise takes a cuboid extent, while dwidenoise2 defaults to a spherical + # kernel, which is sized with a radius instead + ('dwidenoise', {'extent': (5, 5, 5)}), + ('dwidenoise2', {'radius': 5.0}), + ], +) +def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method, kernel_input): """Resolve the default ``auto`` window size for every dwidenoise variant.""" monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 'auto') monkeypatch.setattr(config.workflow, 'unringing_method', 'none') monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) workflow = init_dwi_denoising_wf( @@ -52,7 +72,8 @@ def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method): denoiser = workflow.get_node('denoiser') # cbrt(30) rounded up to the closest odd integer - assert denoiser.inputs.extent == (5, 5, 5) + for name, value in kernel_input.items(): + assert getattr(denoiser.inputs, name) == value def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): @@ -65,6 +86,7 @@ def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) monkeypatch.setattr(config.workflow, 'unringing_method', 'none') monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) workflow = init_dwi_denoising_wf( @@ -72,12 +94,342 @@ def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): partial_fourier=1.0, phase_encoding_direction='j', n_volumes=30, - use_phase=False, + # demodulation is only valid for complex-valued data + use_phase=True, do_biascorr=False, ) denoiser = workflow.get_node('denoiser') assert denoiser.inputs.demodulate == 'nonlinear' assert denoiser.inputs.decomposition == 'bdcsvd' - assert denoiser.inputs.onepass is True - assert denoiser.inputs.subsample == 1 + # Parameters that weren't requested are left at the dwidenoise2 defaults + assert not isdefined(denoiser.inputs.onepass) + assert not isdefined(denoiser.inputs.subsample) + + +@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2']) +def test_denoising_wf_builds_one_mask_for_denoising_and_biascorr(monkeypatch, denoise_method): + """Build the brain mask once and hand it to both the denoiser and bias correction.""" + monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=False, + do_biascorr=True, + ) + + node_names = [node.name for node in workflow._get_all_nodes()] + assert node_names.count('quick_mask') == 1 + assert node_names.count('get_b0s') == 1 + + quick_mask = workflow.get_node('quick_mask') + consumers = { + (dest.name, dest_field) + for src, dest, data in workflow._graph.edges(data=True) + if src is quick_mask + for _, dest_field in data['connect'] + } + assert consumers == {('denoiser', 'mask'), ('biascorr', 'mask')} + + # The mask has to come from the raw series: denoising runs first, so deriving it from + # any later buffer would be circular + get_b0s = workflow.get_node('get_b0s') + assert {src.name for src, dest, _ in workflow._graph.edges(data=True) if dest is get_b0s} == { + 'inputnode' + } + + +@pytest.mark.parametrize('demodulate', ['linear', 'nonlinear']) +def test_dwidenoise2_rejects_demodulation_without_phase(monkeypatch, demodulate): + """Reject phase demodulation unless phase data are available. + + ``dwidenoise2`` errors out partway through a run when asked to demodulate + magnitude-only data, so the workflow rejects the request up front instead. + """ + monkeypatch.setattr(config.workflow, 'denoise_method', f'dwidenoise2;demodulate:{demodulate}') + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + kwargs = { + 'source_file': 'sub-01_dwi.nii.gz', + 'partial_fourier': 1.0, + 'phase_encoding_direction': 'j', + 'n_volumes': 30, + 'do_biascorr': False, + } + + with pytest.raises(ValueError, match='magnitude-only data'): + init_dwi_denoising_wf(use_phase=False, **kwargs) + + # The same request is fine once phase data are available + workflow = init_dwi_denoising_wf(use_phase=True, **kwargs) + assert workflow.get_node('denoiser').inputs.demodulate == demodulate + + +def _run_denoising_wf( + monkeypatch, + tmp_path, + nibs_dwi, + denoise_method, + use_phase, + dwi_denoise_window='auto', +): + """Build and execute a denoising workflow on the nibs-ci DWI series. + + Unringing, bias correction and b=0 harmonization are all disabled so that only the + denoising step is exercised. + + Returns + ------- + nodes : dict + The executed nodes, keyed by node name. + sink_dir : :obj:`pathlib.Path` + Directory holding the files that reached the workflow's ``outputnode``. + """ + monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', dwi_denoise_window) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + metadata = json.loads(Path(nibs_dwi['json_file']).read_text()) + + denoise_wf = init_dwi_denoising_wf( + source_file=nibs_dwi['dwi_file'], + partial_fourier=metadata['PartialFourier'], + phase_encoding_direction=metadata['PhaseEncodingDirection'].replace('-', ''), + n_volumes=nb.load(nibs_dwi['dwi_file']).shape[3], + use_phase=use_phase, + do_biascorr=False, + ) + denoise_wf.inputs.inputnode.dwi_file = nibs_dwi['dwi_file'] + denoise_wf.inputs.inputnode.bval_file = nibs_dwi['bval_file'] + denoise_wf.inputs.inputnode.bvec_file = nibs_dwi['bvec_file'] + if use_phase: + denoise_wf.inputs.inputnode.dwi_phase_file = nibs_dwi['phase_file'] + + # nipype prunes IdentityInterface nodes out of the execution graph, so ``outputnode`` + # can't be inspected directly. Routing its outputs to a DataSink both survives the + # pruning and checks that the workflow really connects them. + sink_dir = tmp_path / 'sink' + sink = pe.Node( + nio.DataSink(base_directory=str(sink_dir), parameterization=False), + name='sink', + ) + workflow = pe.Workflow(name='denoise_test_wf', base_dir=str(tmp_path)) + workflow.connect([ + (denoise_wf, sink, [ + ('outputnode.dwi_file', 'dwi_file'), + ('outputnode.noise_image', 'noise_image'), + ('outputnode.confounds', 'confounds'), + ]), + ]) # fmt:skip + + # nipype raises if any node fails, so a returned graph means every node ran + graph = workflow.run(plugin='Linear') + + return {node.name: node for node in graph.nodes}, sink_dir + + +def _field_of_view(img): + """Return the spatial extent of an image in mm.""" + return np.array(img.shape[:3]) * np.array(img.header.get_zooms()[:3]) + + +def _sink_output(sink_dir, field): + """Return the single file the DataSink wrote for ``field``.""" + matches = sorted((sink_dir / field).glob('*')) + assert len(matches) == 1, f'expected one {field} file, found {matches}' + return matches[0] + + +def _assert_denoiser_is_masked(nodes, denoise_method, raw_file): + """Check that the dwidenoise variants are handed a brain mask built from the raw data.""" + uses_mask = denoise_method.startswith('dwidenoise') + assert ('quick_mask' in nodes) is uses_mask + if not uses_mask: + assert not isdefined(nodes['denoiser'].inputs.mask) + return + + mask_file = nodes['quick_mask'].result.outputs.out_mask + assert nodes['denoiser'].inputs.mask == mask_file + + raw_img = nb.load(raw_file) + mask_img = nb.load(mask_file) + assert mask_img.shape == raw_img.shape[:3] + assert np.allclose(mask_img.affine, raw_img.affine) + # A mask that selected everything or nothing would silently defeat the point + mask_data = mask_img.get_fdata() + assert set(np.unique(mask_data)) <= {0.0, 1.0} + assert 0 < mask_data.sum() < mask_data.size + + +def _assert_denoising_outputs(nodes, sink_dir, raw_file): + """Check the files produced by an executed denoising workflow.""" + raw_img = nb.load(raw_file) + denoiser_outputs = nodes['denoiser'].result.outputs + + denoised_img = nb.load(_sink_output(sink_dir, 'dwi_file')) + assert denoised_img.shape == raw_img.shape + assert np.allclose(denoised_img.affine, raw_img.affine) + assert np.all(np.isfinite(denoised_img.get_fdata())) + # The workflow always returns magnitude data, even when it denoises complex data + assert not np.issubdtype(denoised_img.header.get_data_dtype(), np.complexfloating) + + noise_img = nb.load(_sink_output(sink_dir, 'noise_image')) + assert noise_img.ndim == 3 + # dwidenoise2 subsamples by default, so its noise map sits on a coarser grid than the + # input. Whatever the grid, it has to cover the same field of view. + assert np.allclose(_field_of_view(noise_img), _field_of_view(raw_img), rtol=0.05) + noise_data = noise_img.get_fdata() + finite = np.isfinite(noise_data) + assert finite.any() + assert np.all(noise_data[finite] >= 0) + + assert len(pd.read_csv(_sink_output(sink_dir, 'confounds'))) == raw_img.shape[3] + + assert os.path.isfile(denoiser_outputs.out_report) + + +@pytest.mark.parametrize( + ('denoise_method', 'dwi_denoise_window', 'interface', 'expected_inputs'), + [ + pytest.param( + 'dwidenoise', 5, mrtrix.DWIDenoise, {'extent': (5, 5, 5)}, id='dwidenoise_window5' + ), + pytest.param( + 'dwidenoise', 'auto', mrtrix.DWIDenoise, {'extent': (5, 5, 5)}, id='dwidenoise_auto' + ), + pytest.param( + 'dwidenoise2', + 'auto', + mrtrix.DWIDenoise2, + {'shape': 'sphere', 'radius': 5.0}, + id='dwidenoise2_sphere', + ), + pytest.param( + 'dwidenoise2;decomposition:bdcsvd', + 'auto', + mrtrix.DWIDenoise2, + {'decomposition': 'bdcsvd'}, + id='dwidenoise2_bdcsvd', + ), + pytest.param( + 'dwidenoise2;filter_method:optthresh', + 'auto', + mrtrix.DWIDenoise2, + {'filter_method': 'optthresh'}, + id='dwidenoise2_optthresh', + ), + pytest.param( + 'dwidenoise2;estimator:MRM2023', + 'auto', + mrtrix.DWIDenoise2, + {'estimator': 'MRM2023'}, + id='dwidenoise2_mrm2023', + marks=pytest.mark.xfail( + strict=True, + reason=( + 'the MRM2023 estimator returns a noise map with negative values ' + '(~41% of voxels, down to -15.7 on this series), while every other ' + 'estimator stays positive' + ), + ), + ), + pytest.param('patch2self', 'auto', Patch2Self, {}, id='patch2self'), + ], +) +def test_denoising_wf_magnitude( + monkeypatch, + tmp_path, + nibs_dwi, + denoise_method, + dwi_denoise_window, + interface, + expected_inputs, +): + """Denoise magnitude-only DWI data with each supported method.""" + nodes, sink_dir = _run_denoising_wf( + monkeypatch, + tmp_path, + nibs_dwi, + denoise_method=denoise_method, + use_phase=False, + dwi_denoise_window=dwi_denoise_window, + ) + + denoiser = nodes['denoiser'] + assert isinstance(denoiser.interface, interface) + for name, value in expected_inputs.items(): + assert getattr(denoiser.inputs, name) == value + + # Magnitude-only data never goes through the complex-valued path + assert 'combine_complex' not in nodes + assert 'split_complex' not in nodes + + _assert_denoiser_is_masked(nodes, denoise_method, nibs_dwi['dwi_file']) + _assert_denoising_outputs(nodes, sink_dir, nibs_dwi['dwi_file']) + + +@pytest.mark.parametrize( + ('denoise_method', 'interface', 'expected_inputs'), + [ + pytest.param('dwidenoise', mrtrix.DWIDenoise, {'extent': (5, 5, 5)}, id='dwidenoise'), + pytest.param('dwidenoise2', mrtrix.DWIDenoise2, {'shape': 'sphere'}, id='dwidenoise2'), + pytest.param( + 'dwidenoise2;demodulate:nonlinear', + mrtrix.DWIDenoise2, + {'demodulate': 'nonlinear'}, + id='dwidenoise2_demodulate', + ), + pytest.param('patch2self', Patch2Self, {}, id='patch2self_ignores_phase'), + ], +) +def test_denoising_wf_complex( + monkeypatch, + tmp_path, + nibs_dwi, + denoise_method, + interface, + expected_inputs, +): + """Denoise DWI data when phase data are available. + + Only the dwidenoise variants combine the magnitude and phase data into a + complex-valued series. ``patch2self`` ignores the phase data and denoises the + magnitude data alone. + """ + nodes, sink_dir = _run_denoising_wf( + monkeypatch, + tmp_path, + nibs_dwi, + denoise_method=denoise_method, + use_phase=True, + ) + + denoiser = nodes['denoiser'] + assert isinstance(denoiser.interface, interface) + for name, value in expected_inputs.items(): + assert getattr(denoiser.inputs, name) == value + + uses_complex = denoise_method.startswith('dwidenoise') + assert ('combine_complex' in nodes) is uses_complex + assert ('split_complex' in nodes) is uses_complex + if uses_complex: + complex_img = nb.load(nodes['combine_complex'].result.outputs.out_file) + assert np.issubdtype(complex_img.header.get_data_dtype(), np.complexfloating) + + _assert_denoiser_is_masked(nodes, denoise_method, nibs_dwi['dwi_file']) + _assert_denoising_outputs(nodes, sink_dir, nibs_dwi['dwi_file']) diff --git a/qsiprep/utils/misc.py b/qsiprep/utils/misc.py index 5c34ca1bb..3490e012d 100644 --- a/qsiprep/utils/misc.py +++ b/qsiprep/utils/misc.py @@ -85,6 +85,10 @@ def parse_denoise_method(spec): choices = _DWIDENOISE_ENUM_PARAMETERS[name] if value not in choices: raise ValueError(f'Invalid value for {name!r}: {value!r}; choose from {choices}') + if name == 'shape' and value == 'cuboid': + # dwidenoise2 rejects an odd -extent, but the denoising window is always + # rounded up to an odd number, so cuboid kernels can't be used yet. + raise ValueError("'shape:cuboid' is not supported yet; use 'shape:sphere' instead") parsed_value = value elif name == 'onepass': bool_values = {'true': True, 'false': False, '1': True, '0': False} diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index eab5f89ca..b4741e961 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -27,6 +27,7 @@ DWIDenoise, DWIDenoise2, MRDeGibbs, + MRTrixGradientTable, PolarToComplex, ) from ...interfaces.nilearn import MaskEPI, Merge @@ -418,6 +419,14 @@ def get_buffernode(): # Which steps to apply? denoise_method, dwidenoise_params = parse_denoise_method(config.workflow.denoise_method) + if denoise_method == 'dwidenoise2' and not use_phase: + demodulation = dwidenoise_params.get('demodulate', 'none') + if demodulation != 'none': + raise ValueError( + f'dwidenoise2 cannot apply {demodulation!r} phase demodulation to ' + 'magnitude-only data. Provide phase data or use "demodulate:none".' + ) + unringing_method = config.workflow.unringing_method do_denoise = denoise_method in ('patch2self', 'dwidenoise', 'dwidenoise2') do_unringing = config.workflow.unringing_method in ('mrdegibbs', 'rpg') @@ -427,6 +436,23 @@ def get_buffernode(): num_steps = sum(map(int, [do_denoise, do_unringing, do_biascorr, harmonize_b0s])) merge_confounds = pe.Node(niu.Merge(num_steps), name='merge_confounds') + # A single brain mask is shared by the denoising and bias correction steps. It is built + # from the raw series because denoising runs first, so the mask cannot be derived from + # the output of any earlier step. + # ``dwidenoise`` restricts the voxels it processes to this mask, whereas ``dwidenoise2`` + # has no -mask option, so there the mask only sets the contour drawn on the report. + mask_denoiser = denoise_method.startswith('dwidenoise') + if mask_denoiser or do_biascorr: + get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') + quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') + workflow.connect([ + (inputnode, get_b0s, [ + ('dwi_file', 'dwi_series'), + ('bval_file', 'bval_file'), + ]), + (get_b0s, quick_mask, [('b0_series', 'in_files')]), + ]) # fmt:skip + # Add the steps step_num = 1 # Merge inputs start at 1 last_step = '' @@ -459,7 +485,11 @@ def get_buffernode(): ) auto_str = 'n automatically-determined' - if denoise_method.startswith('dwidenoise') and use_phase: + # Only the dwidenoise variants can denoise complex-valued data. + # Any other method ignores the phase data and denoises the magnitude data alone. + denoise_complex = denoise_method.startswith('dwidenoise') and use_phase + + if denoise_complex: desc += ( 'Magnitude and phase DWI data were combined into a complex-valued file, ' 'then denoised using the Marchenko-Pastur PCA method implemented in ' @@ -578,13 +608,28 @@ def get_buffernode(): n_procs=omp_nthreads, ) - if denoise_method in ('dwidenoise2', 'patch2self'): + if denoise_method == 'patch2self': workflow.connect([(inputnode, denoiser, [('bval_file', 'bval_file')])]) + elif denoise_method == 'dwidenoise2': + # dwidenoise2 needs the gradient table to demean by shell. The standalone + # dwidenoise2 build misreads the two files given to its -fslgrad option, so + # supply the gradients as a single MRtrix-format table instead. + gradient_table = pe.Node(MRTrixGradientTable(), name='gradient_table') + workflow.connect([ + (inputnode, gradient_table, [ + ('bval_file', 'bval_file'), + ('bvec_file', 'bvec_file'), + ]), + (gradient_table, denoiser, [('gradient_file', 'grad_file')]), + ]) # fmt:skip - if denoise_method == 'dwidenoise2': - workflow.connect([(inputnode, denoiser, [('bvec_file', 'bvec_file')])]) + workflow.connect([ + (quick_mask, denoiser, [('out_mask', 'mask')]) + # The noise image is a derivative, so it always comes straight from the denoiser + (denoiser, outputnode, [('noise_image', 'noise_image')]), + ]) # fmt:skip - if (denoise_method != 'none') and not use_phase: + if not denoise_complex: workflow.connect([ (buffernodes[-2], denoiser, [('dwi_file', 'in_file')]), (denoiser, ds_report_denoising, [('out_report', 'in_file')]), @@ -666,19 +711,14 @@ def get_buffernode(): run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB, ) - get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') - quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') - # Add buffernode for bias-corrected DWI buffernodes.append(get_buffernode()) workflow.connect([ (buffernodes[-2], biascorr, [('dwi_file', 'in_file')]), - (buffernodes[-2], get_b0s, [('dwi_file', 'dwi_series')]), - (inputnode, get_b0s, [('bval_file', 'bval_file')]), - (get_b0s, quick_mask, [('b0_series', 'in_files')]), (quick_mask, biascorr, [('out_mask', 'mask')]), (biascorr, buffernodes[-1], [('out_file', 'dwi_file')]), + (biascorr, outputnode, [('bias_image', 'bias_image')]), (biascorr, ds_report_biascorr, [('out_report', 'in_file')]), (biascorr, merge_confounds, [('nmse_text', f'in{step_num}')]), (inputnode, biascorr, [ From 2f16b182119a3cb95033d4f4ab462ec344fbd620 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 5 Aug 2026 12:05:04 -0400 Subject: [PATCH 11/14] Update dwidenoise2 version. --- .circleci/continue_config.yml | 2 + .circleci/data_versions.txt | 1 + .circleci/get_data.sh | 48 +++++ Dockerfile | 25 ++- qsiprep/cli/parser.py | 44 +++-- qsiprep/data/boilerplate.bib | 138 +++++++++++-- qsiprep/interfaces/mrtrix.py | 86 ++++---- qsiprep/tests/conftest.py | 14 +- qsiprep/tests/test_interfaces_mrtrix.py | 35 ++-- qsiprep/tests/test_utils_misc.py | 164 +++++++++++++-- qsiprep/tests/test_workflows_merge.py | 120 ++++++----- qsiprep/tests/utils.py | 1 + qsiprep/utils/misc.py | 253 ++++++++++++++++++++---- qsiprep/workflows/dwi/merge.py | 240 +++++++++------------- 14 files changed, 817 insertions(+), 354 deletions(-) diff --git a/.circleci/continue_config.yml b/.circleci/continue_config.yml index 4c7fd3f3e..8469428bc 100644 --- a/.circleci/continue_config.yml +++ b/.circleci/continue_config.yml @@ -383,6 +383,7 @@ jobs: echo "export DSCSDSI_FMAP_URL=$DSCSDSI_FMAP_URL" >> "$BASH_ENV" echo "export MATERNAL_BRAIN_PROJECT_URL=$MATERNAL_BRAIN_PROJECT_URL" >> "$BASH_ENV" echo "export FORREST_GUMP_URL=$FORREST_GUMP_URL" >> "$BASH_ENV" + echo "export NIBS_URL=$NIBS_URL" >> "$BASH_ENV" - run: name: Download integration datasets command: | @@ -423,6 +424,7 @@ jobs: download_one DSCSDSI_fmap "$DSCSDSI_FMAP_URL" download_one maternal_brain_project "$MATERNAL_BRAIN_PROJECT_URL" download_one forrest_gump "$FORREST_GUMP_URL" + download_one nibs "$NIBS_URL" - save_cache: key: data-v2-{{ checksum ".circleci/data_versions.txt" }} paths: diff --git a/.circleci/data_versions.txt b/.circleci/data_versions.txt index d50f53354..6907b549c 100644 --- a/.circleci/data_versions.txt +++ b/.circleci/data_versions.txt @@ -14,3 +14,4 @@ DSDTI_FMAP_URL=https://upenn.box.com/shared/static/rxr6qbi6ezku9gw3esfpnvqlcxaw7 DSCSDSI_FMAP_URL=https://upenn.box.com/shared/static/l561psez1ojzi4p3a12eidaw9vbizwdc.gz MATERNAL_BRAIN_PROJECT_URL=https://upenn.box.com/shared/static/tkahg1ctipmfihvpa1gmibvcv0gb721h.xz FORREST_GUMP_URL=https://upenn.box.com/shared/static/qat58an322bzzyixrrsk7cmf52q3bepq.xz +NIBS_URL=https://upenn.box.com/shared/static/bkllff4ik51jy9ju6nben2r5zrq4a5me.xz diff --git a/.circleci/get_data.sh b/.circleci/get_data.sh index 7e1058ce8..bf5cb3a11 100644 --- a/.circleci/get_data.sh +++ b/.circleci/get_data.sh @@ -346,6 +346,42 @@ Contents: - data/singleshell_output/qsiprep/sub-PNC/figures/sub-PNC_t1_2_mni.svg - data/singleshell_output/qsiprep/sub-PNC.html + +nibs: +----- + +Downsampled NIBS data acquired with the HBCD protocol, with both magnitude and +phase parts, used to exercise the complex-valued denoising workflows. + +Contents: +^^^^^^^^^ + + - data/nibs/dataset_description.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_dwi.bval + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_dwi.bvec + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-mag_dwi.bval + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-mag_dwi.bvec + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-mag_dwi.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-mag_dwi.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-mag_sbref.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-mag_sbref.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-phase_dwi.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-phase_dwi.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-phase_sbref.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-AP_part-phase_sbref.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_dwi.bval + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_dwi.bvec + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-mag_dwi.bval + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-mag_dwi.bvec + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-mag_dwi.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-mag_dwi.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-mag_sbref.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-mag_sbref.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-phase_dwi.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-phase_dwi.nii.gz + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-phase_sbref.json + - data/nibs/sub-22449/ses-01/dwi/sub-22449_ses-01_dir-PA_part-phase_sbref.nii.gz + DOC @@ -440,6 +476,18 @@ get_bids_data() { rm tinytensors_epi.tar.xz fi + # Downsampled NIBS (complex-valued HBCD protocol) + # Unlike the other archives, this one has no top-level directory, so it is extracted + # into a directory of its own rather than into data/ directly. + if [[ ${DS} = nibs ]]; then + mkdir -p ${WORKDIR}/data/nibs + ${WGET} \ + -O nibs.tar.xz \ + "https://upenn.box.com/shared/static/bkllff4ik51jy9ju6nben2r5zrq4a5me.xz" + tar xvfJ nibs.tar.xz -C ${WORKDIR}/data/nibs + rm nibs.tar.xz + fi + # name: Get data for fieldmap tests if [[ ${DS} = fmaps ]]; then mkdir -p ${WORKDIR}/data/fmaptests diff --git a/Dockerfile b/Dockerfile index 1f4ddf3f5..ca60dc727 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,6 +1,7 @@ ARG BASE_IMAGE=pennlinc/qsiprep-base:20260804 -ARG DWIDENOISE2_COMMIT=892d5e8dd8f453ce1a561878f7dcb3998ae258ba -ARG MRTRIX3_DWIDENOISE2_COMMIT=fa6ee952913fbc1df79aeca745600852155f533f +ARG DWIDENOISE2_COMMIT=cd08ec1a0f5eb1dbc9962f80c20c2bb3428c4f93 +# MRtrix3 "dev" as at 2026-06-22, the commit dwidenoise2 is developed against +ARG MRTRIX3_DWIDENOISE2_COMMIT=b98b54e9ae8168eeb9af23322a07011d4754456d FROM buildpack-deps:bookworm AS dwidenoise2-build ARG DWIDENOISE2_COMMIT @@ -23,13 +24,21 @@ RUN git init . && \ WORKDIR /src/mrtrix3 RUN git clone --filter=blob:none --no-checkout https://github.com/MRtrix3/mrtrix3.git . && \ git checkout --detach ${MRTRIX3_DWIDENOISE2_COMMIT} +# dwidenoise2 has no external-module build, so its sources are dropped into the MRtrix3 +# tree before configuring. The per-command noise estimation schedules must land in +# share/mrtrix3//, where the built commands look for them relative to the +# executable; "copy-share-data" imports them into the build tree and has to be named +# explicitly, because building named executable targets does not run MRtrix3's ALL targets. RUN cp /src/dwidenoise2/cpp/cmd/dwidenoise2.cpp cpp/cmd/dwidenoise2.cpp && \ + cp /src/dwidenoise2/cpp/cmd/dwi2noise.cpp cpp/cmd/dwi2noise.cpp && \ cp -r /src/dwidenoise2/cpp/core/denoise cpp/core/denoise && \ + cp -r /src/dwidenoise2/share/dwidenoise2/. share/mrtrix3/ && \ cmake -B build -GNinja \ -DMRTRIX_BUILD_GUI=OFF \ + -DMRTRIX_ENABLE_GPU=OFF \ -DCMAKE_COMPILE_WARNING_AS_ERROR=ON \ --preset=release && \ - cmake --build build --target dwidenoise2 + cmake --build build --target dwidenoise2 dwi2noise copy-share-data FROM ghcr.io/prefix-dev/pixi:0.58.0 AS build RUN apt-get update && \ @@ -71,15 +80,23 @@ ENV HOME="/home/qsiprep" COPY --from=dwidenoise2-build \ /src/mrtrix3/build/bin/dwidenoise2 \ /opt/dwidenoise2/bin/dwidenoise2 +COPY --from=dwidenoise2-build \ + /src/mrtrix3/build/bin/dwi2noise \ + /opt/dwidenoise2/bin/dwi2noise COPY --from=dwidenoise2-build \ /src/mrtrix3/build/cpp/core/libmrtrix-core.so \ /opt/dwidenoise2/lib/libmrtrix-core.so +# The bundled schedules are found relative to the executable, at ../share/mrtrix3// +COPY --from=dwidenoise2-build \ + /src/mrtrix3/build/share/mrtrix3 \ + /opt/dwidenoise2/share/mrtrix3 COPY --from=dwidenoise2-build \ /src/dwidenoise2/LICENSE \ /opt/dwidenoise2/LICENSE ENV PATH="/opt/dwidenoise2/bin:$PATH" \ LD_LIBRARY_PATH="/opt/dwidenoise2/lib:$LD_LIBRARY_PATH" -RUN dwidenoise2 -version +RUN dwidenoise2 -version && \ + test -d /opt/dwidenoise2/share/mrtrix3/dwidenoise2 RUN chmod -R go=u $HOME WORKDIR /tmp diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index d53272b5b..020e98f46 100644 --- a/qsiprep/cli/parser.py +++ b/qsiprep/cli/parser.py @@ -393,11 +393,13 @@ def _bids_filter(value, parser): help=( 'Window size in voxels for image-based denoising: odd integer or "auto". ' 'Any non-"auto" value must be an odd, positive integer. ' - 'If using the "dwidenoise" denoising method, ' - 'the "auto" option will calculate a window size ' + 'This argument only applies to the "dwidenoise" denoising method, ' + 'where the "auto" option will calculate a window size ' 'based on the number of volumes according to the method described by the ' 'dwidenoise documentation. ' - 'If using the "patch2self" denoising method, this argument will not be used.' + 'It is not used by the "patch2self" or "dwidenoise2" methods: dwidenoise2 sizes ' + 'its patches per iteration from its multi-resolution schedule, which is selected ' + 'with "dwidenoise2;schedule:" instead.' ), ) g_conf.add_argument( @@ -675,6 +677,32 @@ def _bids_filter(value, parser): return parser +def check_denoise_window(denoise_method, dwi_denoise_window): + """Report a ``--dwi-denoise-window`` that the selected denoising method will ignore. + + Only ``dwidenoise`` takes a window size. Leaving the others to silently ignore it would + hide a request that never took effect. + """ + if dwi_denoise_window == 'auto': + # The default, so an unused value is not a sign that anything was misunderstood + return + + if denoise_method == 'patch2self': + config.loggers.cli.error( + 'The --dwi-denoise-window option is not used when --denoise-method=patch2self' + ) + elif denoise_method == 'dwidenoise2': + config.loggers.cli.warning( + 'The --dwi-denoise-window option is not used when --denoise-method=dwidenoise2. ' + 'dwidenoise2 sizes its patches per iteration from its multi-resolution schedule, ' + 'which can be selected with "dwidenoise2;schedule:" instead.' + ) + elif denoise_method == 'none': + config.loggers.cli.warning( + 'The --dwi-denoise-window option is not used when --denoise-method=none' + ) + + def parse_args(args=None, namespace=None): """Parse args and run further checks on the command line.""" import logging @@ -767,15 +795,7 @@ def parse_args(args=None, namespace=None): # Validate the tricky options here denoise_method, _ = parse_denoise_method(config.workflow.denoise_method) - if config.workflow.dwi_denoise_window != 'auto': - if denoise_method == 'patch2self': - config.loggers.cli.error( - 'The --dwi-denoise-window option is not used when --denoise-method=patch2self' - ) - elif denoise_method == 'none': - config.loggers.cli.warning( - 'The --dwi-denoise-window option is not used when --denoise-method=none' - ) + check_denoise_window(denoise_method, config.workflow.dwi_denoise_window) bids_dir = config.execution.bids_dir output_dir = config.execution.output_dir diff --git a/qsiprep/data/boilerplate.bib b/qsiprep/data/boilerplate.bib index 026d6f278..1bce05646 100644 --- a/qsiprep/data/boilerplate.bib +++ b/qsiprep/data/boilerplate.bib @@ -124,6 +124,14 @@ @article{dwidenoise2 pmid = {26599599}, } +@misc{dwidenoise2software, + title = {dwidenoise2}, + author = {Smith, Robert E.}, + year = {2026}, + publisher = {Zenodo}, + doi = {10.5281/zenodo.21601472}, +} + @article{eddyrepol, title = {Incorporating outlier detection and replacement into a non-parametric framework for movement and distortion correction of diffusion MR images}, author = {Andersson, Jesper L. R and Graham, Mark S. and Zsoldos, Enikő and Sotiropoulos, Stamatios N.}, @@ -205,6 +213,16 @@ @article{flirt pmid = {11516708}, } +@inproceedings{foi2011, + title = {Noise estimation and removal in MR imaging: The variance-stabilization approach}, + shorttitle = {Noise estimation and removal in MR imaging}, + author = {Foi, Alessandro}, + year = {2011}, + booktitle = {2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro}, + pages = {1809--1814}, + doi = {10.1109/ISBI.2011.5872758}, +} + @article{fs_reconall, title = {Cortical Surface-Based Analysis}, shorttitle = {Cortical surface-based analysis}, @@ -257,6 +275,17 @@ @article{fsllsr pmid = {14568458}, } +@article{gavish2014, + title = {The Optimal Hard Threshold for Singular Values is $4/\sqrt{3}$}, + author = {Gavish, Matan and Donoho, David L.}, + year = {2014}, + journal = {IEEE Transactions on Information Theory}, + volume = {60}, + number = {8}, + pages = {5040--5053}, + doi = {10.1109/TIT.2014.2323359}, +} + @article{hcppipelines, title = {The minimal preprocessing pipelines for the Human Connectome Project}, author = {Glasser, Matthew F. and Sotiropoulos, Stamatios N. and Wilson, J. Anthony and Coalson, Timothy S. and Fischl, Bruce and Andersson, Jesper L. and Xu, Junqian and Jbabdi, Saad and Webster, Matthew and Polimeni, Jonathan R. and Van Essen, David C. and Jenkinson, Mark}, @@ -270,6 +299,18 @@ @article{hcppipelines pmid = {23668970}, } +@article{koay2006, + title = {Analytically exact correction scheme for signal extraction from noisy magnitude MR signals}, + author = {Koay, Cheng Guan and Basser, Peter J.}, + year = {2006}, + journal = {Journal of Magnetic Resonance}, + volume = {179}, + number = {2}, + pages = {317--322}, + doi = {10.1016/j.jmr.2006.01.016}, + pmid = {16488635}, +} + @article{lanczos, title = {Evaluation of Noisy Data}, author = {Lanczos, C.}, @@ -281,6 +322,31 @@ @article{lanczos doi = {10.1137/0701007}, } +@article{ma2020, + title = {Denoise magnitude diffusion magnetic resonance images via variance-stabilizing transformation and optimal singular-value manipulation}, + author = {Ma, Xiaodong and U{\u{g}}urbil, Kamil and Wu, Xiaoping}, + year = {2020}, + journal = {NeuroImage}, + volume = {215}, + pages = {116852}, + doi = {10.1016/j.neuroimage.2020.116852}, + pmcid = {PMC7292714}, + pmid = {32407993}, +} + +@article{manjon2013, + title = {Diffusion Weighted Image Denoising Using Overcomplete Local PCA}, + author = {Manj{\'o}n, Jos{\'e} V. and Coup{\'e}, Pierrick and Concha, Luis and Buades, Antonio and Collins, D. Louis and Robles, Montserrat}, + year = {2013}, + journal = {PLoS ONE}, + volume = {8}, + number = {9}, + pages = {e73021}, + doi = {10.1371/journal.pone.0073021}, + pmcid = {PMC3765487}, + pmid = {24019889}, +} + @article{mcflirt, title = {Improved Optimization for the Robust and Accurate Linear Registration and Motion Correction of Brain Images}, author = {Jenkinson, Mark and Bannister, Peter and Brady, Michael and Smith, Stephen}, @@ -406,6 +472,18 @@ @misc{nipype2 year = {2025}, } +@article{olesen2023, + title = {Tensor denoising of multidimensional MRI data}, + author = {Olesen, Jonas L. and Ianus, Andrada and {\O}stergaard, Leif and Shemesh, Noam and Jespersen, Sune N.}, + year = {2023}, + journal = {Magnetic Resonance in Medicine}, + volume = {89}, + number = {3}, + pages = {1160--1172}, + doi = {10.1002/mrm.29478}, + pmid = {36372982}, +} + @article{patch2self, title = {Patch2Self: Denoising diffusion MRI with self-supervised learning​}, shorttitle = {Patch2Self}, @@ -415,6 +493,17 @@ @article{patch2self volume = {33}, } +@article{patron2024, + title = {Denoising diffusion MRI: Considerations and implications for analysis}, + shorttitle = {Denoising diffusion MRI}, + author = {Manzano Patron, Jose Pedro and Moeller, Steen and Andersson, Jesper L. R. and U{\u{g}}urbil, Kamil and Yacoub, Essa and Sotiropoulos, Stamatios N.}, + year = {2024}, + journal = {Imaging Neuroscience}, + volume = {2}, + pages = {1--29}, + doi = {10.1162/imag_a_00060}, +} + @article{pfgibbs, title = {Removal of partial Fourier‐induced Gibbs (RPG) ringing artifacts in MRI}, author = {Lee, Hong‐Hsi and Novikov, Dmitry S. and Fieremans, Els}, @@ -428,6 +517,17 @@ @article{pfgibbs pmid = {34227142}, } +@article{pizzolato2020, + title = {Adaptive phase correction of diffusion-weighted images}, + author = {Pizzolato, Marco and Gilbert, Guillaume and Thiran, Jean-Philippe and Descoteaux, Maxime and Deriche, Rachid}, + year = {2020}, + journal = {NeuroImage}, + volume = {206}, + pages = {116274}, + doi = {10.1016/j.neuroimage.2019.116274}, + pmid = {31629828}, +} + @article{pncprocessing, title = {Neuroimaging of the Philadelphia Neurodevelopmental Cohort}, author = {Satterthwaite, Theodore D. and Elliott, Mark A. and Ruparel, Kosha and Loughead, James and Prabhakaran, Karthik and Calkins, Monica E. and Hopson, Ryan and Jackson, Chad and Keefe, Jack and Riley, Marisa and Mentch, Frank D. and Sleiman, Patrick and Verma, Ragini and Davatzikos, Christos and Hakonarson, Hakon and Gur, Ruben C. and Gur, Raquel E.}, @@ -440,6 +540,18 @@ @article{pncprocessing pmid = {23921101}, } +@article{power2017simple, + title = {A simple but useful way to assess fMRI scan qualities}, + author = {Power, Jonathan D.}, + year = {2017}, + journal = {NeuroImage}, + volume = {154}, + pages = {150--158}, + doi = {10.1016/j.neuroimage.2016.08.009}, + pmcid = {PMC5296400}, + pmid = {27510328}, +} + @article{power_fd_dvars, title = {Methods to detect, characterize, and remove motion artifact in resting state fMRI}, author = {Power, Jonathan D. and Mitra, Anish and Laumann, Timothy O. and Snyder, Abraham Z. and Schlaggar, Bradley L. and Petersen, Steven E.}, @@ -453,18 +565,6 @@ @article{power_fd_dvars pmid = {23994314}, } -@article{power2017simple, - title = {A simple but useful way to assess fMRI scan qualities}, - author = {Power, Jonathan D.}, - year = {2017}, - journal = {NeuroImage}, - volume = {154}, - pages = {150--158}, - doi = {10.1016/j.neuroimage.2016.08.009}, - pmcid = {PMC5296400}, - pmid = {27510328}, -} - @article{synthseg1, title = {SynthSeg: Segmentation of brain MRI scans of any contrast and resolution without retraining}, shorttitle = {SynthSeg}, @@ -525,3 +625,17 @@ @article{tortoisev4 journal = {Imaging Neuroscience}, doi = {10.1162/IMAG.a.948}, } + + +@article{zhu2022, + title = {Denoise Functional Magnetic Resonance Imaging With Random Matrix Theory Based Principal Component Analysis}, + author = {Zhu, Wenchao and Ma, Xiaodong and Zhu, Xiao-Hong and U{\u{g}}urbil, Kamil and Chen, Wei and Wu, Xiaoping}, + year = {2022}, + journal = {IEEE Transactions on Biomedical Engineering}, + volume = {69}, + number = {11}, + pages = {3377--3388}, + doi = {10.1109/TBME.2022.3168592}, + pmcid = {PMC9613202}, + pmid = {35439123}, +} diff --git a/qsiprep/interfaces/mrtrix.py b/qsiprep/interfaces/mrtrix.py index 9c3763794..7205f19ac 100644 --- a/qsiprep/interfaces/mrtrix.py +++ b/qsiprep/interfaces/mrtrix.py @@ -184,7 +184,12 @@ def _get_plotting_images(self): class DWIDenoise2InputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): in_file = File(exists=True, argstr='%s', position=-2, mandatory=True, desc='input DWI image') mask = File(exists=True, desc='mask image') - onepass = traits.Bool(argstr='-onepass', desc='estimate noise and denoise in one pass') + # The sliding-window kernel and the subsampling factor are properties of the + # multi-resolution schedule rather than command-line options + schedule = traits.Str( + argstr='-schedule %s', + desc='name of a bundled noise estimation schedule, or a path to a schedule file', + ) datatype = traits.Enum( 'float32', 'float64', @@ -198,10 +203,11 @@ class DWIDenoise2InputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): desc='patch decomposition method', ) estimator = traits.Enum( - 'Exp1', - 'Exp2', - 'Med', - 'MRM2023', + 'exp1', + 'exp2', + 'med', + 'mrm2023', + 'tbme2022', argstr='-estimator %s', desc='noise level estimator', ) @@ -215,39 +221,11 @@ class DWIDenoise2InputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): fixed_rank = traits.Int( argstr='-fixed_rank %d', xor=('noise_in',), desc='fixed input signal rank' ) - shape = traits.Enum( - 'sphere', - 'cuboid', - argstr='-shape %s', - desc='sliding spatial window shape', - ) - radius = traits.Float( - argstr='-radius %g', - xor=('extent',), - desc='absolute spherical kernel radius in mm', - ) - aspect_ratio = traits.Float( - argstr='-aspect_ratio %g', - desc='ratio of kernel voxels to input volumes', - ) - minvoxels = traits.Int(argstr='-minvoxels %d', desc='minimum voxels in a spherical kernel') - extent = traits.Either( - traits.Int, - traits.Tuple(traits.Int, traits.Int, traits.Int), - argstr='-extent %s', - xor=('radius',), - desc='cuboid window size as one integer or a triplet', - ) - subsample = traits.Either( - traits.Int, - traits.Tuple(traits.Int, traits.Int, traits.Int), - argstr='-subsample %s', - desc='PCA kernel subsampling factor as one integer or a triplet', - ) demodulate = traits.Enum( 'none', 'linear', - 'nonlinear', + 'hann', + 'apc', argstr='-demodulate %s', desc='phase demodulation mode', ) @@ -263,10 +241,28 @@ class DWIDenoise2InputSpec(MRTrix3BaseInputSpec, SeriesPreprocReportInputSpec): argstr='-demean %s', desc='demeaning method before PCA', ) - vst = File( - exists=True, - argstr='-vst %s', - desc='noise map for variance-stabilising transformation', + noise_dof = traits.Int( + argstr='-noise_dof %d', + desc='receive channels combined by sum-of-squares reconstruction of magnitude data', + ) + vst_method = traits.Enum( + 'none', + 'linear', + 'foi', + 'koay', + 'mom', + argstr='-vst_method %s', + desc='variance-stabilising transform applied before PCA', + ) + preserve_noise_bias = traits.Bool( + argstr='-preserve_noise_bias', + desc='retain the noise-floor bias in magnitude output instead of removing it', + ) + debias_anchor = traits.Enum( + 'sample', + 'group_mean', + argstr='-debias_anchor %s', + desc='operating point at which the variance-stabilising inverse is evaluated', ) preconditioned_input = File( argstr='-preconditioned_input %s', @@ -413,22 +409,12 @@ class DWIDenoise2(SeriesPreprocReport, MRTrix3Base): output_spec = DWIDenoise2OutputSpec def _format_arg(self, name, spec, value): - if name in ('extent', 'subsample') and not isinstance(value, int): - value = ','.join(str(item) for item in value) - elif name == 'bvec_file': + if name == 'bvec_file': # -fslgrad takes both files, so format them here rather than passing a tuple # to a File trait, which nipype would try to shell-quote as a single value. return spec.argstr % (value, self.inputs.bval_file) return super()._format_arg(name, spec, value) - def _parse_inputs(self, skip=None): - shape = self.inputs.shape if isdefined(self.inputs.shape) else 'sphere' - if shape == 'sphere' and isdefined(self.inputs.extent): - raise ValueError("'extent' cannot be used when 'shape' is 'sphere'") - if shape == 'cuboid' and isdefined(self.inputs.radius): - raise ValueError("'radius' cannot be used when 'shape' is 'cuboid'") - return super()._parse_inputs(skip=skip) - def _get_plotting_images(self): input_dwi = load_img(self.inputs.in_file) outputs = self._list_outputs() diff --git a/qsiprep/tests/conftest.py b/qsiprep/tests/conftest.py index fcf41d88d..1f68eed67 100644 --- a/qsiprep/tests/conftest.py +++ b/qsiprep/tests/conftest.py @@ -41,13 +41,13 @@ def datasets(data_dir): """Locate downloaded datasets.""" dsets = {} dsets['forrest_gump'] = os.path.join(data_dir, 'forrest_gump') - dsets['nibs-ci'] = os.path.join(data_dir, 'nibs-ci') + dsets['nibs'] = os.path.join(data_dir, 'nibs') return dsets @pytest.fixture(scope='session') def nibs_dwi(data_dir): - """Locate the nibs-ci DWI series used to test the denoising workflow. + """Locate the nibs DWI series used to test the denoising workflow. The series is small (48x48x29x76) and has both magnitude and phase data, so it can exercise the complex-valued denoising paths without a long runtime. @@ -58,17 +58,17 @@ def nibs_dwi(data_dir): if not data_dir: pytest.skip('--data_dir was not provided') - dwi_dir = Path(data_dir) / 'nibs-ci' / 'sub-22449' / 'ses-01' / 'dwi' - stem = 'sub-22449_ses-01_acq-HBCD75_rec-norm_dir-AP_run-01' + dwi_dir = Path(data_dir) / 'nibs' / 'sub-22449' / 'ses-01' / 'dwi' + stem = 'sub-22449_ses-01_dir-AP' files = { 'dwi_file': dwi_dir / f'{stem}_part-mag_dwi.nii.gz', 'phase_file': dwi_dir / f'{stem}_part-phase_dwi.nii.gz', - 'bval_file': dwi_dir / f'{stem}_part-mag_dwi.bval', - 'bvec_file': dwi_dir / f'{stem}_part-mag_dwi.bvec', + 'bval_file': dwi_dir / f'{stem}_dwi.bval', + 'bvec_file': dwi_dir / f'{stem}_dwi.bvec', 'json_file': dwi_dir / f'{stem}_part-mag_dwi.json', } missing = sorted(str(f) for f in files.values() if not f.is_file()) if missing: - pytest.skip(f'nibs-ci dataset is unavailable; missing {missing}') + pytest.skip(f'nibs dataset is unavailable; missing {missing}') return {key: str(value) for key, value in files.items()} diff --git a/qsiprep/tests/test_interfaces_mrtrix.py b/qsiprep/tests/test_interfaces_mrtrix.py index 1f3f67bc4..d06471a26 100644 --- a/qsiprep/tests/test_interfaces_mrtrix.py +++ b/qsiprep/tests/test_interfaces_mrtrix.py @@ -4,6 +4,7 @@ import nibabel as nb import pytest +from traits.trait_errors import TraitError from qsiprep.interfaces import mrtrix @@ -45,10 +46,6 @@ def test_dwidenoise2(datasets, tmp_path_factory): in_img = nb.load(in_file) interface = mrtrix.DWIDenoise2( - shape='sphere', - radius=3, - onepass=True, - subsample=1, in_file=in_file, nthreads=1, ) @@ -68,34 +65,26 @@ def test_dwidenoise2(datasets, tmp_path_factory): @pytest.mark.parametrize( - ('shape', 'kernel_option', 'error'), - [ - ('sphere', {'extent': (5, 5, 5)}, "'extent' cannot be used"), - ('cuboid', {'radius': 2.5}, "'radius' cannot be used"), - ], + 'kernel_option', + ['shape', 'radius', 'extent', 'aspect_ratio', 'minvoxels', 'subsample', 'onepass'], ) -def test_dwidenoise2_kernel_shape_validation(tmp_path, shape, kernel_option, error): - """Reject kernel options that do not apply to the selected shape.""" +def test_dwidenoise2_has_no_kernel_options(tmp_path, kernel_option): + """The kernel and subsampling come from the schedule, not from command-line options.""" in_file = tmp_path / 'dwi.nii.gz' in_file.touch() - interface = mrtrix.DWIDenoise2(in_file=in_file, shape=shape, **kernel_option) - with pytest.raises(ValueError, match=error): - _ = interface.cmdline + with pytest.raises(TraitError, match='undefined'): + mrtrix.DWIDenoise2(in_file=in_file, **{kernel_option: 1}) -def test_dwidenoise2_kernel_options_are_mutually_exclusive(tmp_path): - """Reject simultaneous spherical and cuboid kernel size options.""" +def test_dwidenoise2_passes_schedule(tmp_path): + """Select a bundled noise estimation schedule by name.""" in_file = tmp_path / 'dwi.nii.gz' in_file.touch() - with pytest.raises(OSError, match='mutually exclusive'): - mrtrix.DWIDenoise2( - in_file=in_file, - shape='sphere', - radius=2.5, - extent=(5, 5, 5), - ) + interface = mrtrix.DWIDenoise2(in_file=in_file, schedule='vlarge') + + assert '-schedule vlarge' in interface.cmdline def test_dwidenoise2_formats_fslgrad(tmp_path): diff --git a/qsiprep/tests/test_utils_misc.py b/qsiprep/tests/test_utils_misc.py index 73feb2bf8..208e7b9b0 100644 --- a/qsiprep/tests/test_utils_misc.py +++ b/qsiprep/tests/test_utils_misc.py @@ -6,7 +6,7 @@ import pytest from qsiprep.cli.parser import _build_parser -from qsiprep.utils.misc import parse_denoise_method, safe_unit_vector +from qsiprep.utils.misc import describe_dwidenoise2, parse_denoise_method, safe_unit_vector def test_safe_unit_vector_zero_magnitude_substitutes_x_axis(): @@ -51,17 +51,19 @@ def test_angle_between_finite_for_zero_vector(): def test_parse_denoise_method_parameters(): method, parameters = parse_denoise_method( - 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd;' - 'onepass:true;radius:2.5;subsample:2,2,2' + 'dwidenoise2;demodulate:hann;decomposition:bdcsvd;' + 'preserve_noise_bias:true;noise_dof:8;aggregator_fwhm:2.5;schedule:vlarge', + use_phase=True, ) assert method == 'dwidenoise2' assert parameters == { - 'demodulate': 'nonlinear', + 'demodulate': 'hann', 'decomposition': 'bdcsvd', - 'onepass': True, - 'radius': 2.5, - 'subsample': (2, 2, 2), + 'preserve_noise_bias': True, + 'noise_dof': 8, + 'aggregator_fwhm': 2.5, + 'schedule': 'vlarge', } @@ -75,27 +77,42 @@ def test_parse_denoise_method_parameters(): 'dwidenoise2;unknown:value', 'dwidenoise2;decomposition:bdcsvd;decomposition:selfadjoint', 'dwidenoise2;decomposition:invalid', - 'dwidenoise2;onepass:maybe', + 'dwidenoise2;preserve_noise_bias:maybe', + # The kernel and subsampling are set by the schedule, not by command-line options 'dwidenoise2;extent:1,2', + 'dwidenoise2;shape:sphere', + 'dwidenoise2;radius:2.5', + 'dwidenoise2;subsample:2', + 'dwidenoise2;onepass:true', + # dwidenoise2 renamed its demodulation and estimator choices + 'dwidenoise2;demodulate:nonlinear', + 'dwidenoise2;estimator:MRM2023', ], ) def test_parse_denoise_method_rejects_invalid_specs(spec): with pytest.raises(ValueError, match='.'): - parse_denoise_method(spec) + parse_denoise_method(spec, use_phase=True) -def test_parse_denoise_method_rejects_cuboid_shape(): - """Reject cuboid kernels, which need an even extent that QSIPrep never produces.""" - with pytest.raises(ValueError, match='not supported yet'): - parse_denoise_method('dwidenoise2;shape:cuboid') +@pytest.mark.parametrize('demodulate', ['linear', 'hann', 'apc']) +def test_parse_denoise_method_rejects_demodulation_without_phase(demodulate): + """Reject phase demodulation of magnitude-only data, which dwidenoise2 cannot do.""" + spec = f'dwidenoise2;demodulate:{demodulate}' + with pytest.raises(ValueError, match='magnitude-only data'): + parse_denoise_method(spec, use_phase=False) - method, parameters = parse_denoise_method('dwidenoise2;shape:sphere') - assert method == 'dwidenoise2' - assert parameters == {'shape': 'sphere'} + assert parse_denoise_method(spec, use_phase=True) == ( + 'dwidenoise2', + {'demodulate': demodulate}, + ) + + # The CLI validates the specification before it knows whether phase data exist, so an + # unknown phase state skips the check rather than guessing + assert parse_denoise_method(spec) == ('dwidenoise2', {'demodulate': demodulate}) def test_denoise_method_cli_parameter(tmp_path): - spec = 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd' + spec = 'dwidenoise2;demodulate:apc;decomposition:bdcsvd' opts = _build_parser().parse_args( [ str(tmp_path), @@ -124,3 +141,116 @@ def test_denoise_method_cli_rejects_invalid_parameter(tmp_path): 'dwidenoise;decomposition:invalid', ] ) + + +def test_describe_dwidenoise2_covers_defaults(): + """Describe the methods that run by default, not only the requested parameters.""" + description = describe_dwidenoise2({}, complex_data=False) + + # The software, MP-PCA and the noise mapping paper are always applicable + for citation in ( + '@dwidenoise2software', + '@dwidenoise1', + '@dwidenoise2', + '@cordero2019complex', + ): + assert citation in description + + # ...as are the defaults: the mrm2023 estimator, Gaussian aggregation over overlapping + # patches, and the nonlinear variance-stabilizing transform magnitude data require + assert '@olesen2023' in description + assert '@manjon2013' in description + assert '@foi2011' in description + assert '@ma2020' in description + + # Nothing that did not run should be cited + for citation in ('@pizzolato2020', '@patron2024', '@gavish2014', '@zhu2022', '@koay2006'): + assert citation not in description + + +def test_describe_dwidenoise2_demodulation_is_complex_only(): + """Only describe phase demodulation when there are phase data to demodulate.""" + parameters = {'demodulate': 'apc'} + + assert '@pizzolato2020' in describe_dwidenoise2(parameters, complex_data=True) + assert '@pizzolato2020' not in describe_dwidenoise2(parameters, complex_data=False) + + # Complex data are Gaussian, so they need no nonlinear variance-stabilizing transform + # and carry no noise-floor bias + complex_description = describe_dwidenoise2(parameters, complex_data=True) + assert '@foi2011' not in complex_description + assert 'noise-floor bias' not in complex_description + + +@pytest.mark.parametrize( + ('parameters', 'expected', 'unexpected'), + [ + ({'demodulate': 'hann'}, '@patron2024', '@pizzolato2020'), + ({'demodulate': 'linear'}, '@cordero2019complex', '@pizzolato2020'), + ({'estimator': 'tbme2022'}, '@zhu2022', '@olesen2023'), + ({'estimator': 'med'}, '@gavish2014', '@olesen2023'), + ({'aggregator': 'exclusive'}, 'solely from the patch', '@manjon2013'), + ], +) +def test_describe_dwidenoise2_conditional_citations(parameters, expected, unexpected): + """Follow the conditions dwidenoise2 attaches to each citation in its own help.""" + description = describe_dwidenoise2(parameters, complex_data=True) + + assert expected in description + assert unexpected not in description + + +def test_describe_dwidenoise2_filter_follows_fixed_rank(): + """Describe hard truncation when the rank is given rather than estimated.""" + description = describe_dwidenoise2({'fixed_rank': 12}, complex_data=True) + + assert 'hard truncation' in description + assert 'signal rank was fixed at 12' in description + # The rank was not estimated, so no estimator applies + assert '@olesen2023' not in description + + +@pytest.mark.parametrize( + ('denoise_method', 'window', 'expected'), + [ + # dwidenoise2 has no kernel options at all, so a requested window silently does nothing + ('dwidenoise2', 5, 'not used when --denoise-method=dwidenoise2'), + ('none', 5, 'not used when --denoise-method=none'), + # dwidenoise is the only method that takes a window + ('dwidenoise', 5, None), + # 'auto' is the default, so an unused value is not a sign of a misunderstanding + ('dwidenoise2', 'auto', None), + ('patch2self', 'auto', None), + ], +) +def test_check_denoise_window_warns_when_unused(caplog, denoise_method, window, expected): + """Warn when --dwi-denoise-window cannot affect the selected denoising method.""" + from qsiprep.cli.parser import check_denoise_window + + with caplog.at_level(logging.WARNING, logger='cli'): + check_denoise_window(denoise_method, window) + + messages = ' '.join(record.message for record in caplog.records) + if expected is None: + assert not messages + else: + assert expected in messages + + +def test_check_denoise_window_errors_for_patch2self(caplog): + """patch2self never had a window, so an explicit one is reported as an error.""" + from qsiprep.cli.parser import check_denoise_window + + with caplog.at_level(logging.ERROR, logger='cli'): + check_denoise_window('patch2self', 5) + + assert any(record.levelname == 'ERROR' for record in caplog.records) + + +def test_denoise_window_help_mentions_dwidenoise2(): + """Say in the help text that dwidenoise2 ignores the window.""" + parser = _build_parser() + action = next(a for a in parser._actions if '--dwi-denoise-window' in a.option_strings) + + assert 'dwidenoise2' in action.help + assert 'schedule' in action.help diff --git a/qsiprep/tests/test_workflows_merge.py b/qsiprep/tests/test_workflows_merge.py index 854411f0f..e6ff69c62 100644 --- a/qsiprep/tests/test_workflows_merge.py +++ b/qsiprep/tests/test_workflows_merge.py @@ -43,18 +43,9 @@ def test_dwidenoise_workflow_uses_dwidenoise(monkeypatch, use_phase): assert denoiser.inputs.nthreads == 1 -@pytest.mark.parametrize( - ('denoise_method', 'kernel_input'), - [ - # dwidenoise takes a cuboid extent, while dwidenoise2 defaults to a spherical - # kernel, which is sized with a radius instead - ('dwidenoise', {'extent': (5, 5, 5)}), - ('dwidenoise2', {'radius': 5.0}), - ], -) -def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method, kernel_input): - """Resolve the default ``auto`` window size for every dwidenoise variant.""" - monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) +def test_dwidenoise_workflow_resolves_auto_window(monkeypatch): + """Resolve the default ``auto`` window size into a cuboid extent for dwidenoise.""" + monkeypatch.setattr(config.workflow, 'denoise_method', 'dwidenoise') monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 'auto') monkeypatch.setattr(config.workflow, 'unringing_method', 'none') monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) @@ -72,8 +63,36 @@ def test_dwidenoise_workflow_resolves_auto_window(monkeypatch, denoise_method, k denoiser = workflow.get_node('denoiser') # cbrt(30) rounded up to the closest odd integer - for name, value in kernel_input.items(): - assert getattr(denoiser.inputs, name) == value + assert denoiser.inputs.extent == (5, 5, 5) + + +def test_dwidenoise2_workflow_ignores_denoise_window(monkeypatch): + """Leave the kernel to dwidenoise2's schedule rather than the requested window. + + dwidenoise2 sizes its patches per iteration from its multi-resolution schedule and + exposes no kernel options, so ``--dwi-denoise-window`` cannot apply to it. + """ + monkeypatch.setattr(config.workflow, 'denoise_method', 'dwidenoise2') + monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) + monkeypatch.setattr(config.workflow, 'unringing_method', 'none') + monkeypatch.setattr(config.workflow, 'no_b0_harmonization', True) + monkeypatch.setattr(config.workflow, 'b0_threshold', 100) + monkeypatch.setattr(config.nipype, 'omp_nthreads', 1) + + workflow = init_dwi_denoising_wf( + source_file='sub-01_dwi.nii.gz', + partial_fourier=1.0, + phase_encoding_direction='j', + n_volumes=30, + use_phase=False, + do_biascorr=False, + ) + denoiser = workflow.get_node('denoiser') + + for removed in ('shape', 'radius', 'extent', 'subsample'): + assert not hasattr(denoiser.inputs, removed) + # No schedule is requested either, so dwidenoise2 uses its bundled default + assert not isdefined(denoiser.inputs.schedule) def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): @@ -81,7 +100,7 @@ def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): monkeypatch.setattr( config.workflow, 'denoise_method', - 'dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd', + 'dwidenoise2;demodulate:hann;decomposition:bdcsvd', ) monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) monkeypatch.setattr(config.workflow, 'unringing_method', 'none') @@ -100,14 +119,14 @@ def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): ) denoiser = workflow.get_node('denoiser') - assert denoiser.inputs.demodulate == 'nonlinear' + assert denoiser.inputs.demodulate == 'hann' assert denoiser.inputs.decomposition == 'bdcsvd' # Parameters that weren't requested are left at the dwidenoise2 defaults - assert not isdefined(denoiser.inputs.onepass) - assert not isdefined(denoiser.inputs.subsample) + assert not isdefined(denoiser.inputs.estimator) + assert not isdefined(denoiser.inputs.schedule) -@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2']) +@pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2', 'patch2self']) def test_denoising_wf_builds_one_mask_for_denoising_and_biascorr(monkeypatch, denoise_method): """Build the brain mask once and hand it to both the denoiser and bias correction.""" monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) @@ -147,7 +166,7 @@ def test_denoising_wf_builds_one_mask_for_denoising_and_biascorr(monkeypatch, de } -@pytest.mark.parametrize('demodulate', ['linear', 'nonlinear']) +@pytest.mark.parametrize('demodulate', ['linear', 'hann', 'apc']) def test_dwidenoise2_rejects_demodulation_without_phase(monkeypatch, demodulate): """Reject phase demodulation unless phase data are available. @@ -185,7 +204,7 @@ def _run_denoising_wf( use_phase, dwi_denoise_window='auto', ): - """Build and execute a denoising workflow on the nibs-ci DWI series. + """Build and execute a denoising workflow on the nibs DWI series. Unringing, bias correction and b=0 harmonization are all disabled so that only the denoising step is exercised. @@ -255,14 +274,12 @@ def _sink_output(sink_dir, field): return matches[0] -def _assert_denoiser_is_masked(nodes, denoise_method, raw_file): - """Check that the dwidenoise variants are handed a brain mask built from the raw data.""" - uses_mask = denoise_method.startswith('dwidenoise') - assert ('quick_mask' in nodes) is uses_mask - if not uses_mask: - assert not isdefined(nodes['denoiser'].inputs.mask) - return +def _assert_denoiser_is_masked(nodes, raw_file): + """Check that the denoiser is handed a brain mask built from the raw data. + Every method gets the mask. ``dwidenoise`` restricts the voxels it processes to it, + while ``dwidenoise2`` and ``patch2self`` use it only for the report contour. + """ mask_file = nodes['quick_mask'].result.outputs.out_mask assert nodes['denoiser'].inputs.mask == mask_file @@ -312,41 +329,38 @@ def _assert_denoising_outputs(nodes, sink_dir, raw_file): pytest.param( 'dwidenoise', 'auto', mrtrix.DWIDenoise, {'extent': (5, 5, 5)}, id='dwidenoise_auto' ), + # Every option is left at its default, so the bundled 'default' schedule sizes the + # kernel and the mrm2023 estimator is used + pytest.param('dwidenoise2', 'auto', mrtrix.DWIDenoise2, {}, id='dwidenoise2_default'), pytest.param( - 'dwidenoise2', + 'dwidenoise2;decomposition:selfadjoint', 'auto', mrtrix.DWIDenoise2, - {'shape': 'sphere', 'radius': 5.0}, - id='dwidenoise2_sphere', + {'decomposition': 'selfadjoint'}, + id='dwidenoise2_selfadjoint', ), pytest.param( - 'dwidenoise2;decomposition:bdcsvd', + 'dwidenoise2;filter_method:optthresh', 'auto', mrtrix.DWIDenoise2, - {'decomposition': 'bdcsvd'}, - id='dwidenoise2_bdcsvd', + {'filter_method': 'optthresh'}, + id='dwidenoise2_optthresh', ), pytest.param( - 'dwidenoise2;filter_method:optthresh', + 'dwidenoise2;estimator:exp2', 'auto', mrtrix.DWIDenoise2, - {'filter_method': 'optthresh'}, - id='dwidenoise2_optthresh', + {'estimator': 'exp2'}, + id='dwidenoise2_exp2', ), + # A named schedule only resolves if the bundled schedules were installed alongside + # the executable, so this also covers the container build pytest.param( - 'dwidenoise2;estimator:MRM2023', + 'dwidenoise2;schedule:legacy', 'auto', mrtrix.DWIDenoise2, - {'estimator': 'MRM2023'}, - id='dwidenoise2_mrm2023', - marks=pytest.mark.xfail( - strict=True, - reason=( - 'the MRM2023 estimator returns a noise map with negative values ' - '(~41% of voxels, down to -15.7 on this series), while every other ' - 'estimator stays positive' - ), - ), + {'schedule': 'legacy'}, + id='dwidenoise2_legacy_schedule', ), pytest.param('patch2self', 'auto', Patch2Self, {}, id='patch2self'), ], @@ -379,7 +393,7 @@ def test_denoising_wf_magnitude( assert 'combine_complex' not in nodes assert 'split_complex' not in nodes - _assert_denoiser_is_masked(nodes, denoise_method, nibs_dwi['dwi_file']) + _assert_denoiser_is_masked(nodes, nibs_dwi['dwi_file']) _assert_denoising_outputs(nodes, sink_dir, nibs_dwi['dwi_file']) @@ -387,11 +401,11 @@ def test_denoising_wf_magnitude( ('denoise_method', 'interface', 'expected_inputs'), [ pytest.param('dwidenoise', mrtrix.DWIDenoise, {'extent': (5, 5, 5)}, id='dwidenoise'), - pytest.param('dwidenoise2', mrtrix.DWIDenoise2, {'shape': 'sphere'}, id='dwidenoise2'), + pytest.param('dwidenoise2', mrtrix.DWIDenoise2, {}, id='dwidenoise2'), pytest.param( - 'dwidenoise2;demodulate:nonlinear', + 'dwidenoise2;demodulate:hann', mrtrix.DWIDenoise2, - {'demodulate': 'nonlinear'}, + {'demodulate': 'hann'}, id='dwidenoise2_demodulate', ), pytest.param('patch2self', Patch2Self, {}, id='patch2self_ignores_phase'), @@ -431,5 +445,5 @@ def test_denoising_wf_complex( complex_img = nb.load(nodes['combine_complex'].result.outputs.out_file) assert np.issubdtype(complex_img.header.get_data_dtype(), np.complexfloating) - _assert_denoiser_is_masked(nodes, denoise_method, nibs_dwi['dwi_file']) + _assert_denoiser_is_masked(nodes, nibs_dwi['dwi_file']) _assert_denoising_outputs(nodes, sink_dir, nibs_dwi['dwi_file']) diff --git a/qsiprep/tests/utils.py b/qsiprep/tests/utils.py index b705087a9..fb8639d26 100644 --- a/qsiprep/tests/utils.py +++ b/qsiprep/tests/utils.py @@ -39,6 +39,7 @@ def download_test_data(dset, data_dir=None): 'https://upenn.box.com/shared/static/tkahg1ctipmfihvpa1gmibvcv0gb721h.xz' ), 'forrest_gump': 'https://upenn.box.com/shared/static/qat58an322bzzyixrrsk7cmf52q3bepq.xz', + 'nibs': 'https://upenn.box.com/shared/static/bkllff4ik51jy9ju6nben2r5zrq4a5me.xz', } if dset == '*': for k in URLS: diff --git a/qsiprep/utils/misc.py b/qsiprep/utils/misc.py index 3490e012d..f1f0cb7a7 100644 --- a/qsiprep/utils/misc.py +++ b/qsiprep/utils/misc.py @@ -8,33 +8,40 @@ LOGGER = logging.getLogger('nipype.interface') +# The sliding-window kernel and the subsampling factor are no longer command-line options: +# they are properties of the multi-resolution schedule, chosen with 'schedule'. QSIPrep uses +# the bundled 'default' schedule unless one is requested, so --dwi-denoise-window does not +# apply to dwidenoise2. _DWIDENOISE_ENUM_PARAMETERS = { + 'aggregator': ('exclusive', 'gaussian', 'invl0', 'rank', 'uniform'), 'datatype': ('float32', 'float64'), + 'debias_anchor': ('sample', 'group_mean'), 'decomposition': ('bdcsvd', 'selfadjoint'), - 'estimator': ('Exp1', 'Exp2', 'Med', 'MRM2023'), - 'shape': ('sphere', 'cuboid'), - 'demodulate': ('none', 'linear', 'nonlinear'), 'demean': ('none', 'volume_groups', 'shells', 'all'), + 'demodulate': ('none', 'linear', 'hann', 'apc'), + 'estimator': ('exp1', 'exp2', 'med', 'mrm2023', 'tbme2022'), 'filter_method': ('optshrink', 'optthresh', 'truncate'), - 'aggregator': ('exclusive', 'gaussian', 'invl0', 'rank', 'uniform'), + 'vst_method': ('none', 'linear', 'foi', 'koay', 'mom'), } _DWIDENOISE_STRING_PARAMETERS = { 'demod_axes', - 'vst', + 'eigenspectra', + 'lamplus', + 'max_dist', + 'noise_image', + 'patchcount', 'preconditioned_input', 'preconditioned_output', - 'noise_image', - 'lamplus', - 'rank_pcanonzero', 'rank_input', 'rank_output', - 'variance_removed', - 'eigenspectra', - 'max_dist', - 'voxelcount', - 'patchcount', + 'rank_pcanonzero', + 'rankpermm_in', + 'rankpermm_out', + 'schedule', 'sum_aggregation', 'sum_optshrink', + 'variance_removed', + 'voxelcount', 'grad_file', 'bvec_file', 'bval_file', @@ -43,24 +50,30 @@ set(_DWIDENOISE_ENUM_PARAMETERS) | _DWIDENOISE_STRING_PARAMETERS | { - 'onepass', - 'noise_in', + 'aggregator_fwhm', 'fixed_rank', - 'radius', - 'aspect_ratio', - 'minvoxels', - 'extent', - 'subsample', + 'noise_dof', + 'noise_in', + 'preserve_noise_bias', 'residual_statistics', } ) -def parse_denoise_method(spec): +def parse_denoise_method(spec, use_phase=None): """Parse a denoising method and semicolon-delimited parameters. Parameters for dwidenoise2 use ``name:value`` syntax, for example - ``dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd``. + ``dwidenoise2;demodulate:apc;decomposition:bdcsvd``. + + Parameters + ---------- + spec : str + The ``--denoise-method`` specification. + use_phase : bool or None + Whether phase data are available for the series being denoised. ``None`` means + that is not known yet, as when the CLI validates the specification before any + scan has been selected, and skips the checks that depend on it. """ elements = spec.split(';') method = elements[0].strip() @@ -85,31 +98,22 @@ def parse_denoise_method(spec): choices = _DWIDENOISE_ENUM_PARAMETERS[name] if value not in choices: raise ValueError(f'Invalid value for {name!r}: {value!r}; choose from {choices}') - if name == 'shape' and value == 'cuboid': - # dwidenoise2 rejects an odd -extent, but the denoising window is always - # rounded up to an odd number, so cuboid kernels can't be used yet. - raise ValueError("'shape:cuboid' is not supported yet; use 'shape:sphere' instead") parsed_value = value - elif name == 'onepass': + elif name == 'preserve_noise_bias': bool_values = {'true': True, 'false': False, '1': True, '0': False} try: parsed_value = bool_values[value.lower()] except KeyError as exc: raise ValueError(f'Invalid boolean value for {name!r}: {value!r}') from exc - elif name in ('fixed_rank', 'minvoxels'): + elif name in ('fixed_rank', 'noise_dof'): parsed_value = int(value) - elif name in ('radius', 'aspect_ratio'): + elif name == 'aggregator_fwhm': parsed_value = float(value) elif name == 'noise_in': try: parsed_value = float(value) except ValueError: parsed_value = value - elif name in ('extent', 'subsample'): - values = tuple(int(item.strip()) for item in value.split(',')) - if len(values) not in (1, 3): - raise ValueError(f'{name!r} must contain one or three integers') - parsed_value = values[0] if len(values) == 1 else values elif name == 'residual_statistics': parsed_value = tuple(item.strip() for item in value.split(',')) if len(parsed_value) != 3 or not all(parsed_value): @@ -119,9 +123,190 @@ def parse_denoise_method(spec): parameters[name] = parsed_value + if method == 'dwidenoise2' and use_phase is False: + demodulation = parameters.get('demodulate', 'none') + if demodulation != 'none': + raise ValueError( + f'dwidenoise2 cannot apply {demodulation!r} phase demodulation to ' + 'magnitude-only data. Provide phase data or use "demodulate:none".' + ) + return method, parameters +_DWIDENOISE2_DEFAULTS = { + # Defaults of the dwidenoise2 build QSIPrep ships, read from its source. The boilerplate + # describes these too, so that the methods reflect what actually ran rather than only the + # parameters QSIPrep set explicitly. + 'aggregator': 'gaussian', + 'decomposition': 'bdcsvd', + # 'apc' is the default for complex data; magnitude data are never demodulated + 'demodulate': 'apc', + # 'shells' is the default whenever a gradient table is available, which QSIPrep always + # supplies + 'demean': 'shells', + 'estimator': 'mrm2023', + 'schedule': 'default', +} + +_DWIDENOISE2_ESTIMATORS = { + 'exp1': 'the Marchenko-Pastur threshold search of the original `dwidenoise` [@dwidenoise1]', + 'exp2': 'a refined Marchenko-Pastur threshold search [@cordero2019complex]', + 'med': 'the median eigenvalue [@gavish2014]', + 'mrm2023': 'a Marchenko-Pastur fit generalized to multi-dimensional data [@olesen2023]', + 'tbme2022': 'a multiple-moment generalized quarter-circle estimator [@zhu2022]', +} + +_DWIDENOISE2_FILTERS = { + 'optshrink': ( + 'optimal shrinkage of the singular values, which minimizes the Frobenius norm ' + '[@cordero2019complex]' + ), + 'optthresh': 'an optimal hard threshold on the singular values [@gavish2014]', + 'truncate': 'hard truncation, as in the original `dwidenoise` [@dwidenoise1]', +} + +_DWIDENOISE2_DEMODULATION = { + 'apc': ( + 'noise-adaptive phase correction, which re-estimates the background phase at every ' + 'noise level iteration [@pizzolato2020]' + ), + 'hann': 'a fixed nonlinear phase estimate from a Hann-windowed k-space filter [@patron2024]', + 'linear': 'a strictly linear phase term regressed from each k-space [@cordero2019complex]', +} + +_DWIDENOISE2_DEMEAN = { + 'shells': 'the mean signal of each *b*-value shell was regressed out', + 'volume_groups': 'the mean signal of each volume group was regressed out', + 'all': 'the mean signal across all volumes was regressed out', +} + + +def _join_clauses(clauses): + """Join clauses into a comma-separated list with a trailing 'and'.""" + if len(clauses) == 1: + return clauses[0] + + return f'{", ".join(clauses[:-1])} and {clauses[-1]}' + + +def describe_dwidenoise2(parameters, complex_data): + """Describe a ``dwidenoise2`` call for the methods boilerplate. + + ``dwidenoise2`` applies a number of methods beyond the original ``dwidenoise``, most of + them on by default, and each carries its own citation. Describing only the parameters + QSIPrep passed explicitly would therefore both understate what ran and omit references + the authors ask for, so unset options are described using the defaults of the shipped + build. The conditions attached to each citation follow the reference list that + ``dwidenoise2`` prints in its own help. + + Parameters + ---------- + parameters : dict + DWIDenoise2 parameters, as returned by :func:`parse_denoise_method`. + complex_data : bool + Whether ``dwidenoise2`` is run on complex-valued data. Phase demodulation only + applies to complex data, and only magnitude data need a nonlinear + variance-stabilizing transform. + + Returns + ------- + str + Boilerplate text with inline ``[@citation]`` keys, beginning with 'denoised using' + so that the caller can supply its own subject. + """ + used = {**_DWIDENOISE2_DEFAULTS, **parameters} + # The kernel size and the number of PCAs are set per iteration by the schedule rather + # than by a fixed window + schedule = used['schedule'] + schedule_desc = ( + 'its default schedule' if schedule == 'default' else f'the {schedule!r} schedule' + ) + + sentences = [ + 'denoised using the Marchenko-Pastur PCA method [@dwidenoise1; @dwidenoise2] as ' + 'implemented in `dwidenoise2` [@dwidenoise2software; @cordero2019complex], which ' + 'estimates the noise level over a multi-resolution series of iterations following ' + f'{schedule_desc}, sizing the sliding-window patch for noise estimation and for ' + 'denoising separately.' + ] + + preconditioning = [] + if complex_data and used['demodulate'] != 'none': + demodulation = _DWIDENOISE2_DEMODULATION[used['demodulate']] + preconditioning.append( + f'the complex-valued data were phase-demodulated using {demodulation}' + ) + if used['demean'] != 'none': + preconditioning.append(_DWIDENOISE2_DEMEAN[used['demean']]) + + # Complex data are Gaussian, so they always take the linear transform; magnitude data + # get a nonlinear one to account for the non-central chi noise distribution + vst_method = used.get('vst_method', 'linear' if complex_data else 'foi') + if not complex_data and vst_method in ('foi', 'koay', 'mom'): + vst = ( + 'a nonlinear variance-stabilizing transform was applied to render the ' + 'non-central chi distributed magnitude data approximately Gaussian and ' + 'homoscedastic [@foi2011; @ma2020]' + ) + if vst_method == 'koay': + vst += ', inverted with an analytically exact correction scheme [@koay2006]' + if 'noise_dof' in used: + vst += f', assuming {used["noise_dof"]} receive channels' + preconditioning.append(vst) + elif vst_method == 'linear': + preconditioning.append('the data were scaled by the local noise level') + + if preconditioning: + sentences.append(f'Prior to PCA, {_join_clauses(preconditioning)}.') + + decomposition = ( + 'a bidirectional divide-and-conquer SVD' + if used['decomposition'] == 'bdcsvd' + else 'a self-adjoint eigendecomposition' + ) + if 'noise_in' in used: + estimation = 'the noise level was taken from a pre-estimated noise map' + elif 'fixed_rank' in used: + estimation = f'the signal rank was fixed at {used["fixed_rank"]}' + else: + estimation = ( + 'the noise level was estimated from the eigenspectrum using ' + f'{_DWIDENOISE2_ESTIMATORS[used["estimator"]]}' + ) + sentences.append(f'Each patch was decomposed with {decomposition}, and {estimation}.') + + # dwidenoise2 truncates rather than shrinks when the rank is given rather than estimated + default_filter = 'truncate' if 'fixed_rank' in used else 'optshrink' + filter_method = used.get('filter_method', default_filter) + reconstruction = ( + f'Component contributions were filtered by {_DWIDENOISE2_FILTERS[filter_method]}' + ) + if used['aggregator'] == 'exclusive': + reconstruction += ( + ', and each voxel was reconstructed solely from the patch centered on it.' + ) + elif used['aggregator'] == 'gaussian': + reconstruction += ( + ', and each voxel was reconstructed from every overlapping patch, weighted by a ' + 'Gaussian function of its distance to each patch center [@manjon2013].' + ) + else: + reconstruction += ( + ', and each voxel was reconstructed from every overlapping patch, combined with ' + f'{used["aggregator"]} weighting [@manjon2013].' + ) + sentences.append(reconstruction) + + if not complex_data and not used.get('preserve_noise_bias', False): + sentences.append( + 'The inverse transform was evaluated at the exact-unbiased operating point, ' + 'removing the noise-floor bias from the denoised magnitude data.' + ) + + return ' '.join(sentences) + ' ' + + def safe_unit_vector(vector): """Return the unit vector of ``vector``. diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index b4741e961..92588a41c 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -33,7 +33,7 @@ from ...interfaces.nilearn import MaskEPI, Merge from ...interfaces.tortoise import Gibbs from ...utils.bids import IMPORTANT_DWI_FIELDS, update_metadata_from_nifti_header -from ...utils.misc import parse_denoise_method +from ...utils.misc import describe_dwidenoise2, parse_denoise_method from .qc import init_modelfree_qc_wf from .util import _get_wf_name @@ -418,14 +418,10 @@ def get_buffernode(): ]) # fmt:skip # Which steps to apply? - denoise_method, dwidenoise_params = parse_denoise_method(config.workflow.denoise_method) - if denoise_method == 'dwidenoise2' and not use_phase: - demodulation = dwidenoise_params.get('demodulate', 'none') - if demodulation != 'none': - raise ValueError( - f'dwidenoise2 cannot apply {demodulation!r} phase demodulation to ' - 'magnitude-only data. Provide phase data or use "demodulate:none".' - ) + denoise_method, dwidenoise2_params = parse_denoise_method( + config.workflow.denoise_method, + use_phase=use_phase, + ) unringing_method = config.workflow.unringing_method do_denoise = denoise_method in ('patch2self', 'dwidenoise', 'dwidenoise2') @@ -441,8 +437,7 @@ def get_buffernode(): # the output of any earlier step. # ``dwidenoise`` restricts the voxels it processes to this mask, whereas ``dwidenoise2`` # has no -mask option, so there the mask only sets the contour drawn on the report. - mask_denoiser = denoise_method.startswith('dwidenoise') - if mask_denoiser or do_biascorr: + if do_denoise or do_biascorr: get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') workflow.connect([ @@ -471,170 +466,131 @@ def get_buffernode(): mem_gb=DEFAULT_MEMORY_MIN_GB, ) - dwi_denoise_window = config.workflow.dwi_denoise_window - auto_str = '' - if denoise_method.startswith('dwidenoise') and dwi_denoise_window == 'auto': - # Configure the denoising window - import numpy as np - - dwi_denoise_window = closest_odd(int(np.ceil(np.cbrt(n_volumes)))) - dwi_denoise_window = max(dwi_denoise_window, 3) - config.loggers.workflow.info( - f'Automatically using {dwi_denoise_window}, {dwi_denoise_window}, ' - f'{dwi_denoise_window} window for dwidenoise' - ) - auto_str = 'n automatically-determined' - # Only the dwidenoise variants can denoise complex-valued data. # Any other method ignores the phase data and denoises the magnitude data alone. denoise_complex = denoise_method.startswith('dwidenoise') and use_phase - if denoise_complex: - desc += ( - 'Magnitude and phase DWI data were combined into a complex-valued file, ' - 'then denoised using the Marchenko-Pastur PCA method implemented in ' - f'{denoise_method} ' - '[@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, ' - - # If there are phase files available, then we can use dwidenoise - # on the complex-valued data. - phase_to_radians = pe.Node( - PhaseToRad(), - name='phase_to_radians', + # Build the denoiser. The node is the same whether it is handed magnitude-only or + # complex-valued data; only the data feeding it differs, which is wired up below. + if denoise_method == 'dwidenoise2': + # dwidenoise2 sizes its patches per iteration from its multi-resolution schedule, + # so there is no kernel to configure and dwi_denoise_window does not apply here. + denoiser = pe.Node( + DWIDenoise2(nthreads=omp_nthreads, **dwidenoise2_params), + name='denoiser', n_procs=omp_nthreads, ) - workflow.connect([(inputnode, phase_to_radians, [('dwi_phase_file', 'phase_file')])]) - combine_complex = pe.Node( - PolarToComplex(), - name='combine_complex', - n_procs=omp_nthreads, - ) + # dwidenoise2 needs the gradient table to demean by shell. The standalone + # dwidenoise2 build misreads the two files given to its -fslgrad option, so + # supply the gradients as a single MRtrix-format table instead. + gradient_table = pe.Node(MRTrixGradientTable(), name='gradient_table') workflow.connect([ - (buffernodes[-2], combine_complex, [('dwi_file', 'mag_file')]), - (phase_to_radians, combine_complex, [('phase_file', 'phase_file')]), + (inputnode, gradient_table, [ + ('bval_file', 'bval_file'), + ('bvec_file', 'bvec_file'), + ]), + (gradient_table, denoiser, [('gradient_file', 'grad_file')]), ]) # fmt:skip - - if denoise_method == 'dwidenoise2': - dwidenoise_inputs = {'shape': 'sphere', 'nthreads': omp_nthreads} - dwidenoise_inputs.update(dwidenoise_params) - denoiser = pe.Node( - DWIDenoise2(**dwidenoise_inputs), - name='denoiser', - n_procs=omp_nthreads, - ) - else: - denoiser = pe.Node( - DWIDenoise( - extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), - nthreads=omp_nthreads, - ), - name='denoiser', - n_procs=omp_nthreads, + elif denoise_method == 'dwidenoise': + dwi_denoise_window = config.workflow.dwi_denoise_window + auto_str = '' + if dwi_denoise_window == 'auto': + # Configure the denoising window + import numpy as np + + dwi_denoise_window = closest_odd(int(np.ceil(np.cbrt(n_volumes)))) + dwi_denoise_window = max(dwi_denoise_window, 3) + config.loggers.workflow.info( + f'Automatically using {dwi_denoise_window}, {dwi_denoise_window}, ' + f'{dwi_denoise_window} window for dwidenoise' ) + auto_str = 'n automatically-determined' - workflow.connect([ - (combine_complex, denoiser, [('out_file', 'in_file')]), - (denoiser, ds_report_denoising, [('out_report', 'in_file')]), - (denoiser, merge_confounds, [('nmse_text', f'in{step_num}')]), - ]) # fmt:skip - - split_complex = pe.Node( - ComplexToMagnitude(), - name='split_complex', + denoiser = pe.Node( + DWIDenoise( + extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), + nthreads=omp_nthreads, + ), + name='denoiser', 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.startswith('dwidenoise'): - desc += ( - 'DWI data were ' - f'denoised using the Marchenko-Pastur PCA method implemented in {denoise_method} ' - '[@mrtrix3; @dwidenoise1; @dwidenoise2; @cordero2019complex] ' - f'with a{auto_str} window size of {dwi_denoise_window} voxels. ' + else: + denoiser = pe.Node( + Patch2Self(), + name='denoiser', + n_procs=omp_nthreads, ) - last_step = 'After MP-PCA, ' + workflow.connect([(inputnode, denoiser, [('bval_file', 'bval_file')])]) + if denoise_method.startswith('dwidenoise'): if denoise_method == 'dwidenoise2': - dwidenoise_inputs = { - 'shape': 'sphere', - 'radius': dwi_denoise_window, - 'nthreads': omp_nthreads, - } - if dwidenoise_params.get('shape') == 'cuboid': - for parameter in ('radius',): - if parameter not in dwidenoise_params: - dwidenoise_inputs.pop(parameter) - - # cuboid uses extent instead of radius - dwidenoise_inputs['extent'] = ( - dwi_denoise_window, - dwi_denoise_window, - dwi_denoise_window, - ) - - dwidenoise_inputs.update(dwidenoise_params) - denoiser = pe.Node( - DWIDenoise2(**dwidenoise_inputs), - name='denoiser', - n_procs=omp_nthreads, - ) + # dwidenoise2 turns on a number of methods by default, each with its own + # citation, so the description is compiled from the parameters in effect + mppca_desc = describe_dwidenoise2(dwidenoise2_params, complex_data=denoise_complex) else: - denoiser = pe.Node( - DWIDenoise( - extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), - nthreads=omp_nthreads, - ), - name='denoiser', - n_procs=omp_nthreads, + mppca_desc = ( + 'denoised using the Marchenko-Pastur PCA method implemented in dwidenoise ' + '[@mrtrix3; @dwidenoise1; @dwidenoise2] ' + f'with a{auto_str} window size of {dwi_denoise_window} voxels. ' + ) + + if denoise_complex: + desc += ( + 'Magnitude and phase DWI data were combined into a complex-valued file, then ' + f'{mppca_desc}' + 'After denoising, the complex-valued data were split back into magnitude and ' + 'phase, and the denoised magnitude data were retained. ' ) + else: + desc += f'DWI data were {mppca_desc}' + + last_step = 'After MP-PCA, ' else: desc += ( "DWI data were denoised using DiPy's Patch2Self algorithm [@dipy; @patch2self] " 'with an automatically-defined window size. ' ) last_step = 'After `patch2self`, ' - denoiser = pe.Node( - Patch2Self(), - name='denoiser', - n_procs=omp_nthreads, - ) - - if denoise_method == 'patch2self': - workflow.connect([(inputnode, denoiser, [('bval_file', 'bval_file')])]) - elif denoise_method == 'dwidenoise2': - # dwidenoise2 needs the gradient table to demean by shell. The standalone - # dwidenoise2 build misreads the two files given to its -fslgrad option, so - # supply the gradients as a single MRtrix-format table instead. - gradient_table = pe.Node(MRTrixGradientTable(), name='gradient_table') - workflow.connect([ - (inputnode, gradient_table, [ - ('bval_file', 'bval_file'), - ('bvec_file', 'bvec_file'), - ]), - (gradient_table, denoiser, [('gradient_file', 'grad_file')]), - ]) # fmt:skip + # Wiring that is the same for every denoising method workflow.connect([ - (quick_mask, denoiser, [('out_mask', 'mask')]) + (quick_mask, denoiser, [('out_mask', 'mask')]), + (denoiser, ds_report_denoising, [('out_report', 'in_file')]), + (denoiser, merge_confounds, [('nmse_text', f'in{step_num}')]), # The noise image is a derivative, so it always comes straight from the denoiser (denoiser, outputnode, [('noise_image', 'noise_image')]), ]) # fmt:skip - if not denoise_complex: + # The denoiser's input and output are all that the complex-valued path changes + if denoise_complex: + phase_to_radians = pe.Node( + PhaseToRad(), + name='phase_to_radians', + n_procs=omp_nthreads, + ) + combine_complex = pe.Node( + PolarToComplex(), + name='combine_complex', + n_procs=omp_nthreads, + ) + split_complex = pe.Node( + ComplexToMagnitude(), + name='split_complex', + n_procs=omp_nthreads, + ) + workflow.connect([ + (inputnode, phase_to_radians, [('dwi_phase_file', 'phase_file')]), + (buffernodes[-2], combine_complex, [('dwi_file', 'mag_file')]), + (phase_to_radians, combine_complex, [('phase_file', 'phase_file')]), + (combine_complex, denoiser, [('out_file', 'in_file')]), + (denoiser, split_complex, [('out_file', 'complex_file')]), + (split_complex, buffernodes[-1], [('out_file', 'dwi_file')]), + ]) # fmt:skip + else: workflow.connect([ (buffernodes[-2], denoiser, [('dwi_file', 'in_file')]), - (denoiser, ds_report_denoising, [('out_report', 'in_file')]), (denoiser, buffernodes[-1], [('out_file', 'dwi_file')]), - (denoiser, merge_confounds, [('nmse_text', f'in{step_num}')]), ]) # fmt:skip step_num += 1 From 90a6c6eaa518b48275cb0109e1d6033fdb4eeb5d Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 5 Aug 2026 12:56:38 -0400 Subject: [PATCH 12/14] Update test_interfaces_mrtrix.py --- qsiprep/tests/test_interfaces_mrtrix.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/qsiprep/tests/test_interfaces_mrtrix.py b/qsiprep/tests/test_interfaces_mrtrix.py index d06471a26..6ab318677 100644 --- a/qsiprep/tests/test_interfaces_mrtrix.py +++ b/qsiprep/tests/test_interfaces_mrtrix.py @@ -3,12 +3,18 @@ import os import nibabel as nb +import numpy as np import pytest from traits.trait_errors import TraitError from qsiprep.interfaces import mrtrix +def _field_of_view(img): + """Return the spatial extent of an image in mm.""" + return np.array(img.shape[:3]) * np.array(img.header.get_zooms()[:3]) + + def test_dwidenoise(datasets, tmp_path_factory): """Test qsiprep.interfaces.mrtrix.DWIDenoise.""" tmpdir = tmp_path_factory.mktemp('test_dwidenoise') @@ -57,8 +63,17 @@ def test_dwidenoise2(datasets, tmp_path_factory): assert os.path.isfile(results.outputs.noise_image) noise_img = nb.load(results.outputs.noise_image) - assert noise_img.shape == in_img.shape[:3] assert noise_img.ndim == 3 + # dwidenoise2 estimates the noise level on a subsampled grid, so the noise map is + # coarser than the input rather than voxel-for-voxel with it, and the subsampling + # factor is a property of the schedule rather than something QSIPrep sets. + assert all(n <= i for n, i in zip(noise_img.shape, in_img.shape[:3], strict=True)) + # Whatever factor the schedule uses, the coarse grid has to cover the input: its + # extent is the input's rounded up to a whole number of its own (larger) voxels, so + # it can never fall short, nor overshoot by a full voxel. + overshoot = _field_of_view(noise_img) - _field_of_view(in_img) + assert np.all(overshoot > -1e-4) + assert np.all(overshoot < np.array(noise_img.header.get_zooms()[:3])) assert os.path.isfile(results.outputs.out_report) assert os.path.isfile(results.outputs.nmse_text) From 2ae9aa628784b0f861194b59dbdc196a8201b0e2 Mon Sep 17 00:00:00 2001 From: mattcieslak Date: Wed, 12 Aug 2026 10:39:50 -0400 Subject: [PATCH 13/14] PR review --- qsiprep/cli/parser.py | 19 +++++++--- qsiprep/tests/conftest.py | 2 +- qsiprep/tests/test_cli.py | 4 +-- qsiprep/tests/test_interfaces_mrtrix.py | 16 +++------ qsiprep/tests/test_utils_misc.py | 13 +++++-- qsiprep/tests/test_workflows_merge.py | 47 ++++++++----------------- qsiprep/tests/utils.py | 7 ++++ qsiprep/utils/misc.py | 16 +-------- qsiprep/workflows/dwi/merge.py | 28 +++++---------- 9 files changed, 61 insertions(+), 91 deletions(-) diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index 884ed9451..f073d3af7 100644 --- a/qsiprep/cli/parser.py +++ b/qsiprep/cli/parser.py @@ -436,10 +436,7 @@ def _bids_filter(value, parser): '"patch2self" (DIPY), or "none".\n' 'dwidenoise2 parameters may follow the method as semicolon-delimited ' 'name:value pairs, for example ' - '"dwidenoise2;demodulate:nonlinear;decomposition:bdcsvd".\n' - 'To approximate legacy "dwidenoise" with "dwidenoise2", use ' - '"dwidenoise2;shape:cuboid;subsample:1;demodulate:none;demean:none;' - 'filter_method:truncate;aggregator:exclusive".' + '"dwidenoise2;demodulate:linear;decomposition:bdcsvd".' ), ) g_conf.add_argument( @@ -888,8 +885,20 @@ def parse_args(args=None, namespace=None): ) # Validate the tricky options here - denoise_method, _ = parse_denoise_method(config.workflow.denoise_method) + denoise_method, denoise_params = parse_denoise_method(config.workflow.denoise_method) check_denoise_window(denoise_method, config.workflow.dwi_denoise_window) + if ( + config.workflow.denoise_after_combining + and denoise_params.get('demodulate', 'none') != 'none' + ): + # Temporary workaround for a bug in dwidenoise2: the concatenated series + # cannot be denoised with phase data. + parser.error( + '--denoise-after-combining cannot be used with phase demodulation ' + f'("demodulate:{denoise_params["demodulate"]}"). ' + 'Remove the demodulate parameter and use "--ignore phase" to denoise ' + 'the magnitude data only.' + ) bids_dir = config.execution.bids_dir output_dir = config.execution.output_dir diff --git a/qsiprep/tests/conftest.py b/qsiprep/tests/conftest.py index 1f68eed67..a417f56b1 100644 --- a/qsiprep/tests/conftest.py +++ b/qsiprep/tests/conftest.py @@ -49,7 +49,7 @@ def datasets(data_dir): def nibs_dwi(data_dir): """Locate the nibs DWI series used to test the denoising workflow. - The series is small (48x48x29x76) and has both magnitude and phase data, so it can + The series is small and has both magnitude and phase data, so it can exercise the complex-valued denoising paths without a long runtime. Tests using this fixture are skipped when the dataset is unavailable, which keeps diff --git a/qsiprep/tests/test_cli.py b/qsiprep/tests/test_cli.py index 5e1b2ce21..9ff194b4c 100644 --- a/qsiprep/tests/test_cli.py +++ b/qsiprep/tests/test_cli.py @@ -38,7 +38,7 @@ def test_dsdti_fmap(data_dir, output_dir, working_dir): This tests the following features: - Blip-up + Blip-down DWI series for TOPUP/Eddy - Eddy is run on a CPU - - dwidenoise is enabled explicitly. + - dwidenoise is enabled implicitly. Inputs ------ @@ -60,8 +60,6 @@ def test_dsdti_fmap(data_dir, output_dir, working_dir): '--write-graph', '--mem-mb=4096', '--output-resolution=5', - '--denoise-method=dwidenoise2;shape:cuboid;subsample:1;demodulate:none;demean:none;' - 'filter_method:truncate;aggregator:exclusive', ] _run_and_generate(TEST_NAME, parameters, test_main=False) diff --git a/qsiprep/tests/test_interfaces_mrtrix.py b/qsiprep/tests/test_interfaces_mrtrix.py index 6ab318677..6eec1d0a4 100644 --- a/qsiprep/tests/test_interfaces_mrtrix.py +++ b/qsiprep/tests/test_interfaces_mrtrix.py @@ -8,11 +8,7 @@ from traits.trait_errors import TraitError from qsiprep.interfaces import mrtrix - - -def _field_of_view(img): - """Return the spatial extent of an image in mm.""" - return np.array(img.shape[:3]) * np.array(img.header.get_zooms()[:3]) +from qsiprep.tests.utils import field_of_view def test_dwidenoise(datasets, tmp_path_factory): @@ -64,14 +60,10 @@ def test_dwidenoise2(datasets, tmp_path_factory): assert os.path.isfile(results.outputs.noise_image) noise_img = nb.load(results.outputs.noise_image) assert noise_img.ndim == 3 - # dwidenoise2 estimates the noise level on a subsampled grid, so the noise map is - # coarser than the input rather than voxel-for-voxel with it, and the subsampling - # factor is a property of the schedule rather than something QSIPrep sets. + # The schedule estimates noise on a subsampled grid, so the noise map is coarser + # than the input but must cover the input's field of view. assert all(n <= i for n, i in zip(noise_img.shape, in_img.shape[:3], strict=True)) - # Whatever factor the schedule uses, the coarse grid has to cover the input: its - # extent is the input's rounded up to a whole number of its own (larger) voxels, so - # it can never fall short, nor overshoot by a full voxel. - overshoot = _field_of_view(noise_img) - _field_of_view(in_img) + overshoot = field_of_view(noise_img) - field_of_view(in_img) assert np.all(overshoot > -1e-4) assert np.all(overshoot < np.array(noise_img.header.get_zooms()[:3])) diff --git a/qsiprep/tests/test_utils_misc.py b/qsiprep/tests/test_utils_misc.py index 208e7b9b0..1fd2c444d 100644 --- a/qsiprep/tests/test_utils_misc.py +++ b/qsiprep/tests/test_utils_misc.py @@ -52,7 +52,7 @@ def test_angle_between_finite_for_zero_vector(): def test_parse_denoise_method_parameters(): method, parameters = parse_denoise_method( 'dwidenoise2;demodulate:hann;decomposition:bdcsvd;' - 'preserve_noise_bias:true;noise_dof:8;aggregator_fwhm:2.5;schedule:vlarge', + 'preserve_noise_bias:true;noise_dof:8;schedule:vlarge', use_phase=True, ) @@ -62,7 +62,6 @@ def test_parse_denoise_method_parameters(): 'decomposition': 'bdcsvd', 'preserve_noise_bias': True, 'noise_dof': 8, - 'aggregator_fwhm': 2.5, 'schedule': 'vlarge', } @@ -111,6 +110,16 @@ def test_parse_denoise_method_rejects_demodulation_without_phase(demodulate): assert parse_denoise_method(spec) == ('dwidenoise2', {'demodulate': demodulate}) +def test_denoise_parameters_match_interface(): + """Every allowlisted dwidenoise2 parameter must be a trait on DWIDenoise2InputSpec.""" + from qsiprep.interfaces.mrtrix import DWIDenoise2 + from qsiprep.utils.misc import _DWIDENOISE_PARAMETERS + + trait_names = set(DWIDenoise2.input_spec().trait_names()) + missing = sorted(_DWIDENOISE_PARAMETERS - trait_names) + assert not missing + + def test_denoise_method_cli_parameter(tmp_path): spec = 'dwidenoise2;demodulate:apc;decomposition:bdcsvd' opts = _build_parser().parse_args( diff --git a/qsiprep/tests/test_workflows_merge.py b/qsiprep/tests/test_workflows_merge.py index e6ff69c62..10d2ced91 100644 --- a/qsiprep/tests/test_workflows_merge.py +++ b/qsiprep/tests/test_workflows_merge.py @@ -15,6 +15,7 @@ from qsiprep import config from qsiprep.interfaces import mrtrix from qsiprep.interfaces.dipy import Patch2Self +from qsiprep.tests.utils import field_of_view from qsiprep.workflows.dwi.merge import init_dwi_denoising_wf @@ -127,8 +128,8 @@ def test_dwidenoise2_cli_parameters_reach_workflow(monkeypatch): @pytest.mark.parametrize('denoise_method', ['dwidenoise', 'dwidenoise2', 'patch2self']) -def test_denoising_wf_builds_one_mask_for_denoising_and_biascorr(monkeypatch, denoise_method): - """Build the brain mask once and hand it to both the denoiser and bias correction.""" +def test_denoising_wf_masks_only_biascorr(monkeypatch, denoise_method): + """Build the brain mask for bias correction only; the denoisers get no mask.""" monkeypatch.setattr(config.workflow, 'denoise_method', denoise_method) monkeypatch.setattr(config.workflow, 'dwi_denoise_window', 5) monkeypatch.setattr(config.workflow, 'unringing_method', 'none') @@ -156,13 +157,13 @@ def test_denoising_wf_builds_one_mask_for_denoising_and_biascorr(monkeypatch, de if src is quick_mask for _, dest_field in data['connect'] } - assert consumers == {('denoiser', 'mask'), ('biascorr', 'mask')} + assert consumers == {('biascorr', 'mask')} - # The mask has to come from the raw series: denoising runs first, so deriving it from - # any later buffer would be circular + # The mask comes from the series feeding bias correction, not the raw data get_b0s = workflow.get_node('get_b0s') assert {src.name for src, dest, _ in workflow._graph.edges(data=True) if dest is get_b0s} == { - 'inputnode' + 'inputnode', + 'buffer01', } @@ -262,11 +263,6 @@ def _run_denoising_wf( return {node.name: node for node in graph.nodes}, sink_dir -def _field_of_view(img): - """Return the spatial extent of an image in mm.""" - return np.array(img.shape[:3]) * np.array(img.header.get_zooms()[:3]) - - def _sink_output(sink_dir, field): """Return the single file the DataSink wrote for ``field``.""" matches = sorted((sink_dir / field).glob('*')) @@ -274,23 +270,10 @@ def _sink_output(sink_dir, field): return matches[0] -def _assert_denoiser_is_masked(nodes, raw_file): - """Check that the denoiser is handed a brain mask built from the raw data. - - Every method gets the mask. ``dwidenoise`` restricts the voxels it processes to it, - while ``dwidenoise2`` and ``patch2self`` use it only for the report contour. - """ - mask_file = nodes['quick_mask'].result.outputs.out_mask - assert nodes['denoiser'].inputs.mask == mask_file - - raw_img = nb.load(raw_file) - mask_img = nb.load(mask_file) - assert mask_img.shape == raw_img.shape[:3] - assert np.allclose(mask_img.affine, raw_img.affine) - # A mask that selected everything or nothing would silently defeat the point - mask_data = mask_img.get_fdata() - assert set(np.unique(mask_data)) <= {0.0, 1.0} - assert 0 < mask_data.sum() < mask_data.size +def _assert_denoiser_is_not_masked(nodes): + """Check that the denoiser processes the full FOV rather than a masked subset.""" + assert 'quick_mask' not in nodes + assert not isdefined(nodes['denoiser'].inputs.mask) def _assert_denoising_outputs(nodes, sink_dir, raw_file): @@ -307,9 +290,7 @@ def _assert_denoising_outputs(nodes, sink_dir, raw_file): noise_img = nb.load(_sink_output(sink_dir, 'noise_image')) assert noise_img.ndim == 3 - # dwidenoise2 subsamples by default, so its noise map sits on a coarser grid than the - # input. Whatever the grid, it has to cover the same field of view. - assert np.allclose(_field_of_view(noise_img), _field_of_view(raw_img), rtol=0.05) + assert np.allclose(field_of_view(noise_img), field_of_view(raw_img), rtol=0.05) noise_data = noise_img.get_fdata() finite = np.isfinite(noise_data) assert finite.any() @@ -393,7 +374,7 @@ def test_denoising_wf_magnitude( assert 'combine_complex' not in nodes assert 'split_complex' not in nodes - _assert_denoiser_is_masked(nodes, nibs_dwi['dwi_file']) + _assert_denoiser_is_not_masked(nodes) _assert_denoising_outputs(nodes, sink_dir, nibs_dwi['dwi_file']) @@ -445,5 +426,5 @@ def test_denoising_wf_complex( complex_img = nb.load(nodes['combine_complex'].result.outputs.out_file) assert np.issubdtype(complex_img.header.get_data_dtype(), np.complexfloating) - _assert_denoiser_is_masked(nodes, nibs_dwi['dwi_file']) + _assert_denoiser_is_not_masked(nodes) _assert_denoising_outputs(nodes, sink_dir, nibs_dwi['dwi_file']) diff --git a/qsiprep/tests/utils.py b/qsiprep/tests/utils.py index 02b1ea07a..4358554fb 100644 --- a/qsiprep/tests/utils.py +++ b/qsiprep/tests/utils.py @@ -87,6 +87,13 @@ def download_test_data(dset, data_dir=None): return out_dir +def field_of_view(img): + """Return the spatial extent of an image in mm.""" + import numpy as np + + return np.array(img.shape[:3]) * np.array(img.header.get_zooms()[:3]) + + def get_test_data_path(): """Return the path to test datasets, terminated with separator. diff --git a/qsiprep/utils/misc.py b/qsiprep/utils/misc.py index d7dc17825..1a34b0a62 100644 --- a/qsiprep/utils/misc.py +++ b/qsiprep/utils/misc.py @@ -8,10 +8,6 @@ LOGGER = logging.getLogger('nipype.interface') -# The sliding-window kernel and the subsampling factor are no longer command-line options: -# they are properties of the multi-resolution schedule, chosen with 'schedule'. QSIPrep uses -# the bundled 'default' schedule unless one is requested, so --dwi-denoise-window does not -# apply to dwidenoise2. _DWIDENOISE_ENUM_PARAMETERS = { 'aggregator': ('exclusive', 'gaussian', 'invl0', 'rank', 'uniform'), 'datatype': ('float32', 'float64'), @@ -35,8 +31,6 @@ 'rank_input', 'rank_output', 'rank_pcanonzero', - 'rankpermm_in', - 'rankpermm_out', 'schedule', 'sum_aggregation', 'sum_optshrink', @@ -50,7 +44,6 @@ set(_DWIDENOISE_ENUM_PARAMETERS) | _DWIDENOISE_STRING_PARAMETERS | { - 'aggregator_fwhm', 'fixed_rank', 'noise_dof', 'noise_in', @@ -107,8 +100,6 @@ def parse_denoise_method(spec, use_phase=None): raise ValueError(f'Invalid boolean value for {name!r}: {value!r}') from exc elif name in ('fixed_rank', 'noise_dof'): parsed_value = int(value) - elif name == 'aggregator_fwhm': - parsed_value = float(value) elif name == 'noise_in': try: parsed_value = float(value) @@ -134,16 +125,11 @@ def parse_denoise_method(spec, use_phase=None): return method, parameters +# dwidenoise2's own defaults, mirrored here so the boilerplate describes what actually ran _DWIDENOISE2_DEFAULTS = { - # Defaults of the dwidenoise2 build QSIPrep ships, read from its source. The boilerplate - # describes these too, so that the methods reflect what actually ran rather than only the - # parameters QSIPrep set explicitly. 'aggregator': 'gaussian', 'decomposition': 'bdcsvd', - # 'apc' is the default for complex data; magnitude data are never demodulated 'demodulate': 'apc', - # 'shells' is the default whenever a gradient table is available, which QSIPrep always - # supplies 'demean': 'shells', 'estimator': 'mrm2023', 'schedule': 'default', diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 92588a41c..abe36f98b 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -432,22 +432,6 @@ def get_buffernode(): num_steps = sum(map(int, [do_denoise, do_unringing, do_biascorr, harmonize_b0s])) merge_confounds = pe.Node(niu.Merge(num_steps), name='merge_confounds') - # A single brain mask is shared by the denoising and bias correction steps. It is built - # from the raw series because denoising runs first, so the mask cannot be derived from - # the output of any earlier step. - # ``dwidenoise`` restricts the voxels it processes to this mask, whereas ``dwidenoise2`` - # has no -mask option, so there the mask only sets the contour drawn on the report. - if do_denoise or do_biascorr: - get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') - quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') - workflow.connect([ - (inputnode, get_b0s, [ - ('dwi_file', 'dwi_series'), - ('bval_file', 'bval_file'), - ]), - (get_b0s, quick_mask, [('b0_series', 'in_files')]), - ]) # fmt:skip - # Add the steps step_num = 1 # Merge inputs start at 1 last_step = '' @@ -481,9 +465,9 @@ def get_buffernode(): n_procs=omp_nthreads, ) - # dwidenoise2 needs the gradient table to demean by shell. The standalone - # dwidenoise2 build misreads the two files given to its -fslgrad option, so - # supply the gradients as a single MRtrix-format table instead. + # dwidenoise2 needs the gradient table to demean by shell. Temporary + # workaround for a bug in dwidenoise2: supply the gradients as a single + # MRtrix-format table instead of using -fslgrad. gradient_table = pe.Node(MRTrixGradientTable(), name='gradient_table') workflow.connect([ (inputnode, gradient_table, [ @@ -555,7 +539,6 @@ def get_buffernode(): # Wiring that is the same for every denoising method workflow.connect([ - (quick_mask, denoiser, [('out_mask', 'mask')]), (denoiser, ds_report_denoising, [('out_report', 'in_file')]), (denoiser, merge_confounds, [('nmse_text', f'in{step_num}')]), # The noise image is a derivative, so it always comes straight from the denoiser @@ -657,6 +640,8 @@ def get_buffernode(): last_step = True biascorr = pe.Node(DWIBiasCorrect(method='ants'), name='biascorr', n_procs=omp_nthreads) + get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') + quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') ds_report_biascorr = pe.Node( DerivativesDataSink( datatype='figures', @@ -672,6 +657,9 @@ def get_buffernode(): workflow.connect([ (buffernodes[-2], biascorr, [('dwi_file', 'in_file')]), + (buffernodes[-2], get_b0s, [('dwi_file', 'dwi_series')]), + (inputnode, get_b0s, [('bval_file', 'bval_file')]), + (get_b0s, quick_mask, [('b0_series', 'in_files')]), (quick_mask, biascorr, [('out_mask', 'mask')]), (biascorr, buffernodes[-1], [('out_file', 'dwi_file')]), (biascorr, outputnode, [('bias_image', 'bias_image')]), From 1d366cad4c19b859da6393e11378caf41138e2d2 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 12 Aug 2026 10:56:34 -0400 Subject: [PATCH 14/14] Minimize diff. --- qsiprep/workflows/dwi/merge.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index abe36f98b..3f7c1788f 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -640,8 +640,6 @@ def get_buffernode(): last_step = True biascorr = pe.Node(DWIBiasCorrect(method='ants'), name='biascorr', n_procs=omp_nthreads) - get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') - quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') ds_report_biascorr = pe.Node( DerivativesDataSink( datatype='figures', @@ -652,6 +650,9 @@ def get_buffernode(): run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB, ) + get_b0s = pe.Node(ExtractB0s(b0_threshold=config.workflow.b0_threshold), name='get_b0s') + quick_mask = pe.Node(MaskEPI(lower_cutoff=0.02), name='quick_mask') + # Add buffernode for bias-corrected DWI buffernodes.append(get_buffernode())