Skip to content

ENH: Update multi-stage recon-all directives #1991

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 4 commits into from
May 10, 2017
Merged
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
259 changes: 161 additions & 98 deletions nipype/interfaces/freesurfer/preprocess.py
Original file line number Diff line number Diff line change
Expand Up @@ -611,10 +611,18 @@ def _gen_filename(self, name):
class ReconAllInputSpec(CommandLineInputSpec):
subject_id = traits.Str("recon_all", argstr='-subjid %s',
desc='subject name', usedefault=True)
directive = traits.Enum('all', 'autorecon1', 'autorecon2', 'autorecon2-cp',
'autorecon2-wm', 'autorecon2-inflate1',
'autorecon2-perhemi', 'autorecon3', 'localGI',
'qcache', argstr='-%s', desc='process directive',
directive = traits.Enum('all', 'autorecon1',
# autorecon2 variants
'autorecon2', 'autorecon2-volonly',
'autorecon2-perhemi', 'autorecon2-inflate1',
'autorecon2-cp', 'autorecon2-wm',
# autorecon3 variants
'autorecon3', 'autorecon3-T2pial',
# Mix of autorecon2 and autorecon3 steps
'autorecon-pial', 'autorecon-hemi',
# Not "multi-stage flags"
'localGI', 'qcache',
argstr='-%s', desc='process directive',
usedefault=True, position=0)
hemi = traits.Enum('lh', 'rh', desc='hemisphere to process',
argstr="-hemi %s")
Expand Down Expand Up @@ -709,6 +717,38 @@ class ReconAll(CommandLine):
>>> reconall.inputs.flags = ["-cw256", "-qcache"]
>>> reconall.cmdline # doctest: +ALLOW_UNICODE
'recon-all -all -i structural.nii -cw256 -qcache -subjid foo -sd .'

Hemisphere may be specified regardless of directive:

>>> reconall.inputs.flags = []
>>> reconall.inputs.hemi = 'lh'
>>> reconall.cmdline # doctest: +ALLOW_UNICODE
'recon-all -all -i structural.nii -hemi lh -subjid foo -sd .'

``-autorecon-hemi`` uses the ``-hemi`` input to specify the hemisphere
to operate upon:

>>> reconall.inputs.directive = 'autorecon-hemi'
>>> reconall.cmdline # doctest: +ALLOW_UNICODE
'recon-all -autorecon-hemi lh -i structural.nii -subjid foo -sd .'

Hippocampal subfields can accept T1 and T2 images:

>>> reconall_subfields = ReconAll()
>>> reconall_subfields.inputs.subject_id = 'foo'
>>> reconall_subfields.inputs.directive = 'all'
>>> reconall_subfields.inputs.subjects_dir = '.'
>>> reconall_subfields.inputs.T1_files = 'structural.nii'
>>> reconall_subfields.inputs.hippocampal_subfields_T1 = True
>>> reconall_subfields.cmdline # doctest: +ALLOW_UNICODE
'recon-all -all -i structural.nii -hippocampal-subfields-T1 -subjid foo -sd .'
>>> reconall_subfields.inputs.hippocampal_subfields_T2 = (
... 'structural.nii', 'test')
>>> reconall_subfields.cmdline # doctest: +ALLOW_UNICODE
'recon-all -all -i structural.nii -hippocampal-subfields-T1T2 structural.nii test -subjid foo -sd .'
>>> reconall_subfields.inputs.hippocampal_subfields_T1 = False
>>> reconall_subfields.cmdline # doctest: +ALLOW_UNICODE
'recon-all -all -i structural.nii -hippocampal-subfields-T2 structural.nii test -subjid foo -sd .'
"""

_cmd = 'recon-all'
Expand Down Expand Up @@ -744,7 +784,7 @@ class ReconAll(CommandLine):
'mri/brainmask.mgz'], []),
]
if Info.looseversion() < LooseVersion("6.0.0"):
_autorecon2_steps = [
_autorecon2_volonly_steps = [
('gcareg', ['mri/transforms/talairach.lta'], []),
('canorm', ['mri/norm.mgz'], []),
('careg', ['mri/transforms/talairach.m3z'], []),
Expand All @@ -765,65 +805,51 @@ class ReconAll(CommandLine):
('fill', ['mri/filled.mgz',
# 'scripts/ponscc.cut.log',
], []),
('tessellate', ['surf/lh.orig.nofix', 'surf/rh.orig.nofix'], []),
('smooth1', ['surf/lh.smoothwm.nofix', 'surf/rh.smoothwm.nofix'],
[]),
('inflate1', ['surf/lh.inflated.nofix', 'surf/rh.inflated.nofix'],
[]),
('qsphere', ['surf/lh.qsphere.nofix', 'surf/rh.qsphere.nofix'],
[]),
('fix', ['surf/lh.orig', 'surf/rh.orig'], []),
('white', ['surf/lh.white', 'surf/rh.white',
'surf/lh.curv', 'surf/rh.curv',
'surf/lh.area', 'surf/rh.area',
'label/lh.cortex.label', 'label/rh.cortex.label'], []),
('smooth2', ['surf/lh.smoothwm', 'surf/rh.smoothwm'], []),
('inflate2', ['surf/lh.inflated', 'surf/rh.inflated',
'surf/lh.sulc', 'surf/rh.sulc',
'surf/lh.inflated.H', 'surf/rh.inflated.H',
'surf/lh.inflated.K', 'surf/rh.inflated.K'], []),
]
_autorecon2_lh_steps = [
Copy link
Member

Choose a reason for hiding this comment

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

only a nitpick - is there a reason this is called _lh_ as opposed to say _hemi_

otherwise looks reasonable to me.

Copy link
Member Author

Choose a reason for hiding this comment

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

Because I specify lh in the filenames. I then make a rh copy below, just swapping the filenames out.

('tessellate', ['surf/lh.orig.nofix'], []),
('smooth1', ['surf/lh.smoothwm.nofix'], []),
('inflate1', ['surf/lh.inflated.nofix'], []),
('qsphere', ['surf/lh.qsphere.nofix'], []),
('fix', ['surf/lh.orig'], []),
('white', ['surf/lh.white', 'surf/lh.curv', 'surf/lh.area',
'label/lh.cortex.label'], []),
('smooth2', ['surf/lh.smoothwm'], []),
('inflate2', ['surf/lh.inflated', 'surf/lh.sulc',
'surf/lh.inflated.H', 'surf/lh.inflated.K'], []),
# Undocumented in ReconAllTableStableV5.3
('curvstats', ['stats/lh.curv.stats', 'stats/rh.curv.stats'], []),
('curvstats', ['stats/lh.curv.stats'], []),
]
_autorecon3_steps = [
('sphere', ['surf/lh.sphere', 'surf/rh.sphere'], []),
('surfreg', ['surf/lh.sphere.reg', 'surf/rh.sphere.reg'], []),
('jacobian_white', ['surf/lh.jacobian_white',
'surf/rh.jacobian_white'], []),
('avgcurv', ['surf/lh.avg_curv', 'surf/rh.avg_curv'], []),
('cortparc', ['label/lh.aparc.annot', 'label/rh.aparc.annot'], []),
('pial', ['surf/lh.pial', 'surf/rh.pial',
'surf/lh.curv.pial', 'surf/rh.curv.pial',
'surf/lh.area.pial', 'surf/rh.area.pial',
'surf/lh.thickness', 'surf/rh.thickness'], []),
_autorecon3_lh_steps = [
('sphere', ['surf/lh.sphere'], []),
('surfreg', ['surf/lh.sphere.reg'], []),
('jacobian_white', ['surf/lh.jacobian_white'], []),
('avgcurv', ['surf/lh.avg_curv'], []),
('cortparc', ['label/lh.aparc.annot'], []),
('pial', ['surf/lh.pial', 'surf/lh.curv.pial', 'surf/lh.area.pial',
'surf/lh.thickness'], []),
# Misnamed outputs in ReconAllTableStableV5.3: ?h.w-c.pct.mgz
('pctsurfcon', ['surf/lh.w-g.pct.mgh', 'surf/rh.w-g.pct.mgh'], []),
('parcstats', ['stats/lh.aparc.stats', 'stats/rh.aparc.stats',
'label/aparc.annot.a2009s.ctab'], []),
('cortparc2', ['label/lh.aparc.a2009s.annot',
'label/rh.aparc.a2009s.annot'], []),
('parcstats2', ['stats/lh.aparc.a2009s.stats',
'stats/rh.aparc.a2009s.stats',
'label/aparc.annot.a2009s.ctab'], []),
('pctsurfcon', ['surf/lh.w-g.pct.mgh'], []),
('parcstats', ['stats/lh.aparc.stats'], []),
('cortparc2', ['label/lh.aparc.a2009s.annot'], []),
('parcstats2', ['stats/lh.aparc.a2009s.stats'], []),
# Undocumented in ReconAllTableStableV5.3
('cortparc3', ['label/lh.aparc.DKTatlas40.annot',
'label/rh.aparc.DKTatlas40.annot'], []),
('cortparc3', ['label/lh.aparc.DKTatlas40.annot'], []),
# Undocumented in ReconAllTableStableV5.3
('parcstats3', ['stats/lh.aparc.a2009s.stats',
'stats/rh.aparc.a2009s.stats',
'label/aparc.annot.a2009s.ctab'], []),
('parcstats3', ['stats/lh.aparc.a2009s.stats'], []),
('label-exvivo-ec', ['label/lh.entorhinal_exvivo.label'], []),
]
_autorecon3_added_steps = [
('cortribbon', ['mri/lh.ribbon.mgz', 'mri/rh.ribbon.mgz',
'mri/ribbon.mgz'], []),
('segstats', ['stats/aseg.stats'], []),
('aparc2aseg', ['mri/aparc+aseg.mgz',
'mri/aparc.a2009s+aseg.mgz'], []),
('wmparc', ['mri/wmparc.mgz', 'stats/wmparc.stats'], []),
('balabels', ['label/BA.ctab', 'label/BA.thresh.ctab'], []),
('label-exvivo-ec', ['label/lh.entorhinal_exvivo.label',
'label/rh.entorhinal_exvivo.label'], []),
]
else:
_autorecon2_steps = [
_autorecon2_volonly_steps = [
('gcareg', ['mri/transforms/talairach.lta'], []),
('canorm', ['mri/norm.mgz'], []),
('careg', ['mri/transforms/talairach.m3z'], []),
Expand All @@ -838,53 +864,39 @@ class ReconAll(CommandLine):
('fill', ['mri/filled.mgz',
# 'scripts/ponscc.cut.log',
], []),
('tessellate', ['surf/lh.orig.nofix', 'surf/rh.orig.nofix'], []),
('smooth1', ['surf/lh.smoothwm.nofix', 'surf/rh.smoothwm.nofix'],
[]),
('inflate1', ['surf/lh.inflated.nofix', 'surf/rh.inflated.nofix'],
[]),
('qsphere', ['surf/lh.qsphere.nofix', 'surf/rh.qsphere.nofix'],
[]),
('fix', ['surf/lh.orig', 'surf/rh.orig'], []),
('white', ['surf/lh.white.preaparc', 'surf/rh.white.preaparc',
'surf/lh.curv', 'surf/rh.curv',
'surf/lh.area', 'surf/rh.area',
'label/lh.cortex.label', 'label/rh.cortex.label'], []),
('smooth2', ['surf/lh.smoothwm', 'surf/rh.smoothwm'], []),
('inflate2', ['surf/lh.inflated', 'surf/rh.inflated',
'surf/lh.sulc', 'surf/rh.sulc'], []),
('curvHK', ['surf/lh.white.H', 'surf/rh.white.H',
'surf/lh.white.K', 'surf/rh.white.K',
'surf/lh.inflated.H', 'surf/rh.inflated.H',
'surf/lh.inflated.K', 'surf/rh.inflated.K'], []),
('curvstats', ['stats/lh.curv.stats', 'stats/rh.curv.stats'], []),
]
_autorecon3_steps = [
('sphere', ['surf/lh.sphere', 'surf/rh.sphere'], []),
('surfreg', ['surf/lh.sphere.reg', 'surf/rh.sphere.reg'], []),
('jacobian_white', ['surf/lh.jacobian_white',
'surf/rh.jacobian_white'], []),
('avgcurv', ['surf/lh.avg_curv', 'surf/rh.avg_curv'], []),
('cortparc', ['label/lh.aparc.annot', 'label/rh.aparc.annot'], []),
('pial', ['surf/lh.pial', 'surf/rh.pial',
'surf/lh.curv.pial', 'surf/rh.curv.pial',
'surf/lh.area.pial', 'surf/rh.area.pial',
'surf/lh.thickness', 'surf/rh.thickness'], []),
_autorecon2_lh_steps = [
('tessellate', ['surf/lh.orig.nofix'], []),
('smooth1', ['surf/lh.smoothwm.nofix'], []),
('inflate1', ['surf/lh.inflated.nofix'], []),
('qsphere', ['surf/lh.qsphere.nofix'], []),
('fix', ['surf/lh.orig'], []),
('white', ['surf/lh.white.preaparc', 'surf/lh.curv',
'surf/lh.area', 'label/lh.cortex.label'], []),
('smooth2', ['surf/lh.smoothwm'], []),
('inflate2', ['surf/lh.inflated', 'surf/lh.sulc'], []),
('curvHK', ['surf/lh.white.H', 'surf/lh.white.K',
'surf/lh.inflated.H', 'surf/lh.inflated.K'], []),
('curvstats', ['stats/lh.curv.stats'], []),
]
_autorecon3_lh_steps = [
('sphere', ['surf/lh.sphere'], []),
('surfreg', ['surf/lh.sphere.reg'], []),
('jacobian_white', ['surf/lh.jacobian_white'], []),
('avgcurv', ['surf/lh.avg_curv'], []),
('cortparc', ['label/lh.aparc.annot'], []),
('pial', ['surf/lh.pial', 'surf/lh.curv.pial',
'surf/lh.area.pial', 'surf/lh.thickness'], []),
('parcstats', ['stats/lh.aparc.stats'], []),
('cortparc2', ['label/lh.aparc.a2009s.annot'], []),
('parcstats2', ['stats/lh.aparc.a2009s.stats'], []),
('cortparc3', ['label/lh.aparc.DKTatlas.annot'], []),
('parcstats3', ['stats/lh.aparc.DKTatlas.stats'], []),
('pctsurfcon', ['surf/lh.w-g.pct.mgh'], []),
]
_autorecon3_added_steps = [
('cortribbon', ['mri/lh.ribbon.mgz', 'mri/rh.ribbon.mgz',
'mri/ribbon.mgz'], []),
('parcstats', ['stats/lh.aparc.stats', 'stats/rh.aparc.stats',
'label/aparc.annot.ctab'], []),
('cortparc2', ['label/lh.aparc.a2009s.annot',
'label/rh.aparc.a2009s.annot'], []),
('parcstats2', ['stats/lh.aparc.a2009s.stats',
'stats/rh.aparc.a2009s.stats',
'label/aparc.annot.a2009s.ctab'], []),
('cortparc3', ['label/lh.aparc.DKTatlas.annot',
'label/rh.aparc.DKTatlas.annot'], []),
('parcstats3', ['stats/lh.aparc.DKTatlas.stats',
'stats/rh.aparc.DKTatlas.stats',
'label/aparc.annot.DKTatlas.ctab'], []),
('pctsurfcon', ['surf/lh.w-g.pct.mgh', 'surf/rh.w-g.pct.mgh'], []),
('hyporelabel', ['mri/aseg.presurf.hypos.mgz'], []),
('aparc2aseg', ['mri/aparc+aseg.mgz',
'mri/aparc.a2009s+aseg.mgz',
Expand All @@ -900,6 +912,30 @@ class ReconAll(CommandLine):
'label/rh.entorhinal_exvivo.label'], []),
]

# Fill out autorecon2 steps
_autorecon2_rh_steps = [
(step, [out.replace('lh', 'rh') for out in outs], ins)
for step, outs, ins in _autorecon2_lh_steps]
_autorecon2_perhemi_steps = [
(step, [of for out in outs
for of in (out, out.replace('lh', 'rh'))], ins)
for step, outs, ins in _autorecon2_lh_steps]
_autorecon2_steps = _autorecon2_volonly_steps + _autorecon2_perhemi_steps

# Fill out autorecon3 steps
_autorecon3_rh_steps = [
(step, [out.replace('lh', 'rh') for out in outs], ins)
for step, outs, ins in _autorecon3_lh_steps]
_autorecon3_perhemi_steps = [
(step, [of for out in outs
for of in (out, out.replace('lh', 'rh'))], ins)
for step, outs, ins in _autorecon3_lh_steps]
_autorecon3_steps = _autorecon3_perhemi_steps + _autorecon3_added_steps

# Fill out autorecon-hemi lh/rh steps
_autorecon_lh_steps = (_autorecon2_lh_steps + _autorecon3_lh_steps)
_autorecon_rh_steps = (_autorecon2_rh_steps + _autorecon3_rh_steps)

_steps = _autorecon1_steps + _autorecon2_steps + _autorecon3_steps

_binaries = ['talairach', 'mri_normalize', 'mri_watershed',
Expand Down Expand Up @@ -955,14 +991,24 @@ def _is_resuming(self):
def _format_arg(self, name, trait_spec, value):
if name == 'T1_files':
if self._is_resuming():
return ''
return None
Copy link
Member Author

Choose a reason for hiding this comment

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

Switched to returning None just because that avoids the extra spaces in the cmdline string, making the doctests less painful to debug.

if name == 'hippocampal_subfields_T1' and \
isdefined(self.inputs.hippocampal_subfields_T2):
return ''
return None
if all((name == 'hippocampal_subfields_T2',
isdefined(self.inputs.hippocampal_subfields_T1) and
self.inputs.hippocampal_subfields_T1)):
trait_spec.argstr = trait_spec.argstr.replace('T2', 'T1T2')
argstr = trait_spec.argstr.replace('T2', 'T1T2')
return argstr % value
Copy link
Member Author

Choose a reason for hiding this comment

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

This part here was a (minor) bug, in that trait_spec is non-ephemeral, so changing trait_spec.argstr changes it for all future calls. Unlikely to affect anything except the doctests, but still worth cleaning up.

if name == 'directive' and value == 'autorecon-hemi':
if not isdefined(self.inputs.hemi):
raise ValueError("Directive 'autorecon-hemi' requires hemi "
"input to be set")
value += ' ' + self.inputs.hemi
if all((name == 'hemi',
isdefined(self.inputs.directive) and
self.inputs.directive == 'autorecon-hemi')):
return None
return super(ReconAll, self)._format_arg(name, trait_spec, value)

@property
Expand All @@ -985,8 +1031,25 @@ def cmdline(self):
steps = []
elif directive == 'autorecon1':
steps = self._autorecon1_steps
elif directive == 'autorecon2-volonly':
steps = self._autorecon2_volonly_steps
elif directive == 'autorecon2-perhemi':
steps = self._autorecon2_perhemi_steps
elif directive.startswith('autorecon2'):
steps = self._autorecon2_steps
if isdefined(self.inputs.hemi):
if self.inputs.hemi == 'lh':
steps = (self._autorecon2_volonly_steps +
self._autorecon2_lh_steps)
else:
steps = (self._autorecon2_volonly_steps +
self._autorecon2_rh_steps)
else:
steps = self._autorecon2_steps
elif directive == 'autorecon-hemi':
if self.inputs.hemi == 'lh':
steps = self._autorecon_lh_steps
else:
steps = self._autorecon_rh_steps
elif directive == 'autorecon3':
steps = self._autorecon3_steps
else:
Expand Down