MRG, ENH/FIX: Make cHPI calculation code public

[larsoner]

I spent a lot of time revamping our cHPI localization code to match what MaxFilter does. In particular I needed to:

1. Deal with signs properly in correlation check.
2. Check corrs per coil (not jointly) and skip if at least three match.
3. Remove projections before computing whitener (should not include projs).
4. Choose coils based on each one's individual GOF (and later by translated distance).
5. Iteratively subselect of coils from the N good ones.
6. Use a low-order SSS expansion to improve fits by suppressing noise.
7. Use a 1cm grid of "guesses" to come closer to doing a global optimization (avoid getting stuck in local minima)

After this, the 398 recordings I tested (each between 5 and 20 minutes) all match MaxFilter very well. In places that they don't, our code produces a less "noisy" result, which is an improvement as I have seen MaxFilter's estimates oscillate between two position estimates at times.

@bloyl there are a ton of changes and fixes here, it would be good if you could re-check your use cases to see if they still work properly. I did have to change some of the "spot check" values in the tests. Most were subtle, but one changed from [0.054, 0.009, 0.037] with GOF 0.075 to [-0.065, 0.075, 0.000] with GOF 0.9323 -- so way better GOF but quite far from the old value. I'm inclined to assume that the new value is correct and the old one was wrong.

Also contains changes to ProgressBar that motivated #7155. In this PR things are nicer but not as nice (or maintainable) as in #7155. So my preference is to roll back the changes here and merge #7155, but if we're going to kill #7155 these changes are useful.

[codecov]

Codecov Report

Merging #7290 into master will increase coverage by 0.28%.
The diff coverage is 94.29%.

@@            Coverage Diff             @@
##           master    #7290      +/-   ##
==========================================
+ Coverage   89.71%   89.99%   +0.28%     
==========================================
  Files         447      449       +2     
  Lines       80571    81502     +931     
  Branches    12873    13451     +578     
==========================================
+ Hits        72281    73348    +1067     
+ Misses       5441     5312     -129     
+ Partials     2849     2842       -7

[larsoner]

Hit a UTF8 encoding problem on Windows with the modernized ProgressBar that #7155 would allow us to avoid having to fix ourselves...

[larsoner]

New tutorial showing the API: https://17855-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/preprocessing/plot_59_head_positions.html

[larsoner]

That offset appears to be 1 mm. I find that acceptable.

[wmvanvliet]

fair enough

[bloyl]

Api changes look good. I'll hopefully have time today to look closer and to test on some Artemis data.

…

On Fri, Feb 7, 2020, 7:13 AM Marijn van Vliet ***@***.***> wrote: fair enough — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#7290?email_source=notifications&email_token=ABKTXHPBFAFTIFNTNWDHLJTRBVF6PA5CNFSM4KQOFXU2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELCXE2Q#issuecomment-583365226>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABKTXHLYDYO6IBCMJDGFX3DRBVF6PANCNFSM4KQOFXUQ> .

[agramfort]

+1 for MRG assuming the CI errors are unrelated

[larsoner]

Great, thanks for looking and testing @wmvanvliet and @bloyl

[bloyl]

On artemis data I don't have ground truth to compare with so quantitative testing is harder.

Qualitatively, i've looked and computing dev_head transforms a few of our datasets and I'd say that when I have good data the fits are qualitatively similar.

One issue I did have is when I have bad data the new code returns a significantly higher GOF which makes it hard to catch the bad data and the suspect dev_head transforms.

For instance in the following figure shows the topo for the 4 coils. the 1st and 4th have good topos and localize, unfortunately the 2nd and 3rd don't give good localizations but still have GOFs > 0.99. In the old code their GOFs were more like 0.93 and were easily excluded.

I don't know if this is a problem or an improvement of the GOF code but it is a difference.

@larsoner if your interested I can probably share a dataset and code with you.

Also I think i can acquire a CTF dataset that would enable a comparison between CTF fitting code and the MNE code. I hope to get on the scanner early next week to acquire it.

[larsoner]

My guess is that those coils are actually being localized (orientally potentially much) better than they were before (not getting stuck in local minima). The locations might actually be usable now.

If they're not, you might consider quantifying the coil SNR (can use PSD at the court HPI frequencies we as signal and midlines midpoints between them as noise for example), or a pairwise distance metric instead. The payee pairwise distances would give you a quick "how good are these" relative to ground truth, assuming the could coils have not moved on the head the osteosarcoma inter-coil distances should stay fixed.

EDIT: apologies for typos, on mobile

[larsoner]

But yes feel free to share code and I can do some coil quality quantification.

Along the lines of annotate_movement, I've been using pairwise coil distance code to estimate the number of usable coils and could make a PR to add it at some point. We also quantify coil SNR as a function of time and I could add that, too.

[larsoner]

I don't know if this is a problem or an improvement of the GOF code but it is a difference.

Back to a keyboard where I can actually type. In other words:

1. I suspect that it's likely an improvement in the fitting optimization that it found a location that explains 99%+ of the variance (instead of 93% before; probably from point 7 above). I didn't change the GOF calculations themselves in this PR. You can disable the "guesses" in the code pretty easily by modifying the conditional to be something like if guesses is not None and False:, then it won't do the global thing.
2. You could try running with ext_order=0 to see if it changes things (point 6 above) but it probably won't make too much difference. It's possible that you had some homogeneous fields that were getting you only 93% explained before, and those are removed now, and using ext_order=0 should undo this behavior.
3. Is it possible that the coils actually were rotated down below the helmet to the back, and previously they were being localized incorrectly but now they are correctly localized (or at least localized better than before)?

The pairwise distances might give us some insight, so if you can share that potentially coil-2-/coil-3-problematic dataset that would help. We should see that in master and this PR that the pairwise distances for coils 1 and 4 are good in both cases, and I suspect that we'll see that they are bad for coils 2 and 3 (at least to coils 1 and 4) in master but better in this PR.

I'm happy to check these things if you share the data.

[larsoner]

From talking to @bloyl it seems likely that having access to the actual amplitude estimates could also be useful for QA purposes. So the steps are now split just a little bit more:

• chpi_amplitudes = compute_chpi_amplitudes(raw) (for any vendor)
• chpi_locs = compute_chpi_locs(raw.info, chpi_amplitudes) (for any vendor)
• chpi_locs = extract_chpi_locs_ctf(raw) (only for ctf)
• head_pos = compute_head_pos(chpi_locs) (for any vendor, from either of the two chpi_locs returning functions above)

Making it so that chpi_amplitudes is returned will eventually make it easier to tell for example if one coil has much lower amplitudes than the others, suggesting that it might be much farther from the helmet.

@bloyl in changing this API I found ways of unifying things again, can you look in particular at artemis123.py to see if you're okay with the changes there? Basically we fake a hpi_result to make _get_initial_hpi_fit happy in a way that should be backward compatible with what you were doing before. This has the advantage of getting rid of the need for a separate _fit_device_hpi_positions function (now we can just use compute_chpi_amplitudes and compute_chpi_locs directly).

I also added some basic pairwise-distance code back in, as I think eventually @bloyl might want it for Artemis123, we'll probably want it for QA code, and it's another nice sanity check.

[wmvanvliet]

From an regular user point of view: I'm fine with the extra API steps as long as the maxwell_filter function can do continuous head position correction "automatically", meaning I can just set a flag to True and it will do it, just like the maxfilter command line program.

[larsoner]

From an regular user point of view: I'm fine with the extra API steps as long as the maxwell_filter function can do continuous head position correction "automatically", meaning I can just set a flag to True and it will do it, just like the maxfilter command line program.

I'd rather not do this for the following reasons:

1. We are already headed in a separation-of-functionality direction with other aspects of MF. For just two examples (there might be more):
   • filter_chpi -- you have to call this separately before doing maxwell_filter to get what MaxFilter will do automatically and by default for you.
   • MaxFilter always writes files but we do not have out=True by default (or even an option), so users are already doing raw_sss.save manually afterward.
2. head_pos in master is currently a separate step because we don't support it, so I think the barrier for adoption / annoyance level with still needing to make separate call(s) should stay low (at least now they're in Python!).
3. Adding head_pos=True by itself will leave implicit a bunch of choices that you cannot change, or if we add the kwargs to change them to maxwell_filter, it will further blow up the already long API for maxwell_filter.
4. If people want a single monolithic equivalent, we can later add a new mne maxwell_filter, possibly mimicking the MaxFilter API.

TL;DR: I'd argue in the opposite direction, that maxwell_filter is already too monolithic and ideally would have been broken down into separate steps from the start (but actually a number of the steps are not easily separable, such as the destination head position until we add a FIFFV_MULTIPOLE_MOMENT coil type, tSSS from time-varying SSS decomposition, eventually #7156, etc.).

[bloyl]

Just made a first pass through, here are my initial comments. I'll hopefully have more time tomorrow to dig into it.

[bloyl]

How important is having this initial guess as a requirement?

In the case of the artemis data it needs to be faked before we can do a head localization. You do this now in the reader. but in my actual code I don't use the reader's head localization so that I have a more flexability over which time windows i choose to localize.

This is minor since I doubt we'll do head localization in too many other systems. but figured I'd mention it

[larsoner]

With the guesses code in there now using a grid, it probably does not matter much unless you are quite far out of the helmet (and even then might not matter)

[agramfort]

@bloyl i will let you merge if happy

[larsoner]

@bloyl coil moments added. Feel free to look and merge if you're happy

As a side note, I don't know why our calculation and MaxFilter's differ by _MAG_FACTOR (1e-7). I'm guessing that they accidentally double-multiply by it? cc @agramfort in case you know.

[bloyl]

Thanks @larsoner