BUG: Forward code should raise an error if sensors are inside BEM

[jjnurminen]

Hi,

I get strange BEM results with sample using the precomputed BEM. These manifest as erroneous/unexpected field values, mostly in regions below the brain. The problems seem to go away if I recreate the BEM in mne-python.

To demonstrate, I created an array of 1000 pointlike magnetometers uniformly spread around the brain:

I then computed the forward fields (for one particular source; doesn't really matter which one) using the precomputed BEM and a freshly created one:

Note that the sensor index increases with decreasing z coordinate. It looks like the precomputed BEM starts to increasingly deviate from the newly created one as we go lower in the coordinate system. In addition to this, the values of precomputed BEM for some sensors are really strange (the spikes in the plot).

The results from a spherical model, as well as the mne-python freshly created BEM, look fine.

To recreate, please download this sensor description for the spherical array and unzip it: pointlike.zip

Then run the code below:

import mne
import pathlib
import matplotlib.pyplot as plt
import pickle

data_path = pathlib.Path(mne.datasets.sample.data_path())
subjects_dir = data_path / 'subjects'
subject = 'sample'
head_mri_trans_file = data_path / 'MEG/sample/sample_audvis_raw-trans.fif'
bem_file = subjects_dir / 'sample/bem/sample-5120-5120-5120-bem-sol.fif'
raw_file = data_path / 'MEG/sample/sample_audvis_raw.fif'
info = mne.io.read_info(raw_file)

# replace 'info' with a custom one to get pointlike sensors
alt_array_name = 'pointlike.dat'
with open(alt_array_name, 'rb') as f:
    info = pickle.load(f)

# source space
src_cort = mne.setup_source_space(
    subject, spacing='oct6', subjects_dir=subjects_dir, add_dist=False
)

# precomputed BEM
fwd_bem_precomp = mne.make_forward_solution(
    info,
    head_mri_trans_file,
    src_cort,
    str(bem_file),
    eeg=False,
)
lead_bem_precomp = fwd_bem_precomp['sol']['data']

# recompute BEM
model = mne.make_bem_model(
    subject=subject, ico=4, conductivity=[0.3], subjects_dir=subjects_dir
)
bem = mne.make_bem_solution(model)

fwd_bem = mne.make_forward_solution(
    info,
    head_mri_trans_file,
    src_cort,
    bem,
    eeg=False,
)
lead_bem = fwd_bem['sol']['data']

sphere = mne.make_sphere_model(head_radius=None)
fwd_sph = mne.make_forward_solution(
    info,
    head_mri_trans_file,
    src_cort,
    sphere,
    eeg=False,
)
lead_sph = fwd_sph['sol']['data']

# plot
src_ind = 3000  # some random source
plt.figure()
# plt.plot(lead_sph[:, src_ind], label='sphere model forward')
plt.plot(lead_bem[:, src_ind], label='BEM model')
# plt.plot(lead_bem_precomp[:, src_ind], label='BEM model, precomputed')
plt.plot(
    abs(lead_bem[:, src_ind] - lead_bem_precomp[:, src_ind]),
    label='Precomputed vs fresh BEM, absolute difference',
)
plt.ylim([-0.5e-5, 0.5e-5])
plt.legend()
plt.xlabel('Sensor index')
plt.ylabel('Field amplitude')

[agramfort]

is this a pb of numerics that change between systems?

[jjnurminen]

@agramfort maybe, but I find it hard to understand why the precomputed BEM blows up pretty much completely for some field points.

[agramfort]

maybe the precomputed BEM was computed by the C code that uses float32?

[jjnurminen]

To fix sample, we can of course just create and save new BEM models using mne-python. However some mne-python docs still suggest creating BEM models with the command line tools, so we might want to make sure they don't produce erroneous models.

BTW, this issue would probably go unnoticed in many cases due to 1) finite-area sensors causing the field values to be spatially averaged, reducing the error, and 2) the fact that MEG systems usually don't have sensors below the neck line. Still, it was bad enough to completely screw up some simulations I was doing :)

[agramfort]

To fix sample, we can of course just create and save new BEM models using mne-python. However some mne-python docs still suggest creating BEM models with the command line tools,

where?

so we might want to make sure they don't produce erroneous models.

is it really erroneous or just numerics compatibility issues?

BTW, this issue would probably go unnoticed by normal sensor setups due to 1) finite-area sensors cause the field values to be spatially averaged and reduce the error and 2) the fact that MEG systems usually don't have sensors below the neck line.

ok

[jjnurminen]

To fix sample, we can of course just create and save new BEM models using mne-python. However some mne-python docs still suggest creating BEM models with the command line tools, where?

At least https://mne.tools/stable/auto_tutorials/forward/30_forward.html says:

"Computing the BEM surfaces requires FreeSurfer and makes use of the command-line tools mne watershed_bem or mne flash_bem, or the related functions mne.bem.make_watershed_bem() or mne.bem.make_flash_bem()"

is it really erroneous or just numerics compatibility issues?

Well, if the BEM example data gives results that are way off, I'd say it's erroneous. Whether it's due to numerics compability or something else is beside the point (IMO).

[jjnurminen]

Of course, one might argue that the pointlike array is just a corner case. However even with the default 306-channel array, I get 14% mean difference* between the precomputed and mne-python computed BEMs.

*according to following definition:

np.abs(lead_bem - lead_bem_precomp).mean() / np.abs(lead_bem_precomp).mean()

[MattiStenroos]

Hi MNEers!

Before going for actual bugs or smaller differences due to different numerical libraries, let's look at the model... The precomputed MNE model is a 3-shell model, and based on what I understand from Jussi's script / its output, the new one is a single-shell model. For differences between a single- and three-shell model in a normal MEG sensor setting, please see (Stenroos et al, NeuroImage 2014). I haven't looked at MEG fields in this kind of closed sensor layout, but based on the geometry differences in head models the largest differences will be at inferior sensors. So this explains part of the difference.

Part of Jussi's inferior sensors are inside the head in the three-shell model, but outside the single-shell model. As the "helmet" intersects the scalp, some sensor may be very close to the computation points in the scalp mesh, causing potentially major numerical trouble.

So I would first move the sensors outside the scalp (say, at least 1/2 triangle side-length outside the scalp).

Hope this helps! Matti

[jjnurminen]

Thanks @MattiStenroos! I feel dumb now. Somehow the fact that the precomputed model is a three-shell one completely escaped me.

I guess I was somewhat misled by https://mne.tools/stable/auto_tutorials/forward/30_forward.html, which uses the same dataset and says:

"Here we’ll assume it’s already computed. It takes a few minutes per subject. For EEG we use 3 layers (inner skull, outer skull, and skin) while for MEG 1 layer (inner skull) is enough."

So I just assumed that the precomputed model is single shell.

[MattiStenroos]

Don't worry, you are not the first one to get confused with relationships between different head models and sensor or source locations...

I think it would always be a good idea to plot at least the scalp mesh and the sensor locations --- maybe add that to some MNE basic example scripts, if it is not there yet? As coregistration sanity-check, I would also check the distance from scalp to sensors. With head movement correction, I would additionally check that the reconstruction helmet is well outside the scalp.

[agramfort]

this tells me that the forward code should crash if it detects some channels inside the BEM surfaces...

also you say it recommends to use the command line but the command line you refer to is just python code so the only difference could be an IO roundtrip.