Skip to content

ENH/WIP: added Apply VDM functionality to FieldMap SPM interface #2879

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.

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 6 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
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
9 changes: 9 additions & 0 deletions .zenodo.json
Original file line number Diff line number Diff line change
Expand Up @@ -672,6 +672,15 @@
"name": "Junhao WEN",
"orcid": "0000-0003-2077-3070"
},
{
"affiliation": "Technische Universit\u00e4t Dresden, Faculty of Medicine",
"name": "Bernardoni, Fabio",
"orcid": "0000-0002-5112-405X"
},
{
"affiliation": "University of illinois urbana champaign",
"name": "Sharp, Paul"
},
{
"affiliation": "TIB \u2013 Leibniz Information Centre for Science and Technology and University Library, Hannover, Germany",
"name": "Leinweber, Katrin",
Expand Down
297 changes: 246 additions & 51 deletions nipype/interfaces/spm/preprocess.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,138 @@


class FieldMapInputSpec(SPMCommandInputSpec):

jobtype = traits.Enum('calculatevdm', 'applyvdm', usedefault=True,
desc='one of: calculatevdm, applyvdm')
phase_file = File(exists=True, copyfile=False,
field='subj.data.presubphasemag.phase',
desc='presubstracted phase file')
magnitude_file = File(exists=True, copyfile=False,
field='subj.data.presubphasemag.magnitude',
desc='presubstracted magnitude file')
echo_times = traits.Tuple(traits.Float, traits.Float,
field='subj.defaults.defaultsval.et',
desc='short and long echo times')
maskbrain = traits.Bool(True, usedefault=True,
field='subj.defaults.defaultsval.maskbrain',
desc='masking or no masking of the brain')
blip_direction = traits.Enum(1, -1,
field='subj.defaults.defaultsval.blipdir',
desc='polarity of the phase-encode blips')
total_readout_time = traits.Float(field='subj.defaults.defaultsval.tert',
desc='total EPI readout time')
epifm = traits.Bool(False, usedefault=True,
field='subj.defaults.defaultsval.epifm',
desc='epi-based field map');
jacobian_modulation = traits.Bool(False, usedefault=True,
field='subj.defaults.defaultsval.ajm',
desc='jacobian modulation');
# Unwarping defaults parameters
method = traits.Enum('Mark3D', 'Mark2D', 'Huttonish', usedefault=True,
desc='One of: Mark3D, Mark2D, Huttonish',
field='subj.defaults.defaultsval.uflags.method');
unwarp_fwhm = traits.Range(low=0, value=10, usedefault=True,
field='subj.defaults.defaultsval.uflags.fwhm',
desc='gaussian smoothing kernel width');
pad = traits.Range(low=0, value=0, usedefault=True,
field='subj.defaults.defaultsval.uflags.pad',
desc='padding kernel width');
ws = traits.Bool(True, usedefault=True,
field='subj.defaults.defaultsval.uflags.ws',
desc='weighted smoothing');
# Brain mask defaults parameters
template = File(copyfile=False, exists=True,
field='subj.defaults.defaultsval.mflags.template',
desc='template image for brain masking');
mask_fwhm = traits.Range(low=0, value=5, usedefault=True,
field='subj.defaults.defaultsval.mflags.fwhm',
desc='gaussian smoothing kernel width');
nerode = traits.Range(low=0, value=2, usedefault=True,
field='subj.defaults.defaultsval.mflags.nerode',
desc='number of erosions');
ndilate = traits.Range(low=0, value=4, usedefault=True,
field='subj.defaults.defaultsval.mflags.ndilate',
desc='number of erosions');
thresh = traits.Float(0.5, usedefault=True,
field='subj.defaults.defaultsval.mflags.thresh',
desc='threshold used to create brain mask from segmented data');
reg = traits.Float(0.02, usedefault=True,
field='subj.defaults.defaultsval.mflags.reg',
desc='regularization value used in the segmentation');
# EPI unwarping for quality check
epi_file = File(copyfile=False, exists=True,
field='subj.session.epi',
desc='EPI to unwarp for quality check');
matchvdm = traits.Bool(True, usedefault=True,
field='subj.matchvdm',
desc='match VDM to EPI');
sessname = Str('_run-', usedefault=True,
field='subj.sessname',
desc='VDM filename extension');
writeunwarped = traits.Bool(False, usedefault=True,
field='subj.writeunwarped',
desc='write unwarped EPI');
anat_file = File(copyfile=False, exists=True,
field='subj.anat',
desc='anatomical image for comparison');
matchanat = traits.Bool(True, usedefault=True,
field='subj.matchanat',
desc='match anatomical image to EPI');

in_files = InputMultiObject(
traits.Either(ImageFileSPM(exists=True),
traits.List(ImageFileSPM(exists=True))),
field='data.scans',mandatory=True,
copyfile=True,
desc='list of filenames to apply the vdm to')
vdmfile = File(
field='data.vdmfile',
desc='Voxel displacement map to use',mandatory=True,
copyfile=True)
distortion_direction = traits.Int(
2, field='roptions.pedir', desc='phase encode direction input data have been acquired with',
usedefault=True)
write_which = traits.ListInt(
[2, 1],
field='roptions.which',
minlen=2,
maxlen=2,
usedefault=True,
desc='determines which images to apply vdm to')
interpolation = traits.Int(
4, field='roptions.rinterp', desc='phase encode direction input data have been acquired with',
usedefault=True)
reslice_interp = traits.Range(
low=0,
high=7,
field='roptions.rinterp',
desc='degree of b-spline used for interpolation')
write_wrap = traits.List(
traits.Int(),
minlen=3,
maxlen=3,
field='roptions.wrap',
desc=('Check if interpolation should wrap in [x,y,z]'))
write_mask = traits.Bool(
field='roptions.mask', desc='True/False mask time series images')
out_prefix = traits.String(
'u',
field='roptions.prefix',
usedefault=True,
desc='fieldmap corrected output prefix')


class FieldMapOutputSpec(TraitedSpec):
vdm = File(exists=True, desc='voxel difference map')
out_files = OutputMultiPath(
traits.Either(traits.List(File(exists=True)), File(exists=True)),
desc=('If jobtype is applyvdm, '
'these will be the fieldmap corrected files.'
' Otherwise, they will be copies '
'of in_files that have had their '
'headers rewritten.'))
mean_image = File(exists=True, desc='Mean image')

jobtype = traits.Enum(
"calculatevdm",
"applyvdm",
Expand Down Expand Up @@ -196,8 +328,7 @@ class FieldMapInputSpec(SPMCommandInputSpec):
)


class FieldMapOutputSpec(TraitedSpec):
vdm = File(exists=True, desc="voxel difference map")



class FieldMap(SPMCommand):
Expand Down Expand Up @@ -230,22 +361,83 @@ class FieldMap(SPMCommand):
_jobname = "fieldmap"

def _format_arg(self, opt, spec, val):
"""Convert input to appropriate format for spm"""
if opt in ["phase_file", "magnitude_file", "anat_file", "epi_file"]:

"""Convert input to appropriate format for spm
"""
if ((self.inputs.jobtype == "calculatevdm") and (opt in ['phase_file', 'magnitude_file', 'anat_file', 'epi_file'])):

return scans_for_fname(ensure_list(val))

if ((self.inputs.jobtype == "applyvdm") and (opt =='in_files')):
return scans_for_fnames(ensure_list(val))
if ((self.inputs.jobtype == "applyvdm") and (opt =='vdmfile')):
return scans_for_fname(ensure_list(val))
return super(FieldMap, self)._format_arg(opt, spec, val)



def _parse_inputs(self):
"""validate spm fieldmap options if set to None ignore"""
einputs = super(FieldMap, self)._parse_inputs()
return [{self.inputs.jobtype: einputs[0]}]

"""validate spm fieldmap options if set to None ignore
"""
if self.inputs.jobtype == "applyvdm":
einputs = (super(FieldMap, self)
._parse_inputs(skip=('jobtype','phase_file', 'magnitude_file',
'echo_times', 'blip_direction',
'total_readout_time','maskbrain',
'epifm','jacobian_modulation',
'method','unwarp_fwhm','pad','ws',
'template','mask_fwhm','nerode','ndilate',
'thresh','reg','epi_file','matchvdm',
'sessname','writeunwarped',
'anat_file','matchanat')))

else:
einputs = (super(FieldMap, self)
._parse_inputs(skip=('jobtype','in_files', 'vdmfile')))
jobtype = self.inputs.jobtype
# einputs = super(FieldMap, self)._parse_inputs()
# return [{self.inputs.jobtype: einputs[0]}]

return [{'%s' % (jobtype): einputs[0]}]



def _list_outputs(self):
outputs = self._outputs().get()
jobtype = self.inputs.jobtype
resliced_all = self.inputs.write_which[0] > 0
resliced_mean = self.inputs.write_which[1] > 0
if jobtype == "calculatevdm":
outputs["vdm"] = fname_presuffix(self.inputs.phase_file, prefix="vdm5_sc")

outputs['vdm'] = fname_presuffix(self.inputs.phase_file, prefix='vdm5_sc')
elif jobtype == "applyvdm":
# outputs['out_files'] = fname_presuffix(
# self.inputs.in_files, prefix=self.inputs.out_prefix)
if resliced_mean:
if isinstance(self.inputs.in_files[0], list):
first_image = self.inputs.in_files[0][0]
else:
first_image = self.inputs.in_files[0]
outputs['mean_image'] = fname_presuffix(
first_image, prefix='meanu')

if resliced_all:
outputs['out_files'] = []
for idx, imgf in enumerate(ensure_list(self.inputs.in_files)):
appliedvdm_run = []
if isinstance(imgf, list):
for i, inner_imgf in enumerate(ensure_list(imgf)):
newfile = fname_presuffix(inner_imgf,
prefix=self.inputs.out_prefix)
appliedvdm_run.append(newfile)
else:
appliedvdm_run = fname_presuffix(imgf,
prefix=self.inputs.out_prefix)
outputs['out_files'].append(appliedvdm_run)
return outputs



return outputs

Expand Down Expand Up @@ -282,12 +474,18 @@ class SliceTimingInputSpec(SPMCommandInputSpec):
ref_slice = traits.Either(
traits.Int(),
traits.Float(),
field="refslice",
desc="1-based Number of the reference slice or "
"reference time point if slice_order is in "
"onsets (ms)",
mandatory=True,
)

field='so',
desc=('1-based order or onset (in ms) in which '
'slices are acquired'),
mandatory=True)
ref_slice = traits.Int(
field='refslice',
desc='1-based Number of the reference slice or '
'reference time point if slice_order is in '
'onsets (ms)',
mandatory=True)

out_prefix = traits.String(
"a", field="prefix", usedefault=True, desc="slicetimed output prefix"
)
Expand Down Expand Up @@ -543,7 +741,6 @@ def _list_outputs(self):


class RealignUnwarpInputSpec(SPMCommandInputSpec):

in_files = InputMultiObject(
traits.Either(
ImageFileSPM(exists=True), traits.List(ImageFileSPM(exists=True))
Expand Down Expand Up @@ -613,12 +810,11 @@ class RealignUnwarpInputSpec(SPMCommandInputSpec):
desc="Regularisation factor. Default: 100000 (medium).",
)
est_jacobian_deformations = traits.Bool(
field="uweoptions.jm",
desc=(
"Jacobian deformations. In theory a good idea to include them, "
" in practice a bad idea. Default: No."
),
)

field='uweoptions.jm',
desc=('Jacobian deformations. In theory a good idea to include them, '
' in practice a bad idea. Default: No.'))

est_first_order_effects = traits.List(
traits.Int(),
minlen=1,
Expand Down Expand Up @@ -758,11 +954,12 @@ def _parse_inputs(self, skip=()):

if isdefined(self.inputs.in_files):
if isinstance(self.inputs.in_files, list):
data = [
dict(scans=sess, pmscan=pmscan) for sess in spmdict["data"]["scans"]
]

data = [dict(scans=sess, pmscan=pmscan)
for sess in spmdict['data']['scans']]
else:
data = [dict(scans=spmdict["data"]["scans"], pmscan=pmscan)]
data = [dict(scans=spmdict['data']['scans'], pmscan=pmscan)]


spmdict["data"] = data

Expand Down Expand Up @@ -2550,7 +2747,6 @@ def _list_outputs(self):


class VBMSegmentInputSpec(SPMCommandInputSpec):

in_files = InputMultiPath(
ImageFileSPM(exists=True),
desc="A list of files to be segmented",
Expand Down Expand Up @@ -2720,7 +2916,6 @@ class VBMSegmentInputSpec(SPMCommandInputSpec):


class VBMSegmentOuputSpec(TraitedSpec):

native_class_images = traits.List(
traits.List(File(exists=True)), desc="native space probability maps"
)
Expand Down Expand Up @@ -2815,33 +3010,33 @@ def _list_outputs(self):
for i, tis in enumerate(["gm", "wm", "csf"]):
# native space

if getattr(self.inputs, "%s_native" % tis):
outputs["native_class_images"][i].append(
os.path.join(pth, "p%d%s.nii" % (i + 1, base))
)
if getattr(self.inputs, "%s_dartel" % tis) == 1:
outputs["dartel_input_images"][i].append(
os.path.join(pth, "rp%d%s.nii" % (i + 1, base))
)
elif getattr(self.inputs, "%s_dartel" % tis) == 2:
outputs["dartel_input_images"][i].append(
os.path.join(pth, "rp%d%s_affine.nii" % (i + 1, base))
)

if getattr(self.inputs, '%s_native' % tis):
outputs['native_class_images'][i].append(
os.path.join(pth, "p%d%s.nii" % (i + 1, base)))
if getattr(self.inputs, '%s_dartel' % tis) == 1:
outputs['dartel_input_images'][i].append(
os.path.join(pth, "rp%d%s.nii" % (i + 1, base)))
elif getattr(self.inputs, '%s_dartel' % tis) == 2:
outputs['dartel_input_images'][i].append(
os.path.join(pth, "rp%d%s_affine.nii" % (i + 1, base)))

# normalized space
if getattr(self.inputs, "%s_normalized" % tis):
outputs["normalized_class_images"][i].append(
os.path.join(pth, "w%sp%d%s.nii" % (dartel_px, i + 1, base))
)
if getattr(self.inputs, '%s_normalized' % tis):
outputs['normalized_class_images'][i].append(
os.path.join(pth, "w%sp%d%s.nii" % (dartel_px, i + 1,
base)))

if getattr(self.inputs, '%s_modulated_normalized' % tis) == 1:
outputs['modulated_class_images'][i].append(
os.path.join(pth, "mw%sp%d%s.nii" % (dartel_px, i + 1,
base)))
elif getattr(self.inputs,
'%s_modulated_normalized' % tis) == 2:
outputs['normalized_class_images'][i].append(
os.path.join(pth, "m0w%sp%d%s.nii" % (dartel_px, i + 1,
base)))

if getattr(self.inputs, "%s_modulated_normalized" % tis) == 1:
outputs["modulated_class_images"][i].append(
os.path.join(pth, "mw%sp%d%s.nii" % (dartel_px, i + 1, base))
)
elif getattr(self.inputs, "%s_modulated_normalized" % tis) == 2:
outputs["normalized_class_images"][i].append(
os.path.join(pth, "m0w%sp%d%s.nii" % (dartel_px, i + 1, base))
)

if self.inputs.pve_label_native:
outputs["pve_label_native_images"].append(
Expand Down
Loading