diff --git a/CHANGES b/CHANGES index 73c13a2815..34c2dcfead 100644 --- a/CHANGES +++ b/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) ============================ diff --git a/nipype/interfaces/fsl/preprocess.py b/nipype/interfaces/fsl/preprocess.py index d777042952..c16649ac61 100644 --- a/nipype/interfaces/fsl/preprocess.py +++ b/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,6 +1186,17 @@ 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 @@ -1181,11 +1204,34 @@ def _list_outputs(self): 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 diff --git a/nipype/workflows/dmri/fsl/dti.py b/nipype/workflows/dmri/fsl/dti.py index 6fd5c2fe45..1cce5a5a52 100644 --- a/nipype/workflows/dmri/fsl/dti.py +++ b/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 \ No newline at end of file + 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