From 819406336d234c90780c8be824d610025fdb2e56 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 19 Feb 2025 10:19:52 -0500 Subject: [PATCH 01/14] fix: Find boldref2anat, regardless of name of target --- fmriprep/data/io_spec.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fmriprep/data/io_spec.json b/fmriprep/data/io_spec.json index 364cff7e3..d4c54c666 100644 --- a/fmriprep/data/io_spec.json +++ b/fmriprep/data/io_spec.json @@ -34,7 +34,7 @@ "boldref2anat": { "datatype": "func", "from": "boldref", - "to": "anat", + "to": ["anat", "T1w", "T2w"], "mode": "image", "suffix": "xfm", "extension": ".txt" From fabab8b8369f3ffe718c4931a88e6cdcf9fd9106 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 5 Mar 2025 14:31:06 -0500 Subject: [PATCH 02/14] fix: Report pre-computed bold2t1w registration --- fmriprep/interfaces/reports.py | 33 ++++++++++++++++++++------------- fmriprep/workflows/bold/fit.py | 8 +++++++- 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/fmriprep/interfaces/reports.py b/fmriprep/interfaces/reports.py index 6d199ccee..49ac59054 100644 --- a/fmriprep/interfaces/reports.py +++ b/fmriprep/interfaces/reports.py @@ -204,7 +204,11 @@ class FunctionalSummaryInputSpec(TraitedSpec): desc='Phase-encoding direction detected', ) registration = traits.Enum( - 'FSL', 'FreeSurfer', mandatory=True, desc='Functional/anatomical registration method' + 'FSL', + 'FreeSurfer', + 'Precomputed', + mandatory=True, + desc='Functional/anatomical registration method', ) fallback = traits.Bool(desc='Boundary-based registration rejected') registration_dof = traits.Enum( @@ -239,18 +243,21 @@ def _generate_segment(self): else: stc = 'n/a' # TODO: Add a note about registration_init below? - reg = { - 'FSL': [ - 'FSL flirt with boundary-based registration' - f' (BBR) metric - {dof} dof', - 'FSL flirt rigid registration - 6 dof', - ], - 'FreeSurfer': [ - 'FreeSurfer bbregister ' - f'(boundary-based registration, BBR) - {dof} dof', - f'FreeSurfer mri_coreg - {dof} dof', - ], - }[self.inputs.registration][self.inputs.fallback] + if self.inputs.registration == 'Precomputed': + reg = 'Precomputed affine transformation' + else: + reg = { + 'FSL': [ + 'FSL flirt with boundary-based registration' + f' (BBR) metric - {dof} dof', + 'FSL flirt rigid registration - 6 dof', + ], + 'FreeSurfer': [ + 'FreeSurfer bbregister ' + f'(boundary-based registration, BBR) - {dof} dof', + f'FreeSurfer mri_coreg - {dof} dof', + ], + }[self.inputs.registration][self.inputs.fallback] pedir = get_world_pedir(self.inputs.orientation, self.inputs.pe_direction) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index a585d5646..9bca49abb 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -307,7 +307,13 @@ def init_bold_fit_wf( summary = pe.Node( FunctionalSummary( distortion_correction='None', # Can override with connection - registration=('FSL', 'FreeSurfer')[config.workflow.run_reconall], + registration=( + 'Precomputed' + if boldref2anat_xform + else 'FreeSurfer' + if config.workflow.run_reconall + else 'FSL' + ), registration_dof=config.workflow.bold2anat_dof, registration_init=config.workflow.bold2anat_init, pe_direction=metadata.get('PhaseEncodingDirection'), From bdc3564fb86ae980b0695b06e89977a736214c13 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 5 Mar 2025 14:31:39 -0500 Subject: [PATCH 03/14] data: Add fmap_spec.json for discovering precomputed fieldmaps --- fmriprep/data/fmap_spec.json | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 fmriprep/data/fmap_spec.json diff --git a/fmriprep/data/fmap_spec.json b/fmriprep/data/fmap_spec.json new file mode 100644 index 000000000..3091f33b3 --- /dev/null +++ b/fmriprep/data/fmap_spec.json @@ -0,0 +1,33 @@ +{ + "queries": { + "fieldmaps": { + "fieldmap": { + "datatype": "fmap", + "desc": "preproc", + "suffix": "fieldmap", + "extension": [ + ".nii.gz", + ".nii" + ] + }, + "coeffs": { + "datatype": "fmap", + "desc": ["coeff", "coeff0", "coeff1"], + "suffix": "fieldmap", + "extension": [ + ".nii.gz", + ".nii" + ] + }, + "magnitude": { + "datatype": "fmap", + "desc": "magnitude", + "suffix": "fieldmap", + "extension": [ + ".nii.gz", + ".nii" + ] + } + } + } +} From 31c772dd9d1345f763fd845d499ca5c589d6ec1c Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 5 Mar 2025 14:33:41 -0500 Subject: [PATCH 04/14] rf: Cache layout for precomputed derivatives --- fmriprep/utils/bids.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py index db18d8c01..0ee2e3331 100644 --- a/fmriprep/utils/bids.py +++ b/fmriprep/utils/bids.py @@ -28,6 +28,7 @@ import os import sys from collections import defaultdict +from functools import cache from pathlib import Path from bids.layout import BIDSLayout @@ -38,6 +39,15 @@ from ..data import load as load_data +@cache +def _get_layout(derivatives_dir: Path) -> BIDSLayout: + import niworkflows.data + + return BIDSLayout( + derivatives_dir, config=[niworkflows.data.load('nipreps.json')], validate=False + ) + + def collect_derivatives( derivatives_dir: Path, entities: dict, @@ -57,8 +67,7 @@ def collect_derivatives( patterns = _patterns derivs_cache = defaultdict(list, {}) - layout = BIDSLayout(derivatives_dir, config=['bids', 'derivatives'], validate=False) - derivatives_dir = Path(derivatives_dir) + layout = _get_layout(derivatives_dir) # search for both boldrefs for k, q in spec['baseline'].items(): From 70b716c3c946c83b24c9fea820c23e31538288e9 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 5 Mar 2025 14:34:24 -0500 Subject: [PATCH 05/14] feat: Discover pre-computed fieldmaps --- fmriprep/utils/bids.py | 25 +++++++++++++++++++++++++ fmriprep/workflows/base.py | 12 ++++++++++++ 2 files changed, 37 insertions(+) diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py index 0ee2e3331..3f26c71ce 100644 --- a/fmriprep/utils/bids.py +++ b/fmriprep/utils/bids.py @@ -95,6 +95,31 @@ def collect_derivatives( return derivs_cache +def collect_fieldmaps( + derivatives_dir: Path, + entities: dict, + spec: dict | None = None, +): + """Gather existing derivatives and compose a cache.""" + if spec is None: + spec = json.loads(load_data.readable('fmap_spec.json').read_text())['queries'] + + fmap_cache = defaultdict(dict, {}) + layout = _get_layout(derivatives_dir) + + fmapids = layout.get_fmapids(**entities) + + for fmapid in fmapids: + for k, q in spec['fieldmaps'].items(): + query = {**entities, **q} + item = layout.get(return_type='filename', fmapid=fmapid, **query) + if not item: + continue + fmap_cache[fmapid][k] = item[0] if len(item) == 1 else item + + return fmap_cache + + def write_bidsignore(deriv_dir): bids_ignore = ( '*.html', diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index c8a2fb8a7..f9229afca 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -533,6 +533,18 @@ def init_single_subject_wf(subject_id: str): if config.workflow.anat_only: return clean_datasinks(workflow) + fmap_cache = {} + if config.execution.derivatives: + from fmriprep.utils.bids import collect_fieldmaps + + for deriv_dir in config.execution.derivatives.values(): + fmap_cache.update( + collect_fieldmaps( + derivatives_dir=deriv_dir, + entities={'subject': subject_id}, + ) + ) + fmap_estimators, estimator_map = map_fieldmap_estimation( layout=config.execution.layout, subject_id=subject_id, From 80942f0a59c3afbcdcc9488adef01cf25b8e9692 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 5 Mar 2025 14:42:18 -0500 Subject: [PATCH 06/14] rf: Populate buffers with precomputed derivatives early --- fmriprep/workflows/bold/fit.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 9bca49abb..a58dba40f 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -304,6 +304,13 @@ def init_bold_fit_wf( niu.IdentityInterface(fields=['boldref', 'boldmask']), name='regref_buffer' ) + # Set default buffer values; can be overridden by workflows + hmcref_buffer.inputs.boldref = precomputed.get('hmc_boldref') + fmapref_buffer.inputs.sbref_files = sbref_files + hmc_buffer.inputs.hmc_xforms = hmc_xforms + fmapreg_buffer.inputs.boldref2fmap_xfm = boldref2fmap_xform + regref_buffer.inputs.boldref = precomputed.get('coreg_boldref') + summary = pe.Node( FunctionalSummary( distortion_correction='None', # Can override with connection @@ -342,7 +349,6 @@ def init_bold_fit_wf( output_dir=config.execution.fmriprep_dir, ) - # fmt:off workflow.connect([ (hmcref_buffer, outputnode, [ ('boldref', 'hmc_boldref'), @@ -371,8 +377,7 @@ def init_bold_fit_wf( ('boldref2anat_xfm', 'inputnode.boldref2anat_xfm'), ]), (summary, func_fit_reports_wf, [('out_report', 'inputnode.summary_report')]), - ]) - # fmt:on + ]) # fmt:skip # Stage 1: Generate motion correction boldref hmc_boldref_source_buffer = pe.Node( @@ -413,7 +418,6 @@ def init_bold_fit_wf( ]) # fmt:skip else: config.loggers.workflow.info('Found HMC boldref - skipping Stage 1') - hmcref_buffer.inputs.boldref = precomputed['hmc_boldref'] validation_and_dummies_wf = init_validation_and_dummies_wf(bold_file=bold_file) @@ -451,7 +455,6 @@ def init_bold_fit_wf( ]) # fmt:skip else: config.loggers.workflow.info('Found motion correction transforms - skipping Stage 2') - hmc_buffer.inputs.hmc_xforms = hmc_xforms # Stage 3: Create coregistration reference # Fieldmap correction only happens during fit if this stage is needed @@ -459,7 +462,6 @@ def init_bold_fit_wf( config.loggers.workflow.info('Stage 3: Adding coregistration boldref workflow') # Select initial boldref, enhance contrast, and generate mask - fmapref_buffer.inputs.sbref_files = sbref_files if sbref_files and nb.load(sbref_files[0]).ndim > 3: raw_sbref_wf = init_raw_boldref_wf( name='raw_sbref_wf', @@ -523,7 +525,6 @@ def init_bold_fit_wf( ) ds_fmapreg_wf.inputs.inputnode.source_files = [bold_file] - # fmt:off workflow.connect([ (enhance_boldref_wf, fmapreg_wf, [ ('outputnode.bias_corrected_file', 'inputnode.target_ref'), @@ -536,10 +537,7 @@ def init_bold_fit_wf( (fmapreg_wf, itk_mat2txt, [('outputnode.target2fmap_xfm', 'in_xfms')]), (itk_mat2txt, ds_fmapreg_wf, [('out_xfm', 'inputnode.xform')]), (ds_fmapreg_wf, fmapreg_buffer, [('outputnode.xform', 'boldref2fmap_xfm')]), - ]) - # fmt:on - else: - fmapreg_buffer.inputs.boldref2fmap_xfm = boldref2fmap_xform + ]) # fmt:skip unwarp_wf = init_unwarp_wf( free_mem=config.environment.free_mem, @@ -598,7 +596,6 @@ def init_bold_fit_wf( ]) # fmt:skip else: config.loggers.workflow.info('Found coregistration reference - skipping Stage 3') - regref_buffer.inputs.boldref = precomputed['coreg_boldref'] # TODO: Allow precomputed bold masks to be passed # Also needs consideration for how it interacts above From f41f77ceb654f9de0f5b424537a5fbaf0ff98a9d Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Wed, 5 Mar 2025 14:51:10 -0500 Subject: [PATCH 07/14] rf: Collate pre-computed and to-compute fieldmaps --- fmriprep/workflows/base.py | 89 +++++++++++++++++++++++++++++--------- 1 file changed, 69 insertions(+), 20 deletions(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index f9229afca..c0c1b1526 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -30,6 +30,7 @@ """ import os +import re import sys import warnings from copy import deepcopy @@ -555,11 +556,61 @@ def init_single_subject_wf(subject_id: str): filters=config.execution.get().get('bids_filters', {}).get('fmap'), ) + fmap_merge_nodes = { + field: pe.Node(niu.Merge(2), name=f'{field}_merge', run_without_submitting=True) + for field in ['fmap', 'fmap_ref', 'fmap_coeff', 'fmap_mask', 'fmap_id', 'sdc_method'] + } + + fmap_buffer = pe.Node( + niu.IdentityInterface( + fields=['fmap', 'fmap_ref', 'fmap_coeff', 'fmap_mask', 'fmap_id', 'sdc_method'] + ), + name='fmap_buffer', + ) + + workflow.connect( + [(fmap_merge_nodes[field], fmap_buffer, [('out', field)]) for field in fmap_merge_nodes] + ) + + missing_estimators = [] if fmap_estimators: + # Map fmapid entity to internal bids_id + fmapid_map = {re.sub(r'[^a-zA-Z0-9]', '', e.bids_id): e.bids_id for e in fmap_estimators} + pared_cache = { + fmapid_map[fmapid]: fmap_cache[fmapid] for fmapid in fmap_cache if fmapid in fmapid_map + } + + missing_estimators = [e for e in fmap_estimators if e.bids_id not in pared_cache] + + if pared_cache: + config.loggers.workflow.info( + 'Precomputed B0 field inhomogeneity maps found for the following ' + f'{len(pared_cache)} estimator(s): {list(pared_cache)}.' + ) + + fmap_merge_nodes['fmap'].inputs.in1 = [ + pared_cache[fmapid]['fieldmap'] for fmapid in pared_cache + ] + fmap_merge_nodes['fmap_ref'].inputs.in1 = [ + pared_cache[fmapid]['magnitude'] for fmapid in pared_cache + ] + fmap_merge_nodes['fmap_coeff'].inputs.in1 = [ + pared_cache[fmapid]['coeffs'] for fmapid in pared_cache + ] + # Note that masks are not emitted. The BOLD-fmap transforms cannot be + # computed with precomputed fieldmaps until we either start emitting masks + # or start skull-stripping references on the fly. + fmap_merge_nodes['fmap_mask'].inputs.in1 = [ + pared_cache[fmapid].get('mask', 'MISSING') for fmapid in pared_cache + ] + fmap_merge_nodes['fmap_id'].inputs.in1 = list(pared_cache) + fmap_merge_nodes['sdc_method'].inputs.in1 = ['precomputed'] * len(pared_cache) + + if missing_estimators: config.loggers.workflow.info( 'B0 field inhomogeneity map will be estimated with the following ' - f'{len(fmap_estimators)} estimator(s): ' - f'{[e.method for e in fmap_estimators]}.' + f'{len(missing_estimators)} estimator(s): ' + f'{[e.method for e in missing_estimators]}.' ) from sdcflows import fieldmaps as fm @@ -585,24 +636,31 @@ def init_single_subject_wf(subject_id: str): if node.split('.')[-1].startswith('ds_'): fmap_wf.get_node(node).interface.out_path_base = '' + workflow.connect([ + (fmap_wf, fmap_merge_nodes[field], [ + # We get "sdc_method" as "method" from estimator workflows + # All else stays the same, and no other sdc_ prefixes are used + (f'outputnode.{field.removeprefix("sdc_")}', 'in2'), + ]) + for field in fmap_merge_nodes + ]) # fmt:skip + fmap_select_std = pe.Node( KeySelect(fields=['std2anat_xfm'], key='MNI152NLin2009cAsym'), name='fmap_select_std', run_without_submitting=True, ) if any(estimator.method == fm.EstimatorType.ANAT for estimator in fmap_estimators): - # fmt:off workflow.connect([ (anat_fit_wf, fmap_select_std, [ ('outputnode.std2anat_xfm', 'std2anat_xfm'), ('outputnode.template', 'keys')]), - ]) - # fmt:on + ]) # fmt:skip for estimator in fmap_estimators: config.loggers.workflow.info( f"""\ -Setting-up fieldmap "{estimator.bids_id}" ({estimator.method}) with \ +Setting up fieldmap "{estimator.bids_id}" ({estimator.method}) with \ <{', '.join(s.path.name for s in estimator.sources)}>""" ) @@ -645,7 +703,6 @@ def init_single_subject_wf(subject_id: str): syn_preprocessing_wf.inputs.inputnode.in_epis = sources syn_preprocessing_wf.inputs.inputnode.in_meta = source_meta - # fmt:off workflow.connect([ (anat_fit_wf, syn_preprocessing_wf, [ ('outputnode.t1w_preproc', 'inputnode.in_anat'), @@ -661,8 +718,7 @@ def init_single_subject_wf(subject_id: str): ('outputnode.anat_mask', f'in_{estimator.bids_id}.anat_mask'), ('outputnode.sd_prior', f'in_{estimator.bids_id}.sd_prior'), ]), - ]) - # fmt:on + ]) # fmt:skip # Append the functional section to the existing anatomical excerpt # That way we do not need to stream down the number of bold datasets @@ -729,18 +785,11 @@ def init_single_subject_wf(subject_id: str): 'inputnode.sphere_reg_fsLR', ), ]), + (fmap_buffer, bold_wf, [ + (field, f'inputnode.{field}') + for field in fmap_merge_nodes + ]), ]) # fmt:skip - if fieldmap_id: - workflow.connect([ - (fmap_wf, bold_wf, [ - ('outputnode.fmap', 'inputnode.fmap'), - ('outputnode.fmap_ref', 'inputnode.fmap_ref'), - ('outputnode.fmap_coeff', 'inputnode.fmap_coeff'), - ('outputnode.fmap_mask', 'inputnode.fmap_mask'), - ('outputnode.fmap_id', 'inputnode.fmap_id'), - ('outputnode.method', 'inputnode.sdc_method'), - ]), - ]) # fmt:skip if config.workflow.level == 'full': if template_iterator_wf is not None: From eb1bb73c99a48c0c9d988319d92093be8f9919c3 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 18 Mar 2025 16:02:27 -0400 Subject: [PATCH 08/14] fix: Check fieldmap reference orientation when interpolating --- fmriprep/interfaces/resampling.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py index de5c7fcf7..603b901fc 100644 --- a/fmriprep/interfaces/resampling.py +++ b/fmriprep/interfaces/resampling.py @@ -671,6 +671,15 @@ def reconstruct_fieldmap( target.__class__(target.dataobj, projected_affine, target.header), ) else: + # Hack. Sometimes the reference array is rotated relative to the fieldmap + # and coefficient grids. As far as I know, coefficients are always RAS, + # but good to check before doing this. + if ( + nb.aff2axcodes(coefficients[-1].affine) + == ('R', 'A', 'S') + != nb.aff2axcodes(fmap_reference.affine) + ): + fmap_reference = nb.as_closest_canonical(fmap_reference) if not aligned(fmap_reference.affine, coefficients[-1].affine): raise ValueError('Reference passed is not aligned with spline grids') reference, _ = ensure_positive_cosines(fmap_reference) From a47b5a0b5891efa5cdfba4030927e45bfaf34178 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 18 Mar 2025 18:58:03 -0400 Subject: [PATCH 09/14] rf: Drop unnecessary IdentityInterface node --- fmriprep/workflows/base.py | 38 ++++++++++++++------------------------ 1 file changed, 14 insertions(+), 24 deletions(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index c0c1b1526..9c79bd307 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -556,22 +556,11 @@ def init_single_subject_wf(subject_id: str): filters=config.execution.get().get('bids_filters', {}).get('fmap'), ) - fmap_merge_nodes = { + fmap_buffers = { field: pe.Node(niu.Merge(2), name=f'{field}_merge', run_without_submitting=True) for field in ['fmap', 'fmap_ref', 'fmap_coeff', 'fmap_mask', 'fmap_id', 'sdc_method'] } - fmap_buffer = pe.Node( - niu.IdentityInterface( - fields=['fmap', 'fmap_ref', 'fmap_coeff', 'fmap_mask', 'fmap_id', 'sdc_method'] - ), - name='fmap_buffer', - ) - - workflow.connect( - [(fmap_merge_nodes[field], fmap_buffer, [('out', field)]) for field in fmap_merge_nodes] - ) - missing_estimators = [] if fmap_estimators: # Map fmapid entity to internal bids_id @@ -588,23 +577,23 @@ def init_single_subject_wf(subject_id: str): f'{len(pared_cache)} estimator(s): {list(pared_cache)}.' ) - fmap_merge_nodes['fmap'].inputs.in1 = [ + fmap_buffers['fmap'].inputs.in1 = [ pared_cache[fmapid]['fieldmap'] for fmapid in pared_cache ] - fmap_merge_nodes['fmap_ref'].inputs.in1 = [ + fmap_buffers['fmap_ref'].inputs.in1 = [ pared_cache[fmapid]['magnitude'] for fmapid in pared_cache ] - fmap_merge_nodes['fmap_coeff'].inputs.in1 = [ + fmap_buffers['fmap_coeff'].inputs.in1 = [ pared_cache[fmapid]['coeffs'] for fmapid in pared_cache ] # Note that masks are not emitted. The BOLD-fmap transforms cannot be # computed with precomputed fieldmaps until we either start emitting masks # or start skull-stripping references on the fly. - fmap_merge_nodes['fmap_mask'].inputs.in1 = [ + fmap_buffers['fmap_mask'].inputs.in1 = [ pared_cache[fmapid].get('mask', 'MISSING') for fmapid in pared_cache ] - fmap_merge_nodes['fmap_id'].inputs.in1 = list(pared_cache) - fmap_merge_nodes['sdc_method'].inputs.in1 = ['precomputed'] * len(pared_cache) + fmap_buffers['fmap_id'].inputs.in1 = list(pared_cache) + fmap_buffers['sdc_method'].inputs.in1 = ['precomputed'] * len(pared_cache) if missing_estimators: config.loggers.workflow.info( @@ -637,12 +626,12 @@ def init_single_subject_wf(subject_id: str): fmap_wf.get_node(node).interface.out_path_base = '' workflow.connect([ - (fmap_wf, fmap_merge_nodes[field], [ + (fmap_wf, fmap_buffers[field], [ # We get "sdc_method" as "method" from estimator workflows # All else stays the same, and no other sdc_ prefixes are used (f'outputnode.{field.removeprefix("sdc_")}', 'in2'), ]) - for field in fmap_merge_nodes + for field in fmap_buffers ]) # fmt:skip fmap_select_std = pe.Node( @@ -785,10 +774,11 @@ def init_single_subject_wf(subject_id: str): 'inputnode.sphere_reg_fsLR', ), ]), - (fmap_buffer, bold_wf, [ - (field, f'inputnode.{field}') - for field in fmap_merge_nodes - ]), + ]) # fmt:skip + + workflow.connect([ + (buffer, bold_wf, [('out', f'inputnode.{field}')]) + for field, buffer in fmap_buffers.items() ]) # fmt:skip if config.workflow.level == 'full': From 68f20790740843dba76deb0beafb776e4504ae48 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 20 Mar 2025 11:36:38 -0400 Subject: [PATCH 10/14] rf: Rework pared_cache/fmap_estimator logic --- fmriprep/workflows/base.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 9c79bd307..10ec80801 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -546,7 +546,7 @@ def init_single_subject_wf(subject_id: str): ) ) - fmap_estimators, estimator_map = map_fieldmap_estimation( + all_estimators, estimator_map = map_fieldmap_estimation( layout=config.execution.layout, subject_id=subject_id, bold_data=bold_runs, @@ -561,15 +561,16 @@ def init_single_subject_wf(subject_id: str): for field in ['fmap', 'fmap_ref', 'fmap_coeff', 'fmap_mask', 'fmap_id', 'sdc_method'] } - missing_estimators = [] - if fmap_estimators: - # Map fmapid entity to internal bids_id - fmapid_map = {re.sub(r'[^a-zA-Z0-9]', '', e.bids_id): e.bids_id for e in fmap_estimators} - pared_cache = { - fmapid_map[fmapid]: fmap_cache[fmapid] for fmapid in fmap_cache if fmapid in fmapid_map - } - - missing_estimators = [e for e in fmap_estimators if e.bids_id not in pared_cache] + fmap_estimators = [] + if all_estimators: + # Find precomputed fieldmaps that apply to this workflow + pared_cache = {} + for est in all_estimators: + fmapid = re.sub(r'[^a-zA-Z0-9]', '', est.bids_id) + if fmapid in fmap_cache: + pared_cache[fmapid] = fmap_cache[fmapid] + else: + fmap_estimators.append(est) if pared_cache: config.loggers.workflow.info( @@ -595,11 +596,12 @@ def init_single_subject_wf(subject_id: str): fmap_buffers['fmap_id'].inputs.in1 = list(pared_cache) fmap_buffers['sdc_method'].inputs.in1 = ['precomputed'] * len(pared_cache) - if missing_estimators: + # Estimators without precomputed fieldmaps + if fmap_estimators: config.loggers.workflow.info( 'B0 field inhomogeneity map will be estimated with the following ' - f'{len(missing_estimators)} estimator(s): ' - f'{[e.method for e in missing_estimators]}.' + f'{len(fmap_estimators)} estimator(s): ' + f'{[e.method for e in fmap_estimators]}.' ) from sdcflows import fieldmaps as fm From 24b04e370e3e13383725c1826b25e534840b2f7a Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 20 Mar 2025 12:03:42 -0400 Subject: [PATCH 11/14] feat: Add debugging output to indicate what files are reused --- fmriprep/workflows/base.py | 35 ++++++++++++++++++++-------------- fmriprep/workflows/bold/fit.py | 29 +++++++++++++++++----------- 2 files changed, 39 insertions(+), 25 deletions(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 10ec80801..22cc86b17 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -338,6 +338,9 @@ def init_single_subject_wf(subject_id: str): # allow to run with anat-fast-track on fMRI-only dataset if 't1w_preproc' in anatomical_cache and not subject_data['t1w']: + workflow.debug( + 'No T1w image found; using precomputed T1w image: %s', anatomical_cache['t1w_preproc'] + ) workflow.connect([ (bidssrc, bids_info, [(('bold', fix_multi_T1w_source_name), 'in_file')]), (anat_fit_wf, summary, [('outputnode.t1w_preproc', 't1w')]), @@ -539,12 +542,14 @@ def init_single_subject_wf(subject_id: str): from fmriprep.utils.bids import collect_fieldmaps for deriv_dir in config.execution.derivatives.values(): - fmap_cache.update( - collect_fieldmaps( - derivatives_dir=deriv_dir, - entities={'subject': subject_id}, - ) + fmaps = collect_fieldmaps( + derivatives_dir=deriv_dir, + entities={'subject': subject_id}, ) + config.loggers.workflow.debug( + 'Detected precomputed fieldmaps in %d for fieldmap IDs: %s', deriv_dir, list(fmaps) + ) + fmap_cache.update(fmaps) all_estimators, estimator_map = map_fieldmap_estimation( layout=config.execution.layout, @@ -578,15 +583,17 @@ def init_single_subject_wf(subject_id: str): f'{len(pared_cache)} estimator(s): {list(pared_cache)}.' ) - fmap_buffers['fmap'].inputs.in1 = [ - pared_cache[fmapid]['fieldmap'] for fmapid in pared_cache - ] - fmap_buffers['fmap_ref'].inputs.in1 = [ - pared_cache[fmapid]['magnitude'] for fmapid in pared_cache - ] - fmap_buffers['fmap_coeff'].inputs.in1 = [ - pared_cache[fmapid]['coeffs'] for fmapid in pared_cache - ] + fieldmaps = [fmap['fieldmap'] for fmap in pared_cache.values()] + refs = [fmap['magnitude'] for fmap in pared_cache.values()] + coeffs = [fmap['coeffs'] for fmap in pared_cache.values()] + config.loggers.workflow.debug('Reusing fieldmaps: %s', fieldmaps) + config.loggers.workflow.debug('Reusing references: %s', refs) + config.loggers.workflow.debug('Reusing coefficients: %s', coeffs) + + fmap_buffers['fmap'].inputs.in1 = fieldmaps + fmap_buffers['fmap_ref'].inputs.in1 = refs + fmap_buffers['fmap_coeff'].inputs.in1 = coeffs + # Note that masks are not emitted. The BOLD-fmap transforms cannot be # computed with precomputed fieldmaps until we either start emitting masks # or start skull-stripping references on the fly. diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index a58dba40f..0445e4de3 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -236,8 +236,8 @@ def init_bold_fit_wf( # Boolean used to update workflow self-descriptions multiecho = len(bold_series) > 1 - have_hmcref = 'hmc_boldref' in precomputed - have_coregref = 'coreg_boldref' in precomputed + hmc_boldref = precomputed.get('hmc_boldref') + coreg_boldref = precomputed.get('coreg_boldref') # Can contain # 1) boldref2fmap # 2) boldref2anat @@ -304,12 +304,19 @@ def init_bold_fit_wf( niu.IdentityInterface(fields=['boldref', 'boldmask']), name='regref_buffer' ) - # Set default buffer values; can be overridden by workflows - hmcref_buffer.inputs.boldref = precomputed.get('hmc_boldref') + if hmc_boldref: + hmcref_buffer.inputs.boldref = hmc_boldref + config.loggers.workflow.debug('Reusing motion correction reference: %s', hmc_boldref) + if hmc_xforms: + hmc_buffer.inputs.hmc_xforms = hmc_xforms + config.loggers.workflow.debug('Reusing motion correction transforms: %s', hmc_xforms) + if boldref2fmap_xform: + fmapreg_buffer.inputs.boldref2fmap_xfm = boldref2fmap_xform + config.loggers.workflow.debug('Reusing BOLD-to-fieldmap transform: %s', boldref2fmap_xform) + if coreg_boldref: + regref_buffer.inputs.boldref = coreg_boldref + config.loggers.workflow.debug('Reusing coregistration reference: %s', coreg_boldref) fmapref_buffer.inputs.sbref_files = sbref_files - hmc_buffer.inputs.hmc_xforms = hmc_xforms - fmapreg_buffer.inputs.boldref2fmap_xfm = boldref2fmap_xform - regref_buffer.inputs.boldref = precomputed.get('coreg_boldref') summary = pe.Node( FunctionalSummary( @@ -344,7 +351,7 @@ def init_bold_fit_wf( func_fit_reports_wf = init_func_fit_reports_wf( # TODO: Enable sdc report even if we find coregref - sdc_correction=not (have_coregref or fieldmap_id is None), + sdc_correction=not (coreg_boldref or fieldmap_id is None), freesurfer=config.workflow.run_reconall, output_dir=config.execution.fmriprep_dir, ) @@ -384,7 +391,7 @@ def init_bold_fit_wf( niu.IdentityInterface(fields=['in_file']), name='hmc_boldref_source_buffer', ) - if not have_hmcref: + if not hmc_boldref: config.loggers.workflow.info('Stage 1: Adding HMC boldref workflow') hmc_boldref_wf = init_raw_boldref_wf( name='hmc_boldref_wf', @@ -458,7 +465,7 @@ def init_bold_fit_wf( # Stage 3: Create coregistration reference # Fieldmap correction only happens during fit if this stage is needed - if not have_coregref: + if not coreg_boldref: config.loggers.workflow.info('Stage 3: Adding coregistration boldref workflow') # Select initial boldref, enhance contrast, and generate mask @@ -600,7 +607,7 @@ def init_bold_fit_wf( # TODO: Allow precomputed bold masks to be passed # Also needs consideration for how it interacts above skullstrip_precomp_ref_wf = init_skullstrip_bold_wf(name='skullstrip_precomp_ref_wf') - skullstrip_precomp_ref_wf.inputs.inputnode.in_file = precomputed['coreg_boldref'] + skullstrip_precomp_ref_wf.inputs.inputnode.in_file = coreg_boldref workflow.connect([ (skullstrip_precomp_ref_wf, regref_buffer, [('outputnode.mask_file', 'boldmask')]) ]) # fmt:skip From 1765d4ac3a2a8986e2ef0a8c30bbc09aa22193cb Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Thu, 20 Mar 2025 12:20:50 -0400 Subject: [PATCH 12/14] fix: Key pared_cache on bids_id --- fmriprep/workflows/base.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index 22cc86b17..bad23bcd0 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -571,9 +571,8 @@ def init_single_subject_wf(subject_id: str): # Find precomputed fieldmaps that apply to this workflow pared_cache = {} for est in all_estimators: - fmapid = re.sub(r'[^a-zA-Z0-9]', '', est.bids_id) - if fmapid in fmap_cache: - pared_cache[fmapid] = fmap_cache[fmapid] + if found := fmap_cache.get(re.sub(r'[^a-zA-Z0-9]', '', est.bids_id)): + pared_cache[est.bids_id] = found else: fmap_estimators.append(est) From a43ea0dca19eaf755276018ba937ff0f6b7be494 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 20 Mar 2025 12:42:55 -0400 Subject: [PATCH 13/14] Update fmriprep/workflows/base.py Co-authored-by: Mathias Goncalves --- fmriprep/workflows/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index bad23bcd0..e02959a5d 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -547,7 +547,7 @@ def init_single_subject_wf(subject_id: str): entities={'subject': subject_id}, ) config.loggers.workflow.debug( - 'Detected precomputed fieldmaps in %d for fieldmap IDs: %s', deriv_dir, list(fmaps) + 'Detected precomputed fieldmaps in %s for fieldmap IDs: %s', deriv_dir, list(fmaps) ) fmap_cache.update(fmaps) From 00fee362da67eef0b5648dcd6f3eaab4b4ee9e96 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 20 Mar 2025 14:59:21 -0400 Subject: [PATCH 14/14] Update fmriprep/workflows/base.py --- fmriprep/workflows/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py index e02959a5d..5cf07de64 100644 --- a/fmriprep/workflows/base.py +++ b/fmriprep/workflows/base.py @@ -338,7 +338,7 @@ def init_single_subject_wf(subject_id: str): # allow to run with anat-fast-track on fMRI-only dataset if 't1w_preproc' in anatomical_cache and not subject_data['t1w']: - workflow.debug( + config.loggers.workflow.debug( 'No T1w image found; using precomputed T1w image: %s', anatomical_cache['t1w_preproc'] ) workflow.connect([