diff --git a/examples/rsfmri_vol_surface_preprocessing_nipy.py b/examples/rsfmri_vol_surface_preprocessing_nipy.py index f72de82d13..f99118c012 100644 --- a/examples/rsfmri_vol_surface_preprocessing_nipy.py +++ b/examples/rsfmri_vol_surface_preprocessing_nipy.py @@ -558,6 +558,103 @@ def create_reg_workflow(name='registration'): return register +def create_unwarp_workflow(name='unwarping'): + """Creates a workflow for unwarping the EPI data prior to + coregistration with Siemens fieldmaps. + + Parameters + ---------- + + :: + + name : name of workflow (default: 'unwarping') + + Inputs:: + + inputspec.realigned_file : realigned time series to unwarp + inputspec.mean_image : mean or median of time series + inputspec.fmap_mag : fieldmap magnitude image + inputspec.fmap_phase : fieldmap phase image + inputspec.delta_te : echo time difference of the fieldmap sequence [ms] + inputspec.pedir : phase encoding direction [x/y/z/x-/y-/z-] + inputspec.echospacing : effective EPI echo spacing (dwell time) [sec] + + Outputs:: + outputspec.unwarped_mean + outputspec.unwarped_file : unwarped time series + + """ + unwarp = Workflow(name=name) + + inputnode = Node(interface=IdentityInterface(fields=['realigned_file', + 'mean_image', + 'fmap_mag', + 'fmap_phase', + 'delta_te', + 'pedir', + 'echospacing']), + name='inputspec') + + outputnode = Node(interface=IdentityInterface(fields=['unwarped_file', + 'unwarped_mean']), + name='outputspec') + + # split first magnitude image from mag input file + split = Node(fsl.ExtractROI(t_min=0, + t_size=1), + name='split') + unwarp.connect(inputnode, 'fmap_mag', split, 'in_file') + + # skullstrip magnitude image and erode further to get rid of noisy edge voxels + bet = Node(fsl.BET(frac=0.5), + name='bet') + unwarp.connect(split,'roi_file', bet,'in_file') + + erode = Node(fsl.maths.ErodeImage(kernel_shape='sphere', + kernel_size=5), + name='erode') + unwarp.connect(bet,'out_file', erode, 'in_file') + + # prepare fieldmap + prep_fmap = Node(fsl.epi.PrepareFieldmap(), + name='prepare_fieldmap') + unwarp.connect([(erode, prep_fmap, [('out_file', 'in_magnitude')]), + (inputnode, prep_fmap, [('fmap_phase', 'in_phase'), + ('delta_te', 'delta_TE')]) + ]) + + # unmask fieldmap to avoid edge effects + fmap_mask = Node(fsl.maths.MathsCommand(args='-abs -bin -fillh'), + name='fmap_mask') + + unmask = Node(fsl.FUGUE(save_unmasked_fmap=True), + name='unmask') + + unwarp.connect([(prep_fmap, fmap_mask, [('out_fieldmap', 'in_file')]), + (fmap_mask, unmask, [('out_file', 'mask_file')]), + (prep_fmap, unmask,[('out_fieldmap','fmap_in_file')]), + (inputnode, unmask, [('pedir', 'unwarp_direction')]) + ]) + + # unwarp epi mean and time series with fieldmap + unwarp_mean = Node(fsl.FUGUE(save_shift=True),name='unwarp_mean') + unwarp_ts = Node(fsl.FUGUE(save_shift=True),name='unwarp_ts') + + unwarp.connect([(inputnode, unwarp_mean, [('mean_image', 'in_file'), + ('echospacing', 'dwell_time'), + ('pedir', 'unwarp_direction')]), + (unmask, unwarp_mean, [('fmap_out_file', 'fmap_in_file')]), + (fmap_mask, unwarp_mean, [('out_file','mask_file')]), + (inputnode, unwarp_ts, [('realigned_file', 'in_file'), + ('echospacing', 'dwell_time'), + ('pedir', 'unwarp_direction')]), + (unmask, unwarp_ts, [('fmap_out_file', 'fmap_in_file')]), + (fmap_mask, unwarp_ts, [('out_file','mask_file')]), + (unwarp_mean, outputnode, [('unwarped_file', 'unwarped_mean')]), + (unwarp_ts, outputnode, [('unwarped_file', 'unwarped_file')])]) + + return unwarp + """ Creates the main preprocessing workflow """ @@ -577,6 +674,12 @@ def create_workflow(files, subjects_dir=None, sink_directory=os.getcwd(), target_subject=['fsaverage3', 'fsaverage4'], + do_unwarp=False, + fmap_mag=None, + fmap_phase=None, + delta_te=None, + pedir=None, + echospacing=None, name='resting'): wf = Workflow(name=name) @@ -767,7 +870,33 @@ def merge_files(in1, in2): maskts = MapNode(fsl.ApplyMask(), iterfield=['in_file'], name='ts_masker') wf.connect(warpall, 'output_image', maskts, 'in_file') wf.connect(mask_target, 'out_file', maskts, 'mask_file') - + + + """ + Rewire for optional fieldmap unwarping + """ + if do_unwarp: + unwarping = create_unwarp_workflow(name='unwarping') + unwarping.inputs.inputspec.fmap_mag=fmap_mag + unwarping.inputs.inputspec.fmap_phase=fmap_phase + unwarping.inputs.inputspec.pedir=pedir + unwarping.inputs.inputspec.delta_te=delta_te + unwarping.inputs.inputspec.echospacing=echospacing + + wf.disconnect([(calc_median, registration, [('median_file', 'inputspec.mean_image')]), + (realign, art, [('out_file', 'realigned_files')]), + (realign, filter1, [('out_file', 'in_file'), + (('out_file', rename, '_filtermotart'), + 'out_res_name')])]) + wf.connect([(realign, unwarping, [('out_file', 'inputspec.realigned_file')]), + (calc_median, unwarping, [('median_file', 'inputspec.mean_image')]), + (unwarping, registration, [('outputspec.unwarped_mean', 'inputspec.mean_image')]), + (unwarping, art, [('outputspec.unwarped_file', 'realigned_files')]), + (unwarping, filter1, [('outputspec.unwarped_file', 'in_file'), + (('outputspec.unwarped_file', rename, '_filtermotart'), + 'out_res_name')])]) + + # map to surface # extract aparc+aseg ROIs # extract subcortical ROIs @@ -945,6 +1074,12 @@ def create_resting_workflow(args, name=None): if args.dicom_file: TR, slice_times, slice_thickness = get_info(args.dicom_file) slice_times = (np.array(slice_times)/1000.).tolist() + if args.do_unwarp: + fmap_mag=os.path.abspath(args.fmap_mag) + fmap_phase=os.path.abspath(args.fmap_phase) + else: + fmap_mag=None + fmap_phase=None if name is None: name = 'resting_' + args.subject_id @@ -961,6 +1096,12 @@ def create_resting_workflow(args, name=None): lowpass_freq=args.lowpass_freq, highpass_freq=args.highpass_freq, sink_directory=os.path.abspath(args.sink), + do_unwarp=args.do_unwarp, + fmap_mag=fmap_mag, + fmap_phase=fmap_phase, + delta_te=args.delta_te, + pedir=args.pedir, + echospacing=args.echospacing, name=name) wf = create_workflow(**kwargs) return wf @@ -1010,6 +1151,21 @@ def create_resting_workflow(args, name=None): help="Plugin to use") parser.add_argument("--plugin_args", dest="plugin_args", help="Plugin arguments") + parser.add_argument("--do_unwarp", dest="do_unwarp", + default=False, + help="Whether fieldmap unwarping should be used") + parser.add_argument("--fmap_mag", dest="fmap_mag", default=None, + help="fieldmap magnitude image") + parser.add_argument("--fmap_phase", dest="fmap_phase", default=None, + help="fieldmap phase image") + parser.add_argument("--delta_te", dest="delta_te", + default=None, type=float, + help="echo time difference of the fieldmap sequence [ms]") + parser.add_argument("--pedir", dest="pedir", default=None, + help="phase encoding direction [x/y/z/x-/y-/z-]") + parser.add_argument("--echospacing", dest="echospacing", + default=None, type=float, + help="echo time difference of the fieldmap sequence [ms]") args = parser.parse_args() wf = create_resting_workflow(args)