Duplication of the code for EEG channel interpolation based on splines

[rbruna]

I am unsure if this should be considered a bug, so I decided to present it as a request.

Context

I am walking through the methods for the calculation of the Laplacian and I see that the spherical splines method, based in Perrin 1989, also provide bad/missing channel interpolation. It uses the function splint.

Meanwhile, the ft_channelrepair function also allows for interpolation of bad/missing channels using spherical splines. It uses the function sphericalSplineInterpolate, which is also based in Perrin 1989.

The implementation is similar in both cases, but there are two functions doing virtually the same, and I think it should be a good idea to merge them, to avoid future problems.

Differences between implementations

There are two differences between the two functions:

• splint generates, in addition to the potential, an estimate of the surface Laplacian. However, it is cheap to compute, so it should not be a problem.
• sphericalSplineInterpolate calculates the Legendre polynomials up to degree 500, stopping when the change is small.

Extra

splint calculates the Legendre polynomials using the plgndr MEX function, and sphericalSplineInterpolate calculates it explicitly in Matlab. In both cases, the polynomials are calculated for all degrees up to a maximum value.

I have created a MEX function similar to plgndr, but which: a) works for matrix inputs and b) returns all the polynomias up to i-th order, instead of only the i-th order polynomial.

Using this code would improve the efficiency of the code, with no numerical error at all. I attach the C-code here.

I offer to modify the functions for the calculation of the surface Laplacian and the channel interpolation myself, if you consider it useful.

my_plgndr.c.zip

[schoffelen]

Hi Ricardo,

Thanks for your detailed report. I agree that it is not optimal to have two different implementations of doing the same thing. I assume that you are referring to ft_scalpcurrentdensity?

So, in brief, ft_channelrepair uses the splines to estimate 'missing' data based on the distribution of voltages, but it does not return the laplacian transformed data. Ft_scalpcurrentdensity uses the laplacian transformed data from splint (which also outputs some mapping matrices for free)

Our philosophy is also that we to avoid compiled as much as possible, as long as this does not at too huge a cost in terms of computation time spent. The reason is that maintaining mex files on the different supported platforms + versions is quite challenging. Moreover, for robust functional support, we'd need a m-file based overloaded version anyhow.

If I read your comments well, then the most robust functional implementation would be to consider extending sphericalSplineInterpolate (and perhaps renaming it along the way) with an output of the surface Laplacian, and then replace the call to splint with the one-shoe-fits-all version. Perhaps the latter won't be even needed, because the 'slap' option produces the mapping matrix to the surface laplacian, correct?). I grepped the codebase and indeed there are only two occurrences where either one of the functions is called (i.e. ft_scalpcurrentdensity and ft_channelrepair).

Ideally, we'd start with a small test function that shows the overall equivalence of the computations performed in splint and sphericalSplineInterpolate. At the moment it's difficult to evaluate for me, because both function produce different outputs.

However, I noticed that the 'gx' variable from splint is (up to eps) equivalent to 'GdS' from sphericalSplineInterpolate, and that 'hx' is equivalent to -'HdS'. So that would make sense.

[schoffelen]

I added a small test function in a PR. Could you check whether my expectations are correct, in that the call to sphericalSplineInterpolate (with the method 'slap') should yield a mapping matrix that allows for the computation of the surface laplacian? If that were true, the resulting plot should yield a point cloud on the diagonal, right?

[rbruna]

Hello, Jan-Mathijs,

First of all, yes, sorry: the two different implementations are called by ft_scalpcurrentdensity and ft_channelrepair. And yes, sphericalSplineInterpolate can be asked to return the surface Laplacian, and it uses the same method than splint for it.

I have reviewed both functions and yes, as you said, gx equals Gss and Gds (Gss is used for the input electrode definition, and Gds for the requested electrode definition), and hx equals Hds (Hss is not used in sphericalSplineInterpolate).

I completely understand the avoidance of MEX files, as they are strongly dependent on the OS. I have found issues with MEX files working in one version of one OS and not in another, and compilation is not always easy. I can provide a m-file for the calculation of the Legendre polynomials which produces the same output than the previous c-file (it will be slower, bit I think it should take only some milliseconds, even for large electrode arrays).

Really I would keep the splint function, instead of sphericalSplineInterpolate, as it cleaner and easier to follow.

Last, regarding your test, there are a couple if things:

1. Your test produces the filter to calculate the surface Laplacian in sphericalSplineInterpolate, so W*V1 should be equal to L2, not to V2.
2. sphericalSplineInterpolate uses Legendre polynomials up to degree 500, but can stop before if the change between iterations is below the tolerance. So the result will not be the same, unless the tolerance is 0.
3. Furthermore, splint uses right-side division (line 94), instead of inversion (line 93, commented). For some reason, this slightly modifies the output.

So, right now, both implementations are not doing exactly the same, but almost the same...

Attached you can see the result for the surface Laplacian (comparing W*V1 with L2 in your test) and for the potential (using 'spline' instead of 'slap', and comparing W*V1 with V2). Note that sphericalSplineInterpolate dos not change the sign of the Laplacian, as Perrin does in Eq. 5b, so the diagonal goes down. Plus, there is a difference in the units, as splint uses geometrical distances and sphericalSplineInterpolate uses normalized distances.

[rbruna]

I have also noticed that ft_channelrepair checks for spherical fit in lines 360 to 371, but after that it forgets the center of the sphere. ft_scalpcurrentdensity, on the other hand, never checks for a spherical fit.

Both splint and sphericalSplineInterpolate presume both spherical fit and "centerness" in the electrode definition, and then project the electrodes over a sphere of radius 1 (lines 60 and 61 in splint, lines 33-34 in sphericalSplineInterpolate). If the electrodes have high eccentricity (for example, the 10-20 electrode definition has some degree of eccentricity) the projection over the sphere will not maintain the angles, and will be wrong.

I recommend adding a step, at the beginning of the function (or functions, if you decide to continue with the two separate functions), to center the electrode definition before projecting.

[schoffelen]

Great thanks for the quick feedback! I am totally fine with keeping splint; that would mean that it needs to output the projection matrix 'W', right?

W.r.t. the electrode centering I suggest to do that close to the computations, i.e. in splint.

Would you mind taking a shot at cleaning up the interface etc.?

[rbruna]

Of course, I am already doing it for myself, so I will prepare it in "FieldTrip format" and I will propose it as a change.

I know that, usually, you write "check" functions to check that the output is correct after a modification. In this case, specially after centering the electrodes, the result will be different. How should I proceed?

[schoffelen]

great! well if the different results are more correct than before, then it is no problem.

[rbruna]

OK, I have done it.

I have generated a new function (sphsplint) to replace both splint and sphericalSplineInterpolate. I have created a test function to check that the output is similar for both old and the new functions. I have also replaced the calls in ft_channelrepair and ft_scalpcurrentdensity

The result is here: https://github.com/fieldtrip/fieldtrip/compare/master...rbruna:issue-%231841

[schoffelen]

Thanks. Could you make a PR out of this? That makes it easier to review on our end.

[rbruna]

Of course, I was not sure if you wanted me to create the pull request before you had the opportunity to review it.

I have done it now.

[rbruna]

Well, great news, I have found a (quite big) bug in the code for the surface Laplacian...

Should I open another report, or should I expose it here?

[schoffelen]

You mean in code that is different from the one which will be there once I reviewed, and merged the PR?

[schoffelen]

If it's different, I'd prefer a new issue

[rbruna]

It is in the Hjorth method, so it is in the ft_scalpcurrentdensity file (which is modified in the PR).

[schoffelen]

OK, then I suggest to fix it in the PR and document it there (or in this issue)

[rbruna]

OK, here it goes.

I tested the three different methods for the calculation of the surface Laplacian. One would expect the results for Hjorth and finite differences to be similar, and the result for spherical splines pointing in a similar direction. However, I found that Hjorth was different from the other two. Moreover, it looks like the reconstructed signals are similar, but scrambled around the scalp.

These are the evoked potentials and evoked surface Lapacians for some data:

When I change the order of the channels so the data fits the electrode (and neighbor) definition, this is what I get:

This is as expected, the results for all three methods are somehow similar, and those for the Hjorth and finite differences methods look much alike.

Looking through the code, the problem seems to be that Hjorth method only uses the neighbor definition to build the transformation matrix (ft_scalpcurrentdensity.m, lines 224 to 250), so the output data is sorted as in the neighbor definition.

  % convert the neighbourhood structure into a montage
  labelnew = {};
  labelold = {};
  for i=1:length(cfg.neighbours)
    labelnew = cat(2, labelnew, cfg.neighbours(i).label);
    labelold = cat(2, labelold, cfg.neighbours(i).neighblabel(:)');
  end
  labelold = cat(2, labelnew, labelold);
  labelold = unique(labelold);
  tra = zeros(length(labelnew), length(labelold));
  for i=1:length(cfg.neighbours)
    thischan   = match_str(labelold, cfg.neighbours(i).label);
    thisneighb = match_str(labelold, cfg.neighbours(i).neighblabel);
    tra(i, thischan) = 1;
    tra(i, thisneighb) = -1/length(thisneighb);
  end
  % combine it in a montage
  montage.tra = tra;
  montage.labelold = labelold;
  montage.labelnew = labelnew;

But in 275 the code changes the labels back to the original ones:

% collect the results
scd.elec    = elec;
scd.time    = data.time;
scd.label   = data.label;

After this, the labels get scrambled.

There are two possible fits:

1. Change the code so the output of the montage is sorted according the input data.
2. Removing line 275, so the labels are not overwritten.

The second option would be easier, but it requires to change the code for the spherical splines, as the scd variable is initialized to contain only the trial field, nothing else.

What do you think?

[schoffelen]

Good catch w.r.t. the channel order. I agree with removing ling 275 (and 274), because that does not make sense. the output of ft_apply_montage(montage, data, ...) should contain a consistent data structure (at least with a time and label field that are consistent with the data in trial). Perhaps only the elec needs to be added.

Indeed, I would not go for your proposed solution 1 .

But as far as I can see in the PR, you already made this fix, right? So we're good to go and I can merge?

[rbruna]

Indeed! As you said, the output of ft_applly_montage has all the fields, except maybe for the electrode definition. So I did exactly what you suggested.

I think everything should be correct and we can merge :)

[schoffelen]

OK, merged! thanks a lot!

[rbruna]

A pleasure :)

[schoffelen]

Hi @rbruna, just pinging you here that our suite of test functions caught a few errors due to the latest PR concerning the spherical splines. I looked at it, and think I fixed it, just to let you know. Along the way I did some cosmetics in sphsplint, where the most important change was that I needed to make a change in the initialization of a 'ones(... ,1)' vector, which crashed in case there were fewer 'goodchans' than 'allchans'. I think I did it correctly :)

[rbruna]

Hello, Jan-Mathijs,

First of all, sorry for the spaces before and after the parentheses! I prefer my code to have more spaces, as it is easy for me to spot the variables. But I should pay more attention to FieldTrip style :)

Second, indeed, originally I had two variables, siz1 and siz2, for the two sets of electrodes, but I lost the train of thought and thought that I only needed siz1. If I were to fix it, I would have added a siz2 variable after siz1, and then used it in line 126 (now 124). But your fix should indeed work nicely :)

I am glad it worked out, in the end. Tell me if there is some other problem.

And thank you (all of you) for your great work building this.

[schoffelen]

No worries Ricardo, And indeed I changed the spacing because Robert is quite keen on these kind of things :). I am just his lieutenant. Thank you for your high quality contributions. Keep them coming in the future!

[robertoostenveld]

thanks for the good work and detailled reporting! Much appreciated.