Maxwell filter reconstruction

[harrisonritz]

Right now maxwell_filter always reconstructs bad channels and then removes their bad label.

I think there's going to be a higher prevalence of bad channels for OPM than SQUID (we're often missing a few). I've had some issues with interpolation, especially for channels near the edge of the helmet (see example; problem channel in top-right was not recording during the session and labelled bad. In generally, missing channels are reconstructed to have high broadband power).

There are some SSS methods that might be better suited for OPM (eg AMM or iterative SSS). But I wonder if it would make sense to make bads resetting optional. For example, could easily add an argument that controls whether _reset_meg_bads is called, interpolating the bad channels but keeping them labelled bad (which would allow the user to change their mind later, I guess).

HFC is another option for OPM for OPM, but it would be nice to support something a little stronger for systems with higher channel counts (eg I find that SSS has better noise suppression).

So: is SSS reconstruction a plausible issue for modalities that are frequently missing channels and does it make sense to give users control over whether bad channels are reset after SSS?

Any other recommendations here would be appreciated as well, even if user control over bads resetting makes sense. Maybe I should be dropping these channels instead?

[welcome[bot]]

Hello! 👋 Thanks for opening your first issue here! ❤️ We will try to get back to you soon. 🚴

[larsoner]

So the current workaround would require two lines like:

bads = raw.info["bads"]  # workaround line 1
raw = maxwell_filter(...)  # the call itself
raw.info["bads"] = bads  # workaround line 2

or even one in the case where you don't overwrite the variable

raw_sss = maxwell_filter(...)
raw_sss.info["bads"] = raw.info["bads"]

Have you considered this?

I'm a bit torn because this workaround doesn't seem too onerous. But it's also not much maintenance burden for us to add a kwarg, document it, and nest the reset call in a conditional, either. So if you're motivated to open a PR @harrisonritz I can live with a new argument!

[harrisonritz]

Ya the fix is easy. Nice for mne_bids_pipeline for there to be a built-in option, but I've added options like this in my fork.

I guess it's a bit of a science question too. For someone new to MEG, full reconstruction of missing sensors from SSS seems like a strong assumption. Is it safe enough to be the default?

I guess the point of SSS is that you are already reconstructing the good sensors. In SQUID, you usually won't have tons of missing sensors. But this could be an issue with OPM where there are fewer sensors (so SSS is less stable) and more missing sensors (so more opportunities for interpolation issues).

Eg, I wonder if it's worth tracking the expected quality of the reconstruction, eg the condition number of the bases set S (re: Holmes 2023), and then modifying the reconstruction appropriately (eg don't recon bad sensors, regularize, or throw an error or warning).

Is this something that I should just keep an eye on in my pipeline, or do you think that there should there be more guard rails here? Maybe I was just overconfident about using SSS with OPM.

[larsoner]

I haven't used SSS much with OPMs yet. Condition number would indeed be good to check somehow. Also the quality / appropriateness of the sphere fit based on expansion origin, head geometry (more like an oblate spheroid), and sensor distance. This is in theory something we could check and quantify... like given a source space, sensor alignment, and SSS origin, what are the radii of 1) the biggest sphere with that origin that doesn't collide with any sensors (or really, has the nearest sensor on its surface), and 2) the smallest sphere with that origin that encompasses all brain sources (or really, has the farthest source point from the origin on its surface). For SQUID (1) will be bigger than (2) and we can quantify how much bigger it is, and for OPM there might be some cases where (2) is actually bigger than (1) in which case SSS shouldn't be used until we extend it to a different basis set (e.g., ellipsoidal harmonics or similar).

To get the ball rolling on this, we could in MNE-BIDS-Pipeline for example if SSS and source modeling are both enabled for a given config, we could compute these two radii and log them during the source step and put them in the report somewhere. We could even make a nice visualization out of it but just having these values be computed and reported would be a great start...

You could do the same with logging the condition number in theory. It would be slightly trickier because we'd either have to capture the SSS logging and extract it, or recompute it afterward (and for movement compensation it's time-varying so yikes!). So maybe we should leave that alone for now perhaps.

[harrisonritz]

I like this, let me chew on it for a bit in my 'free time'. Good idea to generate condition where SSS is/isn't appropriate, and see how close we get to the edge.

I've been poking around a little bit on porting AMM or iterative SSS (IIUC, these haven't been compared). But, as you identified, the hardest part has been that maxwell_filter processes in chunks. Let's talk if/when we get there, some pointers might be helpful (I think I just need to work backwards from get_decomp).