New EPI dewarping workflow, in a similar manner to eddy_current correction

CHANGES

@@ -3,6 +3,9 @@ Next release
 * ENH: New interfaces: nipy.Trim
 
 * ENH: Allow control over terminal output for commandline interfaces
+* ENH: New workflows: susceptibility correction for EPI imaging based on epidewarp.fsl
+       updated to support both dMRI and fMRI, and new parameters.
+* ENH: Minor improvements to FSL's FUGUE interface
 
 Release 0.7.0 (Dec 18, 2012)
 ============================

nipype/interfaces/fsl/preprocess.py

@@ -1092,6 +1092,12 @@ class FUGUEInputSpec(FSLCommandInputSpec):
                    desc='filename of input volume')
     unwarped_file = File(argstr='--unwarp=%s', genfile=True,
                          desc='apply unwarping and save as filename', hash_files=False)
+
+    save_warped = traits.Bool( desc='apply forward warp and save' )
+
+    warped_file = File(argstr='--warp=%s', genfile=True,
+                         desc='apply forward warp and save as filename', hash_files=False)
+
     phasemap_file = File(exists=True, argstr='--phasemap=%s',
                          desc='filename for input phase image')
     dwell_to_asym_ratio = traits.Float(argstr='--dwelltoasym=%.10f',
@@ -1104,8 +1110,12 @@ class FUGUEInputSpec(FSLCommandInputSpec):
                      desc='filename for saving fieldmap (rad/s)', hash_files=False)
     fmap_in_file = File(exists=True, argstr='--loadfmap=%s',
                         desc='filename for loading fieldmap (rad/s)')
-    shift_out_file = File(argstr='--saveshift=%s',
-                          desc='filename for saving pixel shift volume', hash_files=False)
+
+    save_shift = traits.Bool( desc='output pixel shift volume' )
+
+    shift_out_file = traits.File( argstr='--saveshift=%s', genfile=True,
+                           desc='filename for saving pixel shift volume', hash_files=False)
+
     shift_in_file = File(exists=True, argstr='--loadshift=%s',
                          desc='filename for reading pixel shift volume')
     median_2dfilter = traits.Bool(argstr='--median',
@@ -1154,6 +1164,8 @@ class FUGUEInputSpec(FSLCommandInputSpec):
 
 class FUGUEOutputSpec(TraitedSpec):
     unwarped_file = File(exists=True, desc='unwarped file')
+    shift_out_file = File( desc='voxel shift map file' )
+    warped_file = File( desc='warped file' )
 
 
 class FUGUE(FSLCommand):
@@ -1174,18 +1186,52 @@ def __init__(self, **kwargs):
         super(FUGUE, self).__init__(**kwargs)
         warn('This interface has not been fully tested. Please report any failures.')
 
+    def _parse_inputs(self, skip=None):
+        skip = []
+        if not isdefined( self.inputs.save_shift ) or not self.inputs.save_shift:
+            skip.append('shift_out_file')
+        if not isdefined( self.inputs.save_warped ) or not self.inputs.save_warped:
+            skip.append('warped_file')
+        else:
+            skip.append('unwarped_file')
+
+        return super(FUGUE,self)._parse_inputs(skip=skip)
+
     def _list_outputs(self):
         outputs = self._outputs().get()
         out_file = self.inputs.unwarped_file
         if not isdefined(out_file):
             out_file = self._gen_fname(self.inputs.in_file,
                                       suffix='_unwarped')
         outputs['unwarped_file'] = os.path.abspath(out_file)
+
+        if isdefined( self.inputs.save_shift ) and self.inputs.save_shift:
+            shift_out = self.inputs.shift_out_file
+            if not isdefined(shift_out):
+                shift_out = self._gen_fname( self.inputs.in_file, suffix='_shift' )
+
+            outputs['shift_out_file'] = os.path.abspath( shift_out )
+
+        if isdefined( self.inputs.save_warped ) and self.inputs.save_warped:
+            warped_out = self.inputs.warped_file
+            if not isdefined(warped_out):
+                warped_out = self._gen_fname( self.inputs.in_file, suffix='_fwdwarp' )
+
+            outputs['warped_file'] = os.path.abspath( warped_out )
+            del outputs['unwarped_file']
+
         return outputs
 
     def _gen_filename(self, name):
         if name == 'unwarped_file':
             return self._list_outputs()['unwarped_file']
+
+        if name == 'warped_file':
+            return self._list_outputs()['warped_file']
+
+        if name == 'shift_out_file':
+            return self._list_outputs()['shift_out_file']
+
         return None
 
 

nipype/workflows/dmri/fsl/dti.py

@@ -1,7 +1,9 @@
+# coding: utf-8
+
 import nipype.pipeline.engine as pe
 import nipype.interfaces.utility as util
 import nipype.interfaces.fsl as fsl
-
+import os
 
 def transpose(samples_over_fibres):
     import numpy as np
@@ -214,4 +216,264 @@ def create_eddy_correct_pipeline(name="eddy_correct"):
 
     pipeline.connect([(merge, outputnode, [("merged_file", "eddy_corrected")])])
 
-    return pipeline
+    return pipeline
+
+def create_epi_correct_pipeline(name="epi_correct", register_to_ref=False ):
+    """ Replaces the epidewarp.fsl script (http://www.nmr.mgh.harvard.edu/~greve/fbirn/b0/epidewarp.fsl)
+        for EPI distortion correction with the fieldmap information in
+        dMRI and fMRI (Jezzard et al., MRM 1995 ) using FSL's FUGUE.
+
+    Example
+    -------
+
+    >>> nipype_epicorrect = create_epi_correct_pipeline("nipype_epicorrect")
+    >>> nipype_epicorrect.inputs.inputnode.in_file = 'epi_data.nii'
+    >>> nipype_epicorrect.inputs.inputnode.fieldmap_mag = 'magnitude.nii'
+    >>> nipype_epicorrect.inputs.inputnode.fieldmap_pha = 'phase.nii'
+    >>> nipype_epicorrect.inputs.inputnode.te_diff = 2.46
+    >>> nipype_epicorrect.inputs.inputnode.epi_echospacing = 0.51
+    >>> nipype_epicorrect.inputs.inputnode.epi_rev_encoding = False
+    >>> nipype_epicorrect.inputs.inputnode.ref_volume = 0
+    >>> nipype_epicorrect.inputs.inputnode.epi_parallel = True
+    >>> nipype_epicorrect.run() # doctest: +SKIP
+
+    Inputs::
+
+        inputnode.in_file - The volume acquired with EPI sequence
+        inputnode.fieldmap_mag - The magnitude of the fieldmap
+        inputnode.fieldmap_pha - The phase difference of the fieldmap
+        inputnode.te_diff - Time difference between TE in ms.
+        inputnode.epi_echospacing - The echo spacing (aka dwell time) in the EPI sequence
+        inputnode.epi_ph_encoding_dir - The phase encoding direction in EPI acquisition (default y)
+        inputnode.epi_rev_encoding - True if it is acquired with reverse encoding
+        inputnode.epi_parallel - True if EPI was acquired in a parallel imaging scheme
+        inputnode.vsm_sigma - Sigma value of the gaussian smoothing filter applied to the vsm (voxel shift map)
+        inputnode.ref_volume - The reference volume (B=0 in dMRI or a central frame in fMRI)
+
+
+    Outputs::
+
+        outputnode.epi_corrected
+    """
+
+    inputnode = pe.Node(interface = util.IdentityInterface(fields=["in_file", 
+                                                                   "fieldmap_mag", 
+                                                                   "fieldmap_pha", 
+                                                                   "te_diff", 
+                                                                   "epi_echospacing",
+                                                                   "epi_ph_encoding_dir", 
+                                                                   "epi_rev_encoding", 
+                                                                   "epi_parallel",
+                                                                   "vsm_sigma",
+                                                                   "ref_volume",
+                                                                   "unwarp_direction"
+                                                                   ]), name="inputnode")
+
+    pipeline = pe.Workflow(name=name)
+
+
+    matrix_file = os.path.abspath( './flirt.txt' )
+
+    # Keep first frame from magnitude
+    select_mag = pe.Node( interface=fsl.utils.ExtractROI(t_size=1,t_min=0), name="select_magnitude" )
+
+    # mask_brain
+    mask_mag = pe.Node( interface=fsl.BET(mask=True),name='mask_magnitude' )
+    mask_mag_dil = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_dilate_mask), name='mask_dilate' )
+
+    # Compute dwell time
+    dwell_time = pe.Node( interface=util.Function( input_names=["dwell_time","is_parallel","is_reverse_encoding"], output_names=["dwell_time"], function=_compute_dwelltime), name='dwell_time')
+
+    # Normalize phase diff to be [-pi, pi)
+    norm_pha = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_prepare_phasediff ), name='normalize_phasediff')
+    # Execute FSL PRELUDE: prelude -p %s -a %s -o %s -f -v -m %s
+    prelude = pe.Node( interface=fsl.PRELUDE(process3d=True), name='phase_unwrap' )
+    fill_phase = pe.Node( interface=util.Function( input_names=["in_file"], output_names=["out_file"], function=_fill_phase ), name='fill_phasediff' )
+
+    # to assure that vsm is same dimension as mag. The input only affects the output dimension.
+    # The content of the input has no effect on the vsm. The de-warped mag volume is
+    # meaningless and will be thrown away
+    # fugue -i %s -u %s -p %s --dwell=%s --asym=%s --mask=%s --saveshift=%s % ( mag_name, magdw_name, ph_name, esp, tediff, mask_name, vsmmag_name)
+    vsm = pe.Node( interface=fsl.FUGUE(save_shift=True), name="generate_vsm" )
+    vsm_mean = pe.Node( interface=util.Function( input_names=["in_file","mask_file","in_unwarped"], output_names=["out_file"], function=_vsm_remove_mean), name="vsm_mean_shift" )
+
+    # fugue_epi
+    dwi_split = pe.Node( interface=util.Function( input_names=['in_file'], output_names=['out_files'], function=_split_dwi ), name='dwi_split' )
+    # 'fugue -i %s -u %s --loadshift=%s --mask=%s' % ( vol_name, out_vol_name, vsm_name, mask_name )
+    dwi_applyxfm = pe.MapNode( interface=fsl.FUGUE(icorr=True,save_shift=False), iterfield=['in_file'], name='dwi_fugue' ) 
+    # Merge back all volumes
+    dwi_merge = pe.Node( interface=fsl.utils.Merge(dimension='t' ), name='dwi_merge')
+
+    outputnode = pe.Node(interface = util.IdentityInterface(fields=["epi_corrected"]),
+                        name="outputnode")
+
+
+    pipeline.connect([ 
+                     (inputnode,    dwell_time, [('epi_echospacing','dwell_time'), ('epi_parallel','is_parallel'),('epi_rev_encoding','is_reverse_encoding') ])
+                    ,(inputnode,    select_mag, [('fieldmap_mag','in_file')] )
+                    ,(inputnode,      norm_pha, [('fieldmap_pha','in_file')] )
+                    ,(select_mag,     mask_mag, [('roi_file', 'in_file')] )
+                    ,(mask_mag,   mask_mag_dil, [('mask_file','in_file')] )
+                    ,(select_mag,      prelude, [('roi_file', 'magnitude_file') ])
+                    ,(norm_pha,        prelude, [('out_file', 'phase_file' ) ])
+                    ,(mask_mag_dil,    prelude, [('out_file', 'mask_file' ) ])
+                    ,(prelude,      fill_phase, [('unwrapped_phase_file','in_file')] )
+                    ,(inputnode,           vsm, [('fieldmap_mag', 'in_file')])
+                    ,(fill_phase,          vsm, [('out_file','phasemap_file')]) #, (('out_file',_genvsmpath),'shift_out_file') ] )
+                    ,(inputnode,           vsm, [(('te_diff',_ms2sec),'asym_se_time'),('vsm_sigma','smooth2d')])
+                    ,(dwell_time,          vsm, [(('dwell_time',_ms2sec),'dwell_time') ])
+                    ,(mask_mag_dil,        vsm, [('out_file','mask_file')] )
+                    ,(mask_mag_dil,   vsm_mean, [('out_file','mask_file')] )
+                    ,(vsm,            vsm_mean, [('unwarped_file','in_unwarped'),('shift_out_file','in_file')])
+                    ,(inputnode,     dwi_split, [('in_file','in_file')])
+                    ,(dwi_split,  dwi_applyxfm, [('out_files','in_file')])
+                    ,(dwi_applyxfm,  dwi_merge, [('unwarped_file','in_files')])
+                    ,(dwi_merge,    outputnode, [('merged_file','epi_corrected')])
+                    ])
+
+    if register_to_ref:
+        # register to ref volume
+        """ Register magfw to example epi. There are some parameters here that may need to be tweaked. Should probably strip the mag
+            Pre-condition: forward warp the mag in order to reg with func. What does mask do here?
+        """
+        # Select reference volume from EPI (B0 in dMRI and a middle frame in fMRI)
+        select_epi = pe.Node( interface=fsl.utils.ExtractROI(t_size=1), name="select_epi" )
+
+
+        # fugue -i %s -w %s --loadshift=%s --mask=%s % ( mag_name, magfw_name, vsmmag_name, mask_name ), log ) # Forward Map
+        vsm_fwd = pe.Node( interface=fsl.FUGUE(save_warped=True), name="vsm_fwd")
+        vsm_reg = pe.Node( interface=fsl.FLIRT( bins=256, cost='corratio', dof=6, interp='trilinear',  searchr_x=[-10,10], searchr_y=[-10,10], searchr_z=[-10,10] ), name="vsm_registration" )
+        # 'flirt -in %s -ref %s -out %s -init %s -applyxfm' % ( vsmmag_name, ref_epi, vsmmag_name, magfw_mat_out )
+        vsm_applyxfm = pe.Node( interface=fsl.ApplyXfm(), name='vsm_apply_xfm')
+        # 'flirt -in %s -ref %s -out %s -init %s -applyxfm' % ( mask_name, ref_epi, mask_name, magfw_mat_out )
+        msk_applyxfm = pe.Node( interface=fsl.ApplyXfm(), name='msk_apply_xfm')
+
+        pipeline.connect([
+                     (inputnode,    select_epi, [('in_file','in_file'), ('ref_volume','t_min')] )
+                    ,(select_epi,        vsm_reg, [('roi_file','reference')])
+                    ,(vsm,               vsm_fwd, [('shift_out_file','shift_in_file')])
+                    ,(mask_mag_dil,      vsm_fwd, [('out_file','mask_file')] )
+                    ,(inputnode,         vsm_fwd, [('fieldmap_mag', 'in_file')])
+                    ,(vsm_fwd,           vsm_reg, [('warped_file','in_file')])
+                    ,(vsm_reg,      msk_applyxfm, [('out_matrix_file','in_matrix_file')])
+                    ,(select_epi,   msk_applyxfm, [('roi_file','reference')])
+                    ,(mask_mag_dil, msk_applyxfm, [('out_file','in_file') ])
+                    ,(vsm_reg,      vsm_applyxfm, [('out_matrix_file','in_matrix_file')])
+                    ,(select_epi,   vsm_applyxfm, [('roi_file','reference')])
+                    ,(vsm_mean,     vsm_applyxfm, [('out_file','in_file')]) 
+                    ,(msk_applyxfm, dwi_applyxfm, [('out_file','mask_file')])
+                    ,(vsm_applyxfm, dwi_applyxfm, [('out_file','shift_in_file')])
+                    ])
+    else:
+        pipeline.connect([
+                     (mask_mag_dil, dwi_applyxfm, [('out_file','mask_file')])
+                    ,(vsm_mean,     dwi_applyxfm, [('out_file','shift_in_file')])
+                    ])
+
+
+    return pipeline
+
+def _compute_dwelltime( dwell_time=0.68, is_parallel=False, is_reverse_encoding=False ):
+    if is_parallel:
+        dwell_time*=0.5
+
+    if is_reverse_encoding:
+        dwell_time*=-1.0
+
+    return dwell_time
+
+def _prepare_phasediff( in_file ):
+    import nibabel as nib
+    import os
+    import numpy as np
+    img = nib.load( in_file )
+    max_diff = np.max( img.get_data().reshape(-1) )
+    min_diff = np.min( img.get_data().reshape(-1) )
+    A = ( 2.0 * np.pi )/( max_diff-min_diff )
+    B = np.pi - ( A * max_diff )
+    diff_norm = img.get_data() * A + B
+
+    name,fext = os.path.splitext( os.path.basename(in_file) )
+    if fext == '.gz': name,_ = os.path.splitext(name)
+    out_file = os.path.abspath( './%s_2pi.nii.gz' % name )
+    nib.save( nib.Nifti1Image( diff_norm, img.get_affine(), img.get_header() ), out_file )
+    return out_file
+
+def _dilate_mask ( in_file, iterations=4 ):
+    import nibabel as nib
+    import scipy.ndimage as ndimage
+    import os
+    img = nib.load( in_file )
+    img._data = ndimage.binary_dilation( img.get_data(), iterations=iterations )
+
+    name,fext = os.path.splitext( os.path.basename(in_file) )
+    if fext == '.gz': name,_ = os.path.splitext(name)
+    out_file = os.path.abspath( './%s_dil.nii.gz' % name )
+    nib.save( img, out_file )
+    return out_file
+
+
+def _fill_phase( in_file ):
+    import nibabel as nib
+    import os
+    import numpy as np
+    img = nib.load( in_file )
+    dumb_img = nib.Nifti1Image( np.zeros( img.get_shape() ), img.get_affine(), img.get_header() )
+    out_nii = nib.funcs.concat_images( ( img, dumb_img ) )
+    name,fext = os.path.splitext( os.path.basename(in_file) )
+    if fext == '.gz': name,_ = os.path.splitext(name)
+    out_file = os.path.abspath( './%s_phase_unwrapped.nii.gz' % name ) 
+    nib.save( out_nii, out_file )
+    return out_file
+
+def _vsm_remove_mean( in_file, mask_file, in_unwarped ):
+    import nibabel as nib
+    import os
+    import numpy as np
+    import numpy.ma as ma
+    img = nib.load( in_file )
+    msk = nib.load( mask_file ).get_data()
+    img_data = img.get_data()
+    img_data[ msk==0 ] = 0
+    vsmmag_masked = ma.masked_values( img_data.reshape(-1), 0.0 )
+    vsmmag_masked = vsmmag_masked - vsmmag_masked.mean()
+    img._data = vsmmag_masked.reshape( img.get_shape() )
+    name,fext = os.path.splitext( os.path.basename(in_file) )
+    if fext == '.gz': name,_ = os.path.splitext(name)
+    out_file = os.path.abspath( './%s_demeaned.nii.gz' % name ) 
+    nib.save( img, out_file )
+    return out_file
+
+#def _generate_fwdmap( in_file, mask_file, shift_in_file, in_unwarped):
+#    import os
+#    out_file = os.path.abspath( './fwdmap.nii.gz' )
+#    cmd = 'fugue -i %s -w %s --loadshift=%s --mask=%s' % ( in_file, out_file, shift_in_file, mask_file )
+#    print "	 Running:",cmd
+#    os.system(cmd)
+#    return out_file
+
+def _genvsmpath( in_file ):
+    import os
+    name,fext = os.path.splitext( os.path.basename(in_file) )
+    if fext == '.gz': name,_ = os.path.splitext(name)
+    out_file = os.path.abspath( os.path.join( os.path.dirname(in_file), '%s_vsm.nii.gz' % name ))
+    return out_file
+
+def _ms2sec( val ):
+    return val*1e-3;
+
+def _yield( in_file, in_block ):
+    return in_file
+
+def _split_dwi( in_file ):
+    import nibabel as nib
+    import os
+    out_files = []
+    frames = nib.funcs.four_to_three( nib.load( in_file ) )
+    name,fext = os.path.splitext( os.path.basename(in_file) )
+    if fext == '.gz': name,_ = os.path.splitext(name)
+    for i,frame in enumerate(frames):
+	out_file = os.path.abspath( './%s_%03d.nii.gz' % (name,i) )
+        nib.save( frame, out_file )
+        out_files.append( out_file )
+    return out_files