[MRG] Mne watershed bem

.travis.yml

@@ -124,7 +124,7 @@ install:
       fi;
 
 script:
-    - nosetests --with-timer --timer-top-n 30 --verbosity=2 $COVERAGE
+    - nosetests -a '!ultra_slow_test' --with-timer --timer-top-n 30 --verbosity=2 $COVERAGE
     - if [ "${DEPS}" == "nodata" ]; then
         make flake;
       fi;

Makefile

@@ -43,6 +43,10 @@ $(CURDIR)/examples/MNE-sample-data/MEG/sample/sample_audvis_raw.fif:
 	ln -sf ${PWD}/examples/MNE-sample-data ${PWD}/MNE-sample-data
 
 test: in
+	rm -f .coverage
+	$(NOSETESTS) -a '!ultra_slow_test' mne
+
+test-full: in
 	rm -f .coverage
 	$(NOSETESTS) mne
 
@@ -53,10 +57,8 @@ test-no-testing-data: in
 	@MNE_SKIP_TESTING_DATASET_TESTS=true \
 	$(NOSETESTS) mne
 
-
 test-no-sample-with-coverage: in testing_data
 	rm -rf coverage .coverage
-	@MNE_SKIP_SAMPLE_DATASET_TESTS=true \
 	$(NOSETESTS) --with-coverage --cover-package=mne --cover-html --cover-html-dir=coverage
 
 test-doc: sample_data testing_data

appveyor.yml

@@ -33,4 +33,4 @@ build: false  # Not a C# project, build stuff at the test step instead.
 
 test_script:
   # Run the project tests, but (sadly) exclude ones that take a long time
-  - "nosetests --verbosity=2 -a !slow_test --with-timer --timer-top-n=20 --timer-ok 5 --timer-warning 15"
+  - "nosetests --verbosity=2 -a '!slow_test,!ultra_slow_test' --with-timer --timer-top-n=20 --timer-ok 5 --timer-warning 15"

doc/source/python_reference.rst

@@ -557,6 +557,14 @@ Functions:
    write_bem_surface
    write_trans
 
+.. currentmodule:: mne.bem
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   make_watershed_bem
+
 .. currentmodule:: mne.forward
 
 .. autosummary::

mne/bem.py

@@ -4,13 +4,18 @@
 #
 # License: BSD (3-clause)
 
+import sys
+import os
+import os.path as op
+import shutil
 import numpy as np
 from scipy import linalg
 
 from .fixes import partial
-from .utils import verbose, logger
+from .utils import (verbose, logger, run_subprocess, get_subjects_dir)
 from .io.constants import FIFF
 from .externals.six import string_types
+from .surface import read_surface, write_bem_surface
 
 
 # ############################################################################
@@ -332,3 +337,116 @@ def cost_fun(x, points):
     origin = x_opt[:3]
     radius = x_opt[3]
     return radius, origin
+
+
+# ############################################################################
+# Create BEM surfaces
+
+@verbose
+def make_watershed_bem(subject, subjects_dir=None, overwrite=False,
+                       volume='T1', atlas=False, gcaatlas=False, preflood=None,
+                       verbose=None):
+    """
+    Create BEM surfaces using the watershed algorithm included with FreeSurfer
+
+    Parameters
+    ----------
+    subject : str
+        Subject name (required)
+    subjects_dir : str
+        Directory containing subjects data. If None use
+        the Freesurfer SUBJECTS_DIR environment variable.
+    overwrite : bool
+        Write over existing files
+    volume : str
+        Defaults to T1
+    atlas : bool
+        Specify the --atlas option for mri_watershed
+    gcaatlas : bool
+        Use the subcortical atlas
+    preflood : int
+        Change the preflood height
+    verbose : bool, str or None
+        If not None, override default verbose level
+
+    .. versionadded:: 0.10
+    """
+    env = os.environ.copy()
+
+    if not os.environ.get('FREESURFER_HOME'):
+        raise RuntimeError('FREESURFER_HOME environment variable not set')
+
+    env['SUBJECT'] = subject
+
+    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
+    env['SUBJECTS_DIR'] = subjects_dir
+
+    subject_dir = op.join(subjects_dir, subject)
+    mri_dir = op.join(subject_dir, 'mri')
+    T1_dir = op.join(mri_dir, volume)
+    T1_mgz = op.join(mri_dir, volume + '.mgz')
+    bem_dir = op.join(subject_dir, 'bem')
+    ws_dir = op.join(subject_dir, 'bem', 'watershed')
+
+    if not op.isdir(subject_dir):
+        raise RuntimeError('Could not find the MRI data directory "%s"'
+                           % subject_dir)
+    if not op.isdir(bem_dir):
+        os.makedirs(bem_dir)
+    if not op.isdir(T1_dir) and not op.isfile(T1_mgz):
+        raise RuntimeError('Could not find the MRI data')
+    if op.isdir(ws_dir):
+        if not overwrite:
+            raise RuntimeError('%s already exists. Use the --overwrite option'
+                               'to recreate it.' % ws_dir)
+        else:
+            shutil.rmtree(ws_dir)
+    # put together the command
+    cmd = ['mri_watershed']
+    if preflood:
+        cmd += ["-h",  "%s" % int(preflood)]
+
+    if gcaatlas:
+        cmd += ['-atlas', '-T1', '-brain_atlas', env['FREESURFER_HOME'] +
+                '/average/RB_all_withskull_2007-08-08.gca',
+                subject_dir + '/mri/transforms/talairach_with_skull.lta']
+    elif atlas:
+        cmd += ['-atlas']
+    if op.exists(T1_mgz):
+        cmd += ['-useSRAS', '-surf', op.join(ws_dir, subject), T1_mgz,
+                op.join(ws_dir, 'ws')]
+    else:
+        cmd += ['-useSRAS', '-surf', op.join(ws_dir, subject), T1_dir,
+                op.join(ws_dir, 'ws')]
+    # report and run
+    logger.info('\nRunning mri_watershed for BEM segmentation with the '
+                'following parameters:\n\n'
+                'SUBJECTS_DIR = %s\n'
+                'SUBJECT = %s\n'
+                'Results dir = %s\n' % (subjects_dir, subject, ws_dir))
+    os.makedirs(op.join(ws_dir, 'ws'))
+    run_subprocess(cmd, env=env, stdout=sys.stdout)
+    #
+    os.chdir(ws_dir)
+    if op.isfile(T1_mgz):
+        # XXX : do this with python code
+        surfaces = [subject + '_brain_surface', subject +
+                    '_inner_skull_surface', subject + '_outer_skull_surface',
+                    subject + '_outer_skin_surface']
+        for s in surfaces:
+            cmd = ['mne_convert_surface', '--surf', s, '--mghmri', T1_mgz,
+                   '--surfout', s, "--replacegeom"]
+            run_subprocess(cmd, env=env, stdout=sys.stdout)
+    os.chdir(bem_dir)
+    if op.isfile(subject + '-head.fif'):
+        os.remove(subject + '-head.fif')
+
+    # run the equivalent of mne_surf2bem
+    points, tris = read_surface(op.join(ws_dir,
+                                        subject + '_outer_skin_surface'))
+    points *= 1e-3
+    surf = dict(coord_frame=5, id=4, nn=None, np=len(points),
+                ntri=len(tris), rr=points, sigma=1, tris=tris)
+    write_bem_surface(subject + '-head.fif', surf)
+
+    logger.info('Created %s/%s-head.fif\n\nComplete.' % (bem_dir, subject))

mne/commands/mne_watershed_bem.py

@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+# Authors: Lorenzo De Santis
+"""
+
+    Create BEM surfaces using the watershed algorithm included with
+        FreeSurfer
+
+"""
+
+from __future__ import print_function
+import sys
+
+from mne.bem import make_watershed_bem
+
+
+def run():
+    from mne.commands.utils import get_optparser
+
+    parser = get_optparser(__file__)
+
+    parser.add_option("-s", "--subject", dest="subject",
+                      help="Subject name (required)", default=None)
+    parser.add_option("-d", "--subjects-dir", dest="subjects_dir",
+                      help="Subjects directory", default=None)
+    parser.add_option("-o", "--overwrite", dest="overwrite",
+                      help="Write over existing files", action="store_true")
+    parser.add_option("-v", "--volume", dest="volume",
+                      help="Defaults to T1", default='T1')
+    parser.add_option("-a", "--atlas", dest="atlas",
+                      help="Specify the --atlas option for mri_watershed",
+                      default=False, action="store_true")
+    parser.add_option("-g", "--gcaatlas", dest="gcaatlas",
+                      help="Use the subcortical atlas", default=False,
+                      action="store_true")
+    parser.add_option("-p", "--preflood", dest="preflood",
+                      help="Change the preflood height", default=None)
+    parser.add_option("--verbose", dest="verbose",
+                      help="If not None, override default verbose level",
+                      default=None)
+
+    options, args = parser.parse_args()
+
+    if options.subject is None:
+        parser.print_help()
+        sys.exit(1)
+
+    subject = options.subject
+    subjects_dir = options.subjects_dir
+    overwrite = options.overwrite
+    volume = options.volume
+    atlas = options.atlas
+    gcaatlas = options.gcaatlas
+    preflood = options.preflood
+    verbose = options.verbose
+
+    make_watershed_bem(subject=subject, subjects_dir=subjects_dir,
+                       overwrite=overwrite, volume=volume, atlas=atlas,
+                       gcaatlas=gcaatlas, preflood=preflood, verbose=verbose)
+
+is_main = (__name__ == '__main__')
+if is_main:
+    run()

mne/commands/tests/test_commands.py

@@ -10,9 +10,10 @@
                           mne_compute_proj_ecg, mne_compute_proj_eog,
                           mne_coreg, mne_flash_bem_model, mne_kit2fiff,
                           mne_make_scalp_surfaces, mne_maxfilter,
-                          mne_report, mne_surf2bem)
+                          mne_report, mne_surf2bem, mne_watershed_bem)
 from mne.utils import (run_tests_if_main, _TempDir, requires_mne, requires_PIL,
-                       requires_mayavi, requires_tvtk, ArgvSetter, slow_test)
+                       requires_mayavi, requires_tvtk, requires_freesurfer,
+                       ArgvSetter, slow_test, ultra_slow_test)
 from mne.io import Raw
 from mne.datasets import testing
 
@@ -177,4 +178,28 @@ def test_surf2bem():
     check_usage(mne_surf2bem)
 
 
+@ultra_slow_test
+@requires_freesurfer
+@testing.requires_testing_data
+def test_watershed_bem():
+    """Test mne watershed bem"""
+    check_usage(mne_watershed_bem)
+    # Copy necessary files to tempdir
+    tempdir = _TempDir()
+    mridata_path = op.join(subjects_dir, 'sample', 'mri')
+    mridata_path_new = op.join(tempdir, 'sample', 'mri')
+    os.mkdir(op.join(tempdir, 'sample'))
+    os.mkdir(mridata_path_new)
+    if op.exists(op.join(mridata_path, 'T1')):
+        shutil.copytree(op.join(mridata_path, 'T1'), op.join(mridata_path_new,
+                        'T1'))
+    if op.exists(op.join(mridata_path, 'T1.mgz')):
+        shutil.copyfile(op.join(mridata_path, 'T1.mgz'),
+                        op.join(mridata_path_new, 'T1.mgz'))
+
+    with ArgvSetter(('-d', tempdir, '-s', 'sample', '-o'),
+                    disable_stdout=False, disable_stderr=False):
+        mne_watershed_bem.run()
+
+
 run_tests_if_main()

mne/utils.py

@@ -547,6 +547,13 @@ def slow_test(f):
     return f
 
 
+@nottest
+def ultra_slow_test(f):
+    """Decorator for ultra slow tests"""
+    f.ultra_slow_test = True
+    return f
+
+
 def has_nibabel(vox2ras_tkr=False):
     """Determine if nibabel is installed
 
@@ -777,7 +784,7 @@ def run_subprocess(command, verbose=None, *args, **kwargs):
                "that you use '$HOME' instead of '~'.")
         warnings.warn(msg)
 
-    logger.info("Running subprocess: %s" % str(command))
+    logger.info("Running subprocess: %s" % ' '.join(command))
     try:
         p = subprocess.Popen(command, *args, **kwargs)
     except Exception: