BUG: Add estimated fiducials when missing / assumed head coords

[sappelhoff]

I'd argue that the electrode positions as read from the Neuroscan CNT file header are in unknown space, see also the warning here: https://mne.tools/dev/auto_tutorials/io/20_reading_eeg_data.html#import-cnt

In the reader code, I see no indication of "sphere projection" and even if there was one, then we wouldn't know that the sphere corresponds to the "head" coord system.

follow up to:

• BUG: Fix bugs with digitization #11146

fixes issue identified in mne-tools/mne-bids#1066 (comment)

Needed for mne-bids, would be nice to have in 1.2 release

CI issues seem unrelated

cc @larsoner

[agramfort]

any chance to add a test ? 🙏

[larsoner]

@sappelhoff I see you have only modified the CNT reader here. In #11146 I messed with tests for at least curry, fieldtrip, and CNT. More importantly, though, I changed _test_raw_reader, which runs for every io type we have (or at least it should!). IIRC this is actually what caused me to change some of those other readers -- I asserted something that I thought should be the case in _test_raw_reader and it turned out not to be the case -- true TDD for once 😄 So I think to check the behavior here, you should add something to _test_raw_reader to ensure we're consistent about this across readers, then your CNT changes should follow from -- or at least be consistent with -- that change.

What should the consistency be? Ideally I guess that the dig coords for EEG should be in head coords if and only if there are fiducials also present in dig. Does that sound right? This would fail on main but pass on this PR for CNT, right?

Then the question is, how many other readers would this break, and should we also have them set to UNKNOWN in some cases...

Then another question -- probably for another PR -- is: are we then consistent between the coord frame we say in dig vs the one we say in info['chs']['coord_frame'].

Sorry for making this more complicated than it initially seemed, but hopefully this helps continue moving us toward being consistent with all this stuff across all reader!

[larsoner]

@sappelhoff I think this is the last blocker for 1.2 release. Do you think you can get to it in the next few days? If not, I can push some commits

[sappelhoff]

I was already afraid that it wouldn't be this simple 😉

My only expertise with this PR was that:

1. this fixes mne-bids CI failures
2. it's clear to me that elec positions read from a CNT header cannot possibly be in head coord system

Having that said, I agree that we should work on a more complete fix.

Ideally I guess that the dig coords for EEG should be in head coords if and only if there are fiducials also present in dig. Does that sound right?

Yes that sounds right to me! Except maybe for data formats that are guaranteed to ship electrode positions in "head" format (but nothing apart from BrainVision comes to mind). The coord system should be "unknown" otherwise.

Then another question -- probably for another PR -- is: are we then consistent between the coord frame we say in dig vs the one we say in info['chs']['coord_frame'].

I am not so clear about the distinction, and why there even needs to be one -- where should I read up on this?

@sappelhoff I think this is the last blocker for 1.2 release. Do you think you can get to it in the next few days? If not, I can push some commits

Unfortunately I can't -- I'll be on vacation starting next Sunday and need to wrap up some work things as well. But I can probably read comments and reply ... just not really dig into and write code (not before October 10th)

[larsoner]

I am not so clear about the distinction, and why there even needs to be one -- where should I read up on this?

Not sure we have proper docs on it. But the TL;DR is that we store electrode information in two places: info['dig'] and in info['chs'][ii]['loc'][:3] (and ref in [3:6] when possible). For FIF these are always consistent, I think if you apply_montage they will be consistent, too. We've tried to keep them consistent in readers but I'm not sure if we have.

Moreover, historically the FIF format only allowed EEG points to have the FIFFV_COORD_HEAD coord frame for info['chs'][ii]['loc'] -- it wasn't even a value stored in the file, it was just assumed that if you had EEG electrodes they were in this frame. We recently (1-2 years ago?) added "extended channel block" support where you can say, for each channel, what its coord frame actually is. But we haven't used it much.

[larsoner]

@agramfort @sappelhoff this keeps getting slightly more and more hacky, I had to work around this:

______________________ test_csd_degenerate[testing_data] _______________________
mne/preprocessing/tests/test_csd.py:156: in test_csd_degenerate
    raw = compute_current_source_density(raw)
mne/preprocessing/_csd.py:132: in compute_current_source_density
    radius, origin_head, origin_device = fit_sphere_to_headshape(inst.info)
mne/bem.py:1083: in _fit_sphere_to_headshape
    hsp = get_fitting_dig(info, dig_kinds)
mne/bem.py:1060: in get_fitting_dig
    raise RuntimeError('Digitization points not in head coordinates, '
E   RuntimeError: Digitization points not in head coordinates, contact mne-python developers

This is because in this PR, the info['dig'] are left in UNKNOWN coords (while the info['chs'] are in HEAD) :(

This all only really matters in the case where we are setting a montage (or equivalently, reading ch pos in a reader) in the UNKNOWN frame that does not have fiducials. I think I managed to convince myself that a better solution in this case is to insert synthetic fiducials when "converting" to the head coordinate frame. Compared to the current approach:

• This requires creating/estimating some fiducials. At first this might seem like a disadvantage. But:
  1. This is fairly easy: estimate the average radius r of all EEG dig points, and put points at (-r, 0, 0), (0, r, 0), and (0, 0, r). (This ensures that the unknown->head is identity, so they won't move when the montage is converted to the head frame)
  2. This makes explicit what we're implicitly doing conceptually anyway in 1.1 by assigning info['dig'] and info['chs'] the HEAD frame. The extent to which this solution is wrong, it always has been wrong, it will just more obvious/clear/explicit rather than implicit now. (LPA/RPA/Nas generally don't perfectly lie on a sphere on real heads, but I think this is negligible for our purposes.)
  3. Increases usability of the data with coreg, MNE-BIDS, etc.
• This ensures that info['dig'][ii]['coord_frame'] and info['chs'][ii]['coord_frame'] stay consistent, i.e., both in HEAD. In the current PR, chs will have HEAD and dig will have UNKNOWN.
• This maintains backward compat with 1.1 better, where info['dig'] were always in HEAD coords for EEG readers, even when we aren't sure they use the head coordinate frame. (This PR's current method is to leave info['dig'] in UNKNOWN.)

In the long run, a solution where we do this synthetic insertion, then allow/improve the montage code to allow people to adjust the coordinate frame and/or these fiducial points seems better than introducing more inconsistencies.

[drammock]

(LPA/RPA/Nas generally don't perfectly lie on a sphere on real heads, but I think this is negligible for our purposes.)

to this point, we could fit them to a spheroid (revolved ellipse) instead of a sphere. But allowing people to adjust the fiducials later may be enough / better.

[larsoner]

Yeah we could -- but I think the simpler the better here

[sappelhoff]

I agree with @larsoner's proposal, the backward compatibility and the part below convinced me specifically:

The extent to which this solution is wrong, it always has been wrong, it will just more obvious/clear/explicit rather than implicit now.

If this fixes our current issues and makes problems more obvious to fix in the future, while not changing behavior, I am fine with it! Thanks!

[larsoner]

Okay, proposal implemented. Ready for review/merge from my end

[larsoner]

I will open a new issue for this once this PR is merged

[agramfort]

is it worth running circle full to make sure the examples look good? I am thinking about our tutorials using EEG template MRI.

🙏 @sappelhoff @larsoner

[agramfort]

a file in testing should be fixed?

[larsoner]

No idea if the test file is broken or our reading of it. The associated test appears to be written by @agramfort so maybe they'll be able to figure it out most quickly :)

[sappelhoff]

Thanks a lot @larsoner