Skip to content

WIP/including optional fieldmap unwarping in rsfmri example workflow #1081

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

WIP/including optional fieldmap unwarping in rsfmri example workflow #1081

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
158 changes: 157 additions & 1 deletion examples/rsfmri_vol_surface_preprocessing_nipy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down