enh:new fmb_sdc workflow

Author: oesteban

nipype/workflows/dmri/preprocess/epi.py

@@ -76,6 +76,8 @@ def hmc_pipeline(name='motion_correct'):
     rot_bvec = pe.Node(niu.Function(input_names=['in_bvec', 'in_matrix'],
                            output_names=['out_file'], function=rotate_bvecs),
                            name='Rotate_Bvec')
+    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_bvec', 'out_xfms']),
@@ -91,7 +93,8 @@ def hmc_pipeline(name='motion_correct'):
         ,(inputnode,  rot_bvec,   [('in_bvec', 'in_bvec')])
         ,(flirt,      rot_bvec,   [('out_matrix_file', 'in_matrix')])
         ,(pick_ref,   flirt,      [('out', 'reference')])
-        ,(flirt,      merge,      [('out_file', 'in_files')])
+        ,(flirt,      thres,      [('out_file', 'in_file')])
+        ,(thres,      merge,      [('out_file', 'in_files')])
         ,(merge,      outputnode, [('merged_file', 'out_file')])
         ,(rot_bvec,   outputnode, [('out_file', 'out_bvec')])
         ,(flirt,      outputnode, [('out_matrix_file', 'out_xfms')])
@@ -207,9 +210,11 @@ def ecc_pipeline(name='eddy_correct'):
     return wf
 
 def sdc_fmb(name='fmb_correction',
-            fugue_params=dict(smooth3d=2.0, despike_2dfilter=True),
+            fugue_params=dict(smooth3d=2.0),
             bmap_params=dict(delta_te=2.46e-3),
-            epi_params=dict()):
+            epi_params=dict(echospacing=0.77e-3,
+                            acc_factor=3,
+                            enc_dir='y-')):
     """
     SDC stands for susceptibility distortion correction. FMB stands for fieldmap-based.
 
@@ -233,6 +238,8 @@ def sdc_fmb(name='fmb_correction',
     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']),
+                         name='outputnode')
 
 
     firstmag = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name='GetFirst')
@@ -265,37 +272,58 @@ def sdc_fmb(name='fmb_correction',
 
     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 = epi_params['echospacing'] / (1.0 * epi_params['acc_factor'])
+
+    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')
+
     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'),
+         (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')])
+        ,(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,       outputnode, [('merged_file', 'out_file')])
+        ,(vsm,         outputnode, [('shift_out_file', 'out_vsm')])
     ])
     return wf
 

nipype/workflows/dmri/preprocess/utils.py

@@ -5,7 +5,7 @@
 # @Author: oesteban
 # @Date:   2014-08-30 10:53:13
 # @Last Modified by:   oesteban
-# @Last Modified time: 2014-08-30 16:06:15
+# @Last Modified time: 2014-08-30 18:54:43
 
 def cleanup_edge_pipeline(name='Cleanup'):
     """
@@ -190,9 +190,9 @@ def siemens2rads(in_file, out_file=None):
     hdr = im.get_header().copy()
 
     if len(in_file) == 2:
-        data = data - nb.load(in_file[1]).get_data().astype(np.float32)
+        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[...,0] - data[...,1])
+        data = np.squeeze(data[...,1] - data[...,0])
         hdr.set_data_shape(data.shape[:3])
 
     imin = data.min()
@@ -253,3 +253,25 @@ def demean_image(in_file, in_mask=None, out_file=None):
     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