Skip to content

[ENH/WIP] Add SPM Fieldmap Tool wrapper #1905

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

Merged
merged 8 commits into from
Feb 20, 2018
Merged
Show file tree
Hide file tree
Changes from 3 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
2 changes: 1 addition & 1 deletion nipype/interfaces/spm/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

from .base import (Info, SPMCommand, logger, no_spm, scans_for_fname,
scans_for_fnames)
from .preprocess import (SliceTiming, Realign, Coregister, Normalize,
from .preprocess import (FieldMap, SliceTiming, Realign, Coregister, Normalize,
Normalize12, Segment, Smooth, NewSegment, DARTEL,
DARTELNorm2MNI, CreateWarped, VBMSegment)
from .model import (Level1Design, EstimateModel, EstimateContrast, Threshold,
Expand Down
147 changes: 143 additions & 4 deletions nipype/interfaces/spm/preprocess.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,14 +22,153 @@
# Local imports
from ...utils.filemanip import (fname_presuffix, filename_to_list,
list_to_filename, split_filename)
from ..base import (OutputMultiPath, TraitedSpec, isdefined, traits,
InputMultiPath, File)
from .base import (SPMCommand, scans_for_fname, func_is_3d, scans_for_fnames,
SPMCommandInputSpec, ImageFileSPM)
from ..base import (OutputMultiPath, TraitedSpec, isdefined,
traits, InputMultiPath, File, Str)
from .base import (SPMCommand, scans_for_fname, func_is_3d,
scans_for_fnames, SPMCommandInputSpec, ImageFileSPM)

__docformat__ = 'restructuredtext'


class FieldMapInputSpec(SPMCommandInputSpec):
jobtype = traits.Enum('calculatevdm', 'applyvdm', usedefault=True,
desc='one of: calculatevdm, applyvdm')
phase = File(mandatory=True, exists=True, copyfile=False,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We often use the suffix _file to indicate inputs that should be files. So consider making this phase_file. Similarly below with magnitude_file.

Not a big deal, if you think that would be disorienting to SPM users.

field='subj.data.presubphasemag.phase',
desc='presubstracted phase file')
magnitude = File(mandatory=True, exists=True, copyfile=False,
field='subj.data.presubphasemag.magnitude',
desc='presubstracted magnitude file')
et = traits.Tuple(traits.Float, traits.Float, mandatory=True,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see this is matches the final part of the field name, so it may be that SPM users will find this intuitive, but you may want to consider echo_times or similar for clarity.

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')
blipdir = traits.Enum(1, -1, mandatory=True,
field='subj.defaults.defaultsval.blipdir',
desc='polarity of the phase-encode blips')
tert = traits.Float(mandatory=True,
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');
ajm = traits.Bool(False, usedefault=True,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, again, I would consider jacobian_modulation or similar.

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');
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');
fwhm = traits.Range(low=0, value=5, usedefault=True,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You've got fwhm twice. Perhaps this one should be mask_fwhm, and the other unwarp_fwhm?

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(copyfile=False, exists=True, mandatory=True,
field='subj.session.epi',
desc='EPI to unwarp');
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(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');


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


class FieldMap(SPMCommand):
"""Use the fieldmap toolbox from spm to calculate the voxel displacement map (VDM).

http://www.fil.ion.ucl.ac.uk/spm/doc/manual.pdf#page=173

To do
-----
Deal with real/imag magnitude images and with the two phase files case.

Examples
--------
>>> from nipype.interfaces.spm import FieldMap
>>> fm = FieldMap()
>>> fm.inputs.phase = 'phasediff.nii'
>>> fm.inputs.magnitude = 'magnitude1.nii'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

phasediff.nii and magnitude1.nii don't exist in nipype/testing/data. Either add them or use existing files (such as phase.nii and magnitude.nii).

>>> fm.inputs.et = [5.19, 7.65]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be a tuple, not a list.

>>> fm.inputs.blipdir = 1
>>> fm.inputs.tert = 15.6
>>> fm.inputs.epi = 'bold.nii'
>>> fm.run() # doctest: +SKIP

"""

input_spec = FieldMapInputSpec
output_spec = FieldMapOutputSpec
_jobtype = 'tools'
_jobname = 'fieldmap'

def _format_arg(self, opt, spec, val):
"""Convert input to appropriate format for spm
"""
if opt in ['phase', 'magnitude', 'anat', 'epi']:
return scans_for_fname(filename_to_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()
jobtype = self.inputs.jobtype
return [{'%s' % (jobtype): einputs[0]}]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the advantage of this over the following?

return [{self.inputs.jobtype: einputs[0]}]

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like better your style, but I originally used the syntax of another wrapper (typically Coregister).


def _list_outputs(self):
outputs = self._outputs().get()
jobtype = self.inputs.jobtype
if jobtype == "calculatevdm":
outputs['vdm'] = []
for phase in filename_to_list(self.inputs.phase):
outputs['vdm'].append(fname_presuffix(phase, prefix='vdm5_sc'))
outputs['vdm'] = list_to_filename(outputs['vdm'])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm a bit confused by this. self.inputs.phase is a single file, so this could just be:

if self.inputs.jobtype == 'calculatevdm':
    outputs['vdm'] = fname_presuffix(self.inputs.phase, prefix='vdm5_sc')

Or am I missing something?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're not, I just adapted another more complicated wrapper's piece of code without thinking about making it clearer. Thanks for noticing!


return outputs


class SliceTimingInputSpec(SPMCommandInputSpec):
in_files = InputMultiPath(
traits.Either(
Expand Down