I've been debugging my code by comparing an EEG forward model that I
constructed using MNE python to one constructed by a former labmate using
MNE C. The lead field matrices are different, even when I use identical
source space, BEM solution, and transformation file.
I've narrowed down the problem to info['chs'][...]['eeg_loc'], which in the
old data is 3x2, with the second column having the location of the
reference electrode. In my data, eeg_loc is 3x1, with no information about
the location of the reference electrode. The information in this second
column appears to be used by _bem_specify_els, which loops over the columns
of the electrode location matrix.
I've played around a bit with add_reference_channels and set_eeg_reference,
but neither of them changes the second column in eeg_loc for the
non-reference electrodes. The location of the reference electrode is
specified in my montage file (hpts), but there doesn't seem to be a way to
apply it to the eeg_loc of all of the other electrodes.
This is what I thought might have worked, but it seems to just add another
electrode, without influencing the chs of the other electrodes:
montage = mne.channels.read_montage(montage_hpts, transform=True, unit='mm')
raw = mne.io.read_raw_fif(rawfile,preload=True)
mne.io.add_reference_channels(raw,'EEG000',copy=False)
raw.set_montage(montage)
mne.io.set_eeg_reference(raw,['EEG000'],copy=False)
Yes. Well, in the old data the average reference projector's active
property is set to 'on', and in my data it is set to 'off'.
I tried running mne.io.proj.activate_proj to make it active in the new
data, and it did not affect the lead field matrix. Is there some other way
to activate the proj?
the `set_eeg_reference` and `add_reference_channels` are functions that transform the EEG data. They allow you to specify different referencing schemes, which is useful for sensor level analysis. They also set the `custom_ref_applied` flag in the info structure to explicitly override the average reference projector. Since you are doing source localisation and therefore want the average reference, don?t use these functions, they don?t do what you wanted them to do.
Let?s move on to your problem:
I was always wondering why `eeg_loc` was a 3x2 matrix? thank you for tracking it down. I think many of our data import routines only set the first row of `eeg_loc`; this should be fixed. However, where in `_bem_specify_els` does it loop over the `eeg_loc` matrix exactly? I can?t find it...
In line 141, _bem_specify_els loops over the columns of el_r, which is a
transformed version of els[...][rmag].
els comes from the coils variable in _prep_field_computation, which come
from fwd_data from _compute_forwards, yada yada yada...
If you trace it all the way to its source, the field rmag is set through
make_forward_solution > _prep_channels > _create_coils > _create_eeg_el,
which sets it to a (potentially transformed) version of ch[eeg_loc].
The reason it matters is that in _bem_specify_els, the sol variable is
constructed as a weighted sum of the bem solution at tri vertices near the
locations in the columns of rmag. The weight w on the second column (from
elw, also set in _create_eeg_el) is -1. So the bem is subtracting off the
influence of the reference electrode.
Apologies if this is vague -- I'm still familiarizing myself with the
theory. At this point I have a better grasp on the code, since I've been
stepping through it debugging.
I think the code works properly if the eeg_loc is set to the correct 3x2
matrix: I reran the forward solution for the old data using MNE python and
I get basically the same answer as MNE C (differences on the order of
10^-3, when the lead field matrix entries are on the order of 100). So MNE
python is loading the data properly and using it properly.
The problem is that I can't create a Raw object using *only* MNE python
that has the supposedly correct 3x2 eeg_loc parameters. I think the
solution would involve adding the reference electrode to the Montage or the
DigMontage, since that's where the information is (at least in my data). In
case it matters, I'm getting my Montage from a hpts file that labels the
reference as eeg electrode 0. Since that electrode doesn't exist in the
recordings, it isn't used by set_montage. That's why I thought it might
work to add the reference electrode using add_reference_channels before
running set_montage.
One more thing: the only place in the documentation that talks about the
reference channel location being stored with eeg channel locations is here: http://martinos.org/mne/stable/manual/c_reference.html
under the entries for "ch_pos" in the matlab structures.
I'm assuming ch_pos in the matlab structures serves the same purpose as
eeg_loc in the python info dicts.
why do you the location of a reference channel?
if you average reference you don't need any extra location.
You just need the location of your EEG channels for which you have data.
I'm still hazy on that, actually. What I know is that the answer is
different if the data has the reference electrode location vs if it
doesn't. Based on where it is used in the code, it seems to be affecting
how the BEM model is used in the forward model computation.
Also, I'm starting from the assumption that the forward model constructed
using MNE C by my former labmate is correct. I'm trying to replicate that
result using MNE python. It's possible that the 3x2 representation of the
eeg locations from the MNE C constructed files was a mistake either in the
code or in my labmate's usage of the software.
It's probably best to use e.g. DropBox if you can so other devs not at
Martinos can also have a look. But please, as usual, make sure your data
are appropriately anonymized first before sharing.
I'll also email you privately so we can try to iterate on this a bit, and
report back to the email list what we figure out.