In [1]:
from __future__ import print_function
from __future__ import division
from builtins import str
from builtins import range
import os # system functions
import nipype.interfaces.io as nio # Data i/o
import nipype.interfaces.fsl as fsl # fsl
import nipype.interfaces.utility as util # utility
import nipype.pipeline.engine as pe # pypeline engine
import nipype.algorithms.modelgen as model # model generation
import nipype.algorithms.rapidart as ra # artifact detection
In [ ]:
"""
Preliminaries
-------------
Setup any package specific configuration. The output file format for FSL
routines is being set to compressed NIFTI.
"""
fsl.FSLCommand.set_default_output_type('NIFTI_GZ')
In [ ]:
"""Setup preprocessing workflow
----------------------------
This is a generic fsl feat preprocessing workflow encompassing skull stripping,
motion correction and smoothing operations.
"""
preproc = pe.Workflow(name='preproc')
In [ ]:
inputnode = pe.Node(interface=util.IdentityInterface(fields=['func',
'struct', ]),
name='inputspec')
In [ ]:
"""
Convert functional images to float representation. Since there can be more than
one functional run we use a MapNode to convert each run.
"""
img2float = pe.MapNode(interface=fsl.ImageMaths(out_data_type='float',
op_string='',
suffix='_dtype'),
iterfield=['in_file'],
name='img2float')
preproc.connect(inputnode, 'func', img2float, 'in_file')
In [ ]:
"""
Extract the middle volume of the first run as the reference
"""
extract_ref = pe.Node(interface=fsl.ExtractROI(t_size=1),
name='extractref')
In [ ]:
"""
Define a function to pick the first file from a list of files
"""
def pickfirst(files):
if isinstance(files, list):
return files[0]
else:
return files
preproc.connect(img2float, ('out_file', pickfirst), extract_ref, 'in_file')
In [ ]:
"""
Define a function to return the 1 based index of the middle volume
"""
def getmiddlevolume(func):
from nibabel import load
from nipype.utils import NUMPY_MMAP
funcfile = func
if isinstance(func, list):
funcfile = func[0]
_, _, _, timepoints = load(funcfile, mmap=NUMPY_MMAP).shape
return int(timepoints / 2) - 1
preproc.connect(inputnode, ('func', getmiddlevolume), extract_ref, 't_min')
In [ ]:
"""
Realign the functional runs to the middle volume of the first run
"""
motion_correct = pe.MapNode(interface=fsl.MCFLIRT(save_mats=True,
save_plots=True),
name='realign',
iterfield=['in_file'])
preproc.connect(img2float, 'out_file', motion_correct, 'in_file')
preproc.connect(extract_ref, 'roi_file', motion_correct, 'ref_file')
In [ ]:
"""
Plot the estimated motion parameters
"""
plot_motion = pe.MapNode(interface=fsl.PlotMotionParams(in_source='fsl'),
name='plot_motion',
iterfield=['in_file'])
plot_motion.iterables = ('plot_type', ['rotations', 'translations'])
preproc.connect(motion_correct, 'par_file', plot_motion, 'in_file')
In [ ]:
"""
Extract the mean volume of the first functional run
"""
meanfunc = pe.Node(interface=fsl.ImageMaths(op_string='-Tmean',
suffix='_mean'),
name='meanfunc')
preproc.connect(motion_correct, ('out_file', pickfirst), meanfunc, 'in_file')
In [ ]:
"""
Strip the skull from the mean functional to generate a mask
"""
meanfuncmask = pe.Node(interface=fsl.BET(mask=True,
no_output=True,
frac=0.3),
name='meanfuncmask')
preproc.connect(meanfunc, 'out_file', meanfuncmask, 'in_file')
In [ ]:
"""
Mask the functional runs with the extracted mask
"""
maskfunc = pe.MapNode(interface=fsl.ImageMaths(suffix='_bet',
op_string='-mas'),
iterfield=['in_file'],
name='maskfunc')
preproc.connect(motion_correct, 'out_file', maskfunc, 'in_file')
preproc.connect(meanfuncmask, 'mask_file', maskfunc, 'in_file2')
In [ ]:
"""
Determine the 2nd and 98th percentile intensities of each functional run
"""
getthresh = pe.MapNode(interface=fsl.ImageStats(op_string='-p 2 -p 98'),
iterfield=['in_file'],
name='getthreshold')
preproc.connect(maskfunc, 'out_file', getthresh, 'in_file')
In [ ]:
"""
Threshold the first run of the functional data at 10% of the 98th percentile
"""
threshold = pe.Node(interface=fsl.ImageMaths(out_data_type='char',
suffix='_thresh'),
name='threshold')
preproc.connect(maskfunc, ('out_file', pickfirst), threshold, 'in_file')
In [ ]:
"""
Define a function to get 10% of the intensity
"""
def getthreshop(thresh):
return '-thr %.10f -Tmin -bin' % (0.1 * thresh[0][1])
preproc.connect(getthresh, ('out_stat', getthreshop), threshold, 'op_string')
In [ ]:
"""
Determine the median value of the functional runs using the mask
"""
medianval = pe.MapNode(interface=fsl.ImageStats(op_string='-k %s -p 50'),
iterfield=['in_file'],
name='medianval')
preproc.connect(motion_correct, 'out_file', medianval, 'in_file')
preproc.connect(threshold, 'out_file', medianval, 'mask_file')
In [ ]:
"""
Dilate the mask
"""
dilatemask = pe.Node(interface=fsl.ImageMaths(suffix='_dil',
op_string='-dilF'),
name='dilatemask')
preproc.connect(threshold, 'out_file', dilatemask, 'in_file')
In [ ]:
"""
Mask the motion corrected functional runs with the dilated mask
"""
maskfunc2 = pe.MapNode(interface=fsl.ImageMaths(suffix='_mask',
op_string='-mas'),
iterfield=['in_file'],
name='maskfunc2')
preproc.connect(motion_correct, 'out_file', maskfunc2, 'in_file')
preproc.connect(dilatemask, 'out_file', maskfunc2, 'in_file2')
In [ ]:
"""
Determine the mean image from each functional run
"""
meanfunc2 = pe.MapNode(interface=fsl.ImageMaths(op_string='-Tmean',
suffix='_mean'),
iterfield=['in_file'],
name='meanfunc2')
preproc.connect(maskfunc2, 'out_file', meanfunc2, 'in_file')
In [ ]:
"""
Merge the median values with the mean functional images into a coupled list
"""
mergenode = pe.Node(interface=util.Merge(2, axis='hstack'),
name='merge')
preproc.connect(meanfunc2, 'out_file', mergenode, 'in1')
preproc.connect(medianval, 'out_stat', mergenode, 'in2')
In [ ]:
"""
Smooth each run using SUSAN with the brightness threshold set to 75% of the
median value for each run and a mask constituting the mean functional
"""
smooth = pe.MapNode(interface=fsl.SUSAN(),
iterfield=['in_file', 'brightness_threshold', 'usans'],
name='smooth')
In [ ]:
"""
Define a function to get the brightness threshold for SUSAN
"""
def getbtthresh(medianvals):
return [0.75 * val for val in medianvals]
def getusans(x):
return [[tuple([val[0], 0.75 * val[1]])] for val in x]
preproc.connect(maskfunc2, 'out_file', smooth, 'in_file')
preproc.connect(medianval, ('out_stat', getbtthresh), smooth, 'brightness_threshold')
preproc.connect(mergenode, ('out', getusans), smooth, 'usans')
In [ ]:
"""
Mask the smoothed data with the dilated mask
"""
maskfunc3 = pe.MapNode(interface=fsl.ImageMaths(suffix='_mask',
op_string='-mas'),
iterfield=['in_file'],
name='maskfunc3')
preproc.connect(smooth, 'smoothed_file', maskfunc3, 'in_file')
preproc.connect(dilatemask, 'out_file', maskfunc3, 'in_file2')
In [ ]:
"""
Scale each volume of the run so that the median value of the run is set to 10000
"""
intnorm = pe.MapNode(interface=fsl.ImageMaths(suffix='_intnorm'),
iterfield=['in_file', 'op_string'],
name='intnorm')
preproc.connect(maskfunc3, 'out_file', intnorm, 'in_file')
In [ ]:
"""
Define a function to get the scaling factor for intensity normalization
"""
def getinormscale(medianvals):
return ['-mul %.10f' % (10000. / val) for val in medianvals]
preproc.connect(medianval, ('out_stat', getinormscale), intnorm, 'op_string')
In [ ]:
"""
Perform temporal highpass filtering on the data
"""
highpass = pe.MapNode(interface=fsl.ImageMaths(suffix='_tempfilt'),
iterfield=['in_file'],
name='highpass')
preproc.connect(intnorm, 'out_file', highpass, 'in_file')
In [ ]:
"""
Generate a mean functional image from the first run
"""
meanfunc3 = pe.MapNode(interface=fsl.ImageMaths(op_string='-Tmean',
suffix='_mean'),
iterfield=['in_file'],
name='meanfunc3')
preproc.connect(highpass, ('out_file', pickfirst), meanfunc3, 'in_file')
In [ ]:
"""
Strip the structural image and coregister the mean functional image to the
structural image
"""
nosestrip = pe.Node(interface=fsl.BET(frac=0.3),
name='nosestrip')
skullstrip = pe.Node(interface=fsl.BET(mask=True),
name='stripstruct')
coregister = pe.Node(interface=fsl.FLIRT(dof=6),
name='coregister')
In [ ]:
"""
Use :class:`nipype.algorithms.rapidart` to determine which of the
images in the functional series are outliers based on deviations in
intensity and/or movement.
"""
art = pe.MapNode(interface=ra.ArtifactDetect(use_differences=[True, False],
use_norm=True,
norm_threshold=1,
zintensity_threshold=3,
parameter_source='FSL',
mask_type='file'),
iterfield=['realigned_files', 'realignment_parameters'],
name="art")
preproc.connect([(inputnode, nosestrip, [('struct', 'in_file')]),
(nosestrip, skullstrip, [('out_file', 'in_file')]),
(skullstrip, coregister, [('out_file', 'in_file')]),
(meanfunc2, coregister, [(('out_file', pickfirst), 'reference')]),
(motion_correct, art, [('par_file', 'realignment_parameters')]),
(maskfunc2, art, [('out_file', 'realigned_files')]),
(dilatemask, art, [('out_file', 'mask_file')]),
])