From 8bb6c0d0d96d9e2c0a37b1c738bdc6d8a2dcfb13 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 2 Dec 2024 14:18:52 -0500 Subject: [PATCH] Calculate patch radius from number of volumes in DWI run (#892) * Work on patch radius. * Determine n_volumes and use it. * Fix type. * Convert to integer. * Update boilerplate. * Fix style issue. * Enforce odd integer value. * Update boilerplate. * Remove unused variables. * Capitalize start of sentence. * Apply suggestions from code review Co-authored-by: Matt Cieslak --------- Co-authored-by: Matt Cieslak --- qsiprep/cli/parser.py | 45 +++++++++-- qsiprep/data/boilerplate.bib | 31 ++++++-- qsiprep/workflows/dwi/hmc.py | 2 +- qsiprep/workflows/dwi/merge.py | 128 +++++++++++++++++------------- qsiprep/workflows/fieldmap/syn.py | 2 +- 5 files changed, 137 insertions(+), 71 deletions(-) diff --git a/qsiprep/cli/parser.py b/qsiprep/cli/parser.py index 432bd946..25f080db 100644 --- a/qsiprep/cli/parser.py +++ b/qsiprep/cli/parser.py @@ -99,6 +99,23 @@ def _min_one(value, parser): raise parser.error("Argument can't be less than one.") return value + def _int_or_auto(value, parser): + """Ensure an argument is an integer or 'auto'.""" + if value.lower() == 'auto': + return value + try: + value = int(value) + except ValueError as exc: + raise parser.error('Argument must be an integer or "auto".') from exc + + if value < 1: + raise parser.error('Argument must be greater than zero.') + + if value % 2 == 0: + raise parser.error('Argument must be an odd integer.') + + 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()]) @@ -154,6 +171,7 @@ def _bids_filter(value, parser): PathExists = partial(_path_exists, parser=parser) IsFile = partial(_is_file, parser=parser) PositiveInt = partial(_min_one, parser=parser) + IntOrAuto = partial(_int_or_auto, parser=parser) BIDSFilter = partial(_bids_filter, parser=parser) # Arguments as specified by BIDS-Apps @@ -349,10 +367,17 @@ def _bids_filter(value, parser): g_conf.add_argument( '--dwi-denoise-window', action='store', + type=IntOrAuto, default='auto', - help='window size in voxels for image-based denoising, integer or "auto".' - 'If "auto", 5 will be used for dwidenoise and auto-configured for ' - 'patch2self based on the number of b>0 images.', + 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 ' + '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.' + ), ) g_conf.add_argument( '--denoise-method', @@ -360,7 +385,7 @@ def _bids_filter(value, parser): choices=['dwidenoise', 'patch2self', 'none'], default='dwidenoise', help='Image-based denoising method. Either "dwidenoise" (MRtrix), ' - '"patch2self" (DIPY) or none. (default: dwidenoise)', + '"patch2self" (DIPY) or "none". (default: dwidenoise)', ) g_conf.add_argument( '--unringing-method', @@ -692,10 +717,14 @@ def parse_args(args=None, namespace=None): # Validate the tricky options here if config.workflow.dwi_denoise_window != 'auto': - try: - _ = int(config.workflow.dwi_denoise_window) - except ValueError: - raise Exception("--dwi-denoise-window must be an integer or 'auto'") + if config.workflow.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': + config.loggers.cli.warning( + 'The --dwi-denoise-window option is not used when --denoise-method=none' + ) 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 c3bbff1e..2ae6528f 100644 --- a/qsiprep/data/boilerplate.bib +++ b/qsiprep/data/boilerplate.bib @@ -394,13 +394,15 @@ @article{msmt5tt } @article{mrtrix3, - title = {MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation}, - author = {J-Donald and Smith, Robert and Raffelt, David and Tabbara, Rami and Dhollander, Thijs and Pietsch, Maximilian and Christiaens, Daan and Jeurissen, Ben and Yeh, Chun-Hung and Connelly, Alan"}, - journal = {NeuroImage}, - volume = {202}, - pages = {116137}, - year = {2019}, - issn={1053-8119}, + title={MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation}, + author={Tournier, J-Donald and Smith, Robert and Raffelt, David and Tabbara, Rami and Dhollander, Thijs and Pietsch, Maximilian and Christiaens, Daan and Jeurissen, Ben and Yeh, Chun-Hung and Connelly, Alan}, + journal={Neuroimage}, + volume={202}, + pages={116137}, + year={2019}, + publisher={Elsevier}, + url={https://doi.org/10.1016/j.neuroimage.2019.116137}, + doi={10.1016/j.neuroimage.2019.116137} } @article{originalcsd, @@ -535,7 +537,8 @@ @article{patch2self author={Fadnavis, Shreyas and Batson, Joshua and Garyfallidis, Eleftherios}, journal={Advances in Neural Information Processing Systems}, volume={33}, - year={2020} + year={2020}, + url={https://proceedings.neurips.cc/paper_files/paper/2020/file/bc047286b224b7bfa73d4cb02de1238d-Paper.pdf}, } @article{noddi, @@ -640,3 +643,15 @@ @article{power2017simple url={https://doi.org/10.1016/j.neuroimage.2016.08.009}, doi={10.1016/j.neuroimage.2016.08.009} } + +@article{cordero2019complex, + title={Complex diffusion-weighted image estimation via matrix recovery under general noise models}, + author={Cordero-Grande, Lucilio and Christiaens, Daan and Hutter, Jana and Price, Anthony N and Hajnal, Jo V}, + journal={Neuroimage}, + volume={200}, + pages={391--404}, + year={2019}, + publisher={Elsevier}, + url={https://doi.org/10.1016/j.neuroimage.2019.06.039}, + doi={10.1016/j.neuroimage.2019.06.039} +} diff --git a/qsiprep/workflows/dwi/hmc.py b/qsiprep/workflows/dwi/hmc.py index 1e7fc76a..8263c259 100644 --- a/qsiprep/workflows/dwi/hmc.py +++ b/qsiprep/workflows/dwi/hmc.py @@ -635,7 +635,7 @@ def init_dwi_model_hmc_wf( name='outputnode', ) workflow.__desc__ = ( - 'The the SHORELine method was used to estimate head motion in b>0 ' + 'The SHORELine method was used to estimate head motion in b>0 ' 'images. This entails leaving out each b>0 image and reconstructing ' 'the others using 3dSHORE [@merlet3dshore]. The signal for the left-' f'out image serves as the registration target. A total of {num_iters} ' diff --git a/qsiprep/workflows/dwi/merge.py b/qsiprep/workflows/dwi/merge.py index 191b8f9c..ebca1bac 100644 --- a/qsiprep/workflows/dwi/merge.py +++ b/qsiprep/workflows/dwi/merge.py @@ -149,13 +149,13 @@ def init_merge_and_denoise_wf( # Get a data frame of the raw_dwi_files and their imaging parameters: dwi_df = get_acq_parameters_df(raw_dwi_files, layout=layout) - for dwi_num, row in dwi_df.iterrows(): - dwi_num += 1 # start at 1 + for i_dwi, row in dwi_df.iterrows(): + dwi_num = i_dwi + 1 # start at 1 dwi_file = row.BIDSFile # Conform each image to the requested orientation conformers.append( - pe.Node(ConformDwi(orientation=orientation), name='conform_dwis%02d' % dwi_num), + pe.Node(ConformDwi(orientation=orientation), name=f'conform_dwis{dwi_num:02d}'), ) conformers[-1].inputs.dwi_file = dwi_file @@ -187,11 +187,13 @@ def init_merge_and_denoise_wf( name=f'conform_phase{dwi_num}', ) + n_volumes = row.NumVolumes denoising_wfs.append( init_dwi_denoising_wf( partial_fourier=row.PartialFourier, phase_encoding_direction=row.PhaseEncodingAxis, source_file=dwi_file, + n_volumes=n_volumes, use_phase=use_phase, do_biascorr=do_biascorr, name=wf_name, @@ -279,10 +281,12 @@ def init_merge_and_denoise_wf( # Send the merged series for denoising merge_confounds = pe.Node(niu.Merge(2), name='merge_confounds') hstack_confounds = pe.Node(StackConfounds(axis=1), name='hstack_confounds') + n_volumes = dwi_df['NumVolumes'].sum() denoising_wf = init_dwi_denoising_wf( partial_fourier=get_merged_parameter(dwi_df, 'PartialFourier', 'all'), phase_encoding_direction=get_merged_parameter(dwi_df, 'PhaseEncodingAxis', 'all'), source_file=source_file, + n_volumes=n_volumes, use_phase=False, # can't use phase with concatenated data do_biascorr=do_biascorr, name='merged_denoise', @@ -312,6 +316,7 @@ def init_dwi_denoising_wf( source_file, partial_fourier, phase_encoding_direction, + n_volumes, use_phase, do_biascorr, name='denoise_wf', @@ -326,6 +331,10 @@ def init_dwi_denoising_wf( fraction of k-space acquired phase_encoding_direction : str direction of phase encoding + n_volumes : int + number of volumes in the DWI series. + Used to determine the window size for denoising if dwidenoise is used + and the 'auto' option is selected. use_phase : bool True if phase data are available for the DWI scan. If True, and ``denoise_method`` is ``dwidenoise``, then ``dwidenoise`` @@ -381,6 +390,7 @@ def init_dwi_denoising_wf( ) workflow = Workflow(name=name) omp_nthreads = config.nipype.omp_nthreads + desc = '\n\n' # Get IdentityInterfaces ready to hold intermediate results buffernodes = [] @@ -410,12 +420,6 @@ def get_buffernode(): do_denoise = denoise_method in ('patch2self', 'dwidenoise') do_unringing = config.workflow.unringing_method in ('mrdegibbs', 'rpg') harmonize_b0s = not config.workflow.no_b0_harmonization - # Configure the denoising window - if ( - config.workflow.denoise_method == 'dwidenoise' - ) and config.workflow.dwi_denoise_window == 'auto': - dwi_denoise_window = 5 - config.loggers.workflow.info('Automatically using 5, 5, 5 window for dwidenoise') # How many steps in the denoising pipeline num_steps = sum(map(int, [do_denoise, do_unringing, do_biascorr, harmonize_b0s])) @@ -423,6 +427,7 @@ def get_buffernode(): # Add the steps step_num = 1 # Merge inputs start at 1 + last_step = '' if do_denoise: # Add buffernode for denoised DWI buffernodes.append(get_buffernode()) @@ -438,7 +443,30 @@ def get_buffernode(): mem_gb=DEFAULT_MEMORY_MIN_GB, ) + dwi_denoise_window = config.workflow.dwi_denoise_window + auto_str = '' + if denoise_method == 'dwidenoise' and config.workflow.dwi_denoise_window == 'auto': + # Configure the denoising window + import numpy as np + + dwi_denoise_window = closest_odd(int(np.cbrt(n_volumes))) + 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' + if (denoise_method == '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 ' + '[@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( @@ -485,6 +513,14 @@ def get_buffernode(): ]) # fmt:skip elif denoise_method == 'dwidenoise': + desc += ( + 'DWI data were ' + 'denoised using the Marchenko-Pastur PCA method implemented in dwidenoise ' + '[@mrtrix3; @dwidenoise1; @dwidenoise2; @cordero2019complex] ' + f'with a{auto_str} window size of {dwi_denoise_window} voxels. ' + ) + last_step = 'After MP-PCA, ' + denoiser = pe.Node( DWIDenoise( extent=(dwi_denoise_window, dwi_denoise_window, dwi_denoise_window), @@ -494,8 +530,13 @@ def get_buffernode(): n_procs=omp_nthreads, ) 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(patch_radius=dwi_denoise_window), + Patch2Self(), name='denoiser', n_procs=omp_nthreads, ) @@ -513,12 +554,15 @@ def get_buffernode(): if do_unringing: if unringing_method == 'mrdegibbs': + desc += f'{last_step}Gibbs ringing was removed using MRtrix3 [@mrtrix3; @mrdegibbs]. ' degibbser = pe.Node( MRDeGibbs(nthreads=omp_nthreads), name='degibbser', n_procs=omp_nthreads, ) elif unringing_method == 'rpg': + desc += f'{last_step}Gibbs ringing was removed using TORTOISE [@pfgibbs]. ' + pe_code = { 'i': 0, 'i-': 0, @@ -538,6 +582,8 @@ def get_buffernode(): n_procs=omp_nthreads, ) + last_step = 'After unringing, ' + ds_report_unringing = pe.Node( DerivativesDataSink( datatype='figures', @@ -561,6 +607,12 @@ def get_buffernode(): step_num += 1 if do_biascorr: + desc += ( + f'{last_step}B1 field inhomogeneity was corrected using ' + '`dwibiascorrect` from MRtrix3 with the N4 algorithm [@n4]. ' + ) + last_step = True + biascorr = pe.Node(DWIBiasCorrect(method='ants'), name='biascorr', n_procs=omp_nthreads) ds_report_biascorr = pe.Node( DerivativesDataSink( @@ -597,6 +649,9 @@ def get_buffernode(): # Connect the final buffernode (the most recent output) to the outputnode workflow.connect([(buffernodes[-1], outputnode, [('dwi_file', 'dwi_file')])]) + if not last_step: + desc = 'No denoising steps were applied to the DWI data.' + # If any denoising operations were run, collect their confounds if step_num > 1: hstack_confounds = pe.Node(StackConfounds(axis=1), name='hstack_confounds') @@ -605,6 +660,8 @@ def get_buffernode(): (hstack_confounds, outputnode, [('confounds_file', 'confounds')]), ]) # fmt:skip + workflow.__desc__ = desc + return workflow @@ -615,61 +672,19 @@ def _as_list(item): def gen_denoising_boilerplate(): """Generate a methods boilerplate for the denoising workflow.""" - denoise_method = config.workflow.denoise_method - dwi_denoise_window = config.workflow.dwi_denoise_window - unringing_method = config.workflow.unringing_method b1_biascorrect_stage = config.workflow.b1_biascorrect_stage no_b0_harmonization = config.workflow.no_b0_harmonization b0_threshold = config.workflow.b0_threshold - use_phase = 'phase' not in config.workflow.ignore desc = [ f'Any images with a b-value less than {b0_threshold} s/mm^2 were treated as a ' '*b*=0 image.' ] - do_denoise = denoise_method in ('dwidenoise', 'patch2self') - do_unringing = unringing_method in ('rpg', 'mrdegibbs') harmonize_b0s = not no_b0_harmonization last_step = '' - if do_denoise: - if denoise_method == 'dwidenoise': - desc.append( - "MP-PCA denoising as implemented in MRtrix3's `dwidenoise`" - f'[@dwidenoise1] was applied with a {dwi_denoise_window}-voxel window.' - ) - last_step = 'After MP-PCA, ' - if use_phase: - desc.append( - 'When phase data were available, this was done on complex-valued data.' - ) - - elif denoise_method == 'patch2self': - desc.append('Denoising using `patch2self` [@patch2self] was applied') - if dwi_denoise_window == 'auto': - desc.append('with settings based on developer recommendations.') - else: - desc.append(f'with a {dwi_denoise_window}-voxel window.') - last_step = 'After `patch2self`, ' - - if do_unringing: - unringing_txt = { - 'mrdegibbs': "MRtrix3's `mrdegibbs` [@mrdegibbs].", - 'dipy': 'Dipy [@dipy].', - 'rpg': "TORTOISE's `Gibbs` [@pfgibbs].", - }[unringing_method] - - desc.append(f'{last_step}Gibbs unringing was performed using {unringing_txt}') - last_step = 'Following unringing, ' - - if b1_biascorrect_stage == 'legacy': - desc.append( - f'{last_step}B1 field inhomogeneity was corrected using ' - '`dwibiascorrect` from MRtrix3 with the N4 algorithm [@n4].' - ) - last_step = 'After B1 bias correction, ' if harmonize_b0s: desc.append( - f'{last_step}the mean intensity of the DWI series was adjusted ' + 'The mean intensity of the DWI series was adjusted ' 'so all the mean intensity of the b=0 images matched across each' 'separate DWI scanning sequence.' ) @@ -734,3 +749,10 @@ def get_merged_parameter(parameter_df, parameter_name, selection_mode='all'): return col.mode()[0] raise Exception("selection_mode must be 'all' or 'mode'") + + +def closest_odd(x): + if x % 2 == 0: + return x - 1 + else: + return x diff --git a/qsiprep/workflows/fieldmap/syn.py b/qsiprep/workflows/fieldmap/syn.py index 84c36723..4a6e8ebf 100644 --- a/qsiprep/workflows/fieldmap/syn.py +++ b/qsiprep/workflows/fieldmap/syn.py @@ -135,7 +135,7 @@ def init_syn_sdc_wf(bold_pe=None, atlas_threshold=2): workflow = Workflow(name='syn_sdc_wf') workflow.__desc__ = f"""\ A deformation field to correct for susceptibility distortions was estimated -based on *fmriprep*'s *fieldmap-less* approach. +based on *fMRIPrep*'s *fieldmap-less* approach. The deformation field is that resulting from co-registering the b0 reference to the same-subject T1w-reference with its intensity inverted [@fieldmapless1; @fieldmapless2].