Incorporating CT into 3 layer BEM

[jstout211 (Jeff Stout)]

Hi MNE group,

When I perform the watershed bem on a T1 scan, I typically get a reasonable inner skull boundary that works for a 1 layer MEG forward model. But the outer skull boundary is pretty bad and doesn’t appear that it would work well for an EEG forward model. I know that the Flash sequence is suggested to better resolve the inner/outer skull segmentation, but our group has not historically performed these scans and we have a lot of retrospective data. The epilepsy subjects that we scan with EEG+MEG all will have a presurgical CT - is it possible to incorporate the CT into the BEM model to estimate the skull boundary?

I can coreg and recut the CT to the freesurfer MRI, then threshold to get the bone mask. But at this point, I don’t know how I would incorporate this into the BEM. Is this better suited to forward model in the openMEEG software and subsequently incorporate into MNE afterwards?

Any help is greatly appreciated.

Thank you,
Jeff Stout

[agramfort (Alexandre Gramfort)]

I’ve never done this before but I think it should be possible.

hopefully someone can help here. Maybe @alexrockhill knows ?

Alex

[alexrockhill]

Hi Jeff, it does look like your outer skull surface is not accurate and it should be possible to use the registered CT to compute a different surface. This might be helpful for coregistering the CT if you have not done so already Locating Intracranial Electrode Contacts — MNE 0.24.dev0 documentation. Then the next step should be to compute a surface, which is done by Freesurfer with the watershed algorithm which I am not an expert on the details. MNE does have a marching cubes algorithm (mne.surface._marching_cubes) that can take an image and generate a surface. I’d have to look more into how Freesurfer enforces the surface-to-depth order of the model and how to reimplement this properly in MNE. I can take a bit of time as soon as I have a chance and see if I can get a reasonable example but a ready-made function doesn’t exist in MNE for this yet AFAIK.

[alexrockhill]

You may want to check the tips and see if anything here can help you so that you don’t have to use the CT: Frequently Asked Questions (FAQ) — MNE 0.23.4 documentation.

Otherwise you can try something like this with the CT

import os.path as op
import numpy as np
import mne
import nibabel as nib

misc_path = mne.datasets.misc.data_path()
sample_path = mne.datasets.sample.data_path()
subjects_dir = op.join(sample_path, 'subjects')

# align CT
T1 = nib.load(op.join(misc_path, 'seeg', 'sample_seeg', 'mri', 'T1.mgz'))
CT_orig = nib.load(op.join(misc_path, 'seeg', 'sample_seeg_CT.mgz'))
reg_affine = np.array([
    [0.99270756, -0.03243313, 0.11610254, -133.094156],
    [0.04374389, 0.99439665, -0.09623816, -97.58320673],
    [-0.11233068, 0.10061512, 0.98856381, -84.45551601],
    [0., 0., 0., 1.]])
CT_aligned = mne.transforms.apply_volume_registration(CT_orig, T1, reg_affine)

CT_data = np.array(CT_aligned.dataobj).astype(np.int64)
CT_data[CT_data < 1000] = 0
CT_data[CT_data > 0] = 1

coords, faces = mne.surface._marching_cubes(CT_data, level=[1], smooth=0.95)[0]
coords = mne.transforms.apply_trans(T1.header.get_vox2ras_tkr(), coords)

brain = mne.viz.Brain('sample_seeg', subjects_dir=op.join(misc_path, 'seeg'), alpha=0.8)
brain._renderer.mesh(*coords.T, faces, 'red', opacity=0.5)

# mne.write_surface(op.join(misc_path, 'seeg', 'sample_seeg', 'bem', 'outer_skull.surf'),
#				   coords, faces)

The code is unfinished however because you’d have to remove coordinates and faces that are on the inner skull surface and, for my example data at least (which is in the MNE sample data), there is a bit of a mess on the surface of the skull where electrodes were implanted. I can try and find a bit more time to see if I can get something more reasonable soon but this is a start at least.

[jstout211 (Jeff Stout)]

@alexrockhill

there is a bit of a mess on the surface of the skull where electrodes were implanted. I can try and find a bit more time to

@alexrockhill - thanks so much. I will take a look at this with our data. Appreciate the help.

-Jeff

[alexrockhill]

I spent a bit more time getting a more reasonable skull surface with this code:

import os.path as op
import numpy as np
import mne
import nibabel as nib
from tqdm import tqdm

misc_path = mne.datasets.misc.data_path()
sample_path = mne.datasets.sample.data_path()
subjects_dir = op.join(sample_path, 'subjects')

# align CT
T1 = nib.load(op.join(misc_path, 'seeg', 'sample_seeg', 'mri', 'T1.mgz'))
CT_orig = nib.load(op.join(misc_path, 'seeg', 'sample_seeg_CT.mgz'))
reg_affine = np.array([
    [0.99270756, -0.03243313, 0.11610254, -133.094156],
    [0.04374389, 0.99439665, -0.09623816, -97.58320673],
    [-0.11233068, 0.10061512, 0.98856381, -84.45551601],
    [0., 0., 0., 1.]])
CT_aligned = mne.transforms.apply_volume_registration(CT_orig, T1, reg_affine)

CT_data = np.array(CT_aligned.dataobj).astype(np.int64)
CT_data[CT_data < 1000] = 0
CT_data[CT_data > 0] = 1

coords, faces = mne.surface._marching_cubes(CT_data, level=[1], smooth=0.95)[0]
coords = mne.transforms.apply_trans(T1.header.get_vox2ras_tkr(), coords)

brain = mne.viz.Brain(
    'sample_seeg', subjects_dir=op.join(misc_path, 'seeg'), alpha=0.8)
brain._renderer.mesh(*coords.T, faces, 'red', opacity=0.5)

# load icosahedron and fit it to the surface
# could change to 5 for smoother but slower
ico = mne.surface._get_ico_surface(grade=2)

# find large radius from coords
coords_sph = mne.transforms._cart_to_sph(coords)
r = np.quantile(coords_sph[0], 0.95)
ico_rr = ico['rr'] * r + coords.mean(axis=0)
ico_sph = mne.transforms._cart_to_sph(ico_rr)

brain._renderer.mesh(*ico_rr.T, ico['tris'], 'blue', opacity=0.5)

ico_rr_shrunk = list()
delta_r = 0.5
# divide by 2 to roughly account for doubling of inner skull surface
n_closest = np.round(coords.size / 2 / ico_rr.size).astype(int)
for sph_pt in tqdm(ico_sph):
    pt = mne.transforms._sph_to_cart(np.array([sph_pt]))[0]
    dists = np.linalg.norm(pt - coords, axis=1)
    closest = np.argsort(dists)[:n_closest]  # could change this 100 too
    error = None  # fencepost
    next_error = np.quantile(dists[closest], 0.05)
    while error is None or next_error < error:
        error = next_error
        sph_pt[0] -= delta_r
        pt = mne.transforms._sph_to_cart(np.array([sph_pt]))[0]
        dists = np.linalg.norm(pt - coords, axis=1)
        next_error = np.quantile(dists[closest], 0.05)
    ico_rr_shrunk.append(pt)


ico_rr_shrunk = np.array(ico_rr_shrunk)
brain = mne.viz.Brain(
    'sample_seeg', subjects_dir=op.join(misc_path, 'seeg'), alpha=0.8)
brain._renderer.mesh(*coords.T, faces, 'red', opacity=0.5)
brain._renderer.mesh(*ico_rr_shrunk.T, ico['tris'], 'green', opacity=0.5)

# mne.write_surface(op.join(misc_path, 'seeg', 'sample_seeg',
#                           'bem', 'outer_skull.surf'),
#                   ico_rr_shrunk, ico['tris'])

which might be helpful to you.

I think it’s maybe a bit hacky to go in the MNE codebase unless the full watershed algorithm were reimplemented which involves a smoothness constraint and a template deformation as described here http://surfer.nmr.mgh.harvard.edu/ftp/articles/segonne05.pdf but this is an okay start and could perhaps be useable for your case. If it were me and I had this data

I would probably go into freeview and quickly edit the CT to remove the very large artifacts and then rerun with a finer icosahedron (it took 2.5 minutes on a macbook laptop to run the code above).

Anyway, hope that helps.

[alexrockhill]

I had a thought that incorporating some flooding might make this work much better

import os.path as op
import numpy as np
import mne
import nibabel as nib
import skimage
from tqdm import tqdm
import matplotlib.pyplot as plt


misc_path = mne.datasets.misc.data_path()
sample_path = mne.datasets.sample.data_path()
subjects_dir = op.join(sample_path, 'subjects')

# align CT
T1 = nib.load(op.join(misc_path, 'seeg', 'sample_seeg', 'mri', 'T1.mgz'))
CT_orig = nib.load(op.join(misc_path, 'seeg', 'sample_seeg_CT.mgz'))
reg_affine = np.array([
    [0.99270756, -0.03243313, 0.11610254, -133.094156],
    [0.04374389, 0.99439665, -0.09623816, -97.58320673],
    [-0.11233068, 0.10061512, 0.98856381, -84.45551601],
    [0., 0., 0., 1.]])
CT_aligned = mne.transforms.apply_volume_registration(CT_orig, T1, reg_affine)

CT_data = np.array(CT_aligned.dataobj).astype(np.int32)

plt.imshow(CT_data[128])
plt.show()

seed_vox = (128, 75, 70)  # found manually via plotting
mask = skimage.segmentation.flood(
    CT_data, seed_vox, tolerance=CT_data[seed_vox] / 2)
bounds = skimage.segmentation.find_boundaries(mask)
bounds_filled = skimage.morphology.binary_closing(bounds)

coords, faces = mne.surface._marching_cubes(bounds_filled, level=[1])[0]
coords = mne.transforms.apply_trans(T1.header.get_vox2ras_tkr(), coords)

brain = mne.viz.Brain(
    'sample_seeg', subjects_dir=op.join(misc_path, 'seeg'), alpha=0.8)
brain._renderer.mesh(*coords.T, faces, 'red', opacity=0.5)


# load icosahedron and fit it to the surface
# could change to 5 for smoother but slower
ico = mne.surface._get_ico_surface(grade=5)

# find large radius from coords
coords_sph = mne.transforms._cart_to_sph(coords)
r = coords_sph[0].max() * 1.5
ico_rr = ico['rr'] * r + coords.mean(axis=0)
ico_sph = mne.transforms._cart_to_sph(ico_rr)

brain._renderer.mesh(*ico_rr.T, ico['tris'], 'blue', opacity=0.5)

ico_rr_shrunk = list()
for pt, sph_pt in zip(tqdm(ico_rr), ico_sph):
    dists = np.linalg.norm(pt - coords, axis=1)
    closest = np.argmin(dists)
    ico_rr_shrunk.append(coords[closest])


ico_rr_shrunk = np.array(ico_rr_shrunk)
brain = mne.viz.Brain(
    'sample_seeg', subjects_dir=op.join(misc_path, 'seeg'), alpha=0.8)
brain._renderer.mesh(*coords.T, faces, 'red', opacity=0.5)
brain._renderer.mesh(*ico_rr_shrunk.T, ico['tris'], 'green', opacity=0.5)

# mne.write_surface(op.join(misc_path, 'seeg', 'sample_seeg',
#                           'bem', 'outer_skull.surf'),
#                   ico_rr_shrunk, ico['tris'])

gives

which would work really well without the electrodes and mounts but I think is much better and was worth sharing (also the red surface which has both the inner and outer skull is better).