diff --git a/.circleci/config.yml b/.circleci/config.yml index d3f687906..819b2189d 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -175,7 +175,7 @@ jobs: /out participant --participant-label 01 \ -w /scrth/wkdir_sub01 --output-spaces asl \ --fs-license-file /scrth/licence.txt \ - --scorescrub --basil \ + --scorescrub --basil --use-syn-sdc \ --anat-derivatives /data/smriprep - run: @@ -264,7 +264,7 @@ jobs: /out participant --participant-label A00086748 \ -w /scrth/wkdir_subA00086748 --output-spaces asl \ --fs-license-file /scrth/licence.txt \ - --scorescrub --basil \ + --scorescrub --basil --use-syn-sdc \ --anat-derivatives /data/smriprep - run: diff --git a/aslprep/interfaces/cbf_computation.py b/aslprep/interfaces/cbf_computation.py index 1cd642020..b003cc924 100644 --- a/aslprep/interfaces/cbf_computation.py +++ b/aslprep/interfaces/cbf_computation.py @@ -319,7 +319,7 @@ def cbfcomputation(metadata, mask, m0file, cbffile, m0scale=1): for i in range(cbf_data.shape[1]): cbf_data_ts[:, i] =np.multiply(cbf1[:, i],permfactor[i]) cbf = np.zeros([cbf_data_ts.shape[0], int(cbf_data.shape[1]/len(perfusion_factor))]) - cbf_xx=np.split(cbf_data_ts,int(cbf_data_ts.shape[1]/len(perfusion_factor)),axis=1) + cbf_xx = np.split(cbf_data_ts,int(cbf_data_ts.shape[1]/len(perfusion_factor)),axis=1) # calculate weighted cbf with multiplds # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3791289/ @@ -344,9 +344,9 @@ def cbfcomputation(metadata, mask, m0file, cbffile, m0scale=1): tcbf=np.zeros(maskx.shape) tcbf[maskx==1]=cbf else: - tcbf=np.zeros([maskx.shape[0],maskx.shape[1],maskx.shape[2],cbf.shape[1]]) + tcbf = np.zeros([maskx.shape[0],maskx.shape[1],maskx.shape[2],cbf.shape[1]]) for i in range(cbf.shape[1]): - tcbfx=np.zeros(maskx.shape) + tcbfx = np.zeros(maskx.shape) tcbfx[maskx==1]=cbf[:,i] tcbf[:,:,:,i]=tcbfx if len(tcbf.shape) < 4: @@ -458,7 +458,7 @@ def _weightfun(x, wfun='huber'): weight = 1/(1+np.power(x, 2)) elif wfun == 'logistic': tuner = 1.205 - weight == np.tanh(x)/x + weight = np.tanh(x)/x elif wfun == 'ols': tuner = 1 weight = np.repeat(1, len(x)) diff --git a/aslprep/niworkflows/viz/plots.py b/aslprep/niworkflows/viz/plots.py index b49a4d2e0..1d22fd5b2 100644 --- a/aslprep/niworkflows/viz/plots.py +++ b/aslprep/niworkflows/viz/plots.py @@ -241,9 +241,9 @@ def plotstatsimg(cbf, ref_vol, plot_params=None, order=('z', 'x', 'y'), plot_params['display_mode'] = mode plot_params['cut_coords'] = cuts[mode] plot_params['draw_cross'] = False - plot_params['symmetric_cbar'] = False - plot_params['threshold'] = 0.05 + plot_params['symmetric_cbar'] = True plot_params['vmax'] = 90 + plot_params['threshold'] = 0.02 plot_params['colorbar'] = False plot_params['cmap'] = 'gray' if i == 0: @@ -254,7 +254,8 @@ def plotstatsimg(cbf, ref_vol, plot_params=None, order=('z', 'x', 'y'), # Generate nilearn figure from nilearn.plotting import plot_stat_map - display = plot_stat_map(stat_map_img=cbf, bg_img=ref_vol, **plot_params) + display = plot_stat_map(stat_map_img=cbf, bg_img=ref_vol, + resampling_interpolation='nearest',**plot_params) svg = extract_svg(display, compress=compress) display.close() from lxml import etree diff --git a/aslprep/sdcflows/__init__.py b/aslprep/sdcflows/__init__.py new file mode 100644 index 000000000..74fe6e44e --- /dev/null +++ b/aslprep/sdcflows/__init__.py @@ -0,0 +1,13 @@ +"""SDCflows - :abbr:`SDC (susceptibility distortion correction)` by DUMMIES, for dummies.""" +__packagename__ = "sdcflows" +__copyright__ = "2020, The NiPreps developers" +try: + from ._version import __version__ +except ModuleNotFoundError: + from pkg_resources import get_distribution, DistributionNotFound + try: + __version__ = get_distribution(__packagename__).version + except DistributionNotFound: + __version__ = "unknown" + del get_distribution + del DistributionNotFound diff --git a/aslprep/sdcflows/cli/__init__.py b/aslprep/sdcflows/cli/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aslprep/sdcflows/cli/run.py b/aslprep/sdcflows/cli/run.py new file mode 100644 index 000000000..f6426e9b9 --- /dev/null +++ b/aslprep/sdcflows/cli/run.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python3 +import logging +from pathlib import Path + +logging.addLevelName(25, 'IMPORTANT') # Add a new level between INFO and WARNING +logging.addLevelName(15, 'VERBOSE') # Add a new level between INFO and DEBUG +logger = logging.getLogger('sdcflows') + + +def get_parser(): + """Define the command line interface""" + from argparse import ArgumentParser + from argparse import RawTextHelpFormatter + from .. import __version__ as _vstr + + parser = ArgumentParser(description='SDC Workflows', + formatter_class=RawTextHelpFormatter) + + parser.add_argument( + 'bids_dir', action='store', type=Path, + help='the root folder of a BIDS dataset') + parser.add_argument('output_dir', action='store', type=Path, + help='the output path for the outcomes of preprocessing and visual ' + 'reports') + parser.add_argument('analysis_level', choices=['participant', 'group'], nargs='+', + help='processing stage to be run, "participant" means individual analysis ' + 'and "group" is second level analysis.') + # optional arguments + parser.add_argument('--version', action='version', version='v{}'.format(_vstr)) + + # Options that affect how pyBIDS is configured + g_bids = parser.add_argument_group('Options for filtering BIDS queries') + g_bids.add_argument('--participant-label', action='store', type=str, + nargs='*', dest='subject', help='process only particular subjects') + g_bids.add_argument('--task', action='store', type=str, nargs='*', + help='select a specific task to be processed') + g_bids.add_argument('--dir', action='store', type=str, nargs='*', + help='select a specific direction entity to be processed') + g_bids.add_argument('--acq', action='store', type=str, nargs='*', dest='acquisition', + help='select a specific acquisition entity to be processed') + g_bids.add_argument('--run', action='store', type=int, nargs='*', + help='select a specific run identifier to be processed') + g_bids.add_argument('--suffix', action='store', type=str, nargs='*', default='bold', + help='select a specific run identifier to be processed') + + g_perfm = parser.add_argument_group('Options to handle performance') + g_perfm.add_argument("-v", "--verbose", dest="verbose_count", action="count", default=0, + help="increases log verbosity for each occurence, debug level is -vvv") + g_perfm.add_argument('--ncpus', '--nprocs', action='store', type=int, + help='maximum number of threads across all processes') + g_perfm.add_argument('--nthreads', '--omp-nthreads', action='store', type=int, + help='maximum number of threads per-process') + + g_other = parser.add_argument_group('Other options') + g_other.add_argument('-w', '--work-dir', action='store', type=Path, + help='path where intermediate results should be stored') + + return parser + + +def main(): + """Entry point""" + from os import cpu_count + from multiprocessing import set_start_method + # from bids.layout import BIDSLayout + from nipype import logging as nlogging + set_start_method('forkserver') + + opts = get_parser().parse_args() + + # Retrieve logging level + log_level = int(max(25 - 5 * opts.verbose_count, logging.DEBUG)) + # Set logging + logger.setLevel(log_level) + nlogging.getLogger('nipype.workflow').setLevel(log_level) + nlogging.getLogger('nipype.interface').setLevel(log_level) + nlogging.getLogger('nipype.utils').setLevel(log_level) + + # Resource management options + plugin_settings = { + 'plugin': 'MultiProc', + 'plugin_args': { + 'n_procs': opts.ncpus, + '' + 'raise_insufficient': False, + 'maxtasksperchild': 1, + } + } + # Permit overriding plugin config with specific CLI options + if not opts.ncpus or opts.ncpus < 1: + plugin_settings['plugin_args']['n_procs'] = cpu_count() + + nthreads = opts.nthreads + if not nthreads or nthreads < 1: + nthreads = cpu_count() + + # output_dir = opts.output_dir.resolve() + # bids_dir = opts.bids_dir or output_dir.parent + + # Get absolute path to BIDS directory + # bids_dir = opts.bids_dir.resolve() + # layout = BIDSLayout(str(bids_dir), validate=False, derivatives=str(output_dir)) + # query = {'suffix': opts.suffix, 'extension': ['.nii', '.nii.gz']} + + # for entity in ('subject', 'task', 'dir', 'acquisition', 'run'): + # arg = getattr(opts, entity, None) + # if arg is not None: + # query[entity] = arg + + +if __name__ == '__main__': + raise RuntimeError("sdcflows/cli/run.py should not be run directly;\n" + "Please `pip install` sdcflows and use the `sdcflows` command") diff --git a/aslprep/sdcflows/conftest.py b/aslprep/sdcflows/conftest.py new file mode 100644 index 000000000..1ffdb109d --- /dev/null +++ b/aslprep/sdcflows/conftest.py @@ -0,0 +1,44 @@ +"""py.test configuration""" +import os +from pathlib import Path +import numpy +import pytest +from bids.layout import BIDSLayout + +test_data_env = os.getenv('TEST_DATA_HOME', str(Path.home() / 'sdcflows-tests')) +test_output_dir = os.getenv('TEST_OUTPUT_DIR') +test_workdir = os.getenv('TEST_WORK_DIR') + +layouts = {p.name: BIDSLayout(str(p), validate=False, derivatives=True) + for p in Path(test_data_env).glob('*') if p.is_dir()} + + +def pytest_report_header(config): + msg = "Datasets found: %s" % ', '.join([v.root for v in layouts.values()]) + if test_output_dir is not None: + msg += '\nOutput folder: %s' % Path(test_output_dir).resolve() + return msg + + +@pytest.fixture(autouse=True) +def add_np(doctest_namespace): + doctest_namespace['np'] = numpy + doctest_namespace['os'] = os + doctest_namespace['Path'] = Path + for key, val in list(layouts.items()): + doctest_namespace[key] = Path(val.root) + + +@pytest.fixture +def workdir(): + return None if test_workdir is None else Path(test_workdir) + + +@pytest.fixture +def output_path(): + return None if test_output_dir is None else Path(test_output_dir) + + +@pytest.fixture +def bids_layouts(): + return layouts diff --git a/aslprep/sdcflows/data/__init__.py b/aslprep/sdcflows/data/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aslprep/sdcflows/data/affine.json b/aslprep/sdcflows/data/affine.json new file mode 100644 index 000000000..1c48065fd --- /dev/null +++ b/aslprep/sdcflows/data/affine.json @@ -0,0 +1,24 @@ +{ + "dimension": 3, + "float": false, + "initial_moving_transform_com": 1, + "winsorize_lower_quantile": 0.005, + "winsorize_upper_quantile": 0.995, + "collapse_output_transforms": true, + "use_histogram_matching": [ false, false ], + "use_estimate_learning_rate_once": [ false, false ], + "transforms": [ "Rigid", "Affine" ], + "number_of_iterations": [ [ 1000, 500, 250, 100 ], [ 1000, 500, 250, 100 ] ], + "transform_parameters": [ [ 0.1 ], [ 0.1 ] ], + "convergence_threshold": [ 1e-06, 1e-06 ], + "convergence_window_size": [ 10, 10 ], + "metric": [ "MI", "MI" ], + "sampling_percentage": [ 0.25, 0.25 ], + "sampling_strategy": [ "Regular", "Regular" ], + "smoothing_sigmas": [ [ 3, 2, 1, 0 ], [ 3, 2, 1, 0 ] ], + "sigma_units": [ "vox", "vox" ], + "metric_weight": [ 1.0, 1.0 ], + "shrink_factors": [[ 8, 4, 2, 1] , [ 8, 4, 2, 1 ]], + "radius_or_number_of_bins": [ 32, 32 ], + "interpolation": "Linear" +} diff --git a/aslprep/sdcflows/data/fmap-any_registration.json b/aslprep/sdcflows/data/fmap-any_registration.json new file mode 100644 index 000000000..b97b002fc --- /dev/null +++ b/aslprep/sdcflows/data/fmap-any_registration.json @@ -0,0 +1,24 @@ +{ + "dimension": 3, + "float": true, + "winsorize_lower_quantile": 0.005, + "winsorize_upper_quantile": 0.998, + "collapse_output_transforms": true, + "write_composite_transform": true, + "use_histogram_matching": [ true, true ], + "use_estimate_learning_rate_once": [ true, true ], + "transforms": [ "Translation", "Affine" ], + "number_of_iterations": [ [ 500 ], [ 200 ] ], + "transform_parameters": [ [ 0.05 ], [ 0.01 ] ], + "convergence_threshold": [ 1e-07, 1e-08 ], + "convergence_window_size": [ 200, 100 ], + "metric": [ "Mattes", "Mattes" ], + "sampling_percentage": [ 0.5, 0.5 ], + "sampling_strategy": [ "Random", "Random" ], + "smoothing_sigmas": [ [ 8.0 ], [ 2.0 ] ], + "sigma_units": [ "mm", "mm" ], + "metric_weight": [ 1.0, 1.0 ], + "shrink_factors": [ [ 2 ], [ 1 ] ], + "radius_or_number_of_bins": [ 64, 64 ], + "interpolation": "LanczosWindowedSinc" +} \ No newline at end of file diff --git a/aslprep/sdcflows/data/fmap-any_registration_testing.json b/aslprep/sdcflows/data/fmap-any_registration_testing.json new file mode 100644 index 000000000..207e997ed --- /dev/null +++ b/aslprep/sdcflows/data/fmap-any_registration_testing.json @@ -0,0 +1,24 @@ +{ + "dimension": 3, + "float": true, + "winsorize_lower_quantile": 0.005, + "winsorize_upper_quantile": 0.998, + "collapse_output_transforms": true, + "write_composite_transform": true, + "use_histogram_matching": [ true, true ], + "use_estimate_learning_rate_once": [ true, true ], + "transforms": [ "Rigid", "Rigid" ], + "number_of_iterations": [ [ 500 ], [ 100 ] ], + "transform_parameters": [ [ 1.0 ], [ 0.5 ] ], + "convergence_threshold": [ 1e-07, 1e-08 ], + "convergence_window_size": [ 20, 10 ], + "metric": [ "Mattes", "Mattes" ], + "sampling_percentage": [ 1.0, 1.0 ], + "sampling_strategy": [ "Regular", "Regular" ], + "smoothing_sigmas": [ [ 8.0 ], [ 2.0 ] ], + "sigma_units": [ "mm", "mm" ], + "metric_weight": [ 1.0, 1.0 ], + "shrink_factors": [ [ 2 ], [ 1 ] ], + "radius_or_number_of_bins": [ 64, 64 ], + "interpolation": "LanczosWindowedSinc" +} \ No newline at end of file diff --git a/aslprep/sdcflows/data/fmap_atlas.nii.gz b/aslprep/sdcflows/data/fmap_atlas.nii.gz new file mode 100644 index 000000000..1e3af1ab5 Binary files /dev/null and b/aslprep/sdcflows/data/fmap_atlas.nii.gz differ diff --git a/aslprep/sdcflows/data/fmap_atlas_2_MNI152NLin2009cAsym_affine.mat b/aslprep/sdcflows/data/fmap_atlas_2_MNI152NLin2009cAsym_affine.mat new file mode 100644 index 000000000..d67baf5a1 Binary files /dev/null and b/aslprep/sdcflows/data/fmap_atlas_2_MNI152NLin2009cAsym_affine.mat differ diff --git a/aslprep/sdcflows/data/susceptibility_syn.json b/aslprep/sdcflows/data/susceptibility_syn.json new file mode 100644 index 000000000..d736a70db --- /dev/null +++ b/aslprep/sdcflows/data/susceptibility_syn.json @@ -0,0 +1,24 @@ +{ + "dimension": 3, + "float": true, + "convergence_threshold": [ 1e-08, 1e-08 ], + "convergence_window_size": [ 20, 10 ], + "metric": [ "Mattes", "CC" ], + "metric_weight": [ 1, 1 ], + "radius_or_number_of_bins": [ 56, 5 ], + "transforms": [ "SyN", "SyN" ], + "transform_parameters": [ [ 0.8, 2, 2 ], [ 0.8, 2, 2 ] ], + "number_of_iterations": [ [ 100, 50 ], [ 20, 10 ] ], + "smoothing_sigmas": [ [ 1, 0 ], [ 1, 0 ] ], + "sigma_units": [ "vox", "vox" ], + "shrink_factors": [ [ 2, 1 ], [ 1, 1 ] ], + "winsorize_upper_quantile": 1.0, + "winsorize_lower_quantile": 0.001, + "use_estimate_learning_rate_once": [ true, true ], + "use_histogram_matching": [ true, true ], + "collapse_output_transforms": true, + "write_composite_transform": false, + "output_transform_prefix": "ants_susceptibility", + "output_warped_image": true, + "interpolation": "Linear" +} diff --git a/aslprep/sdcflows/data/translation_rigid.json b/aslprep/sdcflows/data/translation_rigid.json new file mode 100644 index 000000000..6e14b0225 --- /dev/null +++ b/aslprep/sdcflows/data/translation_rigid.json @@ -0,0 +1,24 @@ +{ + "dimension": 3, + "float": true, + "winsorize_lower_quantile": 0.005, + "winsorize_upper_quantile": 0.998, + "collapse_output_transforms": true, + "write_composite_transform": true, + "use_histogram_matching": [ true, true ], + "use_estimate_learning_rate_once": [ true, true ], + "transforms": [ "Translation", "Rigid" ], + "number_of_iterations": [ [ 500 ], [ 200 ] ], + "transform_parameters": [ [ 0.05 ], [ 0.01 ] ], + "convergence_threshold": [ 1e-07, 1e-08 ], + "convergence_window_size": [ 200, 100 ], + "metric": [ "Mattes", "Mattes" ], + "sampling_percentage": [ 0.5, 0.5 ], + "sampling_strategy": [ "Random", "Random" ], + "smoothing_sigmas": [ [ 8.0 ], [ 2.0 ] ], + "sigma_units": [ "mm", "mm" ], + "metric_weight": [ 1.0, 1.0 ], + "shrink_factors": [ [ 2 ], [ 1 ] ], + "radius_or_number_of_bins": [ 64, 64 ], + "interpolation": "LanczosWindowedSinc" +} \ No newline at end of file diff --git a/aslprep/sdcflows/interfaces/__init__.py b/aslprep/sdcflows/interfaces/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aslprep/sdcflows/interfaces/fmap.py b/aslprep/sdcflows/interfaces/fmap.py new file mode 100644 index 000000000..78b62b12c --- /dev/null +++ b/aslprep/sdcflows/interfaces/fmap.py @@ -0,0 +1,715 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Interfaces to deal with the various types of fieldmap sources. + + .. testsetup:: + + >>> tmpdir = getfixture('tmpdir') + >>> tmp = tmpdir.chdir() # changing to a temporary directory + >>> nb.Nifti1Image(np.zeros((90, 90, 60)), None, None).to_filename( + ... tmpdir.join('epi.nii.gz').strpath) + + +""" + +import numpy as np +import nibabel as nb +from nipype import logging +from nipype.utils.filemanip import fname_presuffix +from nipype.interfaces.base import ( + BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits, + SimpleInterface) + +LOGGER = logging.getLogger('nipype.interface') + + +class _SubtractPhasesInputSpec(BaseInterfaceInputSpec): + in_phases = traits.List(File(exists=True), min=1, max=2, + desc='input phase maps') + in_meta = traits.List(traits.Dict(), min=1, max=2, + desc='metadata corresponding to the inputs') + + +class _SubtractPhasesOutputSpec(TraitedSpec): + phase_diff = File(exists=True, desc='phase difference map') + metadata = traits.Dict(desc='output metadata') + + +class SubtractPhases(SimpleInterface): + """Calculate a phase difference map.""" + + input_spec = _SubtractPhasesInputSpec + output_spec = _SubtractPhasesOutputSpec + + def _run_interface(self, runtime): + if len(self.inputs.in_phases) != len(self.inputs.in_meta): + raise ValueError( + 'Length of input phase-difference maps and metadata files ' + 'should match.') + + if len(self.inputs.in_phases) == 1: + self._results['phase_diff'] = self.inputs.in_phases[0] + self._results['metadata'] = self.inputs.in_meta[0] + return runtime + + self._results['phase_diff'], self._results['metadata'] = \ + _subtract_phases(self.inputs.in_phases, + self.inputs.in_meta, + newpath=runtime.cwd) + + return runtime + + +class _FieldEnhanceInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc='input fieldmap') + in_mask = File(exists=True, desc='brain mask') + in_magnitude = File(exists=True, desc='input magnitude') + unwrap = traits.Bool(False, usedefault=True, desc='run phase unwrap') + despike = traits.Bool(True, usedefault=True, desc='run despike filter') + bspline_smooth = traits.Bool(True, usedefault=True, desc='run 3D bspline smoother') + mask_erode = traits.Int(1, usedefault=True, desc='mask erosion iterations') + despike_threshold = traits.Float(0.2, usedefault=True, desc='mask erosion iterations') + num_threads = traits.Int(1, usedefault=True, nohash=True, desc='number of jobs') + + +class _FieldEnhanceOutputSpec(TraitedSpec): + out_file = File(desc='the output fieldmap') + out_unwrapped = File(desc='unwrapped fieldmap') + + +class FieldEnhance(SimpleInterface): + """Massage the input fieldmap (masking, despiking, etc.).""" + + input_spec = _FieldEnhanceInputSpec + output_spec = _FieldEnhanceOutputSpec + + def _run_interface(self, runtime): + from scipy import ndimage as sim + + fmap_nii = nb.load(self.inputs.in_file) + data = np.squeeze(fmap_nii.get_fdata(dtype='float32')) + + # Despike / denoise (no-mask) + if self.inputs.despike: + data = _despike2d(data, self.inputs.despike_threshold) + + mask = None + if isdefined(self.inputs.in_mask): + masknii = nb.load(self.inputs.in_mask) + mask = np.asanyarray(masknii.dataobj).astype('uint8') + + # Dilate mask + if self.inputs.mask_erode > 0: + struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 1) + mask = sim.binary_erosion( + mask, struc, + iterations=self.inputs.mask_erode + ).astype(np.uint8) # pylint: disable=no-member + + self._results['out_file'] = fname_presuffix( + self.inputs.in_file, suffix='_enh', newpath=runtime.cwd) + datanii = nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header) + + if self.inputs.unwrap: + data = _unwrap(data, self.inputs.in_magnitude, mask) + self._results['out_unwrapped'] = fname_presuffix( + self.inputs.in_file, suffix='_unwrap', newpath=runtime.cwd) + nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header).to_filename( + self._results['out_unwrapped']) + + if not self.inputs.bspline_smooth: + datanii.to_filename(self._results['out_file']) + return runtime + else: + from ..utils import bspline as fbsp + from statsmodels.robust.scale import mad + + # Fit BSplines (coarse) + bspobj = fbsp.BSplineFieldmap(datanii, weights=mask, + njobs=self.inputs.num_threads) + bspobj.fit() + smoothed1 = bspobj.get_smoothed() + + # Manipulate the difference map + diffmap = data - smoothed1.get_fdata(dtype='float32') + sderror = mad(diffmap[mask > 0]) + LOGGER.info('SD of error after B-Spline fitting is %f', sderror) + errormask = np.zeros_like(diffmap) + errormask[np.abs(diffmap) > (10 * sderror)] = 1 + errormask *= mask + + nslices = 0 + try: + errorslice = np.squeeze(np.argwhere(errormask.sum(0).sum(0) > 0)) + nslices = errorslice[-1] - errorslice[0] + except IndexError: # mask is empty, do not refine + pass + + if nslices > 1: + diffmapmsk = mask[..., errorslice[0]:errorslice[-1]] + diffmapnii = nb.Nifti1Image( + diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk, + datanii.affine, datanii.header) + + bspobj2 = fbsp.BSplineFieldmap(diffmapnii, knots_zooms=[24., 24., 4.], + njobs=self.inputs.num_threads) + bspobj2.fit() + smoothed2 = bspobj2.get_smoothed().get_fdata(dtype='float32') + + final = smoothed1.get_fdata(dtype='float32').copy() + final[..., errorslice[0]:errorslice[-1]] += smoothed2 + else: + final = smoothed1.get_fdata(dtype='float32') + + nb.Nifti1Image(final, datanii.affine, datanii.header).to_filename( + self._results['out_file']) + + return runtime + + +class _FieldToRadSInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc='input fieldmap') + fmap_range = traits.Float(desc='range of input field map') + + +class _FieldToRadSOutputSpec(TraitedSpec): + out_file = File(desc='the output fieldmap') + fmap_range = traits.Float(desc='range of input field map') + + +class FieldToRadS(SimpleInterface): + """Convert from arbitrary units to rad/s.""" + + input_spec = _FieldToRadSInputSpec + output_spec = _FieldToRadSOutputSpec + + def _run_interface(self, runtime): + fmap_range = None + if isdefined(self.inputs.fmap_range): + fmap_range = self.inputs.fmap_range + self._results['out_file'], self._results['fmap_range'] = _torads( + self.inputs.in_file, fmap_range, newpath=runtime.cwd) + return runtime + + +class _FieldToHzInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc='input fieldmap') + range_hz = traits.Float(mandatory=True, desc='range of input field map') + + +class _FieldToHzOutputSpec(TraitedSpec): + out_file = File(desc='the output fieldmap') + + +class FieldToHz(SimpleInterface): + """Convert from arbitrary units to Hz.""" + + input_spec = _FieldToHzInputSpec + output_spec = _FieldToHzOutputSpec + + def _run_interface(self, runtime): + self._results['out_file'] = _tohz( + self.inputs.in_file, self.inputs.range_hz, newpath=runtime.cwd) + return runtime + + +class _Phasediff2FieldmapInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc='input fieldmap') + metadata = traits.Dict(mandatory=True, desc='BIDS metadata dictionary') + + +class _Phasediff2FieldmapOutputSpec(TraitedSpec): + out_file = File(desc='the output fieldmap') + + +class Phasediff2Fieldmap(SimpleInterface): + """ + Convert a phase difference map into a fieldmap in Hz. + + This interface is equivalent to running the following steps: + #. Convert from rad to rad/s + (``niflow.nipype1.workflows.dmri.fsl.utils.rads2radsec``) + #. FUGUE execution: fsl.FUGUE(save_fmap=True) + #. Conversion from rad/s to Hz (divide by 2pi, ``rsec2hz``). + + """ + + input_spec = _Phasediff2FieldmapInputSpec + output_spec = _Phasediff2FieldmapOutputSpec + + def _run_interface(self, runtime): + self._results['out_file'] = phdiff2fmap( + self.inputs.in_file, + _delta_te(self.inputs.metadata), + newpath=runtime.cwd) + return runtime + + +class _PhaseMap2radsInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc='input (wrapped) phase map') + + +class _PhaseMap2radsOutputSpec(TraitedSpec): + out_file = File(desc='the phase map in the range 0 - 6.28') + + +class PhaseMap2rads(SimpleInterface): + """Convert a phase map in a.u. to radians.""" + + input_spec = _PhaseMap2radsInputSpec + output_spec = _PhaseMap2radsOutputSpec + + def _run_interface(self, runtime): + self._results['out_file'] = au2rads( + self.inputs.in_file, + newpath=runtime.cwd) + return runtime + + +class _FUGUEvsm2ANTSwarpInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, + desc='input displacements field map') + pe_dir = traits.Enum('i', 'i-', 'j', 'j-', 'k', 'k-', + desc='phase-encoding axis') + + +class _FUGUEvsm2ANTSwarpOutputSpec(TraitedSpec): + out_file = File(desc='the output warp field') + fieldmap = File(desc='field map in mm') + + +class FUGUEvsm2ANTSwarp(SimpleInterface): + """Convert a voxel-shift-map to ants warp.""" + + _dtype = ' 1e-6 and + (abs(thisval - patch_med) / patch_range) > thres): + data[i, j, k] = patch_med + return data + + +def _unwrap(fmap_data, mag_file, mask=None): + from math import pi + from nipype.interfaces.fsl import PRELUDE + magnii = nb.load(mag_file) + + if mask is None: + mask = np.ones_like(fmap_data, dtype=np.uint8) + + fmapmax = max(abs(fmap_data[mask > 0].min()), fmap_data[mask > 0].max()) + fmap_data *= pi / fmapmax + + nb.Nifti1Image(fmap_data, magnii.affine).to_filename('fmap_rad.nii.gz') + nb.Nifti1Image(mask, magnii.affine).to_filename('fmap_mask.nii.gz') + nb.Nifti1Image(magnii.get_fdata(dtype='float32'), + magnii.affine).to_filename('fmap_mag.nii.gz') + + # Run prelude + res = PRELUDE(phase_file='fmap_rad.nii.gz', + magnitude_file='fmap_mag.nii.gz', + mask_file='fmap_mask.nii.gz').run() + + unwrapped = nb.load( + res.outputs.unwrapped_phase_file).get_fdata(dtype='float32') * (fmapmax / pi) + return unwrapped + + +def get_ees(in_meta, in_file=None): + r""" + Extract the *effective echo spacing* :math:`t_\text{ees}` from BIDS. + + Calculate the *effective echo spacing* :math:`t_\text{ees}` + for an input :abbr:`EPI (echo-planar imaging)` scan. + + + There are several procedures to calculate the effective + echo spacing. The basic one is that an ``EffectiveEchoSpacing`` + field is set in the JSON sidecar. The following examples + use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the + j-axis encoding direction. + + >>> meta = {'EffectiveEchoSpacing': 0.00059, + ... 'PhaseEncodingDirection': 'j-'} + >>> get_ees(meta) + 0.00059 + + If the *total readout time* :math:`T_\text{ro}` (``TotalReadoutTime`` + BIDS field) is provided, then the effective echo spacing can be + calculated reading the number of voxels :math:`N_\text{PE}` along the + readout direction and the parallel acceleration + factor of the EPI + + .. math :: + + = T_\text{ro} \, (N_\text{PE} / f_\text{acc} - 1)^{-1} + + where :math:`N_y` is the number of pixels along the phase-encoding direction + :math:`y`, and :math:`f_\text{acc}` is the parallel imaging acceleration factor + (:abbr:`GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition)`, + :abbr:`ARC (Autocalibrating Reconstruction for Cartesian imaging)`, etc.). + + >>> meta = {'TotalReadoutTime': 0.02596, + ... 'PhaseEncodingDirection': 'j-', + ... 'ParallelReductionFactorInPlane': 2} + >>> get_ees(meta, in_file='epi.nii.gz') + 0.00059 + + Some vendors, like Philips, store different parameter names (see + http://dbic.dartmouth.edu/pipermail/mrusers/attachments/20141112/eb1d20e6/attachment.pdf + ): + + >>> meta = {'WaterFatShift': 8.129, + ... 'MagneticFieldStrength': 3, + ... 'PhaseEncodingDirection': 'j-', + ... 'ParallelReductionFactorInPlane': 2} + >>> get_ees(meta, in_file='epi.nii.gz') + 0.00041602630141921826 + + """ + + import nibabel as nb + from sdcflows.interfaces.fmap import _get_pe_index + + # Use case 1: EES is defined + ees = in_meta.get('EffectiveEchoSpacing', None) + if ees is not None: + return ees + + # All other cases require the parallel acc and npe (N vox in PE dir) + acc = float(in_meta.get('ParallelReductionFactorInPlane', 1.0)) + npe = nb.load(in_file).shape[_get_pe_index(in_meta)] + etl = npe // acc + + # Use case 2: TRT is defined + trt = in_meta.get('TotalReadoutTime', None) + if trt is not None: + return trt / (etl - 1) + + # Use case 3 (philips scans) + wfs = in_meta.get('WaterFatShift', None) + if wfs is not None: + fstrength = in_meta['MagneticFieldStrength'] + wfd_ppm = 3.4 # water-fat diff in ppm + g_ratio_mhz_t = 42.57 # gyromagnetic ratio for proton (1H) in MHz/T + wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t + return wfs / (wfs_hz * etl) + + raise ValueError('Unknown effective echo-spacing specification') + + +def get_trt(in_meta, in_file=None): + r""" + Extract the *total readout time* :math:`t_\text{RO}` from BIDS. + + Calculate the *total readout time* for an input + :abbr:`EPI (echo-planar imaging)` scan. + + There are several procedures to calculate the total + readout time. The basic one is that a ``TotalReadoutTime`` + field is set in the JSON sidecar. The following examples + use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the + j-axis encoding direction. + + >>> meta = {'TotalReadoutTime': 0.02596} + >>> get_trt(meta) + 0.02596 + + If the *effective echo spacing* :math:`t_\text{ees}` + (``EffectiveEchoSpacing`` BIDS field) is provided, then the + total readout time can be calculated reading the number + of voxels along the readout direction :math:`T_\text{ro}` + and the parallel acceleration factor of the EPI :math:`f_\text{acc}`. + + .. math :: + + T_\text{ro} = t_\text{ees} \, (N_\text{PE} / f_\text{acc} - 1) + + >>> meta = {'EffectiveEchoSpacing': 0.00059, + ... 'PhaseEncodingDirection': 'j-', + ... 'ParallelReductionFactorInPlane': 2} + >>> get_trt(meta, in_file='epi.nii.gz') + 0.02596 + + Some vendors, like Philips, store different parameter names: + + >>> meta = {'WaterFatShift': 8.129, + ... 'MagneticFieldStrength': 3, + ... 'PhaseEncodingDirection': 'j-', + ... 'ParallelReductionFactorInPlane': 2} + >>> get_trt(meta, in_file='epi.nii.gz') + 0.018721183563864822 + + """ + # Use case 1: TRT is defined + trt = in_meta.get('TotalReadoutTime', None) + if trt is not None: + return trt + + # All other cases require the parallel acc and npe (N vox in PE dir) + acc = float(in_meta.get('ParallelReductionFactorInPlane', 1.0)) + npe = nb.load(in_file).shape[_get_pe_index(in_meta)] + etl = npe // acc + + # Use case 2: TRT is defined + ees = in_meta.get('EffectiveEchoSpacing', None) + if ees is not None: + return ees * (etl - 1) + + # Use case 3 (philips scans) + wfs = in_meta.get('WaterFatShift', None) + if wfs is not None: + fstrength = in_meta['MagneticFieldStrength'] + wfd_ppm = 3.4 # water-fat diff in ppm + g_ratio_mhz_t = 42.57 # gyromagnetic ratio for proton (1H) in MHz/T + wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t + return wfs / wfs_hz + + raise ValueError('Unknown total-readout time specification') + + +def _get_pe_index(meta): + pe = meta['PhaseEncodingDirection'] + try: + return {'i': 0, 'j': 1, 'k': 2}[pe[0]] + except KeyError: + raise RuntimeError('"%s" is an invalid PE string' % pe) + + +def _torads(in_file, fmap_range=None, newpath=None): + """ + Convert a field map to rad/s units. + + If fmap_range is None, the range of the fieldmap + will be automatically calculated. + + Use fmap_range=0.5 to convert from Hz to rad/s + + """ + from math import pi + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + out_file = fname_presuffix(in_file, suffix='_rad', newpath=newpath) + fmapnii = nb.load(in_file) + fmapdata = fmapnii.get_fdata(dtype='float32') + + if fmap_range is None: + fmap_range = max(abs(fmapdata.min()), fmapdata.max()) + fmapdata = fmapdata * (pi / fmap_range) + out_img = nb.Nifti1Image(fmapdata, fmapnii.affine, fmapnii.header) + out_img.set_data_dtype('float32') + out_img.to_filename(out_file) + return out_file, fmap_range + + +def _tohz(in_file, range_hz, newpath=None): + """Convert a field map to Hz units.""" + from math import pi + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + out_file = fname_presuffix(in_file, suffix='_hz', newpath=newpath) + fmapnii = nb.load(in_file) + fmapdata = fmapnii.get_fdata(dtype='float32') + fmapdata = fmapdata * (range_hz / pi) + out_img = nb.Nifti1Image(fmapdata, fmapnii.affine, fmapnii.header) + out_img.set_data_dtype('float32') + out_img.to_filename(out_file) + return out_file + + +def phdiff2fmap(in_file, delta_te, newpath=None): + r""" + Convert the input phase-difference map into a fieldmap in Hz. + + Uses eq. (1) of [Hutton2002]_: + + .. math:: + + \Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}} + + + In this case, we do not take into account the gyromagnetic ratio of the + proton (:math:`\gamma`), since it will be applied inside TOPUP: + + .. math:: + + \Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}} + + References + ---------- + .. [Hutton2002] Hutton et al., Image Distortion Correction in fMRI: A Quantitative + Evaluation, NeuroImage 16(1):217-240, 2002. doi:`10.1006/nimg.2001.1054 + `_. + + + """ + import math + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + # GYROMAG_RATIO_H_PROTON_MHZ = 42.576 + + out_file = fname_presuffix(in_file, suffix='_fmap', newpath=newpath) + image = nb.load(in_file) + data = (image.get_fdata(dtype='float32') / (2. * math.pi * delta_te)) + nii = nb.Nifti1Image(data, image.affine, image.header) + nii.set_data_dtype(np.float32) + nii.to_filename(out_file) + return out_file + + +def _delta_te(in_values, te1=None, te2=None): + r"""Read :math:`\Delta_\text{TE}` from BIDS metadata dict.""" + if isinstance(in_values, float): + te2 = in_values + te1 = 0. + + if isinstance(in_values, dict): + te1 = in_values.get('EchoTime1') + te2 = in_values.get('EchoTime2') + + if not all((te1, te2)): + te2 = in_values.get('EchoTimeDifference') + te1 = 0 + + if isinstance(in_values, list): + te2, te1 = in_values + if isinstance(te1, list): + te1 = te1[1] + if isinstance(te2, list): + te2 = te2[1] + + # For convienience if both are missing we should give one error about them + if te1 is None and te2 is None: + raise RuntimeError('EchoTime1 and EchoTime2 metadata fields not found. ' + 'Please consult the BIDS specification.') + if te1 is None: + raise RuntimeError( + 'EchoTime1 metadata field not found. Please consult the BIDS specification.') + if te2 is None: + raise RuntimeError( + 'EchoTime2 metadata field not found. Please consult the BIDS specification.') + + return abs(float(te2) - float(te1)) + + +def au2rads(in_file, newpath=None): + """Convert the input phase difference map in arbitrary units (a.u.) to rads.""" + im = nb.load(in_file) + data = im.get_fdata(caching='unchanged') # Read as float64 for safety + hdr = im.header.copy() + + # Rescale to [0, 2*pi] + data = (data - data.min()) * (2 * np.pi / (data.max() - data.min())) + + # Round to float32 and clip + data = np.clip(np.float32(data), 0.0, 2 * np.pi) + + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units('mm') + out_file = fname_presuffix(in_file, suffix='_rads', newpath=newpath) + nb.Nifti1Image(data, None, hdr).to_filename(out_file) + return out_file + + +def _subtract_phases(in_phases, in_meta, newpath=None): + # Discard traits with copy(), so that pop() works. + in_meta = (in_meta[0].copy(), in_meta[1].copy()) + echo_times = tuple([m.pop('EchoTime', None) for m in in_meta]) + if not all(echo_times): + raise ValueError( + 'One or more missing EchoTime metadata parameter ' + 'associated to one or more phase map(s).') + + if echo_times[0] > echo_times[1]: + in_phases = (in_phases[1], in_phases[0]) + in_meta = (in_meta[1], in_meta[0]) + echo_times = (echo_times[1], echo_times[0]) + + in_phases_nii = [nb.load(ph) for ph in in_phases] + sub_data = in_phases_nii[1].get_fdata(dtype='float32') - \ + in_phases_nii[0].get_fdata(dtype='float32') + + # wrap negative radians back to [0, 2pi] + sub_data[sub_data < 0] += 2 * np.pi + sub_data = np.clip(sub_data, 0.0, 2 * np.pi) + + new_meta = in_meta[1].copy() + new_meta.update(in_meta[0]) + new_meta['EchoTime1'] = echo_times[0] + new_meta['EchoTime2'] = echo_times[1] + + hdr = in_phases_nii[0].header.copy() + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units('mm') + nii = nb.Nifti1Image(sub_data, in_phases_nii[0].affine, hdr) + out_phdiff = fname_presuffix(in_phases[0], suffix='_phdiff', + newpath=newpath) + nii.to_filename(out_phdiff) + return out_phdiff, new_meta diff --git a/aslprep/sdcflows/interfaces/reportlets.py b/aslprep/sdcflows/interfaces/reportlets.py new file mode 100644 index 000000000..9ed584938 --- /dev/null +++ b/aslprep/sdcflows/interfaces/reportlets.py @@ -0,0 +1,96 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Interfaces to generate speciality reportlets.""" +import numpy as np +from nilearn.image import threshold_img, load_img +from ...niworkflows import NIWORKFLOWS_LOG +from ...niworkflows.viz.utils import cuts_from_bbox, compose_view +from nipype.interfaces.base import File, isdefined, traits +from nipype.interfaces.mixins import reporting + +from ..viz.utils import plot_registration, coolwarm_transparent + + +class _FieldmapReportletInputSpec(reporting.ReportCapableInputSpec): + reference = File(exists=True, mandatory=True, desc='input reference') + moving = File(exists=True, desc='input moving') + fieldmap = File(exists=True, mandatory=True, desc='input fieldmap') + max_alpha = traits.Float(0.7, usedefault=True, desc='maximum alpha channel') + mask = File(exists=True, desc='brain mask') + out_report = File('report.svg', usedefault=True, + desc='filename for the visual report') + show = traits.Enum(1, 0, 'both', usedefault=True, + desc='where the fieldmap should be shown') + reference_label = traits.Str('Reference', usedefault=True, + desc='a label name for the reference mosaic') + moving_label = traits.Str('Fieldmap (Hz)', usedefault=True, + desc='a label name for the reference mosaic') + + +class FieldmapReportlet(reporting.ReportCapableInterface): + """An abstract mixin to registration nipype interfaces.""" + + _n_cuts = 7 + input_spec = _FieldmapReportletInputSpec + output_spec = reporting.ReportCapableOutputSpec + + def __init__(self, **kwargs): + """Instantiate FieldmapReportlet.""" + self._n_cuts = kwargs.pop('n_cuts', self._n_cuts) + super(FieldmapReportlet, self).__init__(generate_report=True, **kwargs) + + def _run_interface(self, runtime): + return runtime + + def _generate_report(self): + """Generate a reportlet.""" + NIWORKFLOWS_LOG.info('Generating visual report') + + movnii = refnii = load_img(self.inputs.reference) + fmapnii = load_img(self.inputs.fieldmap) + + if isdefined(self.inputs.moving): + movnii = load_img(self.inputs.moving) + + contour_nii = mask_nii = None + if isdefined(self.inputs.mask): + contour_nii = load_img(self.inputs.mask) + maskdata = contour_nii.get_fdata() > 0 + else: + mask_nii = threshold_img(refnii, 1e-3) + maskdata = mask_nii.get_fdata() > 0 + cuts = cuts_from_bbox(contour_nii or mask_nii, cuts=self._n_cuts) + fmapdata = fmapnii.get_fdata() + vmax = max(abs(np.percentile(fmapdata[maskdata], 99.8)), + abs(np.percentile(fmapdata[maskdata], 0.2))) + + fmap_overlay = [{ + 'overlay': fmapnii, + 'overlay_params': { + 'cmap': coolwarm_transparent(max_alpha=self.inputs.max_alpha), + 'vmax': vmax, + 'vmin': -vmax, + } + }] * 2 + + if self.inputs.show != 'both': + fmap_overlay[not self.inputs.show] = {} + + # Call composer + compose_view( + plot_registration(movnii, 'moving-image', + estimate_brightness=True, + cuts=cuts, + label=self.inputs.moving_label, + contour=contour_nii, + compress=False, + **fmap_overlay[1]), + plot_registration(refnii, 'fixed-image', + estimate_brightness=True, + cuts=cuts, + label=self.inputs.reference_label, + contour=contour_nii, + compress=False, + **fmap_overlay[0]), + out_file=self._out_report + ) diff --git a/aslprep/sdcflows/tests/__init__.py b/aslprep/sdcflows/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aslprep/sdcflows/tests/test_version.py b/aslprep/sdcflows/tests/test_version.py new file mode 100644 index 000000000..ac8be6cf8 --- /dev/null +++ b/aslprep/sdcflows/tests/test_version.py @@ -0,0 +1,42 @@ +"""Test _version.py.""" +import sys +from collections import namedtuple +from pkg_resources import DistributionNotFound +from importlib import reload +import sdcflows + + +def test_version_scm0(monkeypatch): + """Retrieve the version via setuptools_scm.""" + + class _version: + __version__ = "10.0.0" + + monkeypatch.setitem(sys.modules, "sdcflows._version", _version) + reload(sdcflows) + assert sdcflows.__version__ == "10.0.0" + + +def test_version_scm1(monkeypatch): + """Retrieve the version via pkg_resources.""" + monkeypatch.setitem(sys.modules, "sdcflows._version", None) + + def _dist(name): + Distribution = namedtuple("Distribution", ["name", "version"]) + return Distribution(name, "success") + + monkeypatch.setattr("pkg_resources.get_distribution", _dist) + reload(sdcflows) + assert sdcflows.__version__ == "success" + + +def test_version_scm2(monkeypatch): + """Check version could not be interpolated.""" + monkeypatch.setitem(sys.modules, "sdcflows._version", None) + + def _raise(name): + raise DistributionNotFound("No get_distribution mock") + + monkeypatch.setattr("pkg_resources.get_distribution", _raise) + reload(sdcflows) + assert sdcflows.__version__ == "unknown" diff --git a/aslprep/sdcflows/viz/__init__.py b/aslprep/sdcflows/viz/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aslprep/sdcflows/viz/utils.py b/aslprep/sdcflows/viz/utils.py new file mode 100644 index 000000000..102df7939 --- /dev/null +++ b/aslprep/sdcflows/viz/utils.py @@ -0,0 +1,97 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Visualization tooling.""" + + +def plot_registration(anat_nii, div_id, plot_params=None, + order=('z', 'x', 'y'), cuts=None, + estimate_brightness=False, label=None, contour=None, + compress='auto', overlay=None, overlay_params=None): + """ + Plot the foreground and background views. + + Default order is: axial, coronal, sagittal + """ + from uuid import uuid4 + + from lxml import etree + import matplotlib.pyplot as plt + from nilearn.plotting import plot_anat + from svgutils.transform import SVGFigure + from niworkflows.viz.utils import robust_set_limits, extract_svg, SVGNS + + plot_params = plot_params or {} + + # Use default MNI cuts if none defined + if cuts is None: + raise NotImplementedError # TODO + + out_files = [] + if estimate_brightness: + plot_params = robust_set_limits( + anat_nii.get_fdata(dtype='float32').reshape(-1), + plot_params) + + # Plot each cut axis + for i, mode in enumerate(list(order)): + plot_params['display_mode'] = mode + plot_params['cut_coords'] = cuts[mode] + if i == 0: + plot_params['title'] = label + else: + plot_params['title'] = None + + # Generate nilearn figure + display = plot_anat(anat_nii, **plot_params) + if overlay is not None: + _overlay_params = { + 'vmin': overlay.get_fdata(dtype='float32').min(), + 'vmax': overlay.get_fdata(dtype='float32').max(), + 'cmap': plt.cm.gray, + 'interpolation': 'nearest', + } + _overlay_params.update(overlay_params) + display.add_overlay(overlay, **_overlay_params) + if contour is not None: + display.add_contours(contour, colors='g', levels=[0.5], + linewidths=0.5) + + svg = extract_svg(display, compress=compress) + display.close() + + # Find and replace the figure_1 id. + xml_data = etree.fromstring(svg) + find_text = etree.ETXPath("//{%s}g[@id='figure_1']" % SVGNS) + find_text(xml_data)[0].set('id', '%s-%s-%s' % (div_id, mode, uuid4())) + + svg_fig = SVGFigure() + svg_fig.root = xml_data + out_files.append(svg_fig) + + return out_files + + +def coolwarm_transparent(max_alpha=0.7, opaque_perc=30, transparent_perc=8): + """Modify the coolwarm color scale to have full transparency around the middle.""" + import numpy as np + import matplotlib.pylab as pl + from matplotlib.colors import ListedColormap + + # Choose colormap + cmap = pl.cm.coolwarm + + # Get the colormap colors + my_cmap = cmap(np.arange(cmap.N)) + + _20perc = (cmap.N * opaque_perc) // 100 + midpoint = cmap.N // 2 + 1 + _10perc = (cmap.N * transparent_perc) // 100 + # Set alpha + alpha = np.ones(cmap.N) * max_alpha + alpha[midpoint - _10perc:midpoint + _10perc] = 0 + alpha[_20perc:midpoint - _10perc - 1] = np.linspace( + max_alpha, 0, len(alpha[_20perc:midpoint - _10perc - 1])) + alpha[midpoint + _10perc:-_20perc] = np.linspace( + 0, max_alpha, len(alpha[midpoint + _10perc:-_20perc])) + my_cmap[:, -1] = alpha + return ListedColormap(my_cmap) diff --git a/aslprep/sdcflows/workflows/__init__.py b/aslprep/sdcflows/workflows/__init__.py new file mode 100644 index 000000000..80f930255 --- /dev/null +++ b/aslprep/sdcflows/workflows/__init__.py @@ -0,0 +1,37 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Fieldmap estimation and unwarping workflows. + +.. _sdc_base : + +Automatic selection of the appropriate susceptibility-distortion correction (SDC) method +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +If the dataset metadata indicate that more than one field map acquisition is +``IntendedFor`` (see the `BIDS Specification +`__), +the following priority will be used: + + 1. :ref:`sdc_pepolar` (or **blip-up/blip-down**) + + 2. :ref:`sdc_direct_b0` + + 3. :ref:`sdc_phasediff` + + 4. :ref:`sdc_fieldmapless` + + +Table of behavior (fieldmap use-cases): + +=============== =========== ============= =============== +Fieldmaps found ``use_syn`` ``force_syn`` Action +=============== =========== ============= =============== +True * True Fieldmaps + SyN +True * False Fieldmaps +False * True SyN +False True False SyN +False False False HMC only +=============== =========== ============= =============== + +""" diff --git a/aslprep/sdcflows/workflows/base.py b/aslprep/sdcflows/workflows/base.py new file mode 100644 index 000000000..c6e4cd043 --- /dev/null +++ b/aslprep/sdcflows/workflows/base.py @@ -0,0 +1,295 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""SDC workflows coordination.""" +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu +from nipype import logging + +from ...niworkflows.engine.workflows import LiterateWorkflow as Workflow + + +LOGGER = logging.getLogger('nipype.workflow') +FMAP_PRIORITY = { + 'epi': 0, + 'fieldmap': 1, + 'phasediff': 2, + 'syn': 3, +} + +DEFAULT_MEMORY_MIN_GB = 0.01 + + +def init_sdc_estimate_wf(fmaps, epi_meta, omp_nthreads=1, debug=False): + """ + Build a :abbr:`SDC (susceptibility distortion correction)` workflow. + + This workflow implements the heuristics to choose an estimation + methodology for :abbr:`SDC (susceptibility distortion correction)`. + When no field map information is present within the BIDS inputs, + the EXPERIMENTAL "fieldmap-less SyN" can be performed, using + the ``--use-syn`` argument. When ``--force-syn`` is specified, + then the "fieldmap-less SyN" is always executed and reported + despite of other fieldmaps available with higher priority. + In the latter case (some sort of fieldmap(s) is available and + ``--force-syn`` is requested), then the :abbr:`SDC (susceptibility + distortion correction)` method applied is that with the + highest priority. + + Parameters + ---------- + fmaps : list of pybids dicts + A list of dictionaries with the available fieldmaps + (and their metadata using the key ``'metadata'`` for the + case of :abbr:`PEPOLAR (Phase-Encoding POLARity)` fieldmaps). + epi_meta : dict + BIDS metadata dictionary corresponding to the + :abbr:`EPI (echo-planar imaging)` run (i.e., suffix ``bold``, + ``sbref``, or ``dwi``) for which the fieldmap is being estimated. + omp_nthreads : int + Maximum number of threads an individual process may use + debug : bool + Enable debugging outputs + + Inputs + ------ + epi_file + A reference image calculated at a previous stage + epi_brain + Same as above, but brain-masked + epi_mask + Brain mask for the run + t1w_brain + T1w image, brain-masked, for the fieldmap-less SyN method + std2anat_xfm + Standard-to-T1w transform generated during spatial + normalization (only for the fieldmap-less SyN method). + + Outputs + ------- + epi_corrected + The EPI scan reference after unwarping. + epi_mask + The corresponding new mask after unwarping + epi_brain + Brain-extracted, unwarped EPI scan reference + out_warp + The deformation field to unwarp the susceptibility distortions + syn_ref + If ``--force-syn``, an unwarped EPI scan reference with this + method (for reporting purposes) + method : str + Short description of the estimation method that was run. + + """ + workflow = Workflow(name='sdc_estimate_wf' if fmaps else 'sdc_bypass_wf') + inputnode = pe.Node(niu.IdentityInterface( + fields=['epi_file', 'epi_brain', 'epi_mask', 't1w_brain', 'std2anat_xfm']), + name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface( + fields=['epi_corrected', 'epi_mask', 'epi_brain', + 'out_warp', 'syn_ref', 'method']), + name='outputnode') + + # No fieldmaps - forward inputs to outputs + if not fmaps: + workflow.__postdesc__ = """\ +Susceptibility distortion correction (SDC) was omitted. +""" + outputnode.inputs.method = 'None' + outputnode.inputs.out_warp = 'identity' + workflow.connect([ + (inputnode, outputnode, [('epi_file', 'epi_corrected'), + ('epi_mask', 'epi_mask'), + ('epi_brain', 'epi_brain')]), + ]) + return workflow + + workflow.__postdesc__ = """\ +Based on the estimated susceptibility distortion, a corrected +EPI (echo-planar imaging) reference was calculated for a more +accurate co-registration with the anatomical reference. +""" + + only_syn = 'syn' in fmaps and len(fmaps) == 1 + + # PEPOLAR path + if 'epi' in fmaps: + from .pepolar import init_pepolar_unwarp_wf, check_pes + + # SyN works without this metadata + if epi_meta.get('PhaseEncodingDirection') is None: + raise ValueError( + 'PhaseEncodingDirection is not defined within the metadata retrieved ' + 'for the intended EPI (DWI, BOLD, or SBRef) run.') + outputnode.inputs.method = 'PEB/PEPOLAR (phase-encoding based / PE-POLARity)' + + fmaps_epi = [(v[0], v[1].get('PhaseEncodingDirection')) + for v in fmaps['epi']] + + if not all(list(zip(*fmaps_epi))[1]): + raise ValueError( + 'At least one of the EPI runs with alternative phase-encoding ' + 'blips is missing the required "PhaseEncodingDirection" metadata entry.') + + # Find matched PE directions + matched_pe = check_pes(fmaps_epi, epi_meta['PhaseEncodingDirection']) + + # Get EPI polarities and their metadata + sdc_unwarp_wf = init_pepolar_unwarp_wf( + matched_pe=matched_pe, + omp_nthreads=omp_nthreads) + sdc_unwarp_wf.inputs.inputnode.epi_pe_dir = epi_meta['PhaseEncodingDirection'] + sdc_unwarp_wf.inputs.inputnode.fmaps_epi = fmaps_epi + + workflow.connect([ + (inputnode, sdc_unwarp_wf, [ + ('epi_file', 'inputnode.in_reference'), + ('epi_brain', 'inputnode.in_reference_brain'), + ('epi_mask', 'inputnode.in_mask')]), + ]) + + # FIELDMAP path + elif 'fieldmap' in fmaps or 'phasediff' in fmaps: + from .fmap import init_fmap2field_wf + from .unwarp import init_sdc_unwarp_wf + + # SyN works without this metadata + if epi_meta.get('PhaseEncodingDirection') is None: + raise ValueError( + 'PhaseEncodingDirection is not defined within the metadata retrieved ' + 'for the intended EPI (DWI, BOLD, or SBRef) run.') + + if 'fieldmap' in fmaps: + from .fmap import init_fmap_wf + try: + fmap, = fmaps['fieldmap'] + except ValueError: + LOGGER.warning('Several B0 fieldmaps found for the given target, using ' + 'the first one.') + fmap = fmaps['fieldmap'][0] + + outputnode.inputs.method = 'FMB (fieldmap-based) - directly measured B0 map' + fmap_wf = init_fmap_wf( + omp_nthreads=omp_nthreads, + fmap_bspline=False) + # set inputs + fmap_wf.inputs.inputnode.magnitude = [ + m for m, _ in fmap['magnitude']] + fmap_wf.inputs.inputnode.fieldmap = [ + m for m, _ in fmap['fieldmap']] + elif 'phasediff' in fmaps: + from .phdiff import init_phdiff_wf + try: + fmap, = fmaps['phasediff'] + except ValueError: + LOGGER.warning('Several phase-difference maps found for the given target, using ' + 'the first one.') + fmap = fmaps['phasediff'][0] + + outputnode.inputs.method = 'FMB (fieldmap-based) - phase-difference map' + fmap_wf = init_phdiff_wf(omp_nthreads=omp_nthreads) + # set inputs + fmap_wf.inputs.inputnode.magnitude = [ + m for m, _ in fmap['magnitude']] + fmap_wf.inputs.inputnode.phasediff = fmap['phases'] + + fmap2field_wf = init_fmap2field_wf(omp_nthreads=omp_nthreads, debug=debug) + fmap2field_wf.inputs.inputnode.metadata = epi_meta + + sdc_unwarp_wf = init_sdc_unwarp_wf( + omp_nthreads=omp_nthreads, + debug=debug, + name='sdc_unwarp_wf') + + workflow.connect([ + (inputnode, fmap2field_wf, [ + ('epi_file', 'inputnode.in_reference'), + ('epi_brain', 'inputnode.in_reference_brain')]), + (inputnode, sdc_unwarp_wf, [ + ('epi_file', 'inputnode.in_reference'), + ('epi_mask', 'inputnode.in_reference_mask')]), + (fmap_wf, fmap2field_wf, [ + ('outputnode.fmap', 'inputnode.fmap'), + ('outputnode.fmap_ref', 'inputnode.fmap_ref'), + ('outputnode.fmap_mask', 'inputnode.fmap_mask')]), + (fmap2field_wf, sdc_unwarp_wf, [ + ('outputnode.out_warp', 'inputnode.in_warp')]), + + ]) + elif not only_syn: + raise ValueError('Fieldmaps of types %s are not supported' % + ', '.join(['"%s"' % f for f in fmaps])) + + # FIELDMAP-less path + if 'syn' in fmaps: + from .syn import init_syn_sdc_wf + syn_sdc_wf = init_syn_sdc_wf( + epi_pe=epi_meta.get('PhaseEncodingDirection', None), + omp_nthreads=omp_nthreads) + + workflow.connect([ + (inputnode, syn_sdc_wf, [ + ('epi_file', 'inputnode.in_reference'), + ('epi_brain', 'inputnode.in_reference_brain'), + ('t1w_brain', 'inputnode.t1w_brain'), + ('std2anat_xfm', 'inputnode.std2anat_xfm')]), + (syn_sdc_wf, outputnode, [ + ('outputnode.out_reference', 'syn_ref')]), + ]) + + # XXX Eliminate branch when forcing isn't an option + if only_syn: # No fieldmaps, but --use-syn + outputnode.inputs.method = 'FLB ("fieldmap-less", SyN-based)' + sdc_unwarp_wf = syn_sdc_wf + else: # --force-syn was called when other fieldmap was present + sdc_unwarp_wf.__desc__ = None + + workflow.connect([ + (sdc_unwarp_wf, outputnode, [ + ('outputnode.out_warp', 'out_warp'), + ('outputnode.out_reference', 'epi_corrected'), + ('outputnode.out_reference_brain', 'epi_brain'), + ('outputnode.out_mask', 'epi_mask')]), + ]) + + return workflow + + +def fieldmap_wrangler(layout, target_image, use_syn=False, force_syn=False): + """Query the BIDSLayout for fieldmaps, and arrange them for the orchestration workflow.""" + from collections import defaultdict + fmap_bids = layout.get_fieldmap(target_image, return_list=True) + fieldmaps = defaultdict(list) + for fmap in fmap_bids: + if fmap['suffix'] == 'epi': + fieldmaps['epi'].append((fmap['epi'], layout.get_metadata(fmap['epi']))) + + if fmap['suffix'] == 'fieldmap': + fieldmaps['fieldmap'].append({ + 'magnitude': [(fmap['magnitude'], layout.get_metadata(fmap['magnitude']))], + 'fieldmap': [(fmap['fieldmap'], layout.get_metadata(fmap['fieldmap']))], + }) + + if fmap['suffix'] == 'phasediff': + fieldmaps['phasediff'].append({ + 'magnitude': [(fmap[k], layout.get_metadata(fmap[k])) + for k in sorted(fmap.keys()) if k.startswith('magnitude')], + 'phases': [(fmap['phasediff'], layout.get_metadata(fmap['phasediff']))], + }) + + if fmap['suffix'] == 'phase': + fieldmaps['phasediff'].append({ + 'magnitude': [(fmap[k], layout.get_metadata(fmap[k])) + for k in sorted(fmap.keys()) if k.startswith('magnitude')], + 'phases': [(fmap[k], layout.get_metadata(fmap[k])) + for k in sorted(fmap.keys()) if k.startswith('phase')], + }) + + if fieldmaps and force_syn: + # syn: True -> Run SyN in addition to fieldmap-based SDC + fieldmaps['syn'] = True + elif not fieldmaps and (force_syn or use_syn): + # syn: False -> Run SyN as only SDC + fieldmaps['syn'] = False + return fieldmaps diff --git a/aslprep/sdcflows/workflows/fmap.py b/aslprep/sdcflows/workflows/fmap.py new file mode 100644 index 000000000..a7a3c01b2 --- /dev/null +++ b/aslprep/sdcflows/workflows/fmap.py @@ -0,0 +1,267 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Fieldmap-based estimation of susceptibility distortions. + +.. _sdc_direct_b0 : + +Direct B0 mapping sequences +~~~~~~~~~~~~~~~~~~~~~~~~~~~ +When the fieldmap is directly measured with a prescribed sequence (such as +:abbr:`SE (spiral echo)`), we only need to calculate the corresponding +displacements field that accounts for the distortions. +This procedure is described with more detail `here +`__. + +This corresponds to `this section of the BIDS specification +`__. + +""" +import pkg_resources as pkgr + +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu, fsl +from ...niworkflows.engine.workflows import LiterateWorkflow as Workflow +from ...niworkflows.interfaces.images import IntraModalMerge +from ...niworkflows.interfaces.registration import ANTSApplyTransformsRPT, ANTSRegistrationRPT + +from ..interfaces.fmap import get_ees as _get_ees, FieldToRadS, FUGUEvsm2ANTSwarp +from .gre import init_fmap_postproc_wf, init_magnitude_wf + + +def init_fmap_wf(omp_nthreads, fmap_bspline, name='fmap_wf'): + """ + Estimate the fieldmap based on a field-mapping MRI acquisition. + + When we have a sequence that directly measures the fieldmap, + we just need to mask it (using the corresponding magnitude image) + to remove the noise in the surrounding air region, and ensure that + units are Hz. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.fmap import init_fmap_wf + wf = init_fmap_wf(omp_nthreads=6, fmap_bspline=False) + + Parameters + ---------- + omp_nthreads : int + Maximum number of threads an individual process may use. + fmap_bspline : bool + Whether the fieldmap estimate will be smoothed using BSpline basis. + name : str + Unique name of this workflow. + + Inputs + ------ + magnitude : str + Path to the corresponding magnitude image for anatomical reference. + fieldmap : str + Path to the fieldmap acquisition (``*_fieldmap.nii[.gz]`` of BIDS). + + Outputs + ------- + fmap : str + Path to the estimated fieldmap. + fmap_ref : str + Path to a preprocessed magnitude image reference. + fmap_mask : str + Path to a binary brain mask corresponding to the ``fmap`` and ``fmap_ref`` + pair. + + """ + workflow = Workflow(name=name) + workflow.__desc__ = """\ +A B0-nonuniformity map (or *fieldmap*) was directly measured with an MRI scheme +designed with that purpose (typically, a spiral pulse sequence). +""" + inputnode = pe.Node(niu.IdentityInterface( + fields=['magnitude', 'fieldmap']), name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['fmap', 'fmap_ref', 'fmap_mask']), + name='outputnode') + + magnitude_wf = init_magnitude_wf(omp_nthreads=omp_nthreads) + workflow.connect([ + (inputnode, magnitude_wf, [('magnitude', 'inputnode.magnitude')]), + (magnitude_wf, outputnode, [('outputnode.fmap_mask', 'fmap_mask'), + ('outputnode.fmap_ref', 'fmap_ref')]), + ]) + + # Merge input fieldmap images + fmapmrg = pe.Node(IntraModalMerge(zero_based_avg=False, hmc=False), + name='fmapmrg') + applymsk = pe.Node(fsl.ApplyMask(), name='applymsk') + fmap_postproc_wf = init_fmap_postproc_wf(omp_nthreads=omp_nthreads, + fmap_bspline=fmap_bspline) + + workflow.connect([ + (inputnode, fmapmrg, [('fieldmap', 'in_files')]), + (fmapmrg, applymsk, [('out_avg', 'in_file')]), + (magnitude_wf, applymsk, [('outputnode.fmap_mask', 'mask_file')]), + (applymsk, fmap_postproc_wf, [('out_file', 'inputnode.fmap')]), + (magnitude_wf, fmap_postproc_wf, [ + ('outputnode.fmap_mask', 'inputnode.fmap_mask'), + ('outputnode.fmap_ref', 'inputnode.fmap_ref')]), + (fmap_postproc_wf, outputnode, [('outputnode.out_fmap', 'fmap')]), + ]) + return workflow + + +def init_fmap2field_wf(omp_nthreads, debug, name='fmap2field_wf', + generate_report=True): + """ + Convert the estimated fieldmap in Hz into a displacements field. + + This workflow takes in a fieldmap and calculates the corresponding + displacements field (in other words, an ANTs-compatible warp file). + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.fmap import init_fmap2field_wf + wf = init_fmap2field_wf(omp_nthreads=8, + debug=False) + + Parameters + ---------- + omp_nthreads : int + Maximum number of threads an individual process may use. + debug : bool + Run fast configurations of registrations. + name : str + Unique name of this workflow. + + Inputs + ------ + in_reference + the reference image + in_reference_brain + the reference image (skull-stripped) + metadata + metadata associated to the ``in_reference`` EPI input + fmap + the fieldmap in Hz + fmap_ref + the reference (anatomical) image corresponding to ``fmap`` + fmap_mask + a brain mask corresponding to ``fmap`` + + + Outputs + ------- + out_reference + the ``in_reference`` after unwarping + out_reference_brain + the ``in_reference`` after unwarping and skullstripping + out_warp + the corresponding :abbr:`DFM (displacements field map)` compatible with + ANTs + out_jacobian + the jacobian of the field (for drop-out alleviation) + out_mask + mask of the unwarped input file + + """ + workflow = Workflow(name=name) + workflow.__desc__ = """\ +The *fieldmap* was then co-registered to the target EPI (echo-planar imaging) +reference run and converted to a displacements field map (amenable to registration +tools such as ANTs) with FSL's `fugue` and other *SDCflows* tools. +""" + inputnode = pe.Node(niu.IdentityInterface( + fields=['in_reference', 'in_reference_brain', 'metadata', + 'fmap_ref', 'fmap_mask', 'fmap']), name='inputnode') + outputnode = pe.Node(niu.IdentityInterface( + fields=['out_warp']), name='outputnode') + + # Register the reference of the fieldmap to the reference + # of the target image (the one that shall be corrected) + ants_settings = pkgr.resource_filename('sdcflows', 'data/fmap-any_registration.json') + if debug: + ants_settings = pkgr.resource_filename( + 'sdcflows', 'data/fmap-any_registration_testing.json') + + fmap2ref_reg = pe.Node( + ANTSRegistrationRPT(generate_report=False, from_file=ants_settings, + output_inverse_warped_image=True, output_warped_image=True), + name='fmap2ref_reg', n_procs=omp_nthreads) + + # Map the VSM into the EPI space + fmap2ref_apply = pe.Node(ANTSApplyTransformsRPT( + generate_report=False, dimension=3, interpolation='BSpline', float=True), + name='fmap2ref_apply') + + fmap_mask2ref_apply = pe.Node(ANTSApplyTransformsRPT( + generate_report=False, dimension=3, interpolation='MultiLabel', + float=True), + name='fmap_mask2ref_apply') + + # Fieldmap to rads and then to voxels (VSM - voxel shift map) + torads = pe.Node(FieldToRadS(fmap_range=0.5), name='torads') + + get_ees = pe.Node(niu.Function(function=_get_ees, output_names=['ees']), name='get_ees') + + gen_vsm = pe.Node(fsl.FUGUE(save_unmasked_shift=True), name='gen_vsm') + # Convert the VSM into a DFM (displacements field map) + # or: FUGUE shift to ANTS warping. + vsm2dfm = pe.Node(FUGUEvsm2ANTSwarp(), name='vsm2dfm') + + workflow.connect([ + (inputnode, fmap2ref_reg, [('fmap_ref', 'moving_image'), + ('in_reference_brain', 'fixed_image')]), + (inputnode, fmap2ref_apply, [('fmap', 'input_image'), + ('in_reference', 'reference_image')]), + (inputnode, fmap_mask2ref_apply, [('in_reference', 'reference_image'), + ('fmap_mask', 'input_image')]), + (inputnode, get_ees, [('in_reference', 'in_file'), + ('metadata', 'in_meta')]), + (inputnode, gen_vsm, [(('metadata', _get_pedir_fugue), 'unwarp_direction')]), + (inputnode, vsm2dfm, [(('metadata', _get_pedir_bids), 'pe_dir')]), + (fmap2ref_reg, fmap2ref_apply, [ + ('composite_transform', 'transforms')]), + (fmap2ref_reg, fmap_mask2ref_apply, [ + ('composite_transform', 'transforms')]), + (fmap2ref_apply, torads, [('output_image', 'in_file')]), + (fmap_mask2ref_apply, gen_vsm, [('output_image', 'mask_file')]), + (gen_vsm, vsm2dfm, [('shift_out_file', 'in_file')]), + (get_ees, gen_vsm, [('ees', 'dwell_time')]), + (torads, gen_vsm, [('out_file', 'fmap_in_file')]), + (vsm2dfm, outputnode, [('out_file', 'out_warp')]), + ]) + + if generate_report: + from niworkflows.interfaces.bids import DerivativesDataSink + from ..interfaces.reportlets import FieldmapReportlet + + fmap_rpt = pe.Node(FieldmapReportlet( + reference_label='EPI Reference', + moving_label='Magnitude', show='both'), name='fmap_rpt') + ds_report_sdc = pe.Node( + DerivativesDataSink(desc='fieldmap', suffix='bold', datatype='figures'), + name='ds_report_fmap', mem_gb=0.01, run_without_submitting=True + ) + + workflow.connect([ + (inputnode, fmap_rpt, [('in_reference', 'reference')]), + (fmap2ref_reg, fmap_rpt, [('warped_image', 'moving')]), + (fmap_mask2ref_apply, fmap_rpt, [('output_image', 'mask')]), + (vsm2dfm, fmap_rpt, [('fieldmap', 'fieldmap')]), + (fmap_rpt, ds_report_sdc, [('out_report', 'in_file')]), + ]) + + return workflow + + +# Helper functions +# ------------------------------------------------------------ +def _get_pedir_bids(in_dict): + return in_dict['PhaseEncodingDirection'] + + +def _get_pedir_fugue(in_dict): + return in_dict['PhaseEncodingDirection'].replace('i', 'x').replace('j', 'y').replace('k', 'z') diff --git a/aslprep/sdcflows/workflows/gre.py b/aslprep/sdcflows/workflows/gre.py new file mode 100644 index 000000000..8d15a8f0a --- /dev/null +++ b/aslprep/sdcflows/workflows/gre.py @@ -0,0 +1,228 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Processing phase-difference (aka :abbr:`GRE (gradient-recalled echo)`) fieldmaps. + +.. _gre-fieldmaps: + +Workflows for processing :abbr:`GRE (gradient recalled echo)` fieldmaps +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Workflows for preparing the magnitude part of :abbr:`GRE (gradient-recalled echo)` fieldmap +images and cleaning up the fieldmaps created from the phases or phasediff. + +""" + +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu, fsl, ants +from niflow.nipype1.workflows.dmri.fsl.utils import cleanup_edge_pipeline +from ...niworkflows.engine.workflows import LiterateWorkflow as Workflow +from ...niworkflows.interfaces.images import IntraModalMerge +from ...niworkflows.interfaces.masks import BETRPT + + +def init_magnitude_wf(omp_nthreads, name='magnitude_wf'): + """ + Prepare the magnitude part of :abbr:`GRE (gradient-recalled echo)` fieldmaps. + + Average (if not done already) the magnitude part of the + :abbr:`GRE (gradient recalled echo)` images, run N4 to + correct for B1 field nonuniformity, and skull-strip the + preprocessed magnitude. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.fmap import init_magnitude_wf + wf = init_magnitude_wf(omp_nthreads=6) + + Parameters + ---------- + omp_nthreads : int + Maximum number of threads an individual process may use + name : str + Name of workflow (default: ``prepare_magnitude_w``) + + Inputs + ------ + magnitude : pathlike + Path to the corresponding magnitude path(s). + + Outputs + ------- + fmap_ref : pathlike + Path to the fieldmap reference calculated in this workflow. + fmap_mask : pathlike + Path to a binary brain mask corresponding to the reference above. + + """ + workflow = Workflow(name=name) + inputnode = pe.Node( + niu.IdentityInterface(fields=['magnitude']), name='inputnode') + outputnode = pe.Node( + niu.IdentityInterface(fields=['fmap_ref', 'fmap_mask', 'mask_report']), + name='outputnode') + + # Merge input magnitude images + # Do not reorient to RAS to preserve the validity of PhaseEncodingDirection + magmrg = pe.Node(IntraModalMerge(hmc=False, to_ras=False), name='magmrg') + + # de-gradient the fields ("bias/illumination artifact") + n4_correct = pe.Node(ants.N4BiasFieldCorrection(dimension=3, copy_header=True), + name='n4_correct', n_procs=omp_nthreads) + bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True), + name='bet') + + workflow.connect([ + (inputnode, magmrg, [('magnitude', 'in_files')]), + (magmrg, n4_correct, [('out_avg', 'input_image')]), + (n4_correct, bet, [('output_image', 'in_file')]), + (bet, outputnode, [('mask_file', 'fmap_mask'), + ('out_file', 'fmap_ref'), + ('out_report', 'mask_report')]), + ]) + return workflow + + +def init_fmap_postproc_wf(omp_nthreads, fmap_bspline, median_kernel_size=5, + name='fmap_postproc_wf'): + """ + Postprocess a B0 map estimated elsewhere. + + This workflow denoises (mostly via smoothing) a B0 fieldmap. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.fmap import init_fmap_postproc_wf + wf = init_fmap_postproc_wf(omp_nthreads=6, fmap_bspline=False) + + Parameters + ---------- + omp_nthreads : int + Maximum number of threads an individual process may use + fmap_bspline : bool + Whether the fieldmap should be smoothed and extrapolated to off-brain regions + using B-Spline basis. + median_kernel_size : int + Size of the kernel when smoothing is done with a median filter. + name : str + Name of workflow (default: ``fmap_postproc_wf``) + + Inputs + ------ + fmap_mask : pathlike + A brain binary mask corresponding to this fieldmap. + fmap_ref : pathlike + A preprocessed magnitude/reference image for the fieldmap. + fmap : pathlike + A B0-field nonuniformity map (aka fieldmap) estimated elsewhere. + + Outputs + ------- + out_fmap : pathlike + Postprocessed fieldmap. + + """ + workflow = Workflow(name=name) + inputnode = pe.Node(niu.IdentityInterface( + fields=['fmap_mask', 'fmap_ref', 'fmap', 'metadata']), name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_fmap', 'metadata']), + name='outputnode') + if fmap_bspline: + from ..interfaces.fmap import FieldEnhance + # despike_threshold=1.0, mask_erode=1), + fmapenh = pe.Node( + FieldEnhance(unwrap=False, despike=False), + name='fmapenh', mem_gb=4, n_procs=omp_nthreads) + + workflow.connect([ + (inputnode, fmapenh, [('fmap_mask', 'in_mask'), + ('fmap_ref', 'in_magnitude'), + ('fmap_hz', 'in_file')]), + (fmapenh, outputnode, [('out_file', 'out_fmap')]), + ]) + + else: + recenter = pe.Node(niu.Function(function=_recenter), + name='recenter', run_without_submitting=True) + denoise = pe.Node(fsl.SpatialFilter( + operation='median', kernel_shape='sphere', + kernel_size=median_kernel_size), name='denoise') + demean = pe.Node(niu.Function(function=_demean), name='demean') + cleanup_wf = cleanup_edge_pipeline(name="cleanup_wf") + + workflow.connect([ + (inputnode, cleanup_wf, [('fmap_mask', 'inputnode.in_mask')]), + (inputnode, recenter, [(('fmap', _pop), 'in_file')]), + (recenter, denoise, [('out', 'in_file')]), + (denoise, demean, [('out_file', 'in_file')]), + (demean, cleanup_wf, [('out', 'inputnode.in_file')]), + (cleanup_wf, outputnode, [('outputnode.out_file', 'out_fmap')]), + (inputnode, outputnode, [(('metadata', _pop), 'metadata')]), + ]) + + return workflow + + +def _recenter(in_file): + """Recenter the phase-map distribution to the -pi..pi range.""" + from os import getcwd + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + nii = nb.load(in_file) + data = nii.get_fdata(dtype='float32') + msk = data != 0 + msk[data == 0] = False + data[msk] -= np.median(data[msk]) + + out_file = fname_presuffix(in_file, suffix='_recentered', + newpath=getcwd()) + nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file) + return out_file + + +def _demean(in_file, in_mask=None, usemode=True): + """ + Subtract the median (since it is robuster than the mean) from a map. + + Parameters + ---------- + usemode : bool + Use the mode instead of the median (should be even more robust + against outliers). + + """ + from os import getcwd + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + nii = nb.load(in_file) + data = nii.get_fdata(dtype='float32') + + msk = np.ones_like(data, dtype=bool) + if in_mask is not None: + msk[nb.load(in_mask).get_fdata(dtype='float32') < 1e-4] = False + + if usemode: + from scipy.stats import mode + data[msk] -= mode(data[msk], axis=None)[0][0] + else: + data[msk] -= np.median(data[msk], axis=None) + + out_file = fname_presuffix(in_file, suffix='_demean', + newpath=getcwd()) + nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file) + return out_file + + +def _pop(inlist): + if isinstance(inlist, (tuple, list)): + return inlist[0] + return inlist diff --git a/aslprep/sdcflows/workflows/outputs.py b/aslprep/sdcflows/workflows/outputs.py new file mode 100644 index 000000000..b1bc9a582 --- /dev/null +++ b/aslprep/sdcflows/workflows/outputs.py @@ -0,0 +1,82 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Writing out outputs.""" +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu +from ...niworkflows.interfaces.bids import DerivativesDataSink + + +def init_sdc_unwarp_report_wf(name='sdc_unwarp_report_wf', forcedsyn=False): + """ + Save a reportlet showing how SDC unwarping performed. + + This workflow generates and saves a reportlet showing the effect of fieldmap + unwarping a BOLD image. + + Workflow Graph + .. workflow:: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.outputs import init_sdc_unwarp_report_wf + wf = init_sdc_unwarp_report_wf() + + Parameters + ---------- + name : str, optional + Workflow name (default: ``sdc_unwarp_report_wf``) + forcedsyn : bool, optional + Whether SyN-SDC was forced. + + Inputs + ------ + in_pre + Reference image, before unwarping + in_post + Reference image, after unwarping + in_seg + Segmentation of preprocessed structural image, including + gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF) + in_xfm + Affine transform from T1 space to BOLD space (ITK format) + + """ + from ...niworkflows.interfaces import SimpleBeforeAfter + from ...niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms + from ...niworkflows.utils.images import dseg_label as _dseg_label + + DEFAULT_MEMORY_MIN_GB = 0.01 + + workflow = pe.Workflow(name=name) + + inputnode = pe.Node(niu.IdentityInterface( + fields=['in_pre', 'in_post', 'in_seg', 'in_xfm']), name='inputnode') + + map_seg = pe.Node(ApplyTransforms( + dimension=3, float=True, interpolation='MultiLabel'), + name='map_seg', mem_gb=0.3) + + sel_wm = pe.Node(niu.Function(function=_dseg_label), name='sel_wm', + mem_gb=DEFAULT_MEMORY_MIN_GB) + sel_wm.inputs.label = 2 + + bold_rpt = pe.Node(SimpleBeforeAfter(), name='bold_rpt', + mem_gb=0.1) + ds_report_sdc = pe.Node( + DerivativesDataSink(desc=('sdc', 'forcedsyn')[forcedsyn], suffix='bold', + datatype='figures'), name='ds_report_sdc', + mem_gb=DEFAULT_MEMORY_MIN_GB, run_without_submitting=True + ) + + workflow.connect([ + (inputnode, bold_rpt, [('in_post', 'after'), + ('in_pre', 'before')]), + (bold_rpt, ds_report_sdc, [('out_report', 'in_file')]), + (inputnode, map_seg, [('in_post', 'reference_image'), + ('in_seg', 'input_image'), + ('in_xfm', 'transforms')]), + (map_seg, sel_wm, [('output_image', 'in_seg')]), + (sel_wm, bold_rpt, [('out', 'wm_seg')]), + ]) + + return workflow diff --git a/aslprep/sdcflows/workflows/pepolar.py b/aslprep/sdcflows/workflows/pepolar.py new file mode 100644 index 000000000..48678cde2 --- /dev/null +++ b/aslprep/sdcflows/workflows/pepolar.py @@ -0,0 +1,387 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Datasets with multiple phase encoded directions. + +.. _sdc_pepolar : + +:abbr:`PEPOLAR (Phase Encoding POLARity)` techniques +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +This corresponds to `this section of the BIDS specification +`__. + +""" + +import pkg_resources as pkgr + +from ...niworkflows.engine.workflows import LiterateWorkflow as Workflow +from ...niworkflows.interfaces import CopyHeader +from ...niworkflows.interfaces.freesurfer import StructuralReference +from ...niworkflows.interfaces.registration import ANTSApplyTransformsRPT +from ...niworkflows.func.util import init_enhance_and_skullstrip_bold_wf + +from nipype.pipeline import engine as pe +from nipype.interfaces import afni, ants, utility as niu + + +def init_pepolar_unwarp_wf(omp_nthreads=1, matched_pe=False, + name="pepolar_unwarp_wf"): + """ + Create the PE-Polar field estimation workflow. + + This workflow takes in a set of EPI files with opposite phase encoding + direction than the target file and calculates a displacements field + (in other words, an ANTs-compatible warp file). + + This procedure works if there is only one ``_epi`` file is present + (as long as it has the opposite phase encoding direction to the target + file). The target file will be used to estimate the field distortion. + However, if there is another ``_epi`` file present with a matching + phase encoding direction to the target it will be used instead. + + Currently, different phase encoding directions in the target file and the + ``_epi`` file(s) (for example, ``i`` and ``j``) is not supported. + + The warp field correcting for the distortions is estimated using AFNI's + ``3dQwarp``, with displacement estimation limited to the target file phase + encoding direction. + + It also calculates a new mask for the input dataset that takes into + account the distortions. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.pepolar import init_pepolar_unwarp_wf + wf = init_pepolar_unwarp_wf() + + Parameters + ---------- + matched_pe : bool + Whether the input ``fmaps_epi`` will contain images with matched + PE blips or not. Please use :func:`sdcflows.workflows.pepolar.check_pes` + to determine whether they exist or not. + name : str + Name for this workflow + omp_nthreads : int + Parallelize internal tasks across the number of CPUs given by this option. + + Inputs + ------ + fmaps_epi : list of tuple(pathlike, str) + The list of EPI images that will be used in PE-Polar correction, and + their corresponding ``PhaseEncodingDirection`` metadata. + The workflow will use the ``epi_pe_dir`` input to separate out those + EPI acquisitions with opposed PE blips and those with matched PE blips + (the latter could be none, and ``in_reference_brain`` would then be + used). The workflow raises a ``ValueError`` when no images with + opposed PE blips are found. + epi_pe_dir : str + The baseline PE direction. + in_reference : pathlike + The baseline reference image (must correspond to ``epi_pe_dir``). + in_reference_brain : pathlike + The reference image above, but skullstripped. + in_mask : pathlike + Not used, present only for consistency across fieldmap estimation + workflows. + + Outputs + ------- + out_reference : pathlike + The ``in_reference`` after unwarping + out_reference_brain : pathlike + The ``in_reference`` after unwarping and skullstripping + out_warp : pathlike + The corresponding :abbr:`DFM (displacements field map)` compatible with + ANTs. + out_mask : pathlike + Mask of the unwarped input file + + """ + workflow = Workflow(name=name) + workflow.__desc__ = """\ +A B0-nonuniformity map (or *fieldmap*) was estimated based on two (or more) +echo-planar imaging (EPI) references with opposing phase-encoding +directions, with `3dQwarp` @afni (AFNI {afni_ver}). +""".format(afni_ver=''.join(['%02d' % v for v in afni.Info().version() or []])) + + inputnode = pe.Node(niu.IdentityInterface( + fields=['fmaps_epi', 'in_reference', 'in_reference_brain', + 'in_mask', 'epi_pe_dir']), name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface( + fields=['out_reference', 'out_reference_brain', 'out_warp', 'out_mask']), + name='outputnode') + + prepare_epi_wf = init_prepare_epi_wf(omp_nthreads=omp_nthreads, + matched_pe=matched_pe, + name="prepare_epi_wf") + + qwarp = pe.Node(afni.QwarpPlusMinus( + pblur=[0.05, 0.05], blur=[-1, -1], noweight=True, minpatch=9, nopadWARP=True, + environ={'OMP_NUM_THREADS': '%d' % omp_nthreads}), + name='qwarp', n_procs=omp_nthreads) + + to_ants = pe.Node(niu.Function(function=_fix_hdr), name='to_ants', + mem_gb=0.01) + + cphdr_warp = pe.Node(CopyHeader(), name='cphdr_warp', mem_gb=0.01) + + unwarp_reference = pe.Node(ANTSApplyTransformsRPT(dimension=3, + generate_report=False, + float=True, + interpolation='LanczosWindowedSinc'), + name='unwarp_reference') + + enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf( + omp_nthreads=omp_nthreads) + + workflow.connect([ + (inputnode, qwarp, [(('epi_pe_dir', _qwarp_args), 'args')]), + (inputnode, cphdr_warp, [('in_reference', 'hdr_file')]), + (inputnode, prepare_epi_wf, [ + ('fmaps_epi', 'inputnode.maps_pe'), + ('epi_pe_dir', 'inputnode.epi_pe'), + ('in_reference_brain', 'inputnode.ref_brain')]), + (prepare_epi_wf, qwarp, [('outputnode.opposed_pe', 'base_file'), + ('outputnode.matched_pe', 'in_file')]), + (qwarp, cphdr_warp, [('source_warp', 'in_file')]), + (cphdr_warp, to_ants, [('out_file', 'in_file')]), + (to_ants, unwarp_reference, [('out', 'transforms')]), + (inputnode, unwarp_reference, [('in_reference', 'reference_image'), + ('in_reference', 'input_image')]), + (unwarp_reference, enhance_and_skullstrip_bold_wf, [ + ('output_image', 'inputnode.in_file')]), + (unwarp_reference, outputnode, [('output_image', 'out_reference')]), + (enhance_and_skullstrip_bold_wf, outputnode, [ + ('outputnode.mask_file', 'out_mask'), + ('outputnode.skull_stripped_file', 'out_reference_brain')]), + (to_ants, outputnode, [('out', 'out_warp')]), + ]) + + return workflow + + +def init_prepare_epi_wf(omp_nthreads, matched_pe=False, + name="prepare_epi_wf"): + """ + Prepare opposed-PE EPI images for PE-POLAR SDC. + + This workflow takes in a set of EPI files and returns two 3D volumes with + matching and opposed PE directions, ready to be used in field distortion + estimation. + + The procedure involves: estimating a robust template using FreeSurfer's + ``mri_robust_template``, bias field correction using ANTs ``N4BiasFieldCorrection`` + and AFNI ``3dUnifize``, skullstripping using FSL BET and AFNI ``3dAutomask``, + and rigid coregistration to the reference using ANTs. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.pepolar import init_prepare_epi_wf + wf = init_prepare_epi_wf(omp_nthreads=8) + + Parameters + ---------- + matched_pe : bool + Whether the input ``fmaps_epi`` will contain images with matched + PE blips or not. Please use :func:`sdcflows.workflows.pepolar.check_pes` + to determine whether they exist or not. + name : str + Name for this workflow + omp_nthreads : int + Parallelize internal tasks across the number of CPUs given by this option. + + Inputs + ------ + epi_pe : str + Phase-encoding direction of the EPI image to be corrected. + maps_pe : list of tuple(pathlike, str) + list of 3D or 4D NIfTI images + ref_brain + coregistration reference (skullstripped and bias field corrected) + + Outputs + ------- + opposed_pe : pathlike + single 3D NIfTI file + matched_pe : pathlike + single 3D NIfTI file + + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['epi_pe', 'maps_pe', 'ref_brain']), + name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface(fields=['opposed_pe', 'matched_pe']), + name='outputnode') + + ants_settings = pkgr.resource_filename('sdcflows', + 'data/translation_rigid.json') + + split = pe.Node(niu.Function(function=_split_epi_lists), name='split') + + merge_op = pe.Node( + StructuralReference(auto_detect_sensitivity=True, + initial_timepoint=1, + fixed_timepoint=True, # Align to first image + intensity_scaling=True, + # 7-DOF (rigid + intensity) + no_iteration=True, + subsample_threshold=200, + out_file='template.nii.gz'), + name='merge_op') + + ref_op_wf = init_enhance_and_skullstrip_bold_wf( + omp_nthreads=omp_nthreads, name='ref_op_wf') + + op2ref_reg = pe.Node(ants.Registration( + from_file=ants_settings, output_warped_image=True), + name='op2ref_reg', n_procs=omp_nthreads) + + workflow = Workflow(name=name) + workflow.connect([ + (inputnode, split, [('maps_pe', 'in_files'), + ('epi_pe', 'pe_dir')]), + (split, merge_op, [(('out', _front), 'in_files')]), + (merge_op, ref_op_wf, [('out_file', 'inputnode.in_file')]), + (ref_op_wf, op2ref_reg, [ + ('outputnode.skull_stripped_file', 'moving_image')]), + (inputnode, op2ref_reg, [('ref_brain', 'fixed_image')]), + (op2ref_reg, outputnode, [('warped_image', 'opposed_pe')]), + ]) + + if not matched_pe: + workflow.connect([ + (inputnode, outputnode, [('ref_brain', 'matched_pe')]), + ]) + return workflow + + merge_ma = pe.Node( + StructuralReference(auto_detect_sensitivity=True, + initial_timepoint=1, + fixed_timepoint=True, # Align to first image + intensity_scaling=True, + # 7-DOF (rigid + intensity) + no_iteration=True, + subsample_threshold=200, + out_file='template.nii.gz'), + name='merge_ma') + + ref_ma_wf = init_enhance_and_skullstrip_bold_wf( + omp_nthreads=omp_nthreads, name='ref_ma_wf') + + ma2ref_reg = pe.Node(ants.Registration( + from_file=ants_settings, output_warped_image=True), + name='ma2ref_reg', n_procs=omp_nthreads) + + workflow.connect([ + (split, merge_ma, [(('out', _last), 'in_files')]), + (merge_ma, ref_ma_wf, [('out_file', 'inputnode.in_file')]), + (ref_ma_wf, ma2ref_reg, [ + ('outputnode.skull_stripped_file', 'moving_image')]), + (inputnode, ma2ref_reg, [('ref_brain', 'fixed_image')]), + (ma2ref_reg, outputnode, [('warped_image', 'matched_pe')]), + ]) + return workflow + + +def _front(inlist): + if isinstance(inlist, (list, tuple)): + return inlist[0] + return inlist + + +def _last(inlist): + if isinstance(inlist, (list, tuple)): + return inlist[-1] + return inlist + + +def _fix_hdr(in_file, newpath=None): + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + nii = nb.load(in_file) + hdr = nii.header.copy() + hdr.set_data_dtype('`__. + + +""" + +from nipype.interfaces import fsl, utility as niu +from nipype.pipeline import engine as pe +from ...niworkflows.engine.workflows import LiterateWorkflow as Workflow + +from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads, SubtractPhases +from .gre import init_fmap_postproc_wf, init_magnitude_wf + + +def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): + r""" + Distortion correction of EPI sequences using phase-difference maps. + + Estimates the fieldmap using a phase-difference image and one or more + magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)` + acquisitions. + The most delicate bit of this workflow is the phase-unwrapping process: phase maps + are clipped in the range :math:`[0 \dotsb 2\pi )`. + To find the integer number of offsets that make a region continously smooth with + its neighbour, FSL PRELUDE is run [Jenkinson2003]_. + FSL PRELUDE takes wrapped maps in the range 0 to 6.28, `as per the user guide + `__. + For the phase-difference maps, recentering back to :math:`[-\pi \dotsb \pi )` is necessary. + After some massaging and the application of the effective echo spacing parameter, + the phase-difference maps can be converted into a *B0 field map* in Hz units. + The `original code was taken from nipype + `_. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.phdiff import init_phdiff_wf + wf = init_phdiff_wf(omp_nthreads=1) + + Parameters + ---------- + omp_nthreads : int + Maximum number of threads an individual process may use + + Inputs + ------ + magnitude : list of os.pathlike + List of path(s) the GRE magnitude maps. + phasediff : list of tuple(os.pathlike, dict) + List containing one GRE phase-difference map with its corresponding metadata + (requires ``EchoTime1`` and ``EchoTime2``), or the phase maps for the two + subsequent echoes, with their metadata (requires ``EchoTime``). + + Outputs + ------- + fmap_ref : pathlike + The average magnitude image, skull-stripped + fmap_mask : pathlike + The brain mask applied to the fieldmap + fmap : pathlike + The estimated fieldmap in Hz + + References + ---------- + .. [Jenkinson2003] Jenkinson, M. (2003) Fast, automated, N-dimensional phase-unwrapping + algorithm. MRM 49(1):193-197. doi:`10.1002/mrm.10354 <10.1002/mrm.10354>`__. + + """ + workflow = Workflow(name=name) + workflow.__desc__ = """\ +A B0-nonuniformity map (or *fieldmap*) was estimated based on a phase-difference map +calculated with a dual-echo GRE (gradient-recall echo) sequence, processed with a +custom workflow of *SDCFlows* inspired by the +[`epidewarp.fsl` script](http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) +and further improvements in HCP Pipelines [@hcppipelines]. +""" + + inputnode = pe.Node(niu.IdentityInterface(fields=['magnitude', 'phasediff']), + name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface( + fields=['fmap', 'fmap_ref', 'fmap_mask']), name='outputnode') + + split = pe.MapNode(niu.Function(function=_split, output_names=['map_file', 'meta']), + iterfield=['phasediff'], run_without_submitting=True, name='split') + + magnitude_wf = init_magnitude_wf(omp_nthreads=omp_nthreads) + + # phase diff -> radians + phmap2rads = pe.MapNode(PhaseMap2rads(), name='phmap2rads', + iterfield=['in_file'], run_without_submitting=True) + # FSL PRELUDE will perform phase-unwrapping + prelude = pe.Node(fsl.PRELUDE(), name='prelude') + + calc_phdiff = pe.Node(SubtractPhases(), name='calc_phdiff', + run_without_submitting=True) + + fmap_postproc_wf = init_fmap_postproc_wf(omp_nthreads=omp_nthreads, + fmap_bspline=False) + compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap') + + workflow.connect([ + (inputnode, split, [('phasediff', 'phasediff')]), + (inputnode, magnitude_wf, [('magnitude', 'inputnode.magnitude')]), + (magnitude_wf, prelude, [('outputnode.fmap_ref', 'magnitude_file'), + ('outputnode.fmap_mask', 'mask_file')]), + (split, phmap2rads, [('map_file', 'in_file')]), + (phmap2rads, calc_phdiff, [('out_file', 'in_phases')]), + (split, calc_phdiff, [('meta', 'in_meta')]), + (calc_phdiff, prelude, [('phase_diff', 'phase_file')]), + (prelude, fmap_postproc_wf, [('unwrapped_phase_file', 'inputnode.fmap')]), + (calc_phdiff, fmap_postproc_wf, [('metadata', 'inputnode.metadata')]), + (magnitude_wf, fmap_postproc_wf, [ + ('outputnode.fmap_mask', 'inputnode.fmap_mask'), + ('outputnode.fmap_ref', 'inputnode.fmap_ref')]), + (fmap_postproc_wf, compfmap, [('outputnode.out_fmap', 'in_file'), + ('outputnode.metadata', 'metadata')]), + (compfmap, outputnode, [('out_file', 'fmap')]), + (magnitude_wf, outputnode, [('outputnode.fmap_ref', 'fmap_ref'), + ('outputnode.fmap_mask', 'fmap_mask')]), + ]) + + return workflow + + +def _split(phasediff): + return phasediff diff --git a/aslprep/sdcflows/workflows/syn.py b/aslprep/sdcflows/workflows/syn.py new file mode 100644 index 000000000..f16c7eae5 --- /dev/null +++ b/aslprep/sdcflows/workflows/syn.py @@ -0,0 +1,216 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Estimating the susceptibility distortions without fieldmaps. + +.. _sdc_fieldmapless : + +Fieldmap-less estimation (experimental) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +In the absence of direct measurements of fieldmap data, we provide an (experimental) +option to estimate the susceptibility distortion based on the ANTs symmetric +normalization (SyN) technique. +This feature may be enabled, using the ``--use-syn-sdc`` flag, and will only be +applied if fieldmaps are unavailable. + +During the evaluation phase, the ``--force-syn`` flag will cause this estimation to +be performed *in addition to* fieldmap-based estimation, to permit the direct +comparison of the results of each technique. +Note that, even if ``--force-syn`` is given, the functional outputs of FMRIPREP will +be corrected using the fieldmap-based estimates. + +Feedback will be enthusiastically received. + + +""" +from pkg_resources import resource_filename + +from nipype import logging +from nipype.pipeline import engine as pe +from nipype.interfaces import fsl, utility as niu +from nipype.interfaces.image import Rescale + +from niworkflows.engine.workflows import LiterateWorkflow as Workflow +from niworkflows.interfaces.fixes import (FixHeaderApplyTransforms as ApplyTransforms, + FixHeaderRegistration as Registration) +from niworkflows.func.util import init_skullstrip_bold_wf + +DEFAULT_MEMORY_MIN_GB = 0.01 +LOGGER = logging.getLogger('nipype.workflow') + + +def init_syn_sdc_wf(omp_nthreads, epi_pe=None, + atlas_threshold=3, name='syn_sdc_wf'): + """ + Build the *fieldmap-less* susceptibility-distortion estimation workflow. + + This workflow takes a skull-stripped T1w image and reference BOLD image and + estimates a susceptibility distortion correction warp, using ANTs symmetric + normalization (SyN) and the average fieldmap atlas described in + [Treiber2016]_. + + SyN deformation is restricted to the phase-encoding (PE) direction. + If no PE direction is specified, anterior-posterior PE is assumed. + + SyN deformation is also restricted to regions that are expected to have a + >3mm (approximately 1 voxel) warp, based on the fieldmap atlas. + + This technique is a variation on those developed in [Huntenburg2014]_ and + [Wang2017]_. + + Workflow Graph + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.syn import init_syn_sdc_wf + wf = init_syn_sdc_wf( + epi_pe='j', + omp_nthreads=8) + + Inputs + ------ + in_reference + reference image + in_reference_brain + skull-stripped reference image + t1w_brain + skull-stripped, bias-corrected structural image + std2anat_xfm + inverse registration transform of T1w image to MNI template + + Outputs + ------- + out_reference + the ``in_reference`` image after unwarping + out_reference_brain + the ``in_reference_brain`` image after unwarping + out_warp + the corresponding :abbr:`DFM (displacements field map)` compatible with + ANTs + out_mask + mask of the unwarped input file + + References + ---------- + .. [Treiber2016] Treiber, J. M. et al. (2016) Characterization and Correction + of Geometric Distortions in 814 Diffusion Weighted Images, + PLoS ONE 11(3): e0152472. doi:`10.1371/journal.pone.0152472 + `_. + .. [Wang2017] Wang S, et al. (2017) Evaluation of Field Map and Nonlinear + Registration Methods for Correction of Susceptibility Artifacts + in Diffusion MRI. Front. Neuroinform. 11:17. + doi:`10.3389/fninf.2017.00017 + `_. + .. [Huntenburg2014] Huntenburg, J. M. (2014) Evaluating Nonlinear + Coregistration of BOLD EPI and T1w Images. Berlin: Master + Thesis, Freie Universität. `PDF + `_. + + """ + if epi_pe is None or epi_pe[0] not in ['i', 'j']: + LOGGER.warning('Incorrect phase-encoding direction, assuming PA (posterior-to-anterior).') + epi_pe = 'j' + + workflow = Workflow(name=name) + workflow.__desc__ = """\ +A deformation field to correct for susceptibility distortions was estimated +based on *fMRIPrep*'s *fieldmap-less* approach. +The deformation field is that resulting from co-registering the BOLD reference +to the same-subject T1w-reference with its intensity inverted [@fieldmapless1; +@fieldmapless2]. +Registration is performed with `antsRegistration` (ANTs {ants_ver}), and +the process regularized by constraining deformation to be nonzero only +along the phase-encoding direction, and modulated with an average fieldmap +template [@fieldmapless3]. +""".format(ants_ver=Registration().version or '') + inputnode = pe.Node( + niu.IdentityInterface(['in_reference', 'in_reference_brain', + 't1w_brain', 'std2anat_xfm']), + name='inputnode') + outputnode = pe.Node( + niu.IdentityInterface(['out_reference', 'out_reference_brain', + 'out_mask', 'out_warp']), + name='outputnode') + + # Collect predefined data + # Atlas image and registration affine + atlas_img = resource_filename('aslprep', 'sdcflows/data/fmap_atlas.nii.gz') + # Registration specifications + affine_transform = resource_filename('aslprep', 'sdcflows/data/affine.json') + syn_transform = resource_filename('aslprep', 'sdcflows/data/susceptibility_syn.json') + + invert_t1w = pe.Node(Rescale(invert=True), name='invert_t1w', + mem_gb=0.3) + + ref_2_t1 = pe.Node(Registration(from_file=affine_transform), + name='ref_2_t1', n_procs=omp_nthreads) + t1_2_ref = pe.Node(ApplyTransforms(invert_transform_flags=[True]), + name='t1_2_ref', n_procs=omp_nthreads) + + # 1) BOLD -> T1; 2) MNI -> T1; 3) ATLAS -> MNI + transform_list = pe.Node(niu.Merge(3), name='transform_list', + mem_gb=DEFAULT_MEMORY_MIN_GB) + transform_list.inputs.in3 = resource_filename( + 'aslprep', 'sdcflows/data/fmap_atlas_2_MNI152NLin2009cAsym_affine.mat') + + # Inverting (1), then applying in reverse order: + # + # ATLAS -> MNI -> T1 -> BOLD + atlas_2_ref = pe.Node( + ApplyTransforms(invert_transform_flags=[True, False, False]), + name='atlas_2_ref', n_procs=omp_nthreads, + mem_gb=0.3) + atlas_2_ref.inputs.input_image = atlas_img + + threshold_atlas = pe.Node( + fsl.maths.MathsCommand(args='-thr {:.8g} -bin'.format(atlas_threshold), + output_datatype='char'), + name='threshold_atlas', mem_gb=0.3) + + fixed_image_masks = pe.Node(niu.Merge(2), name='fixed_image_masks', + mem_gb=DEFAULT_MEMORY_MIN_GB) + fixed_image_masks.inputs.in1 = 'NULL' + + restrict = [[int(epi_pe[0] == 'i'), int(epi_pe[0] == 'j'), 0]] * 2 + syn = pe.Node( + Registration(from_file=syn_transform, restrict_deformation=restrict), + name='syn', n_procs=omp_nthreads) + + unwarp_ref = pe.Node(ApplyTransforms( + dimension=3, float=True, interpolation='LanczosWindowedSinc'), + name='unwarp_ref') + + skullstrip_bold_wf = init_skullstrip_bold_wf() + + workflow.connect([ + (inputnode, invert_t1w, [('t1w_brain', 'in_file'), + ('in_reference', 'ref_file')]), + (inputnode, ref_2_t1, [('in_reference_brain', 'moving_image')]), + (invert_t1w, ref_2_t1, [('out_file', 'fixed_image')]), + (inputnode, t1_2_ref, [('in_reference', 'reference_image')]), + (invert_t1w, t1_2_ref, [('out_file', 'input_image')]), + (ref_2_t1, t1_2_ref, [('forward_transforms', 'transforms')]), + (ref_2_t1, transform_list, [('forward_transforms', 'in1')]), + (inputnode, transform_list, [ + ('std2anat_xfm', 'in2')]), + (inputnode, atlas_2_ref, [('in_reference', 'reference_image')]), + (transform_list, atlas_2_ref, [('out', 'transforms')]), + (atlas_2_ref, threshold_atlas, [('output_image', 'in_file')]), + (threshold_atlas, fixed_image_masks, [('out_file', 'in2')]), + (inputnode, syn, [('in_reference_brain', 'moving_image')]), + (t1_2_ref, syn, [('output_image', 'fixed_image')]), + (fixed_image_masks, syn, [('out', 'fixed_image_masks')]), + (syn, outputnode, [('forward_transforms', 'out_warp')]), + (syn, unwarp_ref, [('forward_transforms', 'transforms')]), + (inputnode, unwarp_ref, [('in_reference', 'reference_image'), + ('in_reference', 'input_image')]), + (unwarp_ref, skullstrip_bold_wf, [ + ('output_image', 'inputnode.in_file')]), + (unwarp_ref, outputnode, [('output_image', 'out_reference')]), + (skullstrip_bold_wf, outputnode, [ + ('outputnode.skull_stripped_file', 'out_reference_brain'), + ('outputnode.mask_file', 'out_mask')]), + ]) + + return workflow diff --git a/aslprep/sdcflows/workflows/tests/__init__.py b/aslprep/sdcflows/workflows/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/aslprep/sdcflows/workflows/tests/test_base.py b/aslprep/sdcflows/workflows/tests/test_base.py new file mode 100644 index 000000000..8fb347693 --- /dev/null +++ b/aslprep/sdcflows/workflows/tests/test_base.py @@ -0,0 +1,137 @@ +"""Test base workflows.""" +import pytest + +from ..base import init_sdc_estimate_wf + +EPI_METADATA = { + 'MultibandAccelerationFactor': 8, + 'PhaseEncodingDirection': 'i', + 'RepetitionTime': 0.72, + 'TaskName': 'Resting-state fMRI', +} +EPI_FMAP_METADATA_1 = { + 'BandwidthPerPixelPhaseEncode': 2290, + 'EPIFactor': 90, + 'EffectiveEchoSpacing': 0.00058, + 'IntendedFor': ['func/sub-HCP101006_task-rest_dir-LR_bold.nii.gz', + 'func/sub-HCP101006_task-rest_dir-LR_sbref.nii.gz'], + 'MultibandAccelerationFactor': 1, + 'PhaseEncodingDirection': 'i-', +} +EPI_FMAP_METADATA_2 = EPI_FMAP_METADATA_1.copy() +EPI_FMAP_METADATA_2['PhaseEncodingDirection'] = 'i' + +PHDIFF_METADATA = { + 'EchoTime1': 0.00492, + 'EchoTime2': 0.00738, +} +PHASE1_METADATA = { + 'EchoTime': 0.00492, +} +PHASE2_METADATA = { + 'EchoTime': 0.00738, +} + + +FMAP_DICT_ELEMENTS = { + 'epi1': [('sub-HCP101006/fmap/sub-HCP101006_dir-RL_epi.nii.gz', + EPI_FMAP_METADATA_1.copy())], + 'epi2': [('sub-HCP101006/fmap/sub-HCP101006_dir-RL_epi.nii.gz', + EPI_FMAP_METADATA_1.copy()), + ('sub-HCP101006/fmap/sub-HCP101006_dir-LR_epi.nii.gz', + EPI_FMAP_METADATA_2.copy())], + 'phdiff1': { + 'magnitude': [('sub-HCP101006/fmap/sub-HCP101006_magnitude1.nii.gz', {}), + ('sub-HCP101006/fmap/sub-HCP101006_magnitude2.nii.gz', {})], + 'phases': [('sub-HCP101006/fmap/sub-HCP101006_phasediff.nii.gz', PHDIFF_METADATA)], + }, + 'phdiff2': { + 'magnitude': [('sub-HCP101006/fmap/sub-HCP101006_magnitude1.nii.gz', {}), + ('sub-HCP101006/fmap/sub-HCP101006_magnitude2.nii.gz', {})], + 'phases': [('sub-HCP101006/fmap/sub-HCP101006_phase1.nii.gz', + PHASE1_METADATA.copy()), + ('sub-HCP101006/fmap/sub-HCP101006_phase2.nii.gz', + PHASE2_METADATA.copy())], + }, + 'fmap1': { + 'magnitude': [('sub-HCP101006/fmap/sub-HCP101006_magnitude.nii.gz', {})], + 'fieldmap': [('sub-HCP101006/fmap/sub-HCP101006_fieldmap.nii.gz', {})], + }, + 'syn': { + 't1w': [('sub-HCP101006/fmap/sub-HCP101006_T1w.nii.gz', {})] + } +} + + +@pytest.mark.parametrize('method', ['skip', 'phasediff', 'pepolar', 'fieldmap', 'syn']) +def test_base(method): + """Check the heuristics are correctly applied.""" + fieldmaps = { + 'epi': FMAP_DICT_ELEMENTS['epi1'].copy(), + 'fieldmap': [FMAP_DICT_ELEMENTS['fmap1'].copy()], + 'phasediff': [FMAP_DICT_ELEMENTS['phdiff1'].copy()], + } + epi_meta = EPI_METADATA.copy() + + if method == 'skip': + wf = init_sdc_estimate_wf(fmaps=None, epi_meta=epi_meta) + assert wf.inputs.outputnode.method == 'None' + + with pytest.raises(ValueError): + wf = init_sdc_estimate_wf(fmaps={'unsupported': None}, epi_meta=epi_meta) + elif method == 'pepolar': + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + assert 'PEPOLAR' in wf.inputs.outputnode.method + + fieldmaps['epi'][0][1].pop('PhaseEncodingDirection') + with pytest.raises(ValueError): + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + + fieldmaps['epi'][0][1]['PhaseEncodingDirection'] = \ + EPI_FMAP_METADATA_1['PhaseEncodingDirection'] + epi_meta.pop('PhaseEncodingDirection') + with pytest.raises(ValueError): + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + epi_meta = EPI_METADATA.copy() + + fieldmaps['epi'] = FMAP_DICT_ELEMENTS['epi2'] + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + assert 'PEPOLAR' in wf.inputs.outputnode.method + + elif method == 'fieldmap': + fieldmaps = { + 'fieldmap': [FMAP_DICT_ELEMENTS['fmap1'].copy()], + 'phasediff': [FMAP_DICT_ELEMENTS['phdiff1'].copy()], + } + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + assert 'directly measured B0 map' in wf.inputs.outputnode.method + + elif method == 'phasediff': + fieldmaps = { + 'phasediff': [FMAP_DICT_ELEMENTS['phdiff1'].copy()], + } + + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + assert 'phase-difference' in wf.inputs.outputnode.method + + epi_meta.pop('PhaseEncodingDirection') + with pytest.raises(ValueError): + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + + elif method == 'syn': + fmaps_onlysyn = {'syn': FMAP_DICT_ELEMENTS['syn']} + wf = init_sdc_estimate_wf(fmaps=fmaps_onlysyn, epi_meta=epi_meta) + assert 'SyN' in wf.inputs.outputnode.method + + epi_pe = epi_meta.pop('PhaseEncodingDirection') + wf = init_sdc_estimate_wf(fmaps=fmaps_onlysyn, epi_meta=epi_meta) + assert 'SyN' in wf.inputs.outputnode.method + + # Emulate --force-syn + fieldmaps.update(fmaps_onlysyn) + with pytest.raises(ValueError): + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + + epi_meta['PhaseEncodingDirection'] = epi_pe + wf = init_sdc_estimate_wf(fmaps=fieldmaps, epi_meta=epi_meta) + assert 'PEPOLAR' in wf.inputs.outputnode.method diff --git a/aslprep/sdcflows/workflows/tests/test_pepolar.py b/aslprep/sdcflows/workflows/tests/test_pepolar.py new file mode 100644 index 000000000..f2d449104 --- /dev/null +++ b/aslprep/sdcflows/workflows/tests/test_pepolar.py @@ -0,0 +1,194 @@ +"""Test pepolar type of fieldmaps.""" +from os import cpu_count +import pytest +from niworkflows.interfaces.bids import DerivativesDataSink +from nipype.pipeline import engine as pe + +from ..pepolar import ( + check_pes, _split_epi_lists, init_prepare_epi_wf, init_pepolar_unwarp_wf +) + + +def test_split_epi_lists(bids_layouts, tmpdir): + """Test preparation workflow.""" + tmpdir.chdir() + + layout = bids_layouts['testdata'] + + bold = layout.get(suffix='bold', direction='LR', + extension=['.nii.gz', '.nii'])[0] + epidata = layout.get(suffix='epi', desc=None, extension=['.nii.gz', '.nii']) + + # EPI fmaps in HCP101006 have 3 volumes each + a, b = _split_epi_lists( + in_files=[(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata], + pe_dir=bold.get_metadata()['PhaseEncodingDirection'] + ) + + assert len(a) == 3 + assert len(b) == 3 + + # If we append the BOLD aligned (not that you would do this), then the + # second array should have 53 volumes (50 from _bold, 3 from _epi) + a, b = _split_epi_lists( + in_files=[(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata + [bold]], + pe_dir=bold.get_metadata()['PhaseEncodingDirection'] + ) + + assert len(a) == 3 + assert len(b) == 53 + + +def test_prepare_epi_wf0(bids_layouts, tmpdir): + """Test preparation workflow.""" + tmpdir.chdir() + + layout = bids_layouts['testdata'] + + bold = layout.get(suffix='bold', direction='LR', + extension=['.nii.gz', '.nii'])[0] + epidata = layout.get(suffix='epi', direction='LR', + desc=None, extension=['.nii.gz', '.nii']) + + with pytest.raises(ValueError): + check_pes([ + (im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata], bold.get_metadata()['PhaseEncodingDirection']) + + +def test_prepare_epi_wf1(bids_layouts, tmpdir): + """Test preparation workflow.""" + tmpdir.chdir() + + layout = bids_layouts['testdata'] + + bold = layout.get(suffix='bold', direction='LR', + extension=['.nii.gz', '.nii'])[0] + boldref = layout.get(suffix='boldref', direction='LR', + extension=['.nii.gz', '.nii'])[0] + epidata = layout.get(suffix='epi', desc=None, extension=['.nii.gz', '.nii']) + + matched_pe = check_pes([(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata], bold.get_metadata()['PhaseEncodingDirection']) + + assert matched_pe is True + + wf = init_prepare_epi_wf(omp_nthreads=1, matched_pe=matched_pe) + wf.inputs.inputnode.maps_pe = [(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata] + wf.inputs.inputnode.ref_brain = boldref.path + wf.inputs.inputnode.epi_pe = bold.get_metadata()['PhaseEncodingDirection'] + + +def test_prepare_epi_wf2(bids_layouts, tmpdir): + """Test preparation workflow.""" + tmpdir.chdir() + + layout = bids_layouts['testdata'] + + bold = layout.get(suffix='bold', direction='LR', + extension=['.nii.gz', '.nii'])[0] + boldref = layout.get(suffix='boldref', direction='LR', + extension=['.nii.gz', '.nii'])[0] + epidata = layout.get(suffix='epi', direction='RL', + desc=None, extension=['.nii.gz', '.nii']) + + matched_pe = check_pes([(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata], bold.get_metadata()['PhaseEncodingDirection']) + + assert matched_pe is False + + wf = init_prepare_epi_wf(omp_nthreads=1, matched_pe=matched_pe) + wf.inputs.inputnode.maps_pe = [(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata] + wf.inputs.inputnode.ref_brain = boldref.path + wf.inputs.inputnode.epi_pe = bold.get_metadata()['PhaseEncodingDirection'] + + +@pytest.mark.parametrize('dataset', [ + 'ds001600', + 'testdata', +]) +def test_pepolar_wf1(bids_layouts, output_path, dataset, workdir): + """Test preparation workflow.""" + layout = bids_layouts[dataset] + + if dataset == 'testdata': + bold = layout.get(suffix='bold', direction='LR', + extension=['.nii.gz', '.nii'])[0] + boldref = layout.get(suffix='boldref', direction='LR', desc='brain', + extension=['.nii.gz', '.nii'])[0] + epidata = layout.get(suffix='epi', desc=None, extension=['.nii.gz', '.nii']) + elif dataset == 'ds001600': + bold = layout.get(suffix='bold', acquisition='AP', + extension=['.nii.gz', '.nii'])[0] + epidata = layout.get(suffix='epi', extension=['.nii.gz', '.nii']) + + matched_pe = check_pes([(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata], bold.get_metadata()['PhaseEncodingDirection']) + + wf = init_pepolar_unwarp_wf(omp_nthreads=cpu_count(), matched_pe=matched_pe) + wf.inputs.inputnode.fmaps_epi = [(im.path, im.get_metadata()['PhaseEncodingDirection']) + for im in epidata] + wf.inputs.inputnode.epi_pe_dir = bold.get_metadata()['PhaseEncodingDirection'] + + if output_path: + from nipype.interfaces import utility as niu + from ..pepolar import Workflow + from ...interfaces.reportlets import FieldmapReportlet + + boiler = Workflow(name='pepolar_%s' % dataset) + + split_field = pe.Node(niu.Function(function=_split_field), name='split_field') + + if dataset == 'ds001600': + from niworkflows.func.util import init_bold_reference_wf + gen_ref = init_bold_reference_wf( + omp_nthreads=cpu_count(), + bold_file=bold.path) + boiler.connect([ + (gen_ref, wf, [ + ('outputnode.ref_image', 'inputnode.in_reference'), + ('outputnode.ref_image_brain', 'inputnode.in_reference_brain')]) + ]) + else: + wf.inputs.inputnode.in_reference_brain = boldref.path + wf.inputs.inputnode.in_reference = boldref.path + + rep = pe.Node(FieldmapReportlet(reference_label='EPI Reference'), 'simple_report') + rep.interface._always_run = True + dsink = pe.Node(DerivativesDataSink( + base_directory=str(output_path), out_path_base='sdcflows', datatype="figures", + suffix='fieldmap', desc='pepolar', dismiss_entities='fmap'), name='dsink') + dsink.inputs.source_file = epidata[0].path + + boiler.connect([ + (wf, split_field, [ + ('inputnode.epi_pe_dir', 'pe_dir'), + ('outputnode.out_warp', 'in_field')]), + (split_field, rep, [ + ('out', 'fieldmap')]), + (wf, rep, [ + # ('outputnode.out_warp', 'fieldmap'), + ('outputnode.out_reference_brain', 'reference'), + ('outputnode.out_mask', 'mask')]), + (rep, dsink, [('out_report', 'in_file')]), + ]) + + if workdir: + boiler.base_dir = str(workdir) + boiler.run(plugin='MultiProc', plugin_args={'n_proc': cpu_count()}) + + +def _split_field(in_field, pe_dir): + from os.path import abspath + import numpy as np + import nibabel as nb + axis = 'ijk'.index(pe_dir[0]) + im = nb.load(in_field) + data = np.squeeze(im.get_fdata())[..., axis] + dirnii = nb.Nifti1Image(data, im.affine, im.header) + dirnii.to_filename('fieldmap.nii.gz') + return abspath('fieldmap.nii.gz') diff --git a/aslprep/sdcflows/workflows/tests/test_phdiff.py b/aslprep/sdcflows/workflows/tests/test_phdiff.py new file mode 100644 index 000000000..ec936fa7c --- /dev/null +++ b/aslprep/sdcflows/workflows/tests/test_phdiff.py @@ -0,0 +1,120 @@ +"""Test phase-difference type of fieldmaps.""" +import pytest +from niworkflows.interfaces.bids import DerivativesDataSink +from nipype.pipeline import engine as pe + +from ..phdiff import init_phdiff_wf, Workflow + + +@pytest.mark.parametrize('dataset', [ + 'ds001600', + 'testdata', +]) +def test_phdiff(bids_layouts, tmpdir, output_path, dataset, workdir): + """Test creation of the workflow.""" + tmpdir.chdir() + + extra_entities = {} + if dataset == 'ds001600': + extra_entities['acquisition'] = 'v4' + + data = bids_layouts[dataset] + wf = Workflow(name='phdiff_%s' % dataset) + phdiff_wf = init_phdiff_wf(omp_nthreads=2) + phdiff_wf.inputs.inputnode.magnitude = data.get( + suffix=['magnitude1', 'magnitude2'], + return_type='file', + extension=['.nii', '.nii.gz'], + **extra_entities) + + phdiff_files = data.get( + suffix='phasediff', + extension=['.nii', '.nii.gz'], + **extra_entities) + + phdiff_wf.inputs.inputnode.phasediff = [ + (ph.path, ph.get_metadata()) for ph in phdiff_files] + + if output_path: + from ...interfaces.reportlets import FieldmapReportlet + rep = pe.Node(FieldmapReportlet(reference_label='Magnitude'), 'simple_report') + rep.interface._always_run = True + + ds_report = pe.Node(DerivativesDataSink( + base_directory=str(output_path), out_path_base='sdcflows', datatype='figures', + suffix='fieldmap', desc='phasediff', dismiss_entities='fmap'), name='ds_report') + ds_report.inputs.source_file = phdiff_files[0].path + + dsink_fmap = pe.Node(DerivativesDataSink( + base_directory=str(output_path), dismiss_entities='fmap', + desc='phasediff', suffix='fieldmap'), name='dsink_fmap') + dsink_fmap.interface.out_path_base = 'sdcflows' + dsink_fmap.inputs.source_file = phdiff_files[0].path + + wf.connect([ + (phdiff_wf, rep, [ + ('outputnode.fmap', 'fieldmap'), + ('outputnode.fmap_ref', 'reference'), + ('outputnode.fmap_mask', 'mask')]), + (rep, ds_report, [('out_report', 'in_file')]), + (phdiff_wf, dsink_fmap, [('outputnode.fmap', 'in_file')]), + ]) + else: + wf.add_nodes([phdiff_wf]) + + if workdir: + wf.base_dir = str(workdir) + + wf.run() + + +def test_phases(bids_layouts, tmpdir, output_path, workdir): + """Test creation of the workflow.""" + tmpdir.chdir() + + data = bids_layouts['ds001600'] + wf = Workflow(name='phases_ds001600') + phdiff_wf = init_phdiff_wf(omp_nthreads=2) + phdiff_wf.inputs.inputnode.magnitude = data.get( + suffix=['magnitude1', 'magnitude2'], + acquisition='v2', + return_type='file', + extension=['.nii', '.nii.gz']) + + phdiff_files = data.get(suffix=['phase1', 'phase2'], acquisition='v2', + extension=['.nii', '.nii.gz']) + + phdiff_wf.inputs.inputnode.phasediff = [ + (ph.path, ph.get_metadata()) for ph in phdiff_files] + + if output_path: + from ...interfaces.reportlets import FieldmapReportlet + rep = pe.Node(FieldmapReportlet(reference_label='Magnitude'), 'simple_report') + rep.interface._always_run = True + + ds_report = pe.Node(DerivativesDataSink( + base_directory=str(output_path), out_path_base='sdcflows', datatype='figures', + suffix='fieldmap', desc='twophases', dismiss_entities='fmap'), name='ds_report') + ds_report.inputs.source_file = phdiff_files[0].path + + dsink_fmap = pe.Node(DerivativesDataSink( + base_directory=str(output_path), suffix='fieldmap', desc='twophases', + dismiss_entities='fmap'), name='dsink_fmap') + dsink_fmap.interface.out_path_base = 'sdcflows' + dsink_fmap.inputs.source_file = phdiff_files[0].path + + wf.connect([ + (phdiff_wf, rep, [ + ('outputnode.fmap', 'fieldmap'), + ('outputnode.fmap_ref', 'reference'), + ('outputnode.fmap_mask', 'mask')]), + (rep, ds_report, [('out_report', 'in_file')]), + (phdiff_wf, dsink_fmap, [('outputnode.fmap', 'in_file')]), + ]) + else: + wf.add_nodes([phdiff_wf]) + + if workdir: + wf.base_dir = str(workdir) + + wf.run() diff --git a/aslprep/sdcflows/workflows/unwarp.py b/aslprep/sdcflows/workflows/unwarp.py new file mode 100644 index 000000000..02d60c2cb --- /dev/null +++ b/aslprep/sdcflows/workflows/unwarp.py @@ -0,0 +1,99 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Apply the estimated fieldmap to perform susceptibility-derived distortion correction.""" +from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu +from niworkflows.engine.workflows import LiterateWorkflow as Workflow +from niworkflows.interfaces.registration import ANTSApplyTransformsRPT +from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf + + +def init_sdc_unwarp_wf(omp_nthreads, debug, name='sdc_unwarp_wf'): + """ + Apply the warping given by a displacements fieldmap. + + This workflow takes in a displacements field through which the + input reference can be corrected for susceptibility-derived distortion. + + It also calculates a new mask for the input dataset, after the distortions + have been accounted for. + + .. workflow :: + :graph2use: orig + :simple_form: yes + + from sdcflows.workflows.unwarp import init_sdc_unwarp_wf + wf = init_sdc_unwarp_wf(omp_nthreads=8, + debug=False) + + Parameters + ---------- + omp_nthreads : int + Maximum number of threads an individual process may use. + debug : bool + Run fast configurations of registrations. + name : str + Unique name of this workflow. + + Inputs + ------ + in_warp : os.pathlike + The :abbr:`DFM (displacements field map)` that corrects for + susceptibility-derived distortions estimated elsewhere. + in_reference : os.pathlike + the reference image to be unwarped. + in_reference_mask : os.pathlike + the reference image mask to be unwarped + + Outputs + ------- + out_reference : str + the ``in_reference`` after unwarping + out_reference_brain : str + the ``in_reference`` after unwarping and skullstripping + out_warp : str + the ``in_warp`` field is forwarded for compatibility + out_mask : str + mask of the unwarped input file + + """ + workflow = Workflow(name=name) + inputnode = pe.Node(niu.IdentityInterface( + fields=['in_warp', 'in_reference', 'in_reference_mask']), + name='inputnode') + outputnode = pe.Node(niu.IdentityInterface( + fields=['out_reference', 'out_reference_brain', 'out_warp', 'out_mask']), + name='outputnode') + + unwarp_reference = pe.Node(ANTSApplyTransformsRPT(dimension=3, + generate_report=False, + float=True, + interpolation='LanczosWindowedSinc'), + name='unwarp_reference') + + unwarp_mask = pe.Node(ANTSApplyTransformsRPT( + dimension=3, generate_report=False, float=True, + interpolation='NearestNeighbor'), name='unwarp_mask') + + enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf(omp_nthreads=omp_nthreads, + pre_mask=True) + workflow.connect([ + (inputnode, unwarp_reference, [ + ('in_warp', 'transforms'), + ('in_reference', 'reference_image'), + ('in_reference', 'input_image')]), + (inputnode, unwarp_mask, [ + ('in_warp', 'transforms'), + ('in_reference_mask', 'reference_image'), + ('in_reference_mask', 'input_image')]), + (unwarp_reference, enhance_and_skullstrip_bold_wf, [ + ('output_image', 'inputnode.in_file')]), + (unwarp_mask, enhance_and_skullstrip_bold_wf, [ + ('output_image', 'inputnode.pre_mask')]), + (inputnode, outputnode, [('in_warp', 'out_warp')]), + (unwarp_reference, outputnode, [('output_image', 'out_reference')]), + (enhance_and_skullstrip_bold_wf, outputnode, [ + ('outputnode.mask_file', 'out_mask'), + ('outputnode.skull_stripped_file', 'out_reference_brain')]), + ]) + return workflow diff --git a/aslprep/workflows/asl/base.py b/aslprep/workflows/asl/base.py index ecbcfb0aa..cc9556b37 100644 --- a/aslprep/workflows/asl/base.py +++ b/aslprep/workflows/asl/base.py @@ -148,7 +148,7 @@ def init_asl_preproc_wf(asl_file): from ...niworkflows.func.util import init_asl_reference_wf from ...niworkflows.interfaces.nibabel import ApplyMask from ...niworkflows.interfaces.utility import KeySelect - from sdcflows.workflows.base import init_sdc_estimate_wf, fieldmap_wrangler + from ...sdcflows.workflows.base import init_sdc_estimate_wf, fieldmap_wrangler ref_file = asl_file mem_gb = {'filesize': 1, 'resampled': 1, 'largemem': 1} @@ -561,7 +561,7 @@ def init_asl_preproc_wf(asl_file): ]) if fmaps: - from sdcflows.workflows.outputs import init_sdc_unwarp_report_wf + from ...sdcflows.workflows.outputs import init_sdc_unwarp_report_wf # Report on asl correction fmap_unwarp_report_wf = init_sdc_unwarp_report_wf() workflow.connect([ diff --git a/setup.cfg b/setup.cfg index 4d704444f..0d25706f5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -24,15 +24,16 @@ install_requires = nibabel >= 3.0 nipype >= 1.4 nitime + pyslim nitransforms >= 20.0.0rc3,<20.2 - #niworkflows ~= 1.2.5 + niworkflows ~= 1.3.0 numpy >= 1.16.5 sentry-sdk >=0.6.9 pandas psutil >= 5.4 pybids ~= 0.11.0 pyyaml - sdcflows ~= 1.3.1 + #sdcflows ~= 1.3.1 svgutils == 0.3.1 #smriprep ~= 0.6.1 tedana >= 0.0.9a1, < 0.0.10 @@ -98,6 +99,10 @@ aslprep = niworkflows/data/*json niworkflows/data/*txt niworkflows/reports/*tpl + sdcflows/data/*json + sdcflows/data/*txt + sdcflows/data/*mat + sdcflows/data/*nii.gz smriprep/data/*.json smriprep/data/*.bib smriprep/data/*.txt