ENH: Add function to determine time-resolved chpi status

[eort]

References issue #940 in the mne-bids repository that is about mne-bids incorrectly setting the cHPI field in the meg recording sidecar.

Based on that discussion, we (i.e. @larsoner) thought that it might be useful to have a function in mne that returns the time-resolved number of active HPI coils. With that information it should be possible to infer whether continous HPI was active or not.

This PR is work in progress. I mostly wanted to shift the discussion over here. What still misses:

• implementation for CTF files
• implementation for KIT files(?)
• tests

I don't know anything about CTF files nor any files available (other than those in mne-data), but if I understand the test below correctly, for CTF files it holds that if there are any HPI coils in the data file, they are automatically on. If there are not present, they are obviously off. Can someone who has experience with such files (dis)confirm that?

mne-python/mne/tests/test_chpi.py

Line 664 in 9634536

def test_calculate_head_pos_ctf():

Same holds for KIT files. If it is not clear, we could also implement only for neuromag systems and return some defaults for other systems, together with a info/warning or something like that.

Furthermore, I am not sure whether inference by suffix is robust enough to determine the manufacturer (neuromag, ctf, kit, etc). Is there a better way?

[larsoner]

If it is not clear, we could also implement only for neuromag systems and return some defaults for other systems, together with a info/warning or something like that.

For now I would make it just work for Neuromag systems and raise an error for others

I am not sure whether inference by suffix is robust enough to determine the manufacturer (neuromag, ctf, kit, etc). Is there a better way?

The logic in get_meg_helmet_surf should be good enough to tell neuromag vs non-neuromag, I'd reuse that (in whatever way it does it) in a DRY way

[eort]

Both suggestions sound good.

I meant to add a test for inactive hpis in neuromag systems, too, but I can't find files in mne_data that I can use for that. The ones I tried miss the hpi_subsystem field in their info, so that this line throws a ValueError ("No appropriate cHPI information found in info["hpi_meas"] and info["hpi_subsystem"])

mne-python/mne/chpi.py

Line 1410 in 01cad56

chpi_info = get_chpi_info(raw.info)

In my own test file, those fields are present though, such that the method works. Any ideas?

[larsoner]

Did you check test_move_anon_raw.fif ? It should have it

>>> raw = mne.io.read_raw_fif(mne.datasets.testing.data_path() / 'SSS' / 'test_move_anon_raw.fif', allow_maxshield='yes')
>>> raw.info['hpi_subsystem']
{'ncoil': 5, 'hpi_coils': [{'event_bits': array([256,   0, 256, 256], dtype=int32)}, {'event_bits': array([512,   0, 512, 512], dtype=int32)}, {'event_bits': array([1024,    0, 1024, 1024], dtype=int32)}, {'event_bits': array([2048,    0, 2048, 2048], dtype=int32)}, {'event_bits': array([4096,    0, 4096, 4096], dtype=int32)}], 'event_channel': 'STI201'}

for tests that already use it, see https://github.com/mne-tools/mne-python/blob/main/mne/tests/test_chpi.py#L41

[eort]

Yeah, I use it to test whether the method works for when the cHPIs were active. But I though I would also test a file where cHPI were inactive. And a file for that I can't find. Or is that overkill?

[larsoner]

Ideally we would have such a test file but we don't. I would just fake it in the test with something like:

raw_on = read_raw_fif(...).load_data()
status_ch = raw.ch_names.index('SYST201')  # or whatever it is
raw_on._data[status_ch][:1000] = 0
raw_on._data[status_ch][1000:2000] = ... # some set of status bits that show that 3 of the 5 are active for example
...  # some code that checks that it shows off for the first 1000 samples, then on or partiall on for the next 1000

Then also verify it works properly on some of your actual files, and I think we can trust that the hack of setting the status_ch directly in tests is good enough

[eort]

yeah, sounds reasonable. As far as I am concerned, the PR is ready for review.

[larsoner]

LGTM other than a couple of minor things! Can you also please

• update latest.inc as well
• add it to tutorials/preprocessing/60_maxwell_filtering_sss.py or examples/preprocessing/movement_compensation.py in some suitable place

[eort]

I think this is ready for merging @larsoner

[larsoner]

Thanks @eort !