diff --git a/CHANGES b/CHANGES index 7c071b6cb2..998e65f597 100644 --- a/CHANGES +++ b/CHANGES @@ -1,6 +1,8 @@ Next Release ============ +* ENH: Deep revision of workflows for correction of dMRI artifacts. New dmri_preprocessing + example. * ENH: New FSL interface: EpiReg * ENH: New Freesurfer interface: Tkregister2 (for conversion of fsl style matrices to freesurfer format) * ENH: New FSL interfaces: WarpPoints, WarpPointsToStd. diff --git a/examples/dmri_preprocessing.py b/examples/dmri_preprocessing.py new file mode 100644 index 0000000000..373ae64994 --- /dev/null +++ b/examples/dmri_preprocessing.py @@ -0,0 +1,175 @@ +# coding: utf-8 +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + +""" +=================== +dMRI: Preprocessing +=================== + +Introduction +============ + +This script, dmri_preprocessing.py, demonstrates how to prepare dMRI data +for tractography and connectivity analysis with nipype. + +We perform this analysis using the FSL course data, which can be acquired from +here: http://www.fmrib.ox.ac.uk/fslcourse/fsl_course_data2.tar.gz + +Can be executed in command line using ``python dmri_preprocessing.py`` + + +Import necessary modules from nipype. +""" + +import os # system functions +import nipype.interfaces.io as nio # Data i/o +import nipype.interfaces.utility as niu # utility +import nipype.algorithms.misc as misc + +import nipype.pipeline.engine as pe # pypeline engine + +from nipype.interfaces import fsl +from nipype.interfaces import ants + + +""" +Load specific nipype's workflows for preprocessing of dMRI data: +:class:`nipype.workflows.dmri.preprocess.epi.all_peb_pipeline`, +as data include a *b0* volume with reverse encoding direction +(*P>>>A*, or *y*), in contrast with the general acquisition encoding +that is *A>>>P* or *-y* (in RAS systems). +""" + +from nipype.workflows.dmri.fsl.artifacts import all_fsl_pipeline, remove_bias + +""" +Map field names into individual subject runs +""" + +info = dict(dwi=[['subject_id', 'dwidata']], + bvecs=[['subject_id', 'bvecs']], + bvals=[['subject_id', 'bvals']], + dwi_rev=[['subject_id', 'nodif_PA']]) + +infosource = pe.Node(interface=niu.IdentityInterface(fields=['subject_id']), + name="infosource") + +# Set the subject 1 identifier in subject_list, +# we choose the preproc dataset as it contains uncorrected files. +subject_list = ['subj1_preproc'] + + +"""Here we set up iteration over all the subjects. The following line +is a particular example of the flexibility of the system. The +``datasource`` attribute ``iterables`` tells the pipeline engine that +it should repeat the analysis on each of the items in the +``subject_list``. In the current example, the entire first level +preprocessing and estimation will be repeated for each subject +contained in subject_list. +""" + +infosource.iterables = ('subject_id', subject_list) + + +""" +Now we create a :class:`nipype.interfaces.io.DataGrabber` object and +fill in the information from above about the layout of our data. The +:class:`~nipype.pipeline.engine.Node` module wraps the interface object +and provides additional housekeeping and pipeline specific +functionality. +""" + +datasource = pe.Node(nio.DataGrabber(infields=['subject_id'], + outfields=info.keys()), name='datasource') + +datasource.inputs.template = "%s/%s" + +# This needs to point to the fdt folder you can find after extracting +# http://www.fmrib.ox.ac.uk/fslcourse/fsl_course_data2.tar.gz +datasource.inputs.base_directory = os.path.abspath('fdt1') +datasource.inputs.field_template = dict(dwi='%s/%s.nii.gz', + dwi_rev='%s/%s.nii.gz') +datasource.inputs.template_args = info +datasource.inputs.sort_filelist = True + + +""" +An inputnode is used to pass the data obtained by the data grabber to the +actual processing functions +""" + +inputnode = pe.Node(niu.IdentityInterface(fields=["dwi", "bvecs", "bvals", + "dwi_rev"]), name="inputnode") + + +""" + +Setup for dMRI preprocessing +============================ + +In this section we initialize the appropriate workflow for preprocessing of +diffusion images. + +Artifacts correction +-------------------- + +We will use the combination of ``topup`` and ``eddy`` as suggested by FSL. + +In order to configure the susceptibility distortion correction (SDC), we first +write the specific parameters of our echo-planar imaging (EPI) images. + +Particularly, we look into the ``acqparams.txt`` file of the selected subject +to gather the encoding direction, acceleration factor (in parallel sequences +it is > 1), and readout time or echospacing. + +""" + +epi_AP = {'echospacing': 66.5e-3, 'enc_dir': 'y-'} +epi_PA = {'echospacing': 66.5e-3, 'enc_dir': 'y'} +prep = all_fsl_pipeline(epi_params=epi_AP, altepi_params=epi_PA) + + +""" + +Bias field correction +--------------------- + +Finally, we set up a node to correct for a single multiplicative bias field +from computed on the *b0* image, as suggested in [Jeurissen2014]_. + +""" + +bias = remove_bias() + + +""" +Connect nodes in workflow +========================= + +We create a higher level workflow to connect the nodes. Please excuse the +author for writing the arguments of the ``connect`` function in a not-standard +style with readability aims. +""" + +wf = pe.Workflow(name="dMRI_Preprocessing") +wf.base_dir = os.path.abspath('preprocessing_dmri_tutorial') +wf.connect([ + (infosource, datasource, [('subject_id', 'subject_id')]) + ,(datasource, prep, [('dwi', 'inputnode.in_file'), + ('dwi_rev', 'inputnode.alt_file'), + ('bvals', 'inputnode.in_bval'), + ('bvecs', 'inputnode.in_bvec')]) + ,(prep, bias, [('outputnode.out_file', 'inputnode.in_file'), + ('outputnode.out_mask', 'inputnode.in_mask')]) + ,(datasource, bias, [('bvals', 'inputnode.in_bval')]) +]) + + +""" +Run the workflow as command line executable +""" + +if __name__ == '__main__': + wf.run() + wf.write_graph() diff --git a/nipype/algorithms/tests/test_normalize_tpms.py b/nipype/algorithms/tests/test_normalize_tpms.py index 0f4e3e4d3c..1484874751 100644 --- a/nipype/algorithms/tests/test_normalize_tpms.py +++ b/nipype/algorithms/tests/test_normalize_tpms.py @@ -1,12 +1,6 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- +# coding: utf-8 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -# -# @Author: oesteban - code@oscaresteban.es -# @Date: 2014-05-28 17:57:20 -# @Last Modified by: oesteban -# @Last Modified time: 2014-05-29 13:43:09 import os from shutil import rmtree diff --git a/nipype/interfaces/elastix/__init__.py b/nipype/interfaces/elastix/__init__.py index f8b31c5e50..66819f95db 100644 --- a/nipype/interfaces/elastix/__init__.py +++ b/nipype/interfaces/elastix/__init__.py @@ -1,12 +1,7 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- +# coding: utf-8 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -# -# @Author: oesteban - code@oscaresteban.es -# @Date: 2014-06-02 12:06:07 -# @Last Modified by: oesteban -# @Last Modified time: 2014-06-17 10:59:20 + """Top-level namespace for elastix.""" from registration import Registration, ApplyWarp, AnalyzeWarp, PointsWarp diff --git a/nipype/interfaces/elastix/base.py b/nipype/interfaces/elastix/base.py index c837706159..2cbc5326f5 100644 --- a/nipype/interfaces/elastix/base.py +++ b/nipype/interfaces/elastix/base.py @@ -1,12 +1,7 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- +# coding: utf-8 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -# -# @Author: oesteban - code@oscaresteban.es -# @Date: 2014-06-03 13:42:46 -# @Last Modified by: oesteban -# @Last Modified time: 2014-06-17 10:17:43 + """The :py:mod:`nipype.interfaces.elastix` provides the interface to the elastix registration software. diff --git a/nipype/interfaces/elastix/registration.py b/nipype/interfaces/elastix/registration.py index 8e304935e6..6c4c67aaa1 100644 --- a/nipype/interfaces/elastix/registration.py +++ b/nipype/interfaces/elastix/registration.py @@ -1,12 +1,7 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- +# coding: utf-8 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -# -# @Author: oesteban - code@oscaresteban.es -# @Date: 2014-06-02 12:06:50 -# @Last Modified by: oesteban -# @Last Modified time: 2014-09-01 21:03:57 + """ Interfaces to perform image registrations and to apply the resulting displacement maps to images and points. diff --git a/nipype/interfaces/elastix/utils.py b/nipype/interfaces/elastix/utils.py index 2f361996a0..09d46adaf9 100644 --- a/nipype/interfaces/elastix/utils.py +++ b/nipype/interfaces/elastix/utils.py @@ -1,12 +1,7 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- +# coding: utf-8 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -# -# @Author: oesteban - code@oscaresteban.es -# @Date: 2014-06-17 10:17:07 -# @Last Modified by: oesteban -# @Last Modified time: 2014-09-01 21:05:33 + """ Generic interfaces to manipulate registration parameters files, including transform files (to configure warpings) diff --git a/nipype/workflows/data/__init__.py b/nipype/workflows/data/__init__.py new file mode 100644 index 0000000000..3f0dd77db1 --- /dev/null +++ b/nipype/workflows/data/__init__.py @@ -0,0 +1,15 @@ +# coding: utf-8 +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +import os.path as op + + +def get_flirt_schedule(name): + if name == 'ecc': + return op.abspath(op.join(op.dirname(__file__), + 'ecc.sch')) + elif name == 'hmc': + return op.abspath(op.join(op.dirname(__file__), + 'hmc.sch')) + else: + raise RuntimeError('Requested file does not exist.') diff --git a/nipype/workflows/data/ecc.sch b/nipype/workflows/data/ecc.sch new file mode 100644 index 0000000000..a7de1f2b0b --- /dev/null +++ b/nipype/workflows/data/ecc.sch @@ -0,0 +1,67 @@ +# 4mm scale +setscale 4 +setoption smoothing 6 +setoption paramsubset 1 0 0 0 0 0 0 1 1 1 1 1 1 +clear U +clear UA +clear UB +clear US +clear UP +# try the identity transform as a starting point at this resolution +clear UQ +setrow UQ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 +optimise 7 UQ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 rel 4 +sort U +copy U UA +# select best 4 optimised solutions and try perturbations of these +clear U +copy UA:1-4 U +optimise 7 UA:1-4 1.0 0.0 0.0 0.0 0.0 0.0 0.0 rel 4 +optimise 7 UA:1-4 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 1.0 0.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 0.0 1.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 0.1 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 -0.1 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 0.2 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 -0.2 abs 4 +sort U +copy U UB +# 2mm scale +setscale 2 +setoption smoothing 4 +setoption paramsubset 1 0 0 0 0 0 0 1 1 1 1 1 1 +clear U +clear UC +clear UD +clear UE +clear UF +# remeasure costs at this scale +measurecost 7 UB 0 0 0 0 0 0 rel +sort U +copy U UC +clear U +optimise 7 UC:1-3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 2 +copy U UD +sort U +copy U UF +# also try the identity transform as a starting point at this resolution +sort U +clear U UG +clear U +setrow UG 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 +optimise 7 UG 0.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 2 +sort U +copy U UG +# 1mm scale +setscale 1 +setoption smoothing 2 +setoption boundguess 1 +setoption paramsubset 1 0 0 0 0 0 0 1 1 1 1 1 1 +clear U +#also try the identity transform as a starting point at this resolution +setrow UK 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 +optimise 12 UK:1-2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 1 +sort U + diff --git a/nipype/workflows/data/hmc.sch b/nipype/workflows/data/hmc.sch new file mode 100644 index 0000000000..08f3e76e85 --- /dev/null +++ b/nipype/workflows/data/hmc.sch @@ -0,0 +1,64 @@ +# 4mm scale +setscale 4 +setoption smoothing 6 +clear U +clear UA +clear UB +clear US +clear UP +# try the identity transform as a starting point at this resolution +clear UQ +setrow UQ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 +optimise 7 UQ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 rel 4 +sort U +copy U UA +# select best 4 optimised solutions and try perturbations of these +clear U +copy UA:1-4 U +optimise 7 UA:1-4 1.0 0.0 0.0 0.0 0.0 0.0 0.0 rel 4 +optimise 7 UA:1-4 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 1.0 0.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 0.0 1.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 0.1 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 -0.1 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 0.2 abs 4 +optimise 7 UA:1-4 0.0 0.0 0.0 0.0 0.0 0.0 -0.2 abs 4 +sort U +copy U UB +# 2mm scale +setscale 2 +setoption smoothing 4 +clear U +clear UC +clear UD +clear UE +clear UF +# remeasure costs at this scale +measurecost 7 UB 0 0 0 0 0 0 rel +sort U +copy U UC +clear U +optimise 7 UC:1-3 0.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 2 +copy U UD +sort U +copy U UF +# also try the identity transform as a starting point at this resolution +sort U +clear U UG +clear U +setrow UG 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 +optimise 7 UG 0.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 2 +sort U +copy U UG +# 1mm scale +setscale 1 +setoption smoothing 2 +setoption boundguess 1 +clear U +#also try the identity transform as a starting point at this resolution +setrow UK 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 +optimise 12 UK:1-2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 abs 1 +sort U + diff --git a/nipype/workflows/dmri/fsl/__init__.py b/nipype/workflows/dmri/fsl/__init__.py index 1b9e80b03d..c336329fb5 100644 --- a/nipype/workflows/dmri/fsl/__init__.py +++ b/nipype/workflows/dmri/fsl/__init__.py @@ -1,5 +1,9 @@ from dti import create_bedpostx_pipeline +from artifacts import (all_fmb_pipeline, all_peb_pipeline, all_fsl_pipeline, + hmc_pipeline, ecc_pipeline, sdc_fmb, sdc_peb, + remove_bias) + from epi import (fieldmap_correction, topup_correction, create_eddy_correct_pipeline, create_epidewarp_pipeline, create_dmri_preprocessing) diff --git a/nipype/workflows/dmri/fsl/artifacts.py b/nipype/workflows/dmri/fsl/artifacts.py new file mode 100644 index 0000000000..54be28f914 --- /dev/null +++ b/nipype/workflows/dmri/fsl/artifacts.py @@ -0,0 +1,845 @@ +# coding: utf-8 +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + +import nipype.pipeline.engine as pe +import nipype.interfaces.utility as niu +import nipype.interfaces.fsl as fsl +import nipype.interfaces.freesurfer as fs +import nipype.interfaces.ants as ants +import os + +from .utils import * + + +def all_fmb_pipeline(name='hmc_sdc_ecc', + fugue_params=dict(smooth3d=2.0), + bmap_params=dict(delta_te=2.46e-3), + epi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y-')): + """ + Builds a pipeline including three artifact corrections: head-motion + correction (HMC), susceptibility-derived distortion correction (SDC), + and Eddy currents-derived distortion correction (ECC). + + The displacement fields from each kind of distortions are combined. Thus, + only one interpolation occurs between input data and result. + + .. warning:: this workflow rotates the gradients table (*b*-vectors) + [Leemans09]_. + + + Examples + -------- + + >>> from nipype.workflows.dmri.fsl.artifacts import all_fmb_pipeline + >>> allcorr = all_fmb_pipeline() + >>> allcorr.inputs.inputnode.in_file = 'epi.nii' + >>> allcorr.inputs.inputnode.in_bval = 'diffusion.bval' + >>> allcorr.inputs.inputnode.in_bvec = 'diffusion.bvec' + >>> allcorr.inputs.inputnode.bmap_mag = 'magnitude.nii' + >>> allcorr.inputs.inputnode.bmap_pha = 'phase.nii' + >>> allcorr.run() # doctest: +SKIP + + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bvec', + 'in_bval', 'bmap_pha', 'bmap_mag']), name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_mask', + 'out_bvec']), name='outputnode') + avg_b0_0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg_pre') + avg_b0_1 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg_post') + bet_dwi0 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), + name='bet_dwi_pre') + bet_dwi1 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), + name='bet_dwi_post') + + hmc = hmc_pipeline() + sdc = sdc_fmb(fugue_params=fugue_params, bmap_params=bmap_params, + epi_params=epi_params) + ecc = ecc_pipeline() + unwarp = apply_all_corrections() + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, hmc, [('in_file', 'inputnode.in_file'), + ('in_bvec', 'inputnode.in_bvec'), + ('in_bval', 'inputnode.in_bval')]) + ,(inputnode, avg_b0_0, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(avg_b0_0, bet_dwi0, [('out_file', 'in_file')]) + ,(bet_dwi0, hmc, [('mask_file', 'inputnode.in_mask')]) + ,(hmc, sdc, [('outputnode.out_file', 'inputnode.in_file')]) + ,(bet_dwi0, sdc, [('mask_file', 'inputnode.in_mask')]) + ,(inputnode, sdc, [('in_bval', 'inputnode.in_bval'), + ('bmap_pha', 'inputnode.bmap_pha'), + ('bmap_mag', 'inputnode.bmap_mag')]) + ,(hmc, ecc, [('outputnode.out_xfms', 'inputnode.in_xfms')]) + ,(inputnode, ecc, [('in_file', 'inputnode.in_file'), + ('in_bval', 'inputnode.in_bval')]) + ,(bet_dwi0, ecc, [('mask_file', 'inputnode.in_mask')]) + ,(ecc, avg_b0_1, [('outputnode.out_file', 'in_dwi')]) + ,(inputnode, avg_b0_1, [('in_bval', 'in_bval')]) + ,(avg_b0_1, bet_dwi1, [('out_file', 'in_file')]) + + ,(inputnode, unwarp, [('in_file', 'inputnode.in_dwi')]) + ,(hmc, unwarp, [('outputnode.out_xfms', 'inputnode.in_hmc')]) + ,(ecc, unwarp, [('outputnode.out_xfms', 'inputnode.in_ecc')]) + ,(sdc, unwarp, [('outputnode.out_warp', 'inputnode.in_sdc')]) + + ,(hmc, outputnode, [('outputnode.out_bvec', 'out_bvec')]) + ,(unwarp, outputnode, [('outputnode.out_file', 'out_file')]) + ,(bet_dwi1, outputnode, [('mask_file', 'out_mask')]) + ]) + return wf + + +def all_peb_pipeline(name='hmc_sdc_ecc', + epi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y-', + epi_factor=1), + altepi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y', + epi_factor=1)): + """ + Builds a pipeline including three artifact corrections: head-motion + correction (HMC), susceptibility-derived distortion correction (SDC), + and Eddy currents-derived distortion correction (ECC). + + .. warning:: this workflow rotates the gradients table (*b*-vectors) + [Leemans09]_. + + + Examples + -------- + + >>> from nipype.workflows.dmri.fsl.artifacts import all_peb_pipeline + >>> allcorr = all_peb_pipeline() + >>> allcorr.inputs.inputnode.in_file = 'epi.nii' + >>> allcorr.inputs.inputnode.alt_file = 'epi_rev.nii' + >>> allcorr.inputs.inputnode.in_bval = 'diffusion.bval' + >>> allcorr.inputs.inputnode.in_bvec = 'diffusion.bvec' + >>> allcorr.run() # doctest: +SKIP + + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bvec', + 'in_bval', 'alt_file']), name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_mask', + 'out_bvec']), name='outputnode') + + avg_b0_0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg_pre') + avg_b0_1 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg_post') + bet_dwi0 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), + name='bet_dwi_pre') + bet_dwi1 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), + name='bet_dwi_post') + + hmc = hmc_pipeline() + sdc = sdc_peb(epi_params=epi_params, altepi_params=altepi_params) + ecc = ecc_pipeline() + + unwarp = apply_all_corrections() + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, hmc, [('in_file', 'inputnode.in_file'), + ('in_bvec', 'inputnode.in_bvec'), + ('in_bval', 'inputnode.in_bval')]) + ,(inputnode, avg_b0_0, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(avg_b0_0, bet_dwi0, [('out_file', 'in_file')]) + ,(bet_dwi0, hmc, [('mask_file', 'inputnode.in_mask')]) + ,(hmc, sdc, [('outputnode.out_file', 'inputnode.in_file')]) + ,(bet_dwi0, sdc, [('mask_file', 'inputnode.in_mask')]) + ,(inputnode, sdc, [('in_bval', 'inputnode.in_bval'), + ('alt_file', 'inputnode.alt_file')]) + ,(inputnode, ecc, [('in_file', 'inputnode.in_file'), + ('in_bval', 'inputnode.in_bval')]) + ,(bet_dwi0, ecc, [('mask_file', 'inputnode.in_mask')]) + ,(hmc, ecc, [('outputnode.out_xfms', 'inputnode.in_xfms')]) + ,(ecc, avg_b0_1, [('outputnode.out_file', 'in_dwi')]) + ,(inputnode, avg_b0_1, [('in_bval', 'in_bval')]) + ,(avg_b0_1, bet_dwi1, [('out_file', 'in_file')]) + + + ,(inputnode, unwarp, [('in_file', 'inputnode.in_dwi')]) + ,(hmc, unwarp, [('outputnode.out_xfms', 'inputnode.in_hmc')]) + ,(ecc, unwarp, [('outputnode.out_xfms', 'inputnode.in_ecc')]) + ,(sdc, unwarp, [('outputnode.out_warp', 'inputnode.in_sdc')]) + + ,(hmc, outputnode, [('outputnode.out_bvec', 'out_bvec')]) + ,(unwarp, outputnode, [('outputnode.out_file', 'out_file')]) + ,(bet_dwi1, outputnode, [('mask_file', 'out_mask')]) + ]) + return wf + + +def all_fsl_pipeline(name='fsl_all_correct', + epi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y-'), + altepi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y')): + """ + Workflow that integrates FSL ``topup`` and ``eddy``. + + + .. warning:: this workflow rotates the gradients table (*b*-vectors) [Leemans09]_. + + + .. warning:: this workflow does not perform jacobian modulation of each + *DWI* [Jones10]_. + + + Examples + -------- + + >>> from nipype.workflows.dmri.fsl.artifacts import all_fsl_pipeline + >>> allcorr = all_fsl_pipeline() + >>> allcorr.inputs.inputnode.in_file = 'epi.nii' + >>> allcorr.inputs.inputnode.alt_file = 'epi_rev.nii' + >>> allcorr.inputs.inputnode.in_bval = 'diffusion.bval' + >>> allcorr.inputs.inputnode.in_bvec = 'diffusion.bvec' + >>> allcorr.run() # doctest: +SKIP + + """ + + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bvec', + 'in_bval', 'alt_file']), name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_mask', + 'out_bvec']), name='outputnode') + + def _gen_index(in_file): + import numpy as np + import nibabel as nb + import os + out_file = os.path.abspath('index.txt') + vols = nb.load(in_file).get_data().shape[-1] + np.savetxt(out_file, np.ones((vols,)).T) + return out_file + + avg_b0_0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg_pre') + bet_dwi0 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), + name='bet_dwi_pre') + + sdc = sdc_peb(epi_params=epi_params, altepi_params=altepi_params) + ecc = pe.Node(fsl.Eddy(method='jac'), name='fsl_eddy') + rot_bvec = pe.Node(niu.Function(input_names=['in_bvec', 'eddy_params'], + output_names=['out_file'], function=eddy_rotate_bvecs), + name='Rotate_Bvec') + avg_b0_1 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg_post') + bet_dwi1 = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True), + name='bet_dwi_post') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, avg_b0_0, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(avg_b0_0, bet_dwi0, [('out_file','in_file')]) + ,(bet_dwi0, sdc, [('mask_file', 'inputnode.in_mask')]) + ,(inputnode, sdc, [('in_file', 'inputnode.in_file'), + ('alt_file', 'inputnode.alt_file'), + ('in_bval', 'inputnode.in_bval')]) + + ,(sdc, ecc, [('topup.out_enc_file', 'in_acqp'), + ('topup.out_fieldcoef', 'in_topup_fieldcoef'), + ('topup.out_movpar', 'in_topup_movpar')]) + ,(bet_dwi0, ecc, [('mask_file', 'in_mask')]) + ,(inputnode, ecc, [('in_file', 'in_file'), + (('in_file', _gen_index), 'in_index'), + ('in_bval', 'in_bval'), + ('in_bvec', 'in_bvec')]) + ,(inputnode, rot_bvec, [('in_bvec', 'in_bvec')]) + ,(ecc, rot_bvec, [('out_parameter', 'eddy_params')]) + ,(ecc, avg_b0_1, [('out_corrected', 'in_dwi')]) + ,(inputnode, avg_b0_1, [('in_bval', 'in_bval')]) + ,(avg_b0_1, bet_dwi1, [('out_file','in_file')]) + ,(ecc, outputnode, [('out_corrected', 'out_file')]) + ,(rot_bvec, outputnode, [('out_file', 'out_bvec')]) + ,(bet_dwi1, outputnode, [('mask_file', 'out_mask')]) + ]) + return wf + + +def hmc_pipeline(name='motion_correct'): + """ + HMC stands for head-motion correction. + + Creates a pipeline that corrects for head motion artifacts in dMRI + sequences. + It takes a series of diffusion weighted images and rigidly co-registers + them to one reference image. Finally, the `b`-matrix is rotated accordingly + [Leemans09]_ making use of the rotation matrix obtained by FLIRT. + + Search angles have been limited to 4 degrees, based on results in + [Yendiki13]_. + + A list of rigid transformation matrices is provided, so that transforms + can be chained. + This is useful to correct for artifacts with only one interpolation process + (as previously discussed `here + `_), + and also to compute nuisance regressors as proposed by [Yendiki13]_. + + .. warning:: This workflow rotates the `b`-vectors, so please be advised + that not all the dicom converters ensure the consistency between the + resulting nifti orientation and the gradients table (e.g. dcm2nii + checks it). + + .. admonition:: References + + .. [Leemans09] Leemans A, and Jones DK, `The B-matrix must be rotated + when correcting for subject motion in DTI data + `_, + Magn Reson Med. 61(6):1336-49. 2009. doi: 10.1002/mrm.21890. + + .. [Yendiki13] Yendiki A et al., `Spurious group differences due to head + motion in a diffusion MRI study + `_. + Neuroimage. 21(88C):79-90. 2013. doi: 10.1016/j.neuroimage.2013.11.027 + + Example + ------- + + >>> from nipype.workflows.dmri.fsl.artifacts import hmc_pipeline + >>> hmc = hmc_pipeline() + >>> hmc.inputs.inputnode.in_file = 'diffusion.nii' + >>> hmc.inputs.inputnode.in_bvec = 'diffusion.bvec' + >>> hmc.inputs.inputnode.in_bval = 'diffusion.bval' + >>> hmc.inputs.inputnode.in_mask = 'mask.nii' + >>> hmc.run() # doctest: +SKIP + + Inputs:: + + inputnode.in_file - input dwi file + inputnode.in_mask - weights mask of reference image (a file with data range \ +in [0.0, 1.0], indicating the weight of each voxel when computing the metric. + inputnode.in_bvec - gradients file (b-vectors) + inputnode.ref_num (optional, default=0) index of the b0 volume that should be \ +taken as reference + + Outputs:: + + outputnode.out_file - corrected dwi file + outputnode.out_bvec - rotated gradient vectors table + outputnode.out_xfms - list of transformation matrices + + """ + from nipype.workflows.data import get_flirt_schedule + + params = dict(dof=6, bgvalue=0, save_log=True, no_search=True, + #cost='mutualinfo', cost_func='mutualinfo', bins=64, + schedule=get_flirt_schedule('hmc')) + + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'ref_num', + 'in_bvec', 'in_bval', 'in_mask']), name='inputnode') + refid = pe.Node(niu.IdentityInterface(fields=['ref_num']), + name='RefVolume') + pick_ref = pe.Node(fsl.ExtractROI(t_size=1), name='GetRef') + pick_mov = pe.Node(niu.Function(input_names=['in_file', 'in_bval', 'volid'], + output_names=['out_file', 'out_bval'], + function=remove_comp), name='GetMovingDWs') + flirt = dwi_flirt(flirt_param=params) + insmat = pe.Node(niu.Function(input_names=['inlist', 'volid'], + output_names=['out'], function=insert_mat), + name='InsertRefmat') + rot_bvec = pe.Node(niu.Function(input_names=['in_bvec', 'in_matrix'], + output_names=['out_file'], function=rotate_bvecs), + name='Rotate_Bvec') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', + 'out_bvec', 'out_xfms']), + name='outputnode') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, refid, [(('ref_num', _checkrnum), 'ref_num')]) + ,(inputnode, pick_ref, [('in_file', 'in_file')]) + ,(inputnode, pick_mov, [('in_file', 'in_file'), + ('in_bval', 'in_bval')]) + ,(inputnode, flirt, [('in_mask', 'inputnode.ref_mask')]) + ,(refid, pick_ref, [('ref_num', 't_min')]) + ,(refid, pick_mov, [('ref_num', 'volid')]) + ,(pick_mov, flirt, [('out_file', 'inputnode.in_file'), + ('out_bval', 'inputnode.in_bval')]) + + ,(pick_ref, flirt, [('roi_file', 'inputnode.reference')]) + ,(flirt, insmat, [('outputnode.out_xfms', 'inlist')]) + ,(refid, insmat, [('ref_num', 'volid')]) + ,(inputnode, rot_bvec, [('in_bvec', 'in_bvec')]) + ,(insmat, rot_bvec, [('out', 'in_matrix')]) + ,(rot_bvec, outputnode, [('out_file', 'out_bvec')]) + ,(flirt, outputnode, [('outputnode.out_file', 'out_file')]) + ,(insmat, outputnode, [('out', 'out_xfms')]) + ]) + return wf + + +def ecc_pipeline(name='eddy_correct'): + """ + ECC stands for Eddy currents correction. + + Creates a pipeline that corrects for artifacts induced by Eddy currents in dMRI + sequences. + It takes a series of diffusion weighted images and linearly co-registers + them to one reference image (the average of all b0s in the dataset). + + DWIs are also modulated by the determinant of the Jacobian as indicated by [Jones10]_ and + [Rohde04]_. + + A list of rigid transformation matrices can be provided, sourcing from a + :func:`.hmc_pipeline` workflow, to initialize registrations in a *motion free* + framework. + + A list of affine transformation matrices is available as output, so that transforms + can be chained (discussion + `here `_). + + .. admonition:: References + + .. [Jones10] Jones DK, `The signal intensity must be modulated by the determinant of + the Jacobian when correcting for eddy currents in diffusion MRI + `_, + Proc. ISMRM 18th Annual Meeting, (2010). + + .. [Rohde04] Rohde et al., `Comprehensive Approach for Correction of Motion and Distortion + in Diffusion-Weighted MRI + `_, MRM 51:103-114 (2004). + + Example + ------- + + >>> from nipype.workflows.dmri.fsl.artifacts import ecc_pipeline + >>> ecc = ecc_pipeline() + >>> ecc.inputs.inputnode.in_file = 'diffusion.nii' + >>> ecc.inputs.inputnode.in_bval = 'diffusion.bval' + >>> ecc.inputs.inputnode.in_mask = 'mask.nii' + >>> ecc.run() # doctest: +SKIP + + Inputs:: + + inputnode.in_file - input dwi file + inputnode.in_mask - weights mask of reference image (a file with data range \ +sin [0.0, 1.0], indicating the weight of each voxel when computing the metric. + inputnode.in_bval - b-values table + inputnode.in_xfms - list of matrices to initialize registration (from head-motion correction) + + Outputs:: + + outputnode.out_file - corrected dwi file + outputnode.out_xfms - list of transformation matrices + """ + + from nipype.workflows.data import get_flirt_schedule + params = dict(dof=12, no_search=True, interp='spline', bgvalue=0, + schedule=get_flirt_schedule('ecc')) + # cost='normmi', cost_func='normmi', bins=64, + + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bval', + 'in_mask', 'in_xfms']), name='inputnode') + avg_b0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg') + pick_dws = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval', 'b'], + output_names=['out_file'], function=extract_bval), + name='ExtractDWI') + pick_dws.inputs.b = 'diff' + + flirt = dwi_flirt(flirt_param=params, excl_nodiff=True) + + mult = pe.MapNode(fsl.BinaryMaths(operation='mul'), name='ModulateDWIs', + iterfield=['in_file', 'operand_value']) + thres = pe.MapNode(fsl.Threshold(thresh=0.0), iterfield=['in_file'], + name='RemoveNegative') + + split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs') + get_mat = pe.Node(niu.Function(input_names=['in_bval', 'in_xfms'], + output_names=['out_files'], function=recompose_xfm), + name='GatherMatrices') + merge = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval', 'in_corrected'], + output_names=['out_file'], function=recompose_dwi), name='MergeDWIs') + + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_xfms']), + name='outputnode') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, avg_b0, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(inputnode, pick_dws, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(inputnode, merge, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(inputnode, flirt, [('in_mask', 'inputnode.ref_mask'), + ('in_xfms', 'inputnode.in_xfms'), + ('in_bval', 'inputnode.in_bval')]) + ,(inputnode, get_mat, [('in_bval', 'in_bval')]) + ,(avg_b0, flirt, [('out_file', 'inputnode.reference')]) + ,(pick_dws, flirt, [('out_file', 'inputnode.in_file')]) + ,(flirt, get_mat, [('outputnode.out_xfms', 'in_xfms')]) + ,(flirt, mult, [(('outputnode.out_xfms',_xfm_jacobian), + 'operand_value')]) + ,(flirt, split, [('outputnode.out_file', 'in_file')]) + ,(split, mult, [('out_files', 'in_file')]) + ,(mult, thres, [('out_file', 'in_file')]) + ,(thres, merge, [('out_file', 'in_corrected')]) + ,(get_mat, outputnode, [('out_files', 'out_xfms')]) + ,(merge, outputnode, [('out_file', 'out_file')]) + ]) + return wf + + +def sdc_fmb(name='fmb_correction', + fugue_params=dict(smooth3d=2.0), + bmap_params=dict(delta_te=2.46e-3), + epi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y-')): + """ + SDC stands for susceptibility distortion correction. FMB stands for fieldmap-based. + + The fieldmap based method (FMB) implements SDC by using a mapping of the B0 field + as proposed by [Jezzard95]_. This workflow uses the implementation of FSL + (`FUGUE `_). Phase unwrapping is performed + using `PRELUDE `_ + [Jenkinson03]_. Preparation of the fieldmap is performed reproducing the script in FSL + `fsl_prepare_fieldmap `_. + + Example + ------- + + >>> from nipype.workflows.dmri.fsl.artifacts import sdc_fmb + >>> fmb = sdc_fmb() + >>> fmb.inputs.inputnode.in_file = 'diffusion.nii' + >>> fmb.inputs.inputnode.in_bval = 'diffusion.bval' + >>> fmb.inputs.inputnode.in_mask = 'mask.nii' + >>> fmb.inputs.inputnode.bmap_mag = 'magnitude.nii' + >>> fmb.inputs.inputnode.bmap_pha = 'phase.nii' + >>> fmb.run() # doctest: +SKIP + + .. warning:: Only SIEMENS format fieldmaps are supported. + + .. admonition:: References + + .. [Jezzard95] Jezzard P, and Balaban RS, `Correction for geometric distortion in echo + planar images from B0 field variations `_, + MRM 34(1):65-73. (1995). doi: 10.1002/mrm.1910340111. + + .. [Jenkinson03] Jenkinson M., `Fast, automated, N-dimensional phase-unwrapping algorithm + `_, MRM 49(1):193-197, 2003, doi: 10.1002/mrm.10354. + + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bval', + 'in_mask', 'bmap_pha', 'bmap_mag']), + name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_vsm', + 'out_warp']), + name='outputnode') + + delta_te = epi_params['echospacing'] / (1.0 * epi_params['acc_factor']) + + firstmag = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name='GetFirst') + n4 = pe.Node(ants.N4BiasFieldCorrection(dimension=3), name='Bias') + bet = pe.Node(fsl.BET(frac=0.4, mask=True), name='BrainExtraction') + dilate = pe.Node(fsl.maths.MathsCommand(nan2zeros=True, + args='-kernel sphere 5 -dilM'), name='MskDilate') + pha2rads = pe.Node(niu.Function(input_names=['in_file'], output_names=['out_file'], + function=siemens2rads), name='PreparePhase') + prelude = pe.Node(fsl.PRELUDE(process3d=True), name='PhaseUnwrap') + rad2rsec = pe.Node(niu.Function(input_names=['in_file', 'delta_te'], + output_names=['out_file'], function=rads2radsec), name='ToRadSec') + rad2rsec.inputs.delta_te = bmap_params['delta_te'] + + avg_b0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), name='b0_avg') + + flirt = pe.Node(fsl.FLIRT(interp='spline', cost='normmi', cost_func='normmi', + dof=6, bins=64, save_log=True, padding_size=10, + searchr_x=[-4,4], searchr_y=[-4,4], searchr_z=[-4,4], + fine_search=1, coarse_search=10), + name='BmapMag2B0') + applyxfm = pe.Node(fsl.ApplyXfm(interp='spline', padding_size=10, apply_xfm=True), + name='BmapPha2B0') + + pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='PreliminaryFugue') + demean = pe.Node(niu.Function(input_names=['in_file', 'in_mask'], + output_names=['out_file'], function=demean_image), + name='DemeanFmap') + + cleanup = cleanup_edge_pipeline() + + addvol = pe.Node(niu.Function(input_names=['in_file'], output_names=['out_file'], + function=add_empty_vol), name='AddEmptyVol') + + vsm = pe.Node(fsl.FUGUE(save_shift=True, **fugue_params), + name="ComputeVSM") + vsm.inputs.asym_se_time = bmap_params['delta_te'] + vsm.inputs.dwell_time = delta_te + + split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs') + merge = pe.Node(fsl.Merge(dimension='t'), name='MergeDWIs') + unwarp = pe.MapNode(fsl.FUGUE(icorr=True, forward_warping=False), + iterfield=['in_file'], name='UnwarpDWIs') + unwarp.inputs.unwarp_direction = epi_params['enc_dir'] + thres = pe.MapNode(fsl.Threshold(thresh=0.0), iterfield=['in_file'], + name='RemoveNegative') + vsm2dfm = vsm2warp() + vsm2dfm.inputs.inputnode.scaling = 1.0 + vsm2dfm.inputs.inputnode.enc_dir = epi_params['enc_dir'] + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, pha2rads, [('bmap_pha', 'in_file')]) + ,(inputnode, firstmag, [('bmap_mag', 'in_file')]) + ,(inputnode, avg_b0, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(firstmag, n4, [('roi_file', 'input_image')]) + ,(n4, bet, [('output_image', 'in_file')]) + ,(bet, dilate, [('mask_file', 'in_file')]) + ,(pha2rads, prelude, [('out_file', 'phase_file')]) + ,(n4, prelude, [('output_image', 'magnitude_file')]) + ,(dilate, prelude, [('out_file', 'mask_file')]) + ,(prelude, rad2rsec, [('unwrapped_phase_file', 'in_file')]) + ,(avg_b0, flirt, [('out_file', 'reference')]) + ,(inputnode, flirt, [('in_mask', 'ref_weight')]) + ,(n4, flirt, [('output_image', 'in_file')]) + ,(dilate, flirt, [('out_file', 'in_weight')]) + ,(avg_b0, applyxfm, [('out_file', 'reference')]) + ,(rad2rsec, applyxfm, [('out_file', 'in_file')]) + ,(flirt, applyxfm, [('out_matrix_file', 'in_matrix_file')]) + ,(applyxfm, pre_fugue, [('out_file', 'fmap_in_file')]) + ,(inputnode, pre_fugue, [('in_mask', 'mask_file')]) + ,(pre_fugue, demean, [('fmap_out_file', 'in_file')]) + ,(inputnode, demean, [('in_mask', 'in_mask')]) + ,(demean, cleanup, [('out_file', 'inputnode.in_file')]) + ,(inputnode, cleanup, [('in_mask', 'inputnode.in_mask')]) + ,(cleanup, addvol, [('outputnode.out_file', 'in_file')]) + ,(inputnode, vsm, [('in_mask', 'mask_file')]) + ,(addvol, vsm, [('out_file', 'fmap_in_file')]) + ,(inputnode, split, [('in_file', 'in_file')]) + ,(split, unwarp, [('out_files', 'in_file')]) + ,(vsm, unwarp, [('shift_out_file', 'shift_in_file')]) + ,(unwarp, thres, [('unwarped_file', 'in_file')]) + ,(thres, merge, [('out_file', 'in_files')]) + ,(merge, vsm2dfm, [('merged_file', 'inputnode.in_ref')]) + ,(vsm, vsm2dfm, [('shift_out_file', 'inputnode.in_vsm')]) + ,(merge, outputnode, [('merged_file', 'out_file')]) + ,(vsm, outputnode, [('shift_out_file', 'out_vsm')]) + ,(vsm2dfm, outputnode, [('outputnode.out_warp', 'out_warp')]) + ]) + return wf + + +def sdc_peb(name='peb_correction', + epi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y-', + epi_factor=1), + altepi_params=dict(echospacing=0.77e-3, + acc_factor=3, + enc_dir='y', + epi_factor=1)): + """ + SDC stands for susceptibility distortion correction. PEB stands for + phase-encoding-based. + + The phase-encoding-based (PEB) method implements SDC by acquiring + diffusion images with two different enconding directions [Andersson2003]_. + The most typical case is acquiring with opposed phase-gradient blips + (e.g. *A>>>P* and *P>>>A*, or equivalently, *-y* and *y*) + as in [Chiou2000]_, but it is also possible to use orthogonal + configurations [Cordes2000]_ (e.g. *A>>>P* and *L>>>R*, + or equivalently *-y* and *x*). + This workflow uses the implementation of FSL + (`TOPUP `_). + + Example + ------- + + >>> from nipype.workflows.dmri.fsl.artifacts import sdc_peb + >>> peb = sdc_peb() + >>> peb.inputs.inputnode.in_file = 'epi.nii' + >>> peb.inputs.inputnode.alt_file = 'epi_rev.nii' + >>> peb.inputs.inputnode.in_bval = 'diffusion.bval' + >>> peb.inputs.inputnode.in_mask = 'mask.nii' + >>> peb.run() # doctest: +SKIP + + .. admonition:: References + + .. [Andersson2003] Andersson JL et al., `How to correct susceptibility + distortions in spin-echo echo-planar images: application to diffusion + tensor imaging `_. + Neuroimage. 2003 Oct;20(2):870-88. doi: 10.1016/S1053-8119(03)00336-7 + + .. [Cordes2000] Cordes D et al., Geometric distortion correction in EPI + using two images with orthogonal phase-encoding directions, in Proc. + ISMRM (8), p.1712, Denver, US, 2000. + + .. [Chiou2000] Chiou JY, and Nalcioglu O, A simple method to correct + off-resonance related distortion in echo planar imaging, in Proc. + ISMRM (8), p.1712, Denver, US, 2000. + + """ + + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bval', + 'in_mask', 'alt_file', 'ref_num']), + name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_vsm', + 'out_warp']), name='outputnode') + + b0_ref = pe.Node(fsl.ExtractROI(t_size=1), name='b0_ref') + b0_alt = pe.Node(fsl.ExtractROI(t_size=1), name='b0_alt') + b0_comb = pe.Node(niu.Merge(2), name='b0_list') + b0_merge = pe.Node(fsl.Merge(dimension='t'), name='b0_merged') + + topup = pe.Node(fsl.TOPUP(), name='topup') + topup.inputs.encoding_direction = [epi_params['enc_dir'], + altepi_params['enc_dir']] + + readout = compute_readout(epi_params) + topup.inputs.readout_times = [readout, + compute_readout(altepi_params)] + + unwarp = pe.Node(fsl.ApplyTOPUP(in_index=[1], method='jac'), name='unwarp') + + #scaling = pe.Node(niu.Function(input_names=['in_file', 'enc_dir'], + # output_names=['factor'], function=_get_zoom), + # name='GetZoom') + #scaling.inputs.enc_dir = epi_params['enc_dir'] + vsm2dfm = vsm2warp() + vsm2dfm.inputs.inputnode.enc_dir = epi_params['enc_dir'] + vsm2dfm.inputs.inputnode.scaling = readout + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, b0_ref, [('in_file', 'in_file'), + (('ref_num', _checkrnum), 't_min')]) + ,(inputnode, b0_alt, [('alt_file', 'in_file'), + (('ref_num', _checkrnum), 't_min')]) + ,(b0_ref, b0_comb, [('roi_file', 'in1')]) + ,(b0_alt, b0_comb, [('roi_file', 'in2')]) + ,(b0_comb, b0_merge, [('out', 'in_files')]) + ,(b0_merge, topup, [('merged_file', 'in_file')]) + ,(topup, unwarp, [('out_fieldcoef', 'in_topup_fieldcoef'), + ('out_movpar', 'in_topup_movpar'), + ('out_enc_file', 'encoding_file')]) + ,(inputnode, unwarp, [('in_file', 'in_files')]) + ,(unwarp, outputnode, [('out_corrected', 'out_file')]) + + #,(b0_ref, scaling, [('roi_file', 'in_file')]) + #,(scaling, vsm2dfm, [('factor', 'inputnode.scaling')]) + ,(b0_ref, vsm2dfm, [('roi_file', 'inputnode.in_ref')]) + ,(topup, vsm2dfm, [('out_field', 'inputnode.in_vsm')]) + ,(topup, outputnode, [('out_field', 'out_vsm')]) + ,(vsm2dfm, outputnode, [('outputnode.out_warp', 'out_warp')]) + ]) + return wf + + +def remove_bias(name='bias_correct'): + """ + This workflow estimates a single multiplicative bias field from the averaged *b0* + image, as suggested in [Jeurissen2014]_. + + .. admonition:: References + + .. [Jeurissen2014] Jeurissen B. et al., `Multi-tissue constrained spherical deconvolution + for improved analysis of multi-shell diffusion MRI data + `_. NeuroImage (2014). + doi: 10.1016/j.neuroimage.2014.07.061 + + + Example + ------- + + >>> from nipype.workflows.dmri.fsl.artifacts import remove_bias + >>> bias = remove_bias() + >>> bias.inputs.inputnode.in_file = 'epi.nii' + >>> bias.inputs.inputnode.in_bval = 'diffusion.bval' + >>> bias.inputs.inputnode.in_mask = 'mask.nii' + >>> bias.run() # doctest: +SKIP + + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_bval', + 'in_mask']), name='inputnode') + + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file']), + name='outputnode') + + avg_b0 = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval'], + output_names=['out_file'], function=b0_average), + name='b0_avg') + n4 = pe.Node(ants.N4BiasFieldCorrection(dimension=3, + save_bias=True, bspline_fitting_distance=600), name='Bias_b0') + split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs') + mult = pe.MapNode(fsl.MultiImageMaths(op_string='-div %s'), + iterfield=['in_file'], name='RemoveBiasOfDWIs') + thres = pe.MapNode(fsl.Threshold(thresh=0.0), iterfield=['in_file'], + name='RemoveNegative') + merge = pe.Node(fsl.utils.Merge(dimension='t'), name='MergeDWIs') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, avg_b0, [('in_file', 'in_dwi'), + ('in_bval', 'in_bval')]) + ,(avg_b0, n4, [('out_file', 'input_image')]) + ,(inputnode, n4, [('in_mask', 'mask_image')]) + ,(inputnode, split, [('in_file', 'in_file')]) + ,(n4, mult, [('bias_image', 'operand_files')]) + ,(split, mult, [('out_files', 'in_file')]) + ,(mult, thres, [('out_file', 'in_file')]) + ,(thres, merge, [('out_file', 'in_files')]) + ,(merge, outputnode, [('merged_file', 'out_file')]) + ]) + return wf + + +def _checkrnum(ref_num): + from nipype.interfaces.base import isdefined + if (ref_num is None) or not isdefined(ref_num): + return 0 + return ref_num + + +def _nonb0(in_bval): + import numpy as np + bvals = np.loadtxt(in_bval) + return np.where(bvals != 0)[0].tolist() + + +def _xfm_jacobian(in_xfm): + import numpy as np + from math import fabs + return [fabs(np.linalg.det(np.loadtxt(xfm))) for xfm in in_xfm] + + +def _get_zoom(in_file, enc_dir): + import nibabel as nb + + zooms = nb.load(in_file).get_header().get_zooms() + + if 'y' in enc_dir: + return zooms[1] + elif 'x' in enc_dir: + return zooms[0] + elif 'z' in enc_dir: + return zooms[2] + else: + raise ValueError('Wrong encoding direction string') diff --git a/nipype/workflows/dmri/fsl/epi.py b/nipype/workflows/dmri/fsl/epi.py index 48fb6e8b32..2d1da79ca4 100644 --- a/nipype/workflows/dmri/fsl/epi.py +++ b/nipype/workflows/dmri/fsl/epi.py @@ -4,17 +4,22 @@ import nipype.interfaces.utility as niu import nipype.interfaces.fsl as fsl import os +import warnings def create_dmri_preprocessing(name='dMRI_preprocessing', use_fieldmap=True, fieldmap_registration=False): - """Creates a workflow that chains the necessary pipelines to + """ + Creates a workflow that chains the necessary pipelines to correct for motion, eddy currents, and, if selected, susceptibility artifacts in EPI dMRI sequences. - .. warning:: + .. deprecated:: 0.9.3 + Use :func:`nipype.workflows.dmri.preprocess.epi.all_fmb_pipeline` or + :func:`nipype.workflows.dmri.preprocess.epi.all_peb_pipeline` instead. + - IMPORTANT NOTICE: this workflow rotates the b-vectors, so please be adviced - that not all the dicom converters ensure the consistency between the resulting - nifti orientation and the b matrix table (e.g. dcm2nii checks it). + .. warning:: This workflow rotates the b-vectors, so please be + advised that not all the dicom converters ensure the consistency between the resulting + nifti orientation and the b matrix table (e.g. dcm2nii checks it). Example @@ -54,12 +59,16 @@ def create_dmri_preprocessing(name='dMRI_preprocessing', use_fieldmap=True, fiel Optional arguments:: + use_fieldmap - True if there are fieldmap files that should be used (default True) fieldmap_registration - True if registration to fieldmap should be performed (default False) """ + warnings.warn(('This workflow is deprecated from v.1.0.0, use of available ' + 'nipype.workflows.dmri.preprocess.epi.all_*'), DeprecationWarning) + pipeline = pe.Workflow(name=name) inputnode = pe.Node(niu.IdentityInterface( @@ -118,11 +127,14 @@ def create_motion_correct_pipeline(name='motion_correct'): (Leemans et al. 2009 - http://www.ncbi.nlm.nih.gov/pubmed/19319973), making use of the rotation matrix obtained by FLIRT. - .. warning:: - IMPORTANT NOTICE: this workflow rotates the b-vectors, so please be adviced - that not all the dicom converters ensure the consistency between the resulting - nifti orientation and the b matrix table (e.g. dcm2nii checks it). + .. deprecated:: 0.9.3 + Use :func:`nipype.workflows.dmri.preprocess.epi.hmc_pipeline` instead. + + + .. warning:: This workflow rotates the b-vectors, so please be adviced + that not all the dicom converters ensure the consistency between the resulting + nifti orientation and the b matrix table (e.g. dcm2nii checks it). Example @@ -147,6 +159,10 @@ def create_motion_correct_pipeline(name='motion_correct'): """ + warnings.warn(('This workflow is deprecated from v.1.0.0, use ' + 'nipype.workflows.dmri.preprocess.epi.hmc_pipeline instead'), + DeprecationWarning) + inputnode = pe.Node( niu.IdentityInterface( fields=['in_file', 'ref_num', 'in_bvec']), @@ -183,7 +199,13 @@ def create_motion_correct_pipeline(name='motion_correct'): def create_eddy_correct_pipeline(name='eddy_correct'): - """Creates a pipeline that replaces eddy_correct script in FSL. It takes a + """ + + .. deprecated:: 0.9.3 + Use :func:`nipype.workflows.dmri.preprocess.epi.ecc_pipeline` instead. + + + Creates a pipeline that replaces eddy_correct script in FSL. It takes a series of diffusion weighted images and linearly co-registers them to one reference image. No rotation of the B-matrix is performed, so this pipeline should be executed after the motion correction pipeline. @@ -206,6 +228,10 @@ def create_eddy_correct_pipeline(name='eddy_correct'): outputnode.eddy_corrected """ + warnings.warn(('This workflow is deprecated from v.1.0.0, use ' + 'nipype.workflows.dmri.preprocess.epi.ecc_pipeline instead'), + DeprecationWarning) + inputnode = pe.Node( niu.IdentityInterface(fields=['in_file', 'ref_num']), name='inputnode') @@ -237,6 +263,11 @@ def create_eddy_correct_pipeline(name='eddy_correct'): def fieldmap_correction(name='fieldmap_correction', nocheck=False): """ + + .. deprecated:: 0.9.3 + Use :func:`nipype.workflows.dmri.preprocess.epi.sdc_fmb` instead. + + Fieldmap-based retrospective correction of EPI images for the susceptibility distortion artifact (Jezzard et al., 1995). Fieldmap images are assumed to be already registered to EPI data, and a brain mask is required. @@ -246,7 +277,6 @@ def fieldmap_correction(name='fieldmap_correction', nocheck=False): available as of FSL 5.0. - Example ------- @@ -282,6 +312,10 @@ def fieldmap_correction(name='fieldmap_correction', nocheck=False): """ + warnings.warn(('This workflow is deprecated from v.1.0.0, use ' + 'nipype.workflows.dmri.preprocess.epi.sdc_fmb instead'), + DeprecationWarning) + inputnode = pe.Node(niu.IdentityInterface( fields=['in_file', 'in_mask', @@ -360,8 +394,13 @@ def fieldmap_correction(name='fieldmap_correction', nocheck=False): def topup_correction( name='topup_correction' ): """ - Corrects for susceptibilty distortion of EPI images when one reverse encoding dataset has - been acquired + + .. deprecated:: 0.9.3 + Use :func:`nipype.workflows.dmri.preprocess.epi.sdc_peb` instead. + + + Corrects for susceptibilty distortion of EPI images when one reverse encoding dataset has + been acquired Example @@ -374,18 +413,26 @@ def topup_correction( name='topup_correction' ): >>> nipype_epicorrect.inputs.inputnode.ref_num = 0 >>> nipype_epicorrect.run() # doctest: +SKIP + Inputs:: + inputnode.in_file_dir - EPI volume acquired in 'forward' phase encoding inputnode.in_file_rev - EPI volume acquired in 'reversed' phase encoding inputnode.encoding_direction - Direction encoding of in_file_dir inputnode.ref_num - Identifier of the reference volumes (usually B0 volume) + Outputs:: outputnode.epi_corrected """ + + warnings.warn(('This workflow is deprecated from v.1.0.0, use ' + 'nipype.workflows.dmri.preprocess.epi.sdc_peb instead'), + DeprecationWarning) + pipeline = pe.Workflow(name=name) inputnode = pe.Node(niu.IdentityInterface( @@ -435,11 +482,18 @@ def topup_correction( name='topup_correction' ): def create_epidewarp_pipeline(name='epidewarp', fieldmap_registration=False): - """ Replaces the epidewarp.fsl script (http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) + """ + Replaces the epidewarp.fsl script (http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl) for susceptibility distortion correction of dMRI & fMRI acquired with EPI sequences and the fieldmap information (Jezzard et al., 1995) using FSL's FUGUE. The registration to the (warped) fieldmap (strictly following the original script) is available using fieldmap_registration=True. + + .. warning:: This workflow makes use of ``epidewarp.fsl`` a script of FSL deprecated long + time ago. The use of this workflow is not recommended, use + :func:`nipype.workflows.dmri.preprocess.epi.sdc_fmb` instead. + + Example ------- @@ -479,6 +533,9 @@ def create_epidewarp_pipeline(name='epidewarp', fieldmap_registration=False): """ + warnings.warn(('This workflow reproduces a deprecated FSL script.'), + DeprecationWarning) + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'fieldmap_mag', 'fieldmap_pha', diff --git a/nipype/workflows/dmri/fsl/utils.py b/nipype/workflows/dmri/fsl/utils.py new file mode 100644 index 0000000000..4b5561aa12 --- /dev/null +++ b/nipype/workflows/dmri/fsl/utils.py @@ -0,0 +1,701 @@ +# coding: utf-8 +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + +import nipype.pipeline.engine as pe +import nipype.interfaces.utility as niu +from nipype.interfaces import fsl +from nipype.interfaces import ants + + +def cleanup_edge_pipeline(name='Cleanup'): + """ + Perform some de-spiking filtering to clean up the edge of the fieldmap + (copied from fsl_prepare_fieldmap) + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['in_file', 'in_mask']), + name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file']), + name='outputnode') + + fugue = pe.Node(fsl.FUGUE(save_fmap=True, despike_2dfilter=True, + despike_threshold=2.1), name='Despike') + erode = pe.Node(fsl.maths.MathsCommand(nan2zeros=True, + args='-kernel 2D -ero'), name='MskErode') + newmsk = pe.Node(fsl.MultiImageMaths(op_string='-sub %s -thr 0.5 -bin'), + name='NewMask') + applymsk = pe.Node(fsl.ApplyMask(nan2zeros=True), name='ApplyMask') + join = pe.Node(niu.Merge(2), name='Merge') + addedge = pe.Node(fsl.MultiImageMaths(op_string='-mas %s -add %s'), + name='AddEdge') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, fugue, [('in_file', 'fmap_in_file'), + ('in_mask', 'mask_file')]) + ,(inputnode, erode, [('in_mask', 'in_file')]) + ,(inputnode, newmsk, [('in_mask', 'in_file')]) + ,(erode, newmsk, [('out_file', 'operand_files')]) + ,(fugue, applymsk, [('fmap_out_file', 'in_file')]) + ,(newmsk, applymsk, [('out_file', 'mask_file')]) + ,(erode, join, [('out_file', 'in1')]) + ,(applymsk, join, [('out_file', 'in2')]) + ,(inputnode, addedge, [('in_file', 'in_file')]) + ,(join, addedge, [('out', 'operand_files')]) + ,(addedge, outputnode, [('out_file', 'out_file')]) + ]) + return wf + + +def vsm2warp(name='Shiftmap2Warping'): + """ + Converts a voxel shift map (vsm) to a displacements field (warp). + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['in_vsm', + 'in_ref', 'scaling', 'enc_dir']), name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_warp']), + name='outputnode') + fixhdr = pe.Node(niu.Function(input_names=['in_file', 'in_file_hdr'], + output_names=['out_file'], function=copy_hdr), + name='Fix_hdr') + vsm = pe.Node(fsl.maths.BinaryMaths(operation='mul'), name='ScaleField') + vsm2dfm = pe.Node(fsl.ConvertWarp(relwarp=True, out_relwarp=True), + name='vsm2dfm') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, fixhdr, [('in_vsm', 'in_file'), + ('in_ref', 'in_file_hdr')]) + ,(inputnode, vsm, [('scaling', 'operand_value')]) + ,(fixhdr, vsm, [('out_file', 'in_file')]) + + ,(vsm, vsm2dfm, [('out_file', 'shift_in_file')]) + ,(inputnode, vsm2dfm, [('in_ref', 'reference'), + ('enc_dir', 'shift_direction')]) + ,(vsm2dfm, outputnode, [('out_file', 'out_warp')]) + ]) + return wf + + +def dwi_flirt(name='DWICoregistration', excl_nodiff=False, + flirt_param={}): + """ + Generates a workflow for linear registration of dwi volumes + """ + inputnode = pe.Node(niu.IdentityInterface(fields=['reference', + 'in_file', 'ref_mask', 'in_xfms', 'in_bval']), + name='inputnode') + + initmat = pe.Node(niu.Function(input_names=['in_bval', 'in_xfms', + 'excl_nodiff'], output_names=['init_xfms'], + function=_checkinitxfm), name='InitXforms') + initmat.inputs.excl_nodiff = excl_nodiff + dilate = pe.Node(fsl.maths.MathsCommand(nan2zeros=True, + args='-kernel sphere 5 -dilM'), name='MskDilate') + split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs') + pick_ref = pe.Node(niu.Select(), name='Pick_b0') + n4 = pe.Node(ants.N4BiasFieldCorrection(dimension=3), name='Bias') + enhb0 = pe.Node(niu.Function(input_names=['in_file', 'in_mask', + 'clip_limit'], output_names=['out_file'], + function=enhance), name='B0Equalize') + enhb0.inputs.clip_limit = 0.015 + enhdw = pe.MapNode(niu.Function(input_names=['in_file', 'in_mask'], + output_names=['out_file'], function=enhance), + name='DWEqualize', iterfield=['in_file']) + flirt = pe.MapNode(fsl.FLIRT(**flirt_param), name='CoRegistration', + iterfield=['in_file', 'in_matrix_file']) + thres = pe.MapNode(fsl.Threshold(thresh=0.0), iterfield=['in_file'], + name='RemoveNegative') + merge = pe.Node(fsl.Merge(dimension='t'), name='MergeDWIs') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', + 'out_xfms']), name='outputnode') + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, split, [('in_file', 'in_file')]) + ,(inputnode, dilate, [('ref_mask', 'in_file')]) + ,(inputnode, enhb0, [('ref_mask', 'in_mask')]) + ,(inputnode, initmat, [('in_xfms', 'in_xfms'), + ('in_bval', 'in_bval')]) + ,(inputnode, n4, [('reference', 'input_image'), + ('ref_mask', 'mask_image')]) + ,(dilate, flirt, [('out_file', 'ref_weight'), + ('out_file', 'in_weight')]) + ,(n4, enhb0, [('output_image', 'in_file')]) + ,(split, enhdw, [('out_files', 'in_file')]) + ,(dilate, enhdw, [('out_file', 'in_mask')]) + ,(enhb0, flirt, [('out_file', 'reference')]) + ,(enhdw, flirt, [('out_file', 'in_file')]) + ,(initmat, flirt, [('init_xfms', 'in_matrix_file')]) + ,(flirt, thres, [('out_file', 'in_file')]) + ,(thres, merge, [('out_file', 'in_files')]) + ,(merge, outputnode, [('merged_file', 'out_file')]) + ,(flirt, outputnode, [('out_matrix_file', 'out_xfms')]) + ]) + return wf + + +def apply_all_corrections(name='UnwarpArtifacts'): + """ + Combines two lists of linear transforms with the deformation field + map obtained typically after the SDC process. + Additionally, computes the corresponding bspline coefficients and + the map of determinants of the jacobian. + """ + + inputnode = pe.Node(niu.IdentityInterface(fields=['in_sdc', + 'in_hmc', 'in_ecc', 'in_dwi']), name='inputnode') + outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_warp', + 'out_coeff', 'out_jacobian']), name='outputnode') + warps = pe.MapNode(fsl.ConvertWarp(relwarp=True), + iterfield=['premat', 'postmat'], + name='ConvertWarp') + + selref = pe.Node(niu.Select(index=[0]), name='Reference') + + split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs') + unwarp = pe.MapNode(fsl.ApplyWarp(), iterfield=['in_file', 'field_file'], + name='UnwarpDWIs') + + coeffs = pe.MapNode(fsl.WarpUtils(out_format='spline'), + iterfield=['in_file'], name='CoeffComp') + jacobian = pe.MapNode(fsl.WarpUtils(write_jacobian=True), + iterfield=['in_file'], name='JacobianComp') + jacmult = pe.MapNode(fsl.MultiImageMaths(op_string='-mul %s'), + iterfield=['in_file', 'operand_files'], + name='ModulateDWIs') + + thres = pe.MapNode(fsl.Threshold(thresh=0.0), iterfield=['in_file'], + name='RemoveNegative') + merge = pe.Node(fsl.Merge(dimension='t'), name='MergeDWIs') + + wf = pe.Workflow(name=name) + wf.connect([ + (inputnode, warps, [('in_sdc', 'warp1'), + ('in_hmc', 'premat'), + ('in_ecc', 'postmat'), + ('in_dwi', 'reference')]) + ,(inputnode, split, [('in_dwi', 'in_file')]) + ,(split, selref, [('out_files', 'inlist')]) + ,(warps, unwarp, [('out_file', 'field_file')]) + ,(split, unwarp, [('out_files', 'in_file')]) + ,(selref, unwarp, [('out', 'ref_file')]) + ,(selref, coeffs, [('out', 'reference')]) + ,(warps, coeffs, [('out_file', 'in_file')]) + ,(selref, jacobian, [('out', 'reference')]) + ,(coeffs, jacobian, [('out_file', 'in_file')]) + ,(unwarp, jacmult, [('out_file', 'in_file')]) + ,(jacobian, jacmult, [('out_jacobian', 'operand_files')]) + ,(jacmult, thres, [('out_file', 'in_file')]) + ,(thres, merge, [('out_file', 'in_files')]) + + ,(warps, outputnode, [('out_file', 'out_warp')]) + ,(coeffs, outputnode, [('out_file', 'out_coeff')]) + ,(jacobian, outputnode, [('out_jacobian', 'out_jacobian')]) + ,(merge, outputnode, [('merged_file', 'out_file')]) + ]) + return wf + + +def extract_bval(in_dwi, in_bval, b=0, out_file=None): + """ + Writes an image containing only the volumes with b-value specified at + input + """ + import numpy as np + import nibabel as nb + import os.path as op + + if out_file is None: + fname, ext = op.splitext(op.basename(in_dwi)) + if ext == ".gz": + fname, ext2 = op.splitext(fname) + ext = ext2 + ext + out_file = op.abspath("%s_tsoi%s" % (fname, ext)) + + im = nb.load(in_dwi) + dwidata = im.get_data() + bvals = np.loadtxt(in_bval) + + if b == 'diff': + selection = np.where(bvals != 0) + elif b == 'nodiff': + selection = np.where(bvals == 0) + else: + selection = np.where(bvals == b) + + extdata = np.squeeze(dwidata.take(selection, axis=3)) + hdr = im.get_header().copy() + hdr.set_data_shape(extdata.shape) + nb.Nifti1Image(extdata, im.get_affine(), + hdr).to_filename(out_file) + return out_file + + +def remove_comp(in_file, in_bval, volid=0, out_file=None): + """ + Removes the volume ``volid`` from the 4D nifti file + """ + import numpy as np + import nibabel as nb + import os.path as op + + if out_file is None: + fname, ext = op.splitext(op.basename(in_file)) + if ext == ".gz": + fname, ext2 = op.splitext(fname) + ext = ext2 + ext + out_file = op.abspath("%s_extract%s" % (fname, ext)) + + im = nb.load(in_file) + data = im.get_data() + hdr = im.get_header().copy() + bval = np.loadtxt(in_bval) + + if volid == 0: + data = data[..., 1:] + bval = bval[1:] + elif volid == (data.shape[-1] - 1): + data = data[..., :-1] + bval = bval[:-1] + else: + data = np.concatenate((data[..., :volid], data[..., (volid + 1):]), + axis=3) + bval = np.hstack((bval[:volid], bval[(volid + 1):])) + hdr.set_data_shape(data.shape) + nb.Nifti1Image(data, im.get_affine(), hdr).to_filename(out_file) + + out_bval = op.abspath('bval_extract.txt') + np.savetxt(out_bval, bval) + return out_file, out_bval + + +def insert_mat(inlist, volid=0): + import numpy as np + import os.path as op + idfname = op.abspath('identity.mat') + out = inlist + np.savetxt(idfname, np.eye(4)) + out.insert(volid, idfname) + return out + + +def recompose_dwi(in_dwi, in_bval, in_corrected, out_file=None): + """ + Recompose back the dMRI data accordingly the b-values table after EC correction + """ + import numpy as np + import nibabel as nb + import os.path as op + + if out_file is None: + fname, ext = op.splitext(op.basename(in_dwi)) + if ext == ".gz": + fname, ext2 = op.splitext(fname) + ext = ext2 + ext + out_file = op.abspath("%s_eccorrect%s" % (fname, ext)) + + im = nb.load(in_dwi) + dwidata = im.get_data() + bvals = np.loadtxt(in_bval) + dwis = np.where(bvals != 0)[0].tolist() + + if len(dwis) != len(in_corrected): + raise RuntimeError('Length of DWIs in b-values table and after correction should match') + + for bindex, dwi in zip(dwis, in_corrected): + dwidata[..., bindex] = nb.load(dwi).get_data() + + nb.Nifti1Image(dwidata, im.get_affine(), + im.get_header()).to_filename(out_file) + return out_file + + +def recompose_xfm(in_bval, in_xfms): + """ + Insert identity transformation matrices in b0 volumes to build up a list + """ + import numpy as np + import os.path as op + + bvals = np.loadtxt(in_bval) + out_matrix = np.array([np.eye(4)] * len(bvals)) + xfms = iter([np.loadtxt(xfm) for xfm in in_xfms]) + out_files = [] + + for i, b in enumerate(bvals): + if b == 0.0: + mat = np.eye(4) + else: + mat = xfms.next() + + out_name = op.abspath('eccor_%04d.mat' % i) + out_files.append(out_name) + np.savetxt(out_name, mat) + + return out_files + + +def b0_average(in_dwi, in_bval, out_file=None): + """ + A function that averages the *b0* volumes from a DWI dataset. + + .. warning:: *b0* should be already registered (head motion artifact should + be corrected). + + """ + import numpy as np + import nibabel as nb + import os.path as op + + if out_file is None: + fname,ext = op.splitext(op.basename(in_dwi)) + if ext == ".gz": + fname,ext2 = op.splitext(fname) + ext = ext2 + ext + out_file = op.abspath("%s_avg_b0%s" % (fname, ext)) + + imgs = np.array(nb.four_to_three(nb.load(in_dwi))) + bval = np.loadtxt(in_bval) + b0s = [im.get_data().astype(np.float32) for im in imgs[np.where(bval==0)]] + b0 = np.average(np.array(b0s), axis=0) + + hdr = imgs[0].get_header().copy() + hdr.set_data_shape(b0.shape) + hdr.set_xyzt_units('mm') + hdr.set_data_dtype(np.float32) + nb.Nifti1Image(b0, imgs[0].get_affine(), hdr).to_filename(out_file) + return out_file + + +def rotate_bvecs(in_bvec, in_matrix): + """ + Rotates the input bvec file accordingly with a list of matrices. + + .. note:: the input affine matrix transforms points in the destination + image to their corresponding coordinates in the original image. + Therefore, this matrix should be inverted first, as we want to know + the target position of :math:`\\vec{r}`. + + """ + import os + import numpy as np + + name, fext = os.path.splitext(os.path.basename(in_bvec)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('%s_rotated.bvec' % name) + bvecs = np.loadtxt(in_bvec).T + new_bvecs = [] + + if len(bvecs) != len(in_matrix): + raise RuntimeError(('Number of b-vectors (%d) and rotation ' + 'matrices (%d) should match.') % (len(bvecs), + len(in_matrix))) + + for bvec, mat in zip(bvecs, in_matrix): + if np.all(bvec == 0.0): + new_bvecs.append(bvec) + else: + invrot = np.linalg.inv(np.loadtxt(mat))[:3, :3] + newbvec = invrot.dot(bvec) + new_bvecs.append((newbvec/np.linalg.norm(newbvec))) + + np.savetxt(out_file, np.array(new_bvecs).T, fmt='%0.15f') + return out_file + + +def eddy_rotate_bvecs(in_bvec, eddy_params): + """ + Rotates the input bvec file accordingly with a list of parameters sourced + from ``eddy``, as explained `here + `_. + """ + import os + import numpy as np + from math import sin, cos + + name, fext = os.path.splitext(os.path.basename(in_bvec)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('%s_rotated.bvec' % name) + bvecs = np.loadtxt(in_bvec).T + new_bvecs = [] + + params = np.loadtxt(eddy_params) + + if len(bvecs) != len(params): + raise RuntimeError(('Number of b-vectors and rotation ' + 'matrices should match.')) + + for bvec, row in zip(bvecs, params): + if np.all(bvec == 0.0): + new_bvecs.append(bvec) + else: + ax = row[3] + ay = row[4] + az = row[5] + + Rx = np.array([[1.0, 0.0, 0.0], + [0.0, cos(ax), -sin(ax)], + [0.0, sin(ax), cos(ax)]]) + Ry = np.array([[cos(ay), 0.0, sin(ay)], + [0.0, 1.0, 0.0], + [-sin(ay), 0.0, cos(ay)]]) + Rz = np.array([[cos(az), -sin(az), 0.0], + [sin(az), cos(az), 0.0], + [0.0, 0.0, 1.0]]) + R = Rx.dot(Ry).dot(Rz) + + invrot = np.linalg.inv(R) + newbvec = invrot.dot(bvec) + new_bvecs.append((newbvec/np.linalg.norm(newbvec))) + + np.savetxt(out_file, np.array(new_bvecs).T, fmt='%0.15f') + return out_file + + +def compute_readout(params): + """ + Computes readout time from epi params (see `eddy documentation + `_). + + .. warning:: ``params['echospacing']`` should be in *sec* units. + + + """ + epi_factor = 1.0 + acc_factor = 1.0 + try: + if params['epi_factor'] > 1: + epi_factor = float(params['epi_factor'] - 1) + except: + pass + try: + if params['acc_factor'] > 1: + acc_factor = 1.0 / params['acc_factor'] + except: + pass + return acc_factor * epi_factor * params['echospacing'] + + +def siemens2rads(in_file, out_file=None): + """ + Converts input phase difference map to rads + """ + import numpy as np + import nibabel as nb + import os.path as op + import math + + if out_file is None: + fname, fext = op.splitext(op.basename(in_file)) + if fext == '.gz': + fname, _ = op.splitext(fname) + out_file = op.abspath('./%s_rads.nii.gz' % fname) + + in_file = np.atleast_1d(in_file).tolist() + im = nb.load(in_file[0]) + data = im.get_data().astype(np.float32) + hdr = im.get_header().copy() + + if len(in_file) == 2: + data = nb.load(in_file[1]).get_data().astype(np.float32) - data + elif (data.ndim == 4) and (data.shape[-1] == 2): + data = np.squeeze(data[..., 1] - data[..., 0]) + hdr.set_data_shape(data.shape[:3]) + + imin = data.min() + imax = data.max() + data = (2.0 * math.pi * (data - imin)/(imax-imin)) - math.pi + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units('mm') + hdr['datatype'] = 16 + nb.Nifti1Image(data, im.get_affine(), hdr).to_filename(out_file) + return out_file + + +def rads2radsec(in_file, delta_te, out_file=None): + """ + Converts input phase difference map to rads + """ + import numpy as np + import nibabel as nb + import os.path as op + import math + + if out_file is None: + fname, fext = op.splitext(op.basename(in_file)) + if fext == '.gz': + fname, _ = op.splitext(fname) + out_file = op.abspath('./%s_radsec.nii.gz' % fname) + + im = nb.load(in_file) + data = im.get_data().astype(np.float32) * (1.0/delta_te) + nb.Nifti1Image(data, im.get_affine(), + im.get_header()).to_filename(out_file) + return out_file + + +def demean_image(in_file, in_mask=None, out_file=None): + """ + Demean image data inside mask + """ + import numpy as np + import nibabel as nb + import os.path as op + import math + + if out_file is None: + fname, fext = op.splitext(op.basename(in_file)) + if fext == '.gz': + fname, _ = op.splitext(fname) + out_file = op.abspath('./%s_demean.nii.gz' % fname) + + im = nb.load(in_file) + data = im.get_data().astype(np.float32) + msk = np.ones_like(data) + + if in_mask is not None: + msk = nb.load(in_mask).get_data().astype(np.float32) + msk[msk > 0] = 1.0 + msk[msk < 1] = 0.0 + + mean = np.median(data[msk == 1].reshape(-1)) + data[msk == 1] = data[msk == 1] - mean + nb.Nifti1Image(data, im.get_affine(), + im.get_header()).to_filename(out_file) + return out_file + + +def add_empty_vol(in_file, out_file=None): + """ + Adds an empty vol to the phase difference image + """ + import nibabel as nb + import os.path as op + import numpy as np + import math + + if out_file is None: + fname, fext = op.splitext(op.basename(in_file)) + if fext == '.gz': + fname, _ = op.splitext(fname) + out_file = op.abspath('./%s_4D.nii.gz' % fname) + + im = nb.load(in_file) + zim = nb.Nifti1Image(np.zeros_like(im.get_data()), im.get_affine(), + im.get_header()) + nb.funcs.concat_images([im, zim]).to_filename(out_file) + return out_file + + +def reorient_bvecs(in_dwi, old_dwi, in_bvec): + """ + Checks reorientations of ``in_dwi`` w.r.t. ``old_dwi`` and + reorients the in_bvec table accordingly. + """ + import os + import numpy as np + import nibabel as nb + + name, fext = os.path.splitext(os.path.basename(in_bvec)) + if fext == '.gz': + name, _ = os.path.splitext(name) + out_file = os.path.abspath('%s_reorient.bvec' % name) + bvecs = np.loadtxt(in_bvec).T + new_bvecs = [] + + N = nb.load(in_dwi).get_affine() + O = nb.load(old_dwi).get_affine() + RS = N.dot(np.linalg.inv(O))[:3, :3] + sc_idx = np.where((np.abs(RS) != 1) & (RS != 0)) + S = np.ones_like(RS) + S[sc_idx] = RS[sc_idx] + R = RS/S + + new_bvecs = [R.dot(b) for b in bvecs] + np.savetxt(out_file, np.array(new_bvecs).T, fmt='%0.15f') + return out_file + + +def copy_hdr(in_file, in_file_hdr, out_file=None): + import numpy as np + import nibabel as nb + import os.path as op + + if out_file is None: + fname, fext = op.splitext(op.basename(in_file)) + if fext == '.gz': + fname, _ = op.splitext(fname) + out_file = op.abspath('./%s_fixhdr.nii.gz' % fname) + + imref = nb.load(in_file_hdr) + hdr = imref.get_header().copy() + hdr.set_data_dtype(np.float32) + vsm = nb.load(in_file).get_data().astype(np.float32) + hdr.set_data_shape(vsm.shape) + hdr.set_xyzt_units('mm') + nii = nb.Nifti1Image(vsm, imref.get_affine(), hdr) + nii.to_filename(out_file) + return out_file + + +def enhance(in_file, clip_limit=0.010, in_mask=None, out_file=None): + import numpy as np + import nibabel as nb + import os.path as op + from skimage import exposure, img_as_int + + if out_file is None: + fname, fext = op.splitext(op.basename(in_file)) + if fext == '.gz': + fname, _ = op.splitext(fname) + out_file = op.abspath('./%s_enh.nii.gz' % fname) + + im = nb.load(in_file) + imdata = im.get_data() + imshape = im.get_shape() + + if in_mask is not None: + msk = nb.load(in_mask).get_data() + msk[msk > 0] = 1 + msk[msk < 1] = 0 + imdata = imdata * msk + + immin = imdata.min() + imdata = (imdata - immin).astype(np.uint16) + + adapted = exposure.equalize_adapthist(imdata.reshape(imshape[0], -1), + clip_limit=clip_limit) + + nb.Nifti1Image(adapted.reshape(imshape), im.get_affine(), + im.get_header()).to_filename(out_file) + + return out_file + + +def _checkinitxfm(in_bval, excl_nodiff, in_xfms=None): + from nipype.interfaces.base import isdefined + import numpy as np + import os.path as op + bvals = np.loadtxt(in_bval) + + gen_id = ((in_xfms is None) or + (not isdefined(in_xfms)) or + (len(in_xfms) != len(bvals))) + + init_xfms = [] + if excl_nodiff: + dws = np.where(bvals != 0)[0].tolist() + else: + dws = range(len(bvals)) + + if gen_id: + for i in dws: + xfm_file = op.abspath('init_%04d.mat' % i) + np.savetxt(xfm_file, np.eye(4)) + init_xfms.append(xfm_file) + else: + init_xfms = [in_xfms[i] for i in dws] + + return init_xfms diff --git a/nipype/workflows/dmri/setup.py b/nipype/workflows/dmri/setup.py index 599494d5b0..66d5fce20b 100644 --- a/nipype/workflows/dmri/setup.py +++ b/nipype/workflows/dmri/setup.py @@ -1,6 +1,8 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -def configuration(parent_package='',top_path=None): + + +def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration config = Configuration('dmri', parent_package, top_path)