ENH: Add offline homogeneous field correction for OPMs

[larsoner]

Should be simple enough to implement, and should make data for OPMs cleaner, at least for whole-head arrays / multi-axis measurements

https://www.sciencedirect.com/science/article/pii/S1053811921007576

[larsoner]

Looks like I already had this issue open @neurofractal @georgeoneill , so let's discuss here.

This is great news to hear [about field correction]! As you can see we've been doing some initial tests. I've been trying to work out whats the best approach for MNE-Python. We normally do a hard correction to the data (rather than having an online montage) and then update our 'compensation matrix' to update our forward modelling later. But I realise this is not in the nature of how the SSP filters work in MNE-Python (i.e. only activated when needed). This might need some thought.

Yes we'll have to be a bit careful. I have not looked at the math at all but doubt they are true orthogonal projection vectors, so wouldn't be good to call SSP; and I doubt they work like compensation matrices for CTF (which have some nice mathematical property -- like its inverse is its transpose or something like this, I'd have to check to remind myself). We can probably add a new constant for it -- we've wanted one before for the cross-talk correction for MEG data anyway, which is itself just a matrix multiply (that does not modify/reduce rank).

So I say in a second PR after #11257 is done let's get it doing what it needs to do -- probably as a matrix that can be applied (or not) just like a projection -- and then I can talk to some FIF folks about the right way to store the information.

I think the answer to both of the following is "yes", but just to check -- can the higher-order version also be described as just a matrix multiply? And they are not time-varying (i.e., just depend on sensor geometry, not sensor position as a function of time), right?

[georgeoneill]

Hello @larsoner, starting to think about how to best to tackle this, now we have an importer working.

One potential issue we have going forward (other than requiring new classes to represent arbitraty sensor balancing) is the licensing of the original code which generated the basis sets. A simple relicensing and port of spm_opm_vslm.m (which is what SPM uses to generate the regular vector spherical harmonics) is a no-go unfortunately.

However, there's probably a few of BSD-compatible solutions.

1) Build upon mne.preprocessing.maxwell_filter module basis set (specifically the outer space). 2) Use the bfieldtools module, which has functions specifically for geneating vector spherical harmonics. I think this approach means we don't need to do any cartesian to spherical coordinate transforms too. 3) Simply just implement the original Homogenous Field approach where the basis set are the sensor orientations, lose access to gradient, quadratic terms.

Which of these sounds the most appealing?

[larsoner]

Let's use (1) since we should already have all the code we need in MNE-Python already that way

[georgeoneill]

I thought this might be the answer! I'll do some exploration.