BUG: Problematic surfaces extracted using watershed

[sdeslauriers]

Running make_bem_model on the output of watershed_bem on subject 109123 of the HCP raises the following error:

RuntimeError: Surface inner skull is not completely inside surface outer skull

Checking the surfaces using freeview_bem_surfaces shows that the inner skull intersects the outer skull and skin surfaces. In addition, the inner skull surface has loops. This also happens with other subjects.

Steps and/or code to reproduce

I am assuming the HCP data for subject 109123 was downloaded and extracted to HCP_DIR. The data is converted to MNE format using mne-hcp (import hcp in the code) to SUBJECTS_DIR.

import os
import subprocess

import hcp
import mne

HCP_DIR = '...'
SUBJECTS_DIR = '...'

subject = '109123'

hcp.make_mne_anatomy(
    subject,
    subjects_dir=SUBJECTS_DIR,
    hcp_path=HCP_DIR,
    recordings_path=SUBJECTS_DIR)

os.chdir(os.path.join(SUBJECTS_DIR, subject))
subprocess.run([
    'mne', 'watershed_bem',
    '-s', subject,
    '-d', SUBJECTS_DIR,
    '-o'])

surfaces = mne.make_bem_model(
    subject,
    ico=4,
    subjects_dir=SUBJECTS_DIR)

Expected results

Expected non intersecting surfaces and surfaces with correct topology.

Actual results

Traceback (most recent call last):
  File "/.../scripts/solve-forward-problem.py", line 168, in <module>
    main()
  File "/.../scripts/solve-forward-problem.py", line 73, in main
    subjects_dir=SUBJECTS_DIR
  File "<string>", line 2, in make_bem_model
  File "/.../mne-python/mne/utils.py", line 748, in verbose
    return function(*args, **kwargs)
  File "/.../mne-python/mne/bem.py", line 572, in make_bem_model
    surfaces = _surfaces_to_bem(surfaces, ids, conductivity, ico)
  File "/.../mne-python/mne/bem.py", line 502, in _surfaces_to_bem
    _check_surfaces(surfs, incomplete=incomplete)
  File "/.../mne-python/mne/bem.py", line 443, in _check_surfaces
    _assert_inside(surf_2, surf_1)
  File "/.../mne-python/mne/bem.py", line 432, in _assert_inside
    % (_surf_name[fro['id']], _surf_name[to['id']]))
RuntimeError: Surface inner skull is not completely inside surface outer skull

Additional information

Python:        3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 09:53:17)  [GCC 4.8.2 20140120 (Red Hat 4.8.2-15)]
Executable:    /.../bin/python
CPU:           x86_64: 24 cores
Memory:        Unavailable (requires "psutil" package)
mne:           0.17.dev0
numpy:         1.15.0 {blas=mkl_rt, lapack=mkl_rt}
scipy:         1.1.0
matplotlib:    2.2.3 {backend=Qt5Agg}

sklearn:       0.19.1
nibabel:       2.3.0
mayavi:        Not found
cupy:          Not found
pandas:        0.23.4
dipy:          0.14.0

[agramfort]

2 comments: 1) I would use only 1 layer for this subject 3) but the inner brain layer seems not super smooth (see bottom of the brain) see if using the --atlas option (https://github.com/mne-tools/mne-python/blob/master/mne/commands/mne_watershed_bem.py#L30) helps

[sdeslauriers]

1. That is my fix for these subjects, but would not work for EEG data.
2. Using the --atlas option removed the loop, but the inner skull still intersects the other surfaces.

[agramfort]

1. That is my fix for these subjects, but would not work for EEG data. true but then I don't know how to fix this. 1. Using the --atlas option removed the loop, but the inner skull still intersects the other surfaces.

great. Can you open a PR for the next person using these data? I guess it should be on mne-hcp repo.

[sdeslauriers]

I am not sure what you mean. Open a PR for what?

[agramfort]

yes anything so the next person does not get stuck. Add a meaningful error message or something thanks

[sdeslauriers]

Opened an issue instead of a PR, because I really don't know what that PR would contain.
mne-tools/mne-hcp#60

[agramfort]

thanks. That helps.