diff --git a/.zenodo.json b/.zenodo.json index ae3c1d682..50f5009fb 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -106,6 +106,11 @@ "affiliation": "University College London", "orcid": "0000-0002-9699-9052" }, + { + "name": "Cieslak, Matthew", + "affiliation": "Department of Neuropsychiatry, University of Pennsylvania", + "orcid": "0000-0002-1931-4734" + }, { "name": "Halchenko, Yaroslav O.", "affiliation": "Dartmouth College: Hanover, NH, United States", diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib index 8a391aee4..e9e4c8e16 100644 --- a/fmriprep/data/boilerplate.bib +++ b/fmriprep/data/boilerplate.bib @@ -306,6 +306,16 @@ @article{hcppipelines year = 2013 } +@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 others}, + journal={Neuroimage}, + volume={86}, + pages={544--553}, + year={2014}, + publisher={Elsevier} +} + @article{fs_template, author = {Reuter, Martin and Rosas, Herminia Diana and Fischl, Bruce}, doi = {10.1016/j.neuroimage.2010.07.020}, diff --git a/fmriprep/interfaces/__init__.py b/fmriprep/interfaces/__init__.py index 369d8b346..68f8badaa 100644 --- a/fmriprep/interfaces/__init__.py +++ b/fmriprep/interfaces/__init__.py @@ -6,7 +6,7 @@ bids, cifti, freesurfer, images, itk, surf, utils) from .reports import SubjectSummary, FunctionalSummary, AboutSummary -from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap +from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap, Phases2Fieldmap from .confounds import GatherConfounds, ICAConfounds, FMRISummary from .multiecho import T2SMap diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py index f4286c3ea..71c96cc25 100644 --- a/fmriprep/interfaces/fmap.py +++ b/fmriprep/interfaces/fmap.py @@ -21,7 +21,7 @@ from nipype.utils.filemanip import fname_presuffix from nipype.interfaces.base import ( BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits, - SimpleInterface) + SimpleInterface, InputMultiObject) LOGGER = logging.getLogger('nipype.interface') @@ -207,6 +207,34 @@ def _run_interface(self, runtime): return runtime +class Phases2FieldmapInputSpec(BaseInterfaceInputSpec): + phase_files = InputMultiObject( + File(exists=True), mandatory=True, desc='list of phase1, phase2 files') + metadatas = traits.List( + traits.Dict, mandatory=True, desc='list of phase1, phase2 metadata dicts') + + +class Phases2FieldmapOutputSpec(TraitedSpec): + out_file = File(desc='the output fieldmap') + phasediff_metadata = traits.Dict(desc='the phasediff metadata') + + +class Phases2Fieldmap(SimpleInterface): + """ + Convert a phase1, phase2 into a difference map + """ + input_spec = Phases2FieldmapInputSpec + output_spec = Phases2FieldmapOutputSpec + + def _run_interface(self, runtime): + # Get the echo times + fmap_file, merged_metadata = phases2fmap(self.inputs.phase_files, self.inputs.metadatas, + newpath=runtime.cwd) + self._results['phasediff_metadata'] = merged_metadata + self._results['out_file'] = fmap_file + return runtime + + def _despike2d(data, thres, neigh=None): """ despiking as done in FSL fugue @@ -497,8 +525,69 @@ def phdiff2fmap(in_file, delta_te, newpath=None): return out_file +def phases2fmap(phase_files, metadatas, newpath=None): + """Calculates a phasediff from two phase images. Assumes monopolar + readout. """ + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + from copy import deepcopy + + phasediff_file = fname_presuffix(phase_files[0], suffix='_phasediff', newpath=newpath) + echo_times = [meta.get("EchoTime") for meta in metadatas] + if None in echo_times or echo_times[0] == echo_times[1]: + raise RuntimeError() + # Determine the order of subtraction + short_echo_index = echo_times.index(min(echo_times)) + long_echo_index = echo_times.index(max(echo_times)) + + short_phase_image = phase_files[short_echo_index] + long_phase_image = phase_files[long_echo_index] + + image0 = nb.load(short_phase_image) + phase0 = image0.get_fdata() + image1 = nb.load(long_phase_image) + phase1 = image1.get_fdata() + + def rescale_image(img): + if np.any(img < -128): + # This happens sometimes on 7T fieldmaps + LOGGER.info("Found negative values in phase image: rescaling") + imax = img.max() + imin = img.min() + scaled = 2 * ((img - imin) / (imax - imin) - 0.5) + return np.pi * scaled + mask = img > 0 + imax = img.max() + imin = img.min() + max_check = imax - 4096 + if np.abs(max_check) > 10 or np.abs(imin) > 10: + LOGGER.warning("Phase image may be scaled incorrectly: check results") + return mask * (img / 2048 * np.pi - np.pi) + + # Calculate fieldmaps + rad0 = rescale_image(phase0) + rad1 = rescale_image(phase1) + a = np.cos(rad0) + b = np.sin(rad0) + c = np.cos(rad1) + d = np.sin(rad1) + fmap = -np.arctan2(b * c - a * d, a * c + b * d) + + phasediff_nii = nb.Nifti1Image(fmap, image0.affine) + phasediff_nii.set_data_dtype(np.float32) + phasediff_nii.to_filename(phasediff_file) + + merged_metadata = deepcopy(metadatas[0]) + del merged_metadata['EchoTime'] + merged_metadata['EchoTime1'] = float(echo_times[short_echo_index]) + merged_metadata['EchoTime2'] = float(echo_times[long_echo_index]) + + return phasediff_file, merged_metadata + + def _delta_te(in_values, te1=None, te2=None): - """Read :math:`\Delta_\text{TE}` from BIDS metadata dict""" + r"""Read :math:`\Delta_\text{TE}` from BIDS metadata dict""" if isinstance(in_values, float): te2 = in_values te1 = 0. diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index fe7f3aafa..6ce2d9988 100755 --- a/fmriprep/workflows/bold/base.py +++ b/fmriprep/workflows/bold/base.py @@ -336,7 +336,11 @@ def init_func_preproc_wf( if 'fieldmaps' not in ignore: fmaps = layout.get_fieldmap(ref_file, return_list=True) for fmap in fmaps: - fmap['metadata'] = layout.get_metadata(fmap[fmap['suffix']]) + if fmap['suffix'] == 'phase': + fmap_key = 'phase1' + else: + fmap_key = fmap['suffix'] + fmap['metadata'] = layout.get_metadata(fmap[fmap_key]) # Run SyN if forced or in the absence of fieldmap correction if force_syn or (use_syn and not fmaps): diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py index 2c83561f3..b60449cf2 100644 --- a/fmriprep/workflows/fieldmap/base.py +++ b/fmriprep/workflows/fieldmap/base.py @@ -53,7 +53,8 @@ 'epi': 0, 'fieldmap': 1, 'phasediff': 2, - 'syn': 3 + 'phase': 3, + 'syn': 4 } DEFAULT_MEMORY_MIN_GB = 0.01 @@ -199,7 +200,7 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1, ]) # FIELDMAP path - if fmap['suffix'] in ['fieldmap', 'phasediff']: + if fmap['suffix'] in ['fieldmap', 'phasediff', 'phase']: outputnode.inputs.method = 'FMB (%s-based)' % fmap['suffix'] # Import specific workflows here, so we don't break everything with one # unused workflow. @@ -212,11 +213,30 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1, fmap_estimator_wf.inputs.inputnode.fieldmap = fmap['fieldmap'] fmap_estimator_wf.inputs.inputnode.magnitude = fmap['magnitude'] - if fmap['suffix'] == 'phasediff': + if fmap['suffix'] in ('phasediff', 'phase'): from .phdiff import init_phdiff_wf - fmap_estimator_wf = init_phdiff_wf(omp_nthreads=omp_nthreads) + fmap_estimator_wf = init_phdiff_wf(omp_nthreads=omp_nthreads, + phasetype=fmap['suffix']) # set inputs - fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff'] + if fmap['suffix'] == 'phasediff': + fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff'] + elif fmap['suffix'] == 'phase': + # Check that fieldmap is not bipolar + fmap_polarity = fmap['metadata'].get('DiffusionScheme', None) + if fmap_polarity == 'Bipolar': + LOGGER.warning("Bipolar fieldmaps are not supported. Ignoring") + workflow.__postdesc__ = "" + outputnode.inputs.method = 'None' + workflow.connect([ + (inputnode, outputnode, [('bold_ref', 'bold_ref'), + ('bold_mask', 'bold_mask'), + ('bold_ref_brain', 'bold_ref_brain')]), + ]) + return workflow + if fmap_polarity is None: + LOGGER.warning("Assuming phase images are Monopolar") + + fmap_estimator_wf.inputs.inputnode.phasediff = [fmap['phase1'], fmap['phase2']] fmap_estimator_wf.inputs.inputnode.magnitude = [ fmap_ for key, fmap_ in sorted(fmap.items()) if key.startswith("magnitude") diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py index a2a850400..c233f4059 100644 --- a/fmriprep/workflows/fieldmap/phdiff.py +++ b/fmriprep/workflows/fieldmap/phdiff.py @@ -15,6 +15,7 @@ Fieldmap preprocessing workflow for fieldmap data structure 8.9.1 in BIDS 1.0.0: one phase diff and at least one magnitude image +8.9.2 in BIDS 1.0.0: two phases and at least one magnitude image """ @@ -27,10 +28,10 @@ from niworkflows.interfaces.images import IntraModalMerge from niworkflows.interfaces.masks import BETRPT -from ...interfaces import Phasediff2Fieldmap, DerivativesDataSink +from ...interfaces import Phasediff2Fieldmap, Phases2Fieldmap, DerivativesDataSink -def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): +def init_phdiff_wf(omp_nthreads, phasetype='phasediff', name='phdiff_wf'): """ Estimates the fieldmap using a phase-difference image and one or more magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)` @@ -69,12 +70,6 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): outputnode = pe.Node(niu.IdentityInterface( fields=['fmap', 'fmap_ref', 'fmap_mask']), name='outputnode') - def _pick1st(inlist): - return inlist[0] - - # Read phasediff echo times - meta = pe.Node(ReadSidecarJSON(bids_validate=False), name='meta', mem_gb=0.01) - # Merge input magnitude images magmrg = pe.Node(IntraModalMerge(), name='magmrg') @@ -90,9 +85,6 @@ def _pick1st(inlist): # dilate = pe.Node(fsl.maths.MathsCommand( # nan2zeros=True, args='-kernel sphere 5 -dilM'), name='MskDilate') - # phase diff -> radians - pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads') - # FSL PRELUDE will perform phase-unwrapping prelude = pe.Node(fsl.PRELUDE(), name='prelude') @@ -110,6 +102,40 @@ def _pick1st(inlist): # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE') # rsec2hz (divide by 2pi) + if phasetype == "phasediff": + # Read phasediff echo times + meta = pe.Node(ReadSidecarJSON(bids_validate=False), name='meta', mem_gb=0.01) + + # phase diff -> radians + pha2rads = pe.Node(niu.Function(function=siemens2rads), + name='pha2rads') + # Read phasediff echo times + meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01, + run_without_submitting=True) + workflow.connect([ + (meta, compfmap, [('out_dict', 'metadata')]), + (inputnode, pha2rads, [('phasediff', 'in_file')]), + (pha2rads, prelude, [('out', 'phase_file')]), + (inputnode, ds_report_fmap_mask, [('phasediff', 'source_file')]), + ]) + + elif phasetype == "phase": + workflow.__desc__ += """\ +The phase difference used for unwarping was calculated using two separate phase measurements + [@pncprocessing]. + """ + # Special case for phase1, phase2 images + meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01, + run_without_submitting=True, iterfield=['in_file']) + phases2fmap = pe.Node(Phases2Fieldmap(), name='phases2fmap') + workflow.connect([ + (meta, phases2fmap, [('out_dict', 'metadatas')]), + (inputnode, phases2fmap, [('phasediff', 'phase_files')]), + (phases2fmap, prelude, [('out_file', 'phase_file')]), + (phases2fmap, compfmap, [('phasediff_metadata', 'metadata')]), + (phases2fmap, ds_report_fmap_mask, [('out_file', 'source_file')]) + ]) + workflow.connect([ (inputnode, meta, [('phasediff', 'in_file')]), (inputnode, magmrg, [('magnitude', 'in_files')]), @@ -117,9 +143,6 @@ def _pick1st(inlist): (n4, prelude, [('output_image', 'magnitude_file')]), (n4, bet, [('output_image', 'in_file')]), (bet, prelude, [('mask_file', 'mask_file')]), - (inputnode, pha2rads, [('phasediff', 'in_file')]), - (pha2rads, prelude, [('out', 'phase_file')]), - (meta, compfmap, [('out_dict', 'metadata')]), (prelude, denoise, [('unwrapped_phase_file', 'in_file')]), (denoise, demean, [('out_file', 'in_file')]), (demean, cleanup_wf, [('out', 'inputnode.in_file')]), @@ -128,7 +151,6 @@ def _pick1st(inlist): (compfmap, outputnode, [('out_file', 'fmap')]), (bet, outputnode, [('mask_file', 'fmap_mask'), ('out_file', 'fmap_ref')]), - (inputnode, ds_report_fmap_mask, [('phasediff', 'source_file')]), (bet, ds_report_fmap_mask, [('out_report', 'in_file')]), ])