Mismatch between the BEM surfaces from mri_make_bem_surfaces and mri_watershed

[ezemikulan]

Description of the problem

Hi, please correct me if I'm wrong but I think there is a mismatch between the BEM surfaces created with mri_make_bem_surfaces (FLASH) and the ones created with mri_watershed (MPRAGE).

The issue is that the _outerskull surface created with _mri_make_bemsurfaces (top panel in the image) corresponds to the inner skull and the _innerskull surface corresponds to the the brain. At least with the sample subject (I don't have access to FLASH sequences; the issue can be seen in the tutorial) so I don't know if this is a problem of the sample subject or the algorithm.

The surfaces created with the watershed algorithm seem to have the intuitive labeling (Bottom panel in the image).

The two algorithms are working differently: _mri_make_bemsurfaces creates the _outerskull by "stepping out" from the _innerskull whereas _mriwatershed creates it by "shrinking" the _outerskin (that's why it goes to the neck and what was happening in this question ).

Since make_bem_model is using the _inner_skull, outer_skull and outerskin surfaces, it looks like it is expecting the labeling of the surfaces from _mri_make_bemsurfaces and using the ones from the _mriwatershed would point to the wrong surfaces.

Not sure if I'm missing something and what would be the best way to fix this since the code of _mri_make_bemsurfaces is not being maintained anymore.

Maybe @alexrockhill has some ideas or can test with some FLASH images?

Steps to reproduce

import mne
import os.path as op

data_path = mne.datasets.sample.data_path()
subjects_dir = op.join(data_path, 'subjects')
mne.bem.make_watershed_bem('sample', subjects_dir=subjects_dir)

Link to data

No response

Expected results

Consitent labeling

Actual results

Inconsistent labeling

Additional information

Platform Linux-6.2.0-37-generic-x86_64-with-glibc2.35 Python 3.11.5 | packaged by conda-forge | (main, Aug 27 2023, 03:34:09) [GCC 12.3.0] Executable /home/eze/soft/anaconda3/envs/mne/bin/python CPU x86_64 (8 cores) Memory 31.2 GB Core ├☑ mne 1.5.1 ├☑ numpy 1.24.4 (OpenBLAS 0.3.23 with 8 threads) ├☑ scipy 1.11.2 ├☑ matplotlib 3.7.2 (backend=Qt5Agg) ├☑ pooch 1.7.0 └☑ jinja2 3.1.2 Numerical (optional) ├☑ sklearn 1.3.0 ├☑ numba 0.57.1 ├☑ nibabel 5.1.0 ├☑ nilearn 0.10.1 ├☑ dipy 1.7.0 ├☑ openmeeg 2.5.6 ├☑ pandas 2.0.3 └☐ unavailable cupy Visualization (optional) ├☑ pyvista 0.41.1 (OpenGL 4.5.0 NVIDIA 535.129.03 via NVIDIA GeForce GTX 1070/PCIe/SSE2) ├☑ pyvistaqt 0.0.0 ├☑ vtk 9.2.6 ├☑ qtpy 2.3.1 (PyQt5=5.15.8) ├☑ pyqtgraph 0.13.3 ├☑ mne-qt-browser 0.5.2 ├☑ ipywidgets 8.1.0 └☐ unavailable ipympl, trame_client, trame_server, trame_vtk, trame_vuetify Ecosystem (optional) ├☑ mne-bids 0.13 ├☑ mne-connectivity 0.5.0 └☐ unavailable mne-nirs, mne-features, mne-icalabel, mne-bids-pipeline

[alexrockhill]

Hmm yes does seem like an issue but unfortunately I've never touched that code and like you said it's not maintained very often; my understanding was that it was translated from Matti's MNE C code line for line so not sure anyone at the project understands it that well. I can take a look, maybe @larsoner might have some ideas too.

Relatedly, as you mentioned, it might be worth testing to see how much the accuracy of these models changes results. Specifically, does including the neck like that in the outer skull effect results much (I'd guess probably yes somewhat) and does using inner skull as outer skull and brain as inner skull effect results (that might be the best temporary fix).

As I'm writing this, I remembered these might be shallow wrappers around Freesurfer scripts so a full solution might involve changing those or probably more realistically reimplementing some or all of the functionality.

Ref: https://surfer.nmr.mgh.harvard.edu/fswiki/mri_watershed

[ezemikulan]

I was in fact testing how much results change by renaming the surfaces, it seems that in this case (epileptic spikes) results don't change much (on the left with the old surfaces, on the right with the surfaces renamed), there is a slight change in scale but the overall pattern remains the same (of note I'm using an EGI 256 channels which include the cheeks):

[alexrockhill]

One option for an alternative pipeline would be to do a symmetric diffeomorphic transform from fsaverage with a good BEM and then transform the points.

[alexrockhill]

I made a quick example of how to do this, if you think it's useful and what you are looking for @ezemikulan, perhaps you could make it into an mne example in a PR?

Note, it doesn't really work in this case since fsaverage is so blurry. I think it should work for two normal scans.

from shutil import copytree
import numpy as np
import mne
from mne.datasets import sample
import nibabel as nib

data_path = sample.data_path()
subjects_dir = data_path / "subjects"

# plot fsaverage before
mne.viz.plot_bem(subject='fsaverage', subjects_dir=subjects_dir)

# Uh oh, the BEM surfaces have the same issue with the outer skull
# I'd use sample as the template then (switch fsaverage for your scan)

template = 'sample'
subject = 'fsaverage'

# plot sample BEM, looks like good enough template, maybe could be edited a bit
mne.viz.plot_bem(subject=template, subjects_dir=subjects_dir)

# load data for SDR
t1_static = nib.load(subjects_dir / template / 'mri' / 'T1.mgz')
t1_moving = nib.load(subjects_dir / subject / 'mri' / 'T1.mgz')

reg_affine, sdr = mne.transforms.compute_volume_registration(
    static=t1_static, moving=t1_moving)

"""t1_moved = mne.transforms.apply_volume_registration(
    static=t1_static, moving=t1_moving, reg_affine=reg_affine)"""

# first save original BEM files so as not to overwrite!
if (subjects_dir / subject / 'bem').exists():
    copytree(subjects_dir / subject / 'bem', subjects_dir / subject / 'bem2')
(subjects_dir / subject / 'bem').mkdir(exist_ok=True)

# move to subject's space
for surf in ( 'inner_skull', 'outer_skull', 'outer_skin'):
    surf_rr, surf_tris = mne.read_surface(subjects_dir / template / 'bem' / f'{surf}.surf')
    surf_rr = mne.transforms.apply_trans(np.linalg.inv(reg_affine), surf_rr)  # affine registration
    surf_rr = sdr.transform_points(surf_rr, sdr.domain_grid2world, sdr.domain_world2grid)
    mne.write_surface(subjects_dir / subject / 'bem' / f'{surf}.surf', surf_rr, surf_tris, overwrite=True)

# finally, plot transformed BEM
mne.viz.plot_bem(subject=subject, subjects_dir=subjects_dir)

[larsoner]

Yes the algorithms are different. In theory the FLASH one should be more accurate. The watershed one just uses the outer skin (shrunk inward) and brain (expanded outward) to estimate inner and outer skull surfaces, whereas the FLASH code does something smarter (in theory). As for the magic that mri_make_bem_surfaces actually does under the hood I have no idea :(

See also https://github.com/mne-tools/mne-python/pull/8840 for ideas / code for using betsurf in case it's helpful somehow

[ezemikulan]

Thanks for your inputs. My concern is that, as far as I know, most users don't have FLASH sequences and use the watershed, and since the code expects the naming from FLASH (which is different from the watershed), the wrong surfaces are being used.

I'll try your suggestions and see what I can do.

[jhouck]

This probably will not be of any immediate help to @ezemikulan, but there's a summary of the methods involved in mri_make_bem_surfaces at https://surfer.nmr.mgh.harvard.edu/fswiki/mri_make_bem_surfaces. The C++ code is available in the Freesurfer github repository but the last time it came up I think there was some incompatibility between the Freesurfer license and MNE's BSD license.

[ezemikulan]

Thanks, I haven't had time to work on this yet