diff --git a/src/fmripost_aroma/data/io_spec.json b/src/fmripost_aroma/data/io_spec.json index 8b635a7..d75a6bf 100644 --- a/src/fmripost_aroma/data/io_spec.json +++ b/src/fmripost_aroma/data/io_spec.json @@ -1,4 +1,5 @@ { + "name": "aroma", "queries": { "raw": { "bold_raw": { @@ -151,8 +152,8 @@ }, "entities": [ { - "name": "datatype", - "pattern": "[/\\\\]+(anat|func)[/\\\\]+" + "name": "fmapid", + "pattern": "fmapid-([a-zA-Z0-9]+)" }, { "name": "cohort", @@ -172,7 +173,7 @@ "pattern": "(?:^|_)thresh-([a-zA-Z0-9]+)" } ], - "patterns": [ + "default_path_patterns": [ "sub-{subject}[/ses-{session}]/{datatype|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_part-{part}][_space-{space}][_res-{res}][_desc-{desc}]_{suffix}.{extension|nii.gz}", "sub-{subject}[/ses-{session}]/{datatype|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_part-{part}][_space-{space}][_res-{res}][_stat-{statistic}][_desc-{desc}]_{suffix}.{extension|nii.gz}", "sub-{subject}[/ses-{session}]/{datatype|func}/sub-{subject}[_ses-{session}]_task-{task}[_acq-{acquisition}][_ce-{ceagent}][_rec-{reconstruction}][_run-{run}][_echo-{echo}][_part-{part}][_space-{space}][_res-{res}][_stat-{statistic}][_desc-{desc}]_{suffix}.{extension|tsv}", diff --git a/src/fmripost_aroma/interfaces/bids.py b/src/fmripost_aroma/interfaces/bids.py index d90ead0..e36e0e3 100644 --- a/src/fmripost_aroma/interfaces/bids.py +++ b/src/fmripost_aroma/interfaces/bids.py @@ -31,4 +31,4 @@ class DerivativesDataSink(BaseDerivativesDataSink): _allowed_entities = set(config_entities) _config_entities = config_entities _config_entities_dict = merged_entities - _file_patterns = fmripost_aroma_spec['patterns'] + _file_patterns = fmripost_aroma_spec['default_path_patterns'] diff --git a/src/fmripost_aroma/utils/bids.py b/src/fmripost_aroma/utils/bids.py index bbe0701..3a45ea7 100644 --- a/src/fmripost_aroma/utils/bids.py +++ b/src/fmripost_aroma/utils/bids.py @@ -98,6 +98,7 @@ def collect_derivatives( if not entities: entities = {} + _spec = None if spec is None or patterns is None: _spec = json.loads(load_data.readable('io_spec.json').read_text()) @@ -105,7 +106,13 @@ def collect_derivatives( spec = _spec['queries'] if patterns is None: - patterns = _spec['patterns'] + patterns = _spec['default_path_patterns'] + + _spec.pop('queries') + + config = ['bids', 'derivatives'] + if _spec: + config = ['bids', 'derivatives', _spec] # Search for derivatives data derivs_cache = defaultdict(list, {}) @@ -114,7 +121,7 @@ def collect_derivatives( if isinstance(layout, Path): layout = BIDSLayout( layout, - config=['bids', 'derivatives'], + config=config, validate=False, ) diff --git a/src/fmripost_aroma/workflows/base.py b/src/fmripost_aroma/workflows/base.py index 480db47..9bd34b0 100644 --- a/src/fmripost_aroma/workflows/base.py +++ b/src/fmripost_aroma/workflows/base.py @@ -179,6 +179,7 @@ def init_single_subject_wf(subject_id: str): if config.execution.derivatives: # Raw dataset + derivatives dataset config.loggers.workflow.info('Raw+derivatives workflow mode enabled') + # Just build a list of BOLD files right now subject_data = collect_derivatives( raw_dataset=config.execution.layout, derivatives_dataset=None, @@ -191,6 +192,7 @@ def init_single_subject_wf(subject_id: str): else: # Derivatives dataset only config.loggers.workflow.info('Derivatives-only workflow mode enabled') + # Just build a list of BOLD files right now subject_data = collect_derivatives( raw_dataset=None, derivatives_dataset=config.execution.layout, @@ -304,6 +306,24 @@ def init_single_run_wf(bold_file): entities = extract_entities(bold_file) + # Attempt to extract the associated fmap ID + fmapid = None + all_fmapids = config.execution.layout.get_fmapids( + subject=entities['subject'], + session=entities.get('session', None), + ) + if all_fmapids: + fmap_file = config.execution.layout.get_nearest( + bold_file, + to=all_fmapids, + suffix='xfm', + extension='.txt', + strict=False, + **{'from': 'boldref'}, + ) + if fmap_file: + fmapid = config.execution.layout.get_file(fmap_file).entities['to'] + functional_cache = defaultdict(list, {}) if config.execution.derivatives: # Collect native-space derivatives and transforms @@ -311,7 +331,7 @@ def init_single_run_wf(bold_file): raw_dataset=config.execution.layout, derivatives_dataset=None, entities=entities, - fieldmap_id=None, + fieldmap_id=fmapid, allow_multiple=False, spaces=None, ) @@ -322,7 +342,7 @@ def init_single_run_wf(bold_file): raw_dataset=None, derivatives_dataset=deriv_dir, entities=entities, - fieldmap_id=None, + fieldmap_id=fmapid, allow_multiple=False, spaces=spaces, ), @@ -347,7 +367,7 @@ def init_single_run_wf(bold_file): raw_dataset=None, derivatives_dataset=config.execution.layout, entities=entities, - fieldmap_id=None, + fieldmap_id=fmapid, allow_multiple=False, spaces=spaces, ), @@ -371,30 +391,60 @@ def init_single_run_wf(bold_file): ) skip_vols = get_nss(functional_cache['bold_confounds']) + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + 'bold_raw', + 'bold_confounds', + 'bold_mni152nlin6asym', + 'motion_xfm', + 'boldref2fmap_xfm', + 'boldref2anat_xfm', + 'anat2std_xfm', + 'bold_ref_file', + 'fmap', + 'bold_mask_native', + ], + ), + name='inputnode', + ) + inputnode.inputs.bold_raw = functional_cache['bold_raw'] + inputnode.inputs.bold_confounds = functional_cache['bold_confounds'] + inputnode.inputs.bold_mni152nlin6asym = functional_cache['bold_mni152nlin6asym'] + inputnode.inputs.bold_mask_native = functional_cache['bold_mask_native'] + # Transforms + inputnode.inputs.bold_hmc = functional_cache['bold_hmc'] + inputnode.inputs.boldref2fmap = functional_cache['boldref2fmap'] + inputnode.inputs.boldref2anat = functional_cache['boldref2anat'] + inputnode.inputs.anat2mni152nlin6asym = functional_cache['anat2mni152nlin6asym'] + # Field maps + inputnode.inputs.fmap = functional_cache['fmap'] + # Run ICA-AROMA ica_aroma_wf = init_ica_aroma_wf(bold_file=bold_file, metadata=bold_metadata, mem_gb=mem_gb) - ica_aroma_wf.inputs.inputnode.confounds = functional_cache['bold_confounds'] ica_aroma_wf.inputs.inputnode.skip_vols = skip_vols + workflow.connect([(inputnode, ica_aroma_wf, [('bold_confounds', 'inputnode.confounds')])]) mni6_buffer = pe.Node(niu.IdentityInterface(fields=['bold', 'bold_mask']), name='mni6_buffer') if ('bold_mni152nlin6asym' not in functional_cache) and ('bold_raw' in functional_cache): # Resample to MNI152NLin6Asym:res-2, for ICA-AROMA classification - from fmriprep.workflows.bold.apply import init_bold_volumetric_resample_wf from fmriprep.workflows.bold.stc import init_bold_stc_wf from niworkflows.interfaces.header import ValidateImage from templateflow.api import get as get_template from fmripost_aroma.interfaces.misc import ApplyTransforms + from fmripost_aroma.workflows.resample import init_bold_volumetric_resample_wf workflow.__desc__ += """\ Raw BOLD series were resampled to MNI152NLin6Asym:res-2, for ICA-AROMA classification. """ validate_bold = pe.Node( - ValidateImage(in_file=functional_cache['bold_raw']), + ValidateImage(), name='validate_bold', ) + workflow.connect([(inputnode, validate_bold, [('bold_raw', 'in_file')])]) stc_buffer = pe.Node( niu.IdentityInterface(fields=['bold_file']), @@ -427,30 +477,27 @@ def init_single_run_wf(bold_file): bold_MNI6_wf = init_bold_volumetric_resample_wf( metadata=bold_metadata, - fieldmap_id=None, # XXX: Ignoring the field map for now + fieldmap_id=fmapid, omp_nthreads=omp_nthreads, mem_gb=mem_gb, jacobian='fmap-jacobian' not in config.workflow.ignore, name='bold_MNI6_wf', ) - bold_MNI6_wf.inputs.inputnode.motion_xfm = functional_cache['bold_hmc'] - bold_MNI6_wf.inputs.inputnode.boldref2fmap_xfm = functional_cache['boldref2fmap'] - bold_MNI6_wf.inputs.inputnode.boldref2anat_xfm = functional_cache['boldref2anat'] - bold_MNI6_wf.inputs.inputnode.anat2std_xfm = functional_cache['anat2mni152nlin6asym'] bold_MNI6_wf.inputs.inputnode.resolution = '02' - # use mask as boldref? - bold_MNI6_wf.inputs.inputnode.bold_ref_file = functional_cache['bold_mask_native'] bold_MNI6_wf.inputs.inputnode.target_mask = mni6_mask bold_MNI6_wf.inputs.inputnode.target_ref_file = mni6_mask + # Pass transforms and warped fieldmap to bold_MNI6_wf workflow.connect([ + (inputnode, bold_MNI6_wf, [ + ('bold_hmc', 'inputnode.motion_xfm'), + ('boldref2fmap', 'inputnode.boldref2fmap_xfm'), + ('boldref2anat', 'inputnode.boldref2anat_xfm'), + ('anat2mni152nlin6asym', 'inputnode.anat2std_xfm'), + # use mask as boldref? + ('bold_mask_native', 'inputnode.bold_ref_file'), + ]), # Resample BOLD to MNI152NLin6Asym, may duplicate bold_std_wf above - # XXX: Ignoring the field map for now - # (inputnode, bold_MNI6_wf, [ - # ('fmap_ref', 'inputnode.fmap_ref'), - # ('fmap_coeff', 'inputnode.fmap_coeff'), - # ('fmap_id', 'inputnode.fmap_id'), - # ]), (stc_buffer, bold_MNI6_wf, [('bold_file', 'inputnode.bold_file')]), (bold_MNI6_wf, mni6_buffer, [('outputnode.bold_file', 'bold')]), ]) # fmt:skip @@ -459,7 +506,6 @@ def init_single_run_wf(bold_file): mask_to_mni6 = pe.Node( ApplyTransforms( interpolation='GenericLabel', - input_image=functional_cache['bold_mask_native'], reference_image=mni6_mask, transforms=[ functional_cache['anat2mni152nlin6asym'], @@ -468,15 +514,22 @@ def init_single_run_wf(bold_file): ), name='mask_to_mni6', ) - workflow.connect([(mask_to_mni6, mni6_buffer, [('output_image', 'bold_mask')])]) + workflow.connect([ + (inputnode, mask_to_mni6, [('bold_mask_native', 'input_image')]), + (mask_to_mni6, mni6_buffer, [('output_image', 'bold_mask')]), + ]) # fmt:skip elif 'bold_mni152nlin6asym' in functional_cache: workflow.__desc__ += """\ Preprocessed BOLD series in MNI152NLin6Asym:res-2 space were collected for ICA-AROMA classification. """ - mni6_buffer.inputs.bold = functional_cache['bold_mni152nlin6asym'] - mni6_buffer.inputs.bold_mask = functional_cache['bold_mask_mni152nlin6asym'] + workflow.connect([ + (inputnode, mni6_buffer, [ + ('bold_mni152nlin6asym', 'bold'), + ('bold_mask_native', 'bold_mask'), + ]), + ]) # fmt:skip else: raise ValueError('No valid BOLD series found for ICA-AROMA classification.') @@ -491,9 +544,13 @@ def init_single_run_wf(bold_file): # Generate reportlets func_fit_reports_wf = init_func_fit_reports_wf(output_dir=config.execution.output_dir) func_fit_reports_wf.inputs.inputnode.source_file = bold_file - func_fit_reports_wf.inputs.inputnode.anat2std_xfm = functional_cache['anat2mni152nlin6asym'] - func_fit_reports_wf.inputs.inputnode.anat_dseg = functional_cache['anat_dseg'] - workflow.connect([(mni6_buffer, func_fit_reports_wf, [('bold', 'inputnode.bold_mni6')])]) + workflow.connect([ + (inputnode, func_fit_reports_wf, [ + ('anat2mni152nlin6asym', 'inputnode.anat2std_xfm'), + ('anat_dseg', 'inputnode.anat_dseg'), + ]), + (mni6_buffer, func_fit_reports_wf, [('bold', 'inputnode.bold_mni6')]), + ]) # fmt:skip if config.workflow.denoise_method: # Now denoise the output-space BOLD data using ICA-AROMA @@ -501,9 +558,9 @@ def init_single_run_wf(bold_file): denoise_wf.inputs.inputnode.skip_vols = skip_vols denoise_wf.inputs.inputnode.space = 'MNI152NLin6Asym' denoise_wf.inputs.inputnode.res = '2' - denoise_wf.inputs.inputnode.confounds_file = functional_cache['bold_confounds'] workflow.connect([ + (inputnode, denoise_wf, [('bold_confounds', 'inputnode.confounds_file')]), (mni6_buffer, denoise_wf, [ ('bold', 'inputnode.bold_file'), ('bold_mask', 'inputnode.bold_mask'), diff --git a/src/fmripost_aroma/workflows/resample.py b/src/fmripost_aroma/workflows/resample.py new file mode 100644 index 0000000..ef21296 --- /dev/null +++ b/src/fmripost_aroma/workflows/resample.py @@ -0,0 +1,196 @@ +from __future__ import annotations + +import nipype.interfaces.utility as niu +import nipype.pipeline.engine as pe +from fmriprep.interfaces.resampling import DistortionParameters, ResampleSeries +from niworkflows.interfaces.nibabel import GenerateSamplingReference + + +def init_bold_volumetric_resample_wf( + *, + metadata: dict, + mem_gb: dict[str, float], + jacobian: bool, + fieldmap_id: str | None = None, + omp_nthreads: int = 1, + name: str = 'bold_volumetric_resample_wf', +) -> pe.Workflow: + """Resample a BOLD series to a volumetric target space. + + This workflow collates a sequence of transforms to resample a BOLD series + in a single shot, including motion correction and fieldmap correction, if + requested. + + .. workflow:: + + from fmripost_aroma.workflows.bold.resampling import init_bold_volumetric_resample_wf + + wf = init_bold_volumetric_resample_wf( + metadata={ + 'RepetitionTime': 2.0, + 'PhaseEncodingDirection': 'j-', + 'TotalReadoutTime': 0.03 + }, + mem_gb={'resampled': 1}, + jacobian=True, + fieldmap_id='my_fieldmap', + ) + + Parameters + ---------- + metadata + BIDS metadata for BOLD file. + fieldmap_id + Fieldmap identifier, if fieldmap correction is to be applied. + omp_nthreads + Maximum number of threads an individual process may use. + name + Name of workflow (default: ``bold_volumetric_resample_wf``) + + Inputs + ------ + bold_file + BOLD series to resample. + bold_ref_file + Reference image to which BOLD series is aligned. + target_ref_file + Reference image defining the target space. + target_mask + Brain mask corresponding to ``target_ref_file``. + This is used to define the field of view for the resampled BOLD series. + motion_xfm + List of affine transforms aligning each volume to ``bold_ref_file``. + If undefined, no motion correction is performed. + fmap + Fieldmap image. + fmap_id + Fieldmap identifier, to select correct fieldmap in case there are multiple. + boldref2anat_xfm + Affine transform from ``bold_ref_file`` to the anatomical reference image. + anat2std_xfm + Affine transform from the anatomical reference image to standard space. + Leave undefined to resample to anatomical reference space. + + Outputs + ------- + bold_file + The ``bold_file`` input, resampled to ``target_ref_file`` space. + resampling_reference + An empty reference image with the correct affine and header for resampling + further images into the BOLD series' space. + + """ + workflow = pe.Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface( + fields=[ + 'bold_file', + 'bold_ref_file', + 'target_ref_file', + 'target_mask', + # HMC + 'motion_xfm', + # SDC + 'fmap', + 'fmap_id', + # Anatomical + 'boldref2anat_xfm', + # Template + 'anat2std_xfm', + # Entity for selecting target resolution + 'resolution', + ], + ), + name='inputnode', + ) + + outputnode = pe.Node( + niu.IdentityInterface(fields=['bold_file', 'resampling_reference']), + name='outputnode', + ) + + gen_ref = pe.Node(GenerateSamplingReference(), name='gen_ref', mem_gb=0.3) + + boldref2target = pe.Node(niu.Merge(2), name='boldref2target', run_without_submitting=True) + bold2target = pe.Node(niu.Merge(2), name='bold2target', run_without_submitting=True) + resample = pe.Node( + ResampleSeries(jacobian=jacobian), + name='resample', + n_procs=omp_nthreads, + mem_gb=mem_gb['resampled'], + ) + + workflow.connect([ + (inputnode, gen_ref, [ + ('bold_ref_file', 'moving_image'), + ('target_ref_file', 'fixed_image'), + ('target_mask', 'fov_mask'), + (('resolution', _is_native), 'keep_native'), + ]), + (inputnode, boldref2target, [ + ('boldref2anat_xfm', 'in1'), + ('anat2std_xfm', 'in2'), + ]), + (inputnode, bold2target, [('motion_xfm', 'in1')]), + (inputnode, resample, [('bold_file', 'in_file')]), + (gen_ref, resample, [('out_file', 'ref_file')]), + (boldref2target, bold2target, [('out', 'in2')]), + (bold2target, resample, [('out', 'transforms')]), + (gen_ref, outputnode, [('out_file', 'resampling_reference')]), + (resample, outputnode, [('out_file', 'bold_file')]), + ]) # fmt:skip + + if not fieldmap_id: + return workflow + + distortion_params = pe.Node( + DistortionParameters(metadata=metadata), + name='distortion_params', + run_without_submitting=True, + ) + fmap2target = pe.Node(niu.Merge(2), name='fmap2target', run_without_submitting=True) + inverses = pe.Node( + niu.Function(function=_gen_inverses), + name='inverses', + run_without_submitting=True, + ) + fmap_resample = pe.Node(ResampleSeries(jacobian=False), name='fmap_resample') + + workflow.connect([ + # Resample fieldmap to target space + (inputnode, fmap2target, [('boldref2fmap_xfm', 'in1')]), + (boldref2target, fmap2target, [('out', 'in2')]), + (boldref2target, inverses, [('out', 'inlist')]), + (inputnode, fmap_resample, [('fmap', 'in_file')]), + (gen_ref, fmap_resample, [('out_file', 'ref_file')]), + (fmap2target, fmap_resample, [('out', 'transforms')]), + (inverses, fmap_resample, [('out', 'inverse')]), + + # Inject fieldmap correction into resample node + (inputnode, distortion_params, [('bold_file', 'in_file')]), + (distortion_params, resample, [ + ('readout_time', 'ro_time'), + ('pe_direction', 'pe_dir'), + ]), + (fmap_resample, resample, [('out_file', 'fieldmap')]), + ]) # fmt:skip + + return workflow + + +def _gen_inverses(inlist: list) -> list[bool]: + """Create a list indicating the first transform should be inverted. + + The input list is the collection of transforms that follow the + inverted one. + """ + from niworkflows.utils.connections import listify + + if not inlist: + return [True] + return [True] + [False] * len(listify(inlist)) + + +def _is_native(value): + return value == 'native'