Revision of dMRI artifact correction workflows

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.

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()

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

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

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.
 

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.

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)

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.')

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
+