Colocated topomaps for OPM

[harrisonritz]

Reference issue (if any)

Adds features from #11579

What does this implement/fix?

This PR extends topomaps to accommodate OPM data with colocated sensors.

It modifies the code used for colocated fNIRS, which now handles both fnirs and OPM.
In contrast to fNIRS which averages the colocated sensors, for OPM sensors we just take the most radial sensor.

Additional information

• I've written an opm test, but would like guidance about how to access mne.datasets.ucl_opm_auditory. Test won't work as is (eg can't update UCL dataset, and wanted to confirm proper method), but the test is structured on the fnirs overlap test so probably won't need much more
• I rely on some hard-coded look-ups:
  • lookup OPM sensor coil types (quspin currently included)
  • I need to use different loc triplets for dual-axis vs triaxial sensors, based on comparing ucl_opm_auditory vs Princeton's triaxial opms. At UCL, they use *-RAD and *-TAN to indicate the two axes, and at Princeton we use * X or * [X].
• for the UCL dataset, most of the *-RAD channels are the most radial, but in some cases both RAD and TAN have low radial scores, and TAN is the most radial sensors

Attaching my internal testing code to validate against picking RAD channels.
Happy to share Princeton dataset if that's helpful.

# %%
import numpy as np
import mne


opm_file = (
mne.datasets.ucl_opm_auditory.data_path() / "sub-002" / "ses-001" / "meg" / "sub-002_ses-001_task-aef_run-001_meg.bin"
)
raw = mne.io.read_raw_fil(opm_file, verbose="error")
raw.crop(100, 300).load_data()  # crop for speed
rad_string = "^.*-RAD$"
picks = mne.pick_types(raw.info, meg=True)


# %% preprocess --------

raw.notch_filter(50, picks=picks)
raw.filter(1, 40, picks=picks)

projs = mne.preprocessing.compute_proj_hfc(raw.info, order=2)
raw.add_proj(projs).apply_proj(verbose="error")


# %%  epoch data --------

events = mne.find_events(raw, min_duration=0.1)
epochs = mne.Epochs(
    raw, events, tmin=-0.1, tmax=0.4, baseline=(-0.1, 0.0), verbose="error",event_repeated='merge',
)
evoked = epochs.average()
t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]


# %% plot evoked topomap --------

fig = evoked.plot_topomap(times=t_peak, ch_type="mag", time_unit="s")
print('len axes: ', len(fig.axes))

rad_picks = mne.pick_channels_regexp(evoked.info.ch_names, rad_string)
fig = evoked.copy().pick(rad_picks).plot_topomap(times=t_peak, ch_type="mag", time_unit="s")
print('len axes: ', len(fig.axes))


# %% test ICA --------

ica = mne.preprocessing.ICA(n_components=10)
ica.fit(epochs, picks=picks, decim=10)

fig = ica.plot_components()
print('len axes: ', len(fig.axes))


# %% test spectrum --------

spec = raw.compute_psd()
spec.plot_topomap()

rad_picks = mne.pick_channels_regexp(raw.info.ch_names, rad_string)
spec2 = raw.copy().pick(rad_picks).compute_psd()
spec2.plot_topomap()

[welcome]

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

[larsoner]

I've written an opm test, but would like guidance about how to access mne.datasets.ucl_opm_auditory. Test won't work as is (eg can't update UCL dataset, and wanted to confirm proper method), but the test is structured on the fnirs overlap test so probably won't need much more

The best thing to do is probably to take a file from that dataset, crop it to be very short (or maybe better, just save a clean evoked response as it's already probably ~1MB!) and add it to mne-testing-data. We don't want to make tests depend on another dataset, but we can/should update that one.

for the UCL dataset, most of the *-RAD channels are the most radial, but in some cases both RAD and TAN have low radial scores, and TAN is the most radial sensors

Yikes, so the ones that are nominally radial aren't actually the most radial? Do the sensors get turned a lot? Seems like it would need to be off by 45 degrees at least, which is a lot.

I didn't expect this to be the case in practice. We could in theory triage based on system type and channel name, but that doesn't seem all that satisfying. Does the resulting topomap look reasonable? Does it look better if you use your projection-onto-normal method (i.e., a weighted average of the colocated sensors in the direction of the ideal outward normal)?

[harrisonritz]

The best thing to do is probably to take a file from that dataset, crop it to be very short (or maybe better, just save a clean evoked response as it's already probably ~1MB!) and add it to mne-testing-data. We don't want to make tests depend on another dataset, but we can/should update that one.

Cool, I created a PR to add an evoked OPM dataset to mne_testing_data: mne-tools/mne-testing-data#120

Yikes, so the ones that are nominally radial aren't actually the most radial? Do the sensors get turned a lot? Seems like it would need to be off by 45 degrees at least, which is a lot.

I didn't expect this to be the case in practice. We could in theory triage based on system type and channel name, but that doesn't seem all that satisfying. Does the resulting topomap look reasonable? Does it look better if you use your projection-onto-normal method (i.e., a weighted average of the colocated sensors in the direction of the ideal outward normal)?

Yes, the radial-identified sensors look pretty similar to just picking the sensors with the RAD label (see below). It's only 2 / 43 where RAD-labelled sensors weren't identified as the most radial. They are right next to each other on the side of the head, could have easily been a funny placement or helmet (@georgeoneill, @neurofractal). NBD.

In our data, the max-radial picks the Z sensor 100% of the time, so I think it's a reasonable heuristic. What will be annoying is if different loc triplets should be used in different sensors types (eg fieldline). Guess it'll rely on people raising issues when it doesn't work. Down the line, maybe worth having somewhere in the info structure to store that kind of information (even lists of radial sensors).

It's actually pretty annoying to do the weighted average. Sensor radial scores are computed in a part of topomaps that just has access to info, and the combination of sensors happens in a part of layouts that just had access to data. So might require a bit of coordination to combine them.

[georgeoneill]

Good morning @harrisonritz, @larsoner! Just looking at the UCL dataset, it appears that we've got some false friends in the naming of the channels. The -RAD channels have a tangential orientation, whilst the -TAN have a radial profile.

At this point in the development of the UCL system, we laid the sensors on their side for a lower-profile helmet, but the labelling of the sensors never changed to reflect that. Apologies for the confustion!

[harrisonritz]

No worries, anything to get them to work!

When I used the x-axis orientation (loc[3:6]), the RAD sensors were (mostly) radial. I think you mentioned in the issue that we should use the y-axis orientation (loc[6:9])? I'll see how that looks.

This does highlight the trickiness of inferring the right orientation axis from info. Maybe we should add a 'radial sensors' field?

[harrisonritz]

oh great, picking the most-radial z-axis orientation matches the *-TAN sensors (and alignment is above .85 for all sensors), so can use that orientation by default.

[georgeoneill]

loc[6:9] are the orientations for each channel. 'loc[0:5] are redundant for our needs with the FIL format and are used to make visualisations of the sensors make sense.

[harrisonritz]

@georgeoneill I found that picking radial sensors using loc[9:12] selected the *-TAN sensors, which is the same as our system (for us, picking radial sensors based on loc[9:12] selects the z sensors). Just double checking, because you wrote loc[6:9] above.

[georgeoneill]

Hi @harrisonritz, my mistake, you are correct! I forgot what the order of the 12 elements of the loc array corresponded to.

0-2: Position
3-5: Redundant (except for visualisation)
6-8: Redundant (except for visualisation)
9-11: Orientation

[harrisonritz]

all good, thanks! just double-checking

[harrisonritz]

OK -- looks like it's passing major checks, sorry about the wait.

[larsoner]

Just some quick comments based on a preliminary look. Not 100% sure I followed the rename-then-drop mechanics, but had an idea that might make sense to simplify it and potentially make it more robust given that OPM names could contain . and x already. But if I'm off base on this then feel free to ignore (maybe after adding a test that uses rename_channels to add what I thought could be problematic names and shows they are not problematic 😄 ). I won't be able to do a deeper review for another week or so but wanted to give some feedback/food for thought in the meantime...

[larsoner]

We've been trying to add these places in our code nowadays

Suggested change

	def _merge_ch_data(data, ch_type, names, method="rms", modality="opm"):
	def _merge_ch_data(data, ch_type, names, method="rms", *, modality="opm"):

[larsoner]

Tiny bit more efficient because it will only allocate one new array instead of a bunch of them

Suggested change

	for rem in sorted(to_remove, reverse=True):
	del merged_names[rem]
	data = np.delete(data, rem, 0)
	for rem in sorted(to_remove, reverse=True):
	del merged_names[rem]
	data = np.delete(data, to_remove, axis=0)

[larsoner]

Is this description correct? I don't think it involves x or taking a mean?

[larsoner]

Don't use the int but rather the corresponding FIFF.FIFFV_ instead. And there should be more entries, no?

[larsoner]

If we're going to join based on a . and then check it later, better raise an error here if there is already a . in the name. But maybe I'm missing a check somewhere else. And maybe a moot point if you take my idea above

[larsoner]

Looks like it would be S1.S2, no? And maybe the overlapping ones should just end with OPM-MERGE-REMOVE or something, and the ones that we keep should just have the original channel names. Then to figure out which ones to remove you do a .endswith("OPM-MERGE-REMOVE") (which seems way safer than checking for a x or . in the name, which could happen in real data) and keep the original names.

[harrisonritz]

These look super helpful, thank you!

I'll take a pass over the next week or so (going to a week-long conference on weds).

I discovered draft mode, so hopefully now you won't get pinged for every pull.