mne.preprocessing_current_source_density gives wrong CSD estimates

[mservantufc]

Describe the bug

Attempt to compute CSD on evoked object called "monopolar_av_ref.fif" (see attached). The evoked object contains grand-averaged EEG waveforms (64 EEG channels (averaged reference), 'biosemi64' montage + 1 EOG channel).

Steps to reproduce

GA_comp_corr_stim = mne.read_evokeds('monopolar_av_ref.fif') GA_comp_corr_stim_CSD = mne.preprocessing.compute_current_source_density(GA_comp_corr_stim)#should be similar to the attached CSD.fif file

plots obtained: times = np.arange(.3, 1.,.05) fig1 = GA_comp_corr_stim.plot_topomap(times = times ,ch_type = 'eeg', cmap = 'viridis', nrows = 2, average = .03, time_unit='s') fig3 = GA_comp_corr_stim_CSD.plot_topomap(times = times ,ch_type = 'csd', cmap = 'viridis', nrows = 2, average = .03, time_unit='s')

Actual results

plots "topo_CSD.png" and "topo_monopolar.png" in the attached. It is clear that CSD estimates are wrong, since the spatial selectivity appears lower than the spatial selectivity of the monopolar data, and differ from those obtained with a python implemenbtation of Mike Cohen's matlab code ('topo_CSD_Cohen.png') Maybe it's a problem with the sphere= auto argument that returns: Fitted sphere radius: 95.0 mm Origin head coordinates: 0.0 -0.0 40.1 mm Origin device coordinates: 0.0 -0.0 40.1 mm

Additional information

Platform: Windows-10-10.0.18362-SP0 Python: 3.7.6 (default, Jan 8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)] Executable: C:\Users\mservant\Anaconda3\pythonw.exe CPU: Intel64 Family 6 Model 158 Stepping 10, GenuineIntel: 12 cores Memory: 15.8 GB

mne: 0.20.5 numpy: 1.18.1 {blas=mkl_rt, lapack=mkl_rt} scipy: 1.4.1 matplotlib: 3.1.3 {backend=Qt5Agg}

sklearn: 0.22.1 numba: 0.48.0 nibabel: 3.0.1

cupy: Not found pandas: 1.0.1 dipy: 1.1.1 mayavi: 4.7.1 {qt_api=pyqt5, PyQt5=5.9.2} pyvista: Not found vtk: 8.1.2

mne.zip

[agramfort]

@alexrockhill can you have look?

[alexrockhill]

Can you link to the referenced code converted from Mike Cohen's MATLAB code?

I'm just copying the minimally reproducible example with the imports formatted here for convenience:

import numpy as np
import mne

GA_comp_corr_stim = mne.read_evokeds('monopolar_av_ref.fif')[0]
GA_comp_corr_stim_CSD = mne.preprocessing.compute_current_source_density(GA_comp_corr_stim, sphere=(0, 0, 0, 85))  #should be similar to the attached CSD.fif file

times = np.arange(.3, 1.,.05)
fig1 = GA_comp_corr_stim.plot_topomap(times = times ,ch_type = 'eeg', cmap = 'viridis', nrows = 2, average = .03, time_unit='s')
fig3 = GA_comp_corr_stim_CSD.plot_topomap(times = times ,ch_type = 'csd', cmap = 'viridis', nrows = 2, average = .03, time_unit='s')

[mservantufc]

Hi Alex, It's a Python translation of Mike's code (attached). Please note that this code is only applicable to epochs, so you'll not be able to reproduce. But you can see the computations performed. CSD.zip .

[alexrockhill]

Did you look at whether the same stiffness and lambda2 arguments are used between the scripts? As you can see in https://mne.tools/dev/auto_examples/preprocessing/plot_eeg_csd.html, the parametrization makes quite a difference.

[mservantufc]

De: "Alex Rockhill" notifications@github.com À: "mne-tools/mne-python" mne-python@noreply.github.com Cc: "MATHIEU SERVANT" mathieu.servant@univ-fcomte.fr, "Author" author@noreply.github.com Envoyé: Mercredi 3 Juin 2020 17:05:04 Objet: Re: [mne-tools/mne-python] mne.preprocessing_current_source_density gives wrong CSD estimates (#7862)

Did you look at whether the same stiffness and lambda2 arguments are used between the scripts? As you can see in [ https://mne.tools/dev/auto_examples/preprocessing/plot_eeg_csd.html | https://mne.tools/dev/auto_examples/preprocessing/plot_eeg_csd.html ] , the parametrization makes quite a difference.

— You are receiving this because you authored the thread. Reply to this email directly, [ https://github.com/mne-tools/mne-python/issues/7862#issuecomment-638257531 | view it on GitHub ] , or [ https://github.com/notifications/unsubscribe-auth/AKYSBODXDOU5DBNRF7G6MRLRUZRCBANCNFSM4NRPGX7Q | unsubscribe ] .

[mservantufc]

Yes there are exactly the same. Regardless of differences between codes, a CSD transform with stiffness = 4 should provide better spatial selectivity than an average reference, which is not the case with the mne.preprocessing.current_source_density function. The github depository of the python conversion of Mike 's code is https://github.com/alberto-ara/Surface-Laplacian

[alexrockhill]

I ran the example from the github you linked and there csd from mne fails and gives a runtime warning that the head is too big which I assume you are right. I'll check the steps in between first but I assume it's something we are doing wrong with handling the head size as you suggested.

[alexrockhill]

Hmmm after looking more it seems there is a units issue with the data from https://github.com/alberto-ara/Surface-Laplacian where it is provided in mm instead of m, but that can be solved easily. I think the bigger issue is that G and H are calculated differently. @dengemann, you posted the nice-tools version from here https://github.com/nice-tools/pycsd/blob/master/pycsd/csd.py that I used in the implementation. There isn't any calculation of euclidean distance between the electrodes like in

# compute cousine distance between all pairs of electrodes
cosdist = np.zeros((numelectrodes, numelectrodes))
for i in range(numelectrodes):
    for j in range(i+1,numelectrodes):
        cosdist[i,j] = 1 - (((x[i] - x[j])**2 + (y[i] - y[j])**2 + (z[i] - z[j])**2)/2)

I'll look into how nice-tools does it but any idea what the issue is?

[alexrockhill]

I think maybe the issue is with the nice-tool implementation because this matches that exactly as shown in the tests (mne.preprocessing.tests.test_csd.py) maybe @fraimondo can comment?

[fraimondo]

The implementation of pycsd is ported from the CSD toolbox here.

A few months ago, something changed with the MNE montage reading and somehow the positions in the CSD file were modified after reading.

I ended up undoing the transform: https://github.com/nice-tools/pycsd/blob/fd79ec254f39fcf1d979e3a7d3a7b586147bd381/pycsd/csd.py#L40-L42

Then I guess I figured out it was something to do with the head size and took it back (memory is fuzzy now).

However, the tests accuracy were not the same after MNE 0.20. I went from a tolerance of 1e-16 to 1e-5.

[larsoner]

The way the code is in master, it is not calculating a cosine similarity, which is the dot product of two vectors divided by their norms -- we instead divided by some assumed radius, which was not guaranteed to be correct (if the radius was wrong, or varies vector-by-vector, it was no longer computing a cosine similarity).

One difference between Mike's code and Jurgen's (posted in #7863) is that Jurgen projects to a unit sphere with [X,Y,Z] = sph2cart(ThetaRad,PhiRad,1.0), whereas Mike just divides by the max radius. I think the former is safer than the latter -- especially for realistic head positions -- but shouldn't matter in any case when you use a truly spherical montage. (It also seems closer to the definition of cosine similarity, because you are guaranteed to divide the vectors by their norms). So that's what we do over in #7863. We don't use the same math to compute the cosine similarity, but numerically they should be the same (I've now checked finally!).

In Jurgen's code, the radius only ends up being a scaling factor for the data -- basically as the radius shrinks the result will grow (as the inverse square of the radius) which makes some sense because the units are V/m^2, and as the sphere shrinks you need higher and higher current source density to explain the data. But otherwise the result will not change (anymore), which is nice.

@fraimondo @mservantufc can you try #7862, and also look at the PR vs master doc outputs and let me know if it makes sense and works on your data?

[fraimondo]

that Jurgen projects to a unit sphere with [X,Y,Z] = sph2cart(ThetaRad,PhiRad,1.0)

I think this is not a projection but just the transformation of the coordinate system. The 1.0 there is because the locations are read from a CSD file in which the projection was already done. In other words, we don't use the digitised or standard montage positions as in any EEG analysis, we have to use a specific montage file in which the coordinates are projected to a sphere with r=1.

In any case, the pycsd code should be a port of the CSD Toolbox code, which does not mean it's the best approach so far. At that time, with @dengemann we needed to have something working in python so we created that, respecting the original license.

[larsoner]

I think this is not a projection but just the transformation of the coordinate system. The 1.0 there is because the locations are read from a CSD file in which the projection was already done.

Looking a bit, I see where the CSD has been computed/applied at that point in the code -- it looks like it's done as the CSD transform itself is being computed. Basically the code makes it so that it does not matter what the actual montage radius was -- it sets all electrode radii effectively to 1. (projects them to the unit sphere). Mike's code also does this by dividing by the maximal radius (at least if all radii are fixed / spherical montage it's equivalent). These two implementations plus the idea that the code is supposed to compute / operate on a cosine similarity (based on naming of vars) all point to this as being the right thing to do. That's enough evidence for me to not bother reading the original papers :)

[mservantufc]

@larsoner how can I update my mne in a clean way to test the latest version of mne.preprocessing.current_source_density ?

[larsoner]

If you plan on only doing this once in a while (a few times a year?) you can do the following to test Alex's branch from #7863 with:

pip install --upgrade --no-deps https://github.com/alexrockhill/mne-python/archive/csd2.zip

If you want to be able to switch more easily, you can do for example:

git clone https://github.com/mne-tools/mne-python.git
cd mne-python
git remote rename origin upstream
git remote add alexrockhill https://github.com/alexrockhill/mne-python.git
git fetch alexrockhill
git checkout -b csd2 alexrockhill/csd2
pip install -ve .

Then later when we merge into master you can do in the mne-python directory:

git checkout master
git pull

and you'll be up to date with latest mne-tools/mne-python master branch.

[mservantufc]

@larsoner @alexrockhill Just tested the updated version of mne.preprocessing.current_source_density function on my EEG data posted at the beginning of this issue. It now gives results that very similar though not identical to Alberto's Python translation of Mike's code (see attached). How do you guys explain these small differences ? One remaining (potential) issue concerns the units. There seems to be a discrepancy between units of the topomaps of and units of the data. What are the units of the CSD data returned by mne.preprocessing.current_source_density? microvolt/cm^2? I am confused because the topomap function returns V/m^2

test_csd.zip

[larsoner]

The units should be V/m^2.

Mike's code only normalizes/computes cosine similarity when the radii are all equal for the electrodes, is that the case for your data?

[larsoner]

(Also his code does not normalize by the head radius while ours does, but this will just be a constant scale factor. To get rid of it pass sphere=(0,0,0,1) or so)

[mservantufc]

@larsoner I used mne's 'biosemi64' montage so I think so. Ok so topomap returns mV/m2 as units. Maybe it would be good to have consistent V/m2 everywhere. Thanks guys, great work!

[larsoner]

Ok so topomap returns mV/m2 as units.

he computations and data are in V/m^2. Like everything in MNE (or so we try), data is in SI units, but could be scaled to some non-SI unit (using milli/centi/micro/nano/etc.) for display -- here the topomaps display in mV/m^2 because this is a "nice" scaling. In general we try to pick units for display that are human-readable, like uV for EEG data and fT for magnetometers. We could display V/m^2 here but the scale numbers would be uglier.

I used mne's 'biosemi64' montage so I think so

If I do like Mike's code in @alexrockhill's tests:

    pos /= np.linalg.norm(pos, axis=-1).max()
    from scipy.spatial.distance import squareform, pdist
    cos_dist = 1 - squareform(pdist(pos, 'sqeuclidean')) / 2.

instead of what we do (project all points to sphere):

    pos /= np.linalg.norm(pos, axis=1, keepdims=True)

I get results that match to 5 orders of magnitude closer -- rather than to roughly a few percent, it's like 0.0001%. So indeed some small differences (for spherical montages) are expected but they shouldn't be drastic. For non-spherical / realistic positions, they will diverge quite a bit. But I don't think dividing by the max is the right approach in either case, assuming it is supposed to be the cosine similarity / computations done on the sphere surface.

I'll go back to the original paper to double check if this is what is supposed to be done to be sure...

[mservantufc]

Thanks, I understand the difference now.

[larsoner]

Looking through the papers a bit, these computations depend on spherical splines, so indeed I think it's important the electrodes are (effectively) on a unit sphere, at least the way we have it implemented.

[larsoner]

I think this was fixed by #7863 but let's reopen if not