ENH: new artifact correction workflow on fsl

Uses TOPUP and Eddy for correction. Tested on the new dmri_preprocess.py
example with the fsl course data.

Author: oesteban

examples/dmri_preprocessing.py

@@ -3,7 +3,7 @@
 # @Author: oesteban
 # @Date:   2014-08-31 20:32:22
 # @Last Modified by:   oesteban
-# @Last Modified time: 2014-09-02 09:41:46
+# @Last Modified time: 2014-09-02 13:12:12
 """
 ===================
 dMRI: Preprocessing
@@ -15,17 +15,18 @@
 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
+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 util     # utility
+import nipype.interfaces.utility as niu      # utility
 import nipype.algorithms.misc as misc
 
 import nipype.pipeline.engine as pe          # pypeline engine
@@ -34,7 +35,6 @@
 from nipype.interfaces import ants
 
 
-
 """
 Load specific nipype's workflows for preprocessing of dMRI data:
 :class:`nipype.workflows.dmri.preprocess.epi.all_peb_pipeline`,
@@ -43,8 +43,7 @@
 that is *A>>>P* or *-y* (in RAS systems).
 """
 
-from nipype.workflows.dmri.preprocess.epi import all_peb_pipeline
-
+from nipype.workflows.dmri.preprocess.epi import all_fsl_pipeline, remove_bias
 
 """
 Map field names into individual subject runs
@@ -55,7 +54,7 @@
             bvals=[['subject_id', 'bvals']],
             dwi_rev=[['subject_id', 'nodif_PA']])
 
-infosource = pe.Node(interface=util.IdentityInterface(fields=['subject_id']),
+infosource = pe.Node(interface=niu.IdentityInterface(fields=['subject_id']),
                      name="infosource")
 
 # Set the subject 1 identifier in subject_list,
@@ -102,7 +101,7 @@
 actual processing functions
 """
 
-inputnode = pe.Node(util.IdentityInterface(fields=["dwi", "bvecs", "bvals",
+inputnode = pe.Node(niu.IdentityInterface(fields=["dwi", "bvecs", "bvals",
                     "dwi_rev"]), name="inputnode")
 
 
@@ -117,7 +116,7 @@
 Artifacts correction
 --------------------
 
-We will use the combination of ```topup``` and ```eddy``` as suggested by FSL.
+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.
@@ -128,25 +127,23 @@
 
 """
 
-epi_AP = {'echospacing': 66.5e-3, 'acc_factor': 1, 'enc_dir': 'y-'}
-epi_PA = {'echospacing': 66.5e-3, 'acc_factor': 1, 'enc_dir': 'y'}
-prep = all_peb_pipeline(epi_params=epi_AP, altepi_params=epi_PA)
+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 estimate a single multiplicative bias field from
-the *b0* image, as suggested in [Jeurissen2014]_.
+Finally, we set up a node to correct for a single multiplicative bias field
+from computed on the *b0* image, as suggested in [Jeurissen2014]_.
 
 """
 
-n4 = pe.Node(ants.N4BiasFieldCorrection(dimension=3), name='Bias_b0')
-split = pe.Node(fsl.Split(dimension='t'), name='SplitDWIs')
-merge = pe.Node(niu.Function(input_names=['in_dwi', 'in_bval', 'in_corrected'],
-                output_names=['out_file'], function=recompose_dwi), name='MergeDWIs')
+bias = remove_bias()
+
 
 """
 Connect nodes in workflow
@@ -165,6 +162,9 @@
                                      ('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')])
 ])
 
 
@@ -176,12 +176,3 @@
     wf.run()
     wf.write_graph()
 
-
-"""
-
-.. admonition:: References
-
-  .. [Jeurissen2014] Jeurissen B. et al., `Multi-tissue constrained spherical deconvolution
-    for improved analysis of multi-shell diffusion MRI data
-    <http://dx.doi.org/10.1016/j.neuroimage.2014.07.061>`_. NeuroImage (2014).
-    doi: 10.1016/j.neuroimage.2014.07.061

nipype/workflows/dmri/preprocess/__init__.py

@@ -1,2 +1,3 @@
-from epi import (all_fmb_pipeline, all_peb_pipeline,
-                 hmc_pipeline, ecc_pipeline, sdc_fmb, sdc_peb)
+from epi import (all_fmb_pipeline, all_peb_pipeline, all_fsl_pipeline,
+                 hmc_pipeline, ecc_pipeline, sdc_fmb, sdc_peb,
+                 remove_bias)

nipype/workflows/dmri/preprocess/epi.py

@@ -70,10 +70,12 @@ def all_fmb_pipeline(name='hmc_sdc_ecc'):
 def all_peb_pipeline(name='hmc_sdc_ecc',
                      epi_params=dict(echospacing=0.77e-3,
                                      acc_factor=3,
-                                     enc_dir='y-'),
+                                     enc_dir='y-',
+                                     epi_factor=1),
                      altepi_params=dict(echospacing=0.77e-3,
                                         acc_factor=3,
-                                        enc_dir='y')):
+                                        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
@@ -96,6 +98,10 @@ def all_peb_pipeline(name='hmc_sdc_ecc',
     sdc = sdc_peb(epi_params=epi_params, altepi_params=altepi_params)
     ecc = ecc_pipeline()
 
+    rot_bvec = pe.Node(niu.Function(input_names=['in_bvec', 'eddy_params'],
+                       output_names=['out_file'], function=eddy_rotate_bvecs),
+                       name='Rotate_Bvec')
+
     regrid = pe.Node(fs.MRIConvert(vox_size=(2.0, 2.0, 2.0), out_orientation='RAS'),
                      name='Reslice')
 
@@ -120,6 +126,76 @@ def all_peb_pipeline(name='hmc_sdc_ecc',
         ,(regrid,      avg_b0_1,   [('out_file', 'in_dwi')])
         ,(inputnode,   avg_b0_1,   [('in_bval', 'in_bval')])
         ,(avg_b0_1,    bet_dwi1,   [('out_file','in_file')])
+        ,(inputnode,   rot_bvec,   [('in_bvec', 'in_bvec')])
+        ,(ecc,         rot_bvec,   [('out_parameter', 'eddy_params')])
+        ,(bet_dwi1,    outputnode, [('mask_file', 'out_mask')])
+        ,(rot_bvec,    outputnode, [('out_file', 'out_bvec')])
+    ])
+    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``.
+    """
+
+    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')
+
+    regrid = pe.Node(fs.MRIConvert(vox_size=(2.0, 2.0, 2.0), out_orientation='RAS'),
+                     name='Reslice')
+    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('dMRI_Artifacts_FSL')
+    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')])
+        ,(ecc,         regrid,     [('out_corrected', 'in_file')])
+        ,(regrid,      outputnode, [('out_file', 'out_file')])
+        ,(regrid,      avg_b0_1,   [('out_file', 'in_dwi')])
+        ,(inputnode,   avg_b0_1,   [('in_bval', 'in_bval')])
+        ,(avg_b0_1,    bet_dwi1,   [('out_file','in_file')])
         ,(bet_dwi1,    outputnode, [('mask_file', 'out_mask')])
     ])
     return wf
@@ -218,6 +294,7 @@ def hmc_pipeline(name='motion_correct'):
     ])
     return wf
 
+
 def ecc_pipeline(name='eddy_correct'):
     """
     ECC stands for Eddy currents correction.
@@ -326,6 +403,7 @@ def ecc_pipeline(name='eddy_correct'):
     ])
     return wf
 
+
 def sdc_fmb(name='fmb_correction',
             fugue_params=dict(smooth3d=2.0),
             bmap_params=dict(delta_te=2.46e-3),
@@ -462,10 +540,12 @@ def sdc_fmb(name='fmb_correction',
 def sdc_peb(name='peb_correction',
             epi_params=dict(echospacing=0.77e-3,
                             acc_factor=3,
-                            enc_dir='y-'),
+                            enc_dir='y-',
+                            epi_factor=1),
             altepi_params=dict(echospacing=0.77e-3,
                                acc_factor=3,
-                               enc_dir='y')):
+                               enc_dir='y',
+                               epi_factor=1)):
     """
     SDC stands for susceptibility distortion correction. PEB stands for phase-encoding-based.
 
@@ -517,25 +597,72 @@ def sdc_peb(name='peb_correction',
 
     topup = pe.Node(fsl.TOPUP(), name='topup')
     topup.inputs.encoding_direction = [epi_params['enc_dir'], altepi_params['enc_dir']]
-    topup.inputs.readout_times = [epi_params['echospacing']/(1.0*epi_params['acc_factor']),
-                                  altepi_params['echospacing']/(1.0*altepi_params['acc_factor'])]
+    topup.inputs.readout_times = [compute_readout(epi_params),
+                                  compute_readout(altepi_params)]
+
     unwarp = pe.Node(fsl.ApplyTOPUP(in_index=[1], method='jac'), name='unwarp')
 
     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')])
+         (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_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')])
+    ])
+    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
+        <http://dx.doi.org/10.1016/j.neuroimage.2014.07.061>`_. NeuroImage (2014).
+        doi: 10.1016/j.neuroimage.2014.07.061
+
+    """
+    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')])
     ])
     return wf
 
@@ -546,12 +673,13 @@ def _checkrnum(ref_num):
         return 0
     return ref_num
 
+
 def _checkinitxfm(in_bval, in_xfms=None):
     from nipype.interfaces.base import isdefined
     import numpy as np
     import os.path as op
     bvals = np.loadtxt(in_bval)
-    non_b0 = np.where(bvals!=0)[0].tolist()
+    non_b0 = np.where(bvals != 0)[0].tolist()
 
     init_xfms = []
     if (in_xfms is None) or (not isdefined(in_xfms)) or (len(in_xfms)!=len(bvals)):