unstable estimation of the head position based on the HPI coils

Description of the problem

When we estimated the time-varying head position with the MEG data that we had collected (306-channel Vectorview, Elekta Neuromag), we found that sometimes the estimated head position was unstable consistently “teleporting” between two locations over time, which couldn’t be thought of as actual human movements. 5 HPI coils were installed during the MEG acquisition (261, 278, 294, 311, and 328 Hz) and the head position was computed by applying three chpi functions sequentially to these 5 HPI coils (i.e., mne.chpi.compute_chpi_amplitudes -> mne.chpi.compute_chpi_locs -> mne.chpi.compute_head_pos). I tried to figure out with Dr. Burkhard Maess (@bmaess, at MPI CBS) in which part this weird pattern was produced. Adjusting the size of a window and a step didn’t help to solve the problem, and we observed that this problem could happen even when the estimated cHPI locations (results of mne.chpi.compute_chpi_locs) were actually smooth and showed good gofs (goodness of fit).

Steps to reproduce

import os
import mne
import matplotlib.pyplot as plt

raw_file = fpath + os.sep + 'te01a1.fif' # 'fpath' is the path where the sample file is stored

## 1) compute head positions
raw = mne.io.read_raw_fif(raw_file)
# compute cHPI amplitudes
chpi_amps = mne.chpi.compute_chpi_amplitudes(raw, t_step_min=0.25, t_window=0.5) 
# compute cHPI locations
chpi_locs = mne.chpi.compute_chpi_locs(raw.info, chpi_amps)
# compute head position
head_pos = mne.chpi.compute_head_pos(raw.info, chpi_locs) 
mne.viz.plot_head_positions(head_pos)

## 2) plot the estimated cHPI locations (i.e., chpi_locs)
t = chpi_locs['times']
coil2plot = 1 # from 1 to 5
locs = chpi_locs['rrs'][:,coil2plot-1,:]
moments = chpi_locs['moments'][:,coil2plot-1,:]
 
fig, ax = plt.subplots(3,2)
fig.suptitle(f'Locs and Moments of the HPI coil %d' % coil2plot)
 
ax[0,0].set_title('Locations')
ax[0,0].plot(t, locs[:,0])
ax[0,0].set_ylabel('x')
 
ax[1,0].plot(t, locs[:,1])
ax[1,0].set_ylabel('y')
 
ax[2,0].plot(t, locs[:,2])
ax[2,0].set_ylabel('z')
ax[2,0].set_xlabel('time [s]')
 
ax[0,1].set_title('Moments')
ax[0,1].plot(t, moments[:,0])
ax[0,1].set_ylabel('x')
 
ax[1,1].plot(t, moments[:,1])
ax[1,1].set_ylabel('y')
 
ax[2,1].plot(t, moments[:,2])
ax[2,1].set_ylabel('z')
ax[2,1].set_xlabel('time [s]')
 
# The figure title should not overlap with the subplots.
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

Link to data

Download (363.03 MiB)

Expected results

This is the figure that we got after we had debugged mne.chpi.compute_head_pos to fix the problem.

We found that inside of mne.chpi.compute_head_pos, there is a function mne.chpi._fit_chpi_quat_subset that chooses maximum 3 HPI coils based on the gof (goodness of fit) to estimate the head position at a given time point.

https://github.com/mne-tools/mne-python/blob/0fbd8fcae6e6c3a69c47069da9af8127bd113e8a/mne/chpi.py#L939-L949

This sometimes caused the different members of the HPI coils to be selected for the different time points. It could make the basis for estimating the head position change over time, and thus the estimated head position could be unstable. In fact, we observed that the sudden jumps or teleportation in the estimated head position coincide with the changes in the members of the HPI coils in use.

## a code to replicate the above figure
import os
import numpy as np
import mne
import matplotlib.pyplot as plt

raw_file = fpath + os.sep + 'te01a1.fif' # 'fpath' is the path where the sample file is stored

# compute head positions
raw = mne.io.read_raw_fif(raw_file)
chpi_amps = mne.chpi.compute_chpi_amplitudes(raw, t_step_min=0.25, t_window=0.5) # compute cHPI amplitudes
chpi_locs = mne.chpi.compute_chpi_locs(raw.info, chpi_amps) # compute cHPI locations
head_pos = mne.chpi.compute_head_pos(raw.info, chpi_locs) # compute head position
 
# find out the coil indices used for each time point
times = chpi_locs['times']
coil_idx_mtx = np.empty([3, len(times)])
time_idx = 0
 
hpi_dig_head_rrs = mne.chpi._get_hpi_initial_fit(raw.info, adjust=None, verbose='error')
for fit_time, this_coil_dev_rrs, g_coils in zip(
            *(chpi_locs[key] for key in ('times', 'rrs', 'gofs'))):
    use_idx = np.where(g_coils >= 0.98)[0]
 
    _, _, use_idx = mne.chpi._fit_chpi_quat_subset(this_coil_dev_rrs, hpi_dig_head_rrs, use_idx)
    coil_idx_mtx[:,count] = use_idx[:]
    time_idx += 1
 
# plot the sum of coil indices with head position
t = head_pos[:,0]
coils = np.sum(coil_idx_mtx, axis=0)
 
fig, ax = plt.subplots(3,1)
fig.suptitle(f'head position with the selected coils’')
 
ax[0].plot(t, head_pos[:,4]*-1000, 'b')
ax[0].set_ylabel('x [mm]', color='b')
ax[0].set_xlim((t[0], t[-1]))
ax1 = ax[0].twinx()
ax1.plot(t, coils, 'g', alpha=0.6)
ax1.set_ylim((0,10))
 
ax[1].plot(t, head_pos[:,5]*1000, 'b')
ax[1].set_ylabel('y [mm]', color='b')
ax[1].set_xlim((t[0], t[-1]))
ax2 = ax[1].twinx()
ax2.plot(t, coils, 'g', alpha=0.6)
ax2.set_ylabel('Sum of Channel Indices in Use (among 0 to 4)', color='g')
ax2.set_ylim((0,10))
 
ax[2].plot(t, head_pos[:,6]*-1000, 'b')
ax[2].set_ylabel('z [mm]', color='b')
ax[2].set_xlim((t[0], t[-1]))
ax[2].set_xlabel('time [s]')
ax3 = ax[2].twinx()
ax3.plot(t, coils, 'g', alpha=0.6)
ax3.set_ylim((0,10))
 
# The figure title should not overlap with the subplots.
fig.tight_layout(rect=[0, 0.03, 1, 0.95])

This problem has been solved after we inactivated mne.chpi._fit_chpi_quat_subset called by mne.chpi.compute_head_pos such that the head position at a given time point is estimated by all available HPI coils (i.e., which, at least, satisfy both gof_limit and dist_limit constraints).

Actual results

This is the figure showing the unstable head position estimation when mne.chpi._fit_chpi_quat_subset is activated inside of mne.chpi.compute_head_pos

Additional information

Platform: Linux-5.10.0-19-amd64-x86_64-with-glibc2.31 Python: 3.10.5 | packaged by conda-forge | (main, Jun 14 2022, 07:04:59) [GCC 10.3.0] Executable: /data/pt_01585/miniconda3/envs/mne_v110_baek/bin/python CPU: : 32 cores Memory: 503.5 GB

mne: 1.1.0 numpy: 1.22.4 {} scipy: 1.8.1 matplotlib: 3.5.2 {backend=QtAgg}

sklearn: 1.1.2 numba: 0.56.0 nibabel: 4.0.1 nilearn: 0.9.2 dipy: 1.5.0 cupy: Not found pandas: 1.4.3 pyvista: 0.36.1 {OpenGL 4.5 (Core Profile) Mesa 20.3.5 via llvmpipe (LLVM 11.0.1, 256 bits)} pyvistaqt: 0.9.0 ipyvtklink: 0.2.2 vtk: 9.2.0.rc2 qtpy: 2.2.0 {PyQt5=5.15.2} ipympl: Not found pyqtgraph: 0.12.4 pooch: v1.6.0

mne_bids: Not found mne_nirs: Not found mne_features: Not found mne_qt_browser: 0.3.1 mne_connectivity: Not found mne_icalabel: Not found

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

Thanks for the excellent, in-depth analysis @SeungCheolBaek !

I don't have time to look too deeply today, but will try tomorrow. My response from a quick glance at what you wrote is that fitting should be done using as many coils as possible that satisfy the GOF and distance criteria. Somewhere in the code there is something that should look to see if any coils violate the distance criteria, and if they do, subselects (somehow) 4 or 3 coils, stopping once the distiance criteria is met. It very well might not be doing this.

And I, too, have observed this oscillatory behavior but never tracked it down -- if we could get rid of it, that would be great!

If you don't want to wait for me to look and confirm your suspicions, one way for you to test what (I think) is your idea would be to

add a logger.debug line saying how many coils are being used to fit
in the tests do a fit with verbose='debug' and catch_logging
check that the number varies over time in the resulting log, and/or is pinned at 5, rather than 3

This test should pass if everything is working correctly in our code. It sounds from your analysis like it is not, so this test would fail. You could then fix the bug in the code, and the test would pass. (Thus you'd be doing test-driven development, which we try to follow the spirit of at least!)

Thank you for all of the comments!

I have actually done something similar to what you suggested.

Concerning the number of HPI coils being used for fitting, in this version (v.1.1.0), it's always three as long as mne.chpi._fit_chpi_quat_subset is called inside of mne.chpi.compute_head_pos to fit the head position. So, what I tried was to comment out mne.chpi._fit_chpi_quat_subset and replace it with mne.chpi._fit_chpi_quat so that all available channels could be used for the fitting (except for the coils not meeting the gof_limit criterion). Then, the oscillation in the estimated head position disappeared.
I also estimated the head position by varying the number of the coils at different time points (i.e., from 3 to 5). As long as more than 3 coils were used for fitting, the estimated head position was pretty robust to the varying number of coils. Even in the case where only three coils were used and its members didn't change, the estimates were stable. The oscillating head position only occurred when mne.chpi._fit_chpi_quat_subset forced to use only three coils for the fitting at a given time point and changed the members over time (i.e., based on gof; please refer to the second figure in my first message - green line shows the sum of indices of three coils being used for the estimation).

Hello, I would suggest to generally refrain from switching the number of involved coils. We may use the the inverse of the error covariance as a weight within the quaternion fits. Then a bad coil would have a reduced influence. And no switch is needed. A more correct way would be to include error estimates for the estimated (x,y,z) coordinates of each HPI coil. However, these values need to be a (second) result of the coil fits ...

In this paper the authors apply such a weighting (its for a spacecraft and not for a human head, but the difference is negligible here ;) https://ntrs.nasa.gov/api/citations/20000034107/downloads/20000034107.pdf

I would suggest to generally refrain from switching the number of involved coils. We may use the the inverse of the error covariance as a weight within the quaternion fits.

My first reaction to this is that, to a certain extent, we need to stop using some coils sometimes because they either stop emitting signals (rare), are masked by some artifact in the HPI frequency region (rare), or shift location on the head of the participant during the recording (not so rare, and still good GOF) so shouldn't be trusted. But I think this can be made consistent with what you're saying if we think of "stop using" as a smooth downweighting rather than binary elimination during fitting. In other words, our enemy/problem here is that the weights for the coils are always 0 or 1, and they switch discontinuously -- and this correspondingly creates discontinuities in our estimates. If we instead smoothly vary the weights somehow, this problem should go away.

A more correct way would be to include error estimates for the estimated (x,y,z) coordinates of each HPI coil. However, these values need to be a (second) result of the coil fits ...

One way out of this is the pairwise coil distance differences, i.e., you compute the pairwise distances of all 5 coils to the other 4 coils for the original digitization (treated as ground truth), then at each quat-fitting time you recompute these distances and take the difference in these matrices. This gives you a "pairwise distance error", which can easily then by computed per coil (using a row mean of the 5x4 array). This does not require fitting quats, and I think it's actually calculated already. (If one coil gets really bad for example, all pairwise distances will suffer -- but the average pairwise distance for the really bad coil will suffer the most.)

Maybe the thing to do here is:

Still use the gof and dist criteria as strict/binary criteria to say: do we fit a quat (>= 3 coils meeting both criteria) or not
If we do fit one, use the gof and inter-coil-distance error to compute weights for the coil distances during quat fitting.

For (2), we "just" need to decide the right way to compute these weights. There are a number of options, but these two at least seem reasonable potentially:

Say "0 when gof < gof_limit; and rising linear ramp (or quadratic) from 0-to-1 as gof_limit <= gof <= 1", then take this and multiply it by "the average inter-coil-distance error for the given coil"
Same as above, but for the gof part use something like the 10th power of gof, which for gof=0.5 yields a weight of 0.001, for 0.9 a weight of 0.3, for 0.95 a weight of 0.6, and for 0.99 a weight of 0.9.

I'm not sure (2) is really necessary, even though it stays totally smooth (and differentiable, etc.) -- having "linear smoothness" in the usable region (transitioning to zero below) it might be enough to get well behaved fits already. And it keeps us closer to what gof_limit used to do anyway -- completely turn off some coils below the limit -- while hopefully introducing sufficient smoothness in the usable region to avoid a discontinuity when the coil weight goes to zero.

In this paper the authors apply such a weighting

I did not read this, but let me know if you find something interesting :+1:

mne-tools / mne-python