# Authors: The MNE-Python contributors. # License: BSD-3-Clause # Copyright the MNE-Python contributors. from pathlib import Path import numpy as np import pytest from numpy.testing import assert_allclose, assert_array_equal, assert_array_less from scipy.interpolate import interp1d from scipy.spatial.distance import cdist from mne import pick_info, pick_types from mne._fiff.constants import FIFF from mne.chpi import ( _chpi_locs_to_times_dig, _compute_good_distances, _get_hpi_initial_fit, _setup_ext_proj, compute_chpi_amplitudes, compute_chpi_locs, compute_chpi_snr, compute_head_pos, extract_chpi_locs_ctf, extract_chpi_locs_kit, filter_chpi, get_active_chpi, get_chpi_info, head_pos_to_trans_rot_t, read_head_pos, write_head_pos, ) from mne.datasets import testing from mne.forward._compute_forward import _MAG_FACTOR from mne.io import ( RawArray, read_info, read_raw_artemis123, read_raw_ctf, read_raw_fif, read_raw_kit, ) from mne.simulation import add_chpi from mne.transforms import _angle_between_quats, rot_to_quat from mne.utils import ( _record_warnings, assert_meg_snr, catch_logging, object_diff, verbose, ) from mne.viz import plot_head_positions base_dir = Path(__file__).parents[1] / "io" / "tests" / "data" ctf_fname = base_dir / "test_ctf_raw.fif" hp_fif_fname = base_dir / "test_chpi_raw_sss.fif" raw_fname = base_dir / "test_raw.fif" data_path = testing.data_path(download=False) sample_fname = data_path / "MEG" / "sample" / "sample_audvis_trunc_raw.fif" chpi_fif_fname = data_path / "SSS" / "test_move_anon_raw.fif" pos_fname = data_path / "SSS" / "test_move_anon_raw.pos" sss_fif_fname = data_path / "SSS" / "test_move_anon_raw_sss.fif" sss_hpisubt_fname = data_path / "SSS" / "test_move_anon_hpisubt_raw.fif" chpi5_fif_fname = data_path / "SSS" / "chpi5_raw.fif" chpi5_pos_fname = data_path / "SSS" / "chpi5_raw_mc.pos" ctf_chpi_fname = data_path / "CTF" / "testdata_ctf_mc.ds" ctf_chpi_pos_fname = data_path / "CTF" / "testdata_ctf_mc.pos" art_fname = ( data_path / "ARTEMIS123" / "Artemis_Data_2017-04-04-15h-44m-22s_Motion_Translation-z.bin" ) art_mc_fname = ( data_path / "ARTEMIS123" / "Artemis_Data_2017-04-04-15h-44m-22s_Motion_Translation-z_mc.pos" ) con_fname = data_path / "KIT" / "MQKIT_125_2sec.con" mrk_fname = data_path / "KIT" / "MQKIT_125.mrk" elp_fname = data_path / "KIT" / "MQKIT_125.elp" hsp_fname = data_path / "KIT" / "MQKIT_125.hsp" berlin_fname = data_path / "KIT" / "data_berlin.con" @testing.requires_testing_data def test_chpi_adjust(): """Test cHPI logging and adjustment.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes") with catch_logging() as log: _get_hpi_initial_fit(raw.info, adjust=True, verbose="debug") get_chpi_info(raw.info, on_missing="raise", verbose="debug") # Ran MaxFilter (with -list, -v, -movecomp, etc.), and got: msg = [ "HPIFIT: 5 coils digitized in order 5 1 4 3 2", "HPIFIT: 3 coils accepted: 1 2 4", "Hpi coil moments (3, 5):", "2.08542e-15 -1.52486e-15 -1.53484e-15", "2.14516e-15 2.09608e-15 7.30303e-16", "-3.2318e-16 -4.25666e-16 2.69997e-15", "5.21717e-16 1.28406e-15 1.95335e-15", "1.21199e-15 -1.25801e-19 1.18321e-15", "HPIFIT errors: 0.3, 0.3, 5.3, 0.4, 3.2 mm.", "HPI consistency of isotrak and hpifit is OK.", "HP fitting limits: err = 5.0 mm, gval = 0.980.", "Using 5 HPI coils: 83 143 203 263 323 Hz", # actually came earlier ] log = log.getvalue().splitlines() assert set(log) == set(msg), "\n" + "\n".join(set(msg) - set(log)) # Then took the raw file, did this: raw.info["dig"][5]["r"][2] += 1.0 # And checked the result in MaxFilter, which changed the logging as: msg = ( msg[:8] + [ "HPIFIT errors: 0.3, 0.3, 5.3, 999.7, 3.2 mm.", "Note: HPI coil 3 isotrak is adjusted by 5.3 mm!", "Note: HPI coil 5 isotrak is adjusted by 3.2 mm!", ] + msg[-2:] ) with catch_logging() as log: _get_hpi_initial_fit(raw.info, adjust=True, verbose="debug") get_chpi_info(raw.info, on_missing="raise", verbose="debug") log = log.getvalue().splitlines() assert set(log) == set(msg), "\n" + "\n".join(set(msg) - set(log)) @testing.requires_testing_data def test_read_write_head_pos(tmp_path): """Test reading and writing head position quaternion parameters.""" temp_name = tmp_path / "temp.pos" # This isn't a 100% valid quat matrix but it should be okay for tests head_pos_rand = np.random.RandomState(0).randn(20, 10) # This one is valid head_pos_read = read_head_pos(pos_fname) for head_pos_orig in (head_pos_rand, head_pos_read): write_head_pos(temp_name, head_pos_orig) head_pos = read_head_pos(temp_name) assert_allclose(head_pos_orig, head_pos, atol=1e-3) # Degenerate cases pytest.raises(TypeError, write_head_pos, 0, head_pos_read) # not filename pytest.raises(ValueError, write_head_pos, temp_name, "foo") # not array pytest.raises(ValueError, write_head_pos, temp_name, head_pos_read[:, :9]) pytest.raises(TypeError, read_head_pos, 0) pytest.raises(OSError, read_head_pos, "101") @testing.requires_testing_data def test_hpi_info(tmp_path): """Test getting HPI info.""" temp_name = tmp_path / "temp_raw.fif" for fname in (chpi_fif_fname, sss_fif_fname): raw = read_raw_fif(fname, allow_maxshield="yes").crop(0, 0.1) assert len(raw.info["hpi_subsystem"]) > 0 raw.save(temp_name, overwrite=True) info = read_info(temp_name) assert len(info["hpi_subsystem"]) == len(raw.info["hpi_subsystem"]) # test get_chpi_info() info = read_info(chpi_fif_fname) hpi_freqs, stim_ch_idx, hpi_on_codes = get_chpi_info(info) assert_allclose(hpi_freqs, np.array([83.0, 143.0, 203.0, 263.0, 323.0])) assert stim_ch_idx == 378 assert_allclose(hpi_on_codes, np.array([256, 512, 1024, 2048, 4096])) # test get_chpi_info() if no proper cHPI info is available with info._unlock(): info["hpi_subsystem"] = None info["hpi_meas"] = [] info["hpi_results"] = [] with pytest.raises(ValueError, match="No appropriate cHPI information"): get_chpi_info(info) with pytest.warns(RuntimeWarning, match="No appropriate cHPI information"): get_chpi_info(info, on_missing="warn") hpi_freqs, stim_ch_idx, hpi_on_codes = get_chpi_info(info, on_missing="ignore") assert_array_equal([], hpi_freqs) assert stim_ch_idx is None assert_array_equal([], hpi_on_codes) def _assert_quats( actual, desired, dist_tol=0.003, angle_tol=5.0, err_rtol=0.5, gof_rtol=0.001, vel_atol=2e-3, ): # 2 mm/s """Compare estimated cHPI positions.""" __tracebackhide__ = True trans_est, rot_est, t_est = head_pos_to_trans_rot_t(actual) trans, rot, t = head_pos_to_trans_rot_t(desired) quats_est = rot_to_quat(rot_est) gofs, errs, vels = desired[:, 7:].T gofs_est, errs_est, vels_est = actual[:, 7:].T del actual, desired # maxfilter produces some times that are implausibly large (weird) if not np.isclose(t[0], t_est[0], atol=1e-1): # within 100 ms raise AssertionError( f"Start times not within 100 ms: {t[0]:0.3f} != {t_est[0]:0.3f}" ) use_mask = (t >= t_est[0]) & (t <= t_est[-1]) t = t[use_mask] trans = trans[use_mask] quats = rot_to_quat(rot) quats = quats[use_mask] gofs, errs, vels = gofs[use_mask], errs[use_mask], vels[use_mask] # double-check our angle function for q in (quats, quats_est): angles = _angle_between_quats(q, q) assert_allclose(angles, 0.0, atol=1e-5) # limit translation difference between MF and our estimation trans_est_interp = interp1d(t_est, trans_est, axis=0)(t) distances = np.sqrt(np.sum((trans - trans_est_interp) ** 2, axis=1)) assert np.isfinite(distances).all() arg_worst = np.argmax(distances) assert distances[arg_worst] <= dist_tol, ( f"@ {t[arg_worst]:0.3f} seconds: " f"{1000 * distances[arg_worst]:0.3f} > {1000 * dist_tol:0.3f} mm" ) # limit rotation difference between MF and our estimation # (note that the interpolation will make this slightly worse) quats_est_interp = interp1d(t_est, quats_est, axis=0)(t) angles = 180 * _angle_between_quats(quats_est_interp, quats) / np.pi arg_worst = np.argmax(angles) assert angles[arg_worst] <= angle_tol, ( f"@ {t[arg_worst]:0.3f} seconds: " f"{angles[arg_worst]:0.3f} > {angle_tol:0.3f} deg" ) # error calculation difference errs_est_interp = interp1d(t_est, errs_est)(t) assert_allclose( errs_est_interp, errs, rtol=err_rtol, atol=1e-3, err_msg="err" ) # 1 mm # gof calculation difference gof_est_interp = interp1d(t_est, gofs_est)(t) assert_allclose(gof_est_interp, gofs, rtol=gof_rtol, atol=1e-7, err_msg="gof") # velocity calculation difference vel_est_interp = interp1d(t_est, vels_est)(t) assert_allclose(vel_est_interp, vels, atol=vel_atol, err_msg="velocity") def _decimate_chpi(raw, decim=4): """Decimate raw data (with aliasing) in cHPI-fitting compatible way.""" raw_dec = RawArray( raw._data[:, ::decim], raw.info, first_samp=raw.first_samp // decim ) with raw_dec.info._unlock(): raw_dec.info["sfreq"] /= decim for coil in raw_dec.info["hpi_meas"][0]["hpi_coils"]: if coil["coil_freq"] > raw_dec.info["sfreq"]: coil["coil_freq"] = np.mod(coil["coil_freq"], raw_dec.info["sfreq"]) if coil["coil_freq"] > raw_dec.info["sfreq"] / 2.0: coil["coil_freq"] = raw_dec.info["sfreq"] - coil["coil_freq"] return raw_dec # A shortcut method for testing that does both steps @verbose def _calculate_chpi_positions( raw, t_step_min=0.01, t_step_max=1.0, t_window="auto", too_close="raise", dist_limit=0.005, gof_limit=0.98, ext_order=1, verbose=None, ): chpi_amplitudes = compute_chpi_amplitudes( raw, t_step_min=t_step_min, t_window=t_window, ext_order=ext_order, verbose=verbose, ) chpi_locs = compute_chpi_locs( raw.info, chpi_amplitudes, t_step_max=t_step_max, too_close=too_close, verbose=verbose, ) head_pos = compute_head_pos( raw.info, chpi_locs, dist_limit=dist_limit, gof_limit=gof_limit, verbose=verbose ) return head_pos @testing.requires_testing_data def test_calculate_chpi_positions_preload(): """Test calculation of cHPI positions with and without data loaded.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes").crop(0, 2) kwargs = dict(t_step_min=0.1, t_window="auto", verbose=True) pos = compute_chpi_amplitudes(raw, **kwargs) raw.load_data() pos_preload = compute_chpi_amplitudes(raw, **kwargs) assert object_diff(pos, pos_preload) == "" @pytest.mark.slowtest @testing.requires_testing_data def test_calculate_chpi_positions_vv(): """Test calculation of cHPI positions.""" # Check to make sure our fits match MF decently mf_quats = read_head_pos(pos_fname) raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes") raw.crop(0, 5).load_data() # check "auto" t_window estimation at full sampling rate with catch_logging() as log: compute_chpi_amplitudes( raw, t_step_min=0.1, t_window="auto", tmin=0, tmax=2, verbose=True ) assert "83.3 ms" in log.getvalue() # This is a little hack (aliasing while decimating) to make it much faster # for testing purposes only. We can relax this later if we find it breaks # something. raw_dec = _decimate_chpi(raw, 15) with catch_logging() as log: with pytest.warns(RuntimeWarning, match="cannot determine"): py_quats = _calculate_chpi_positions(raw_dec, t_window=0.2, verbose="debug") log = log.getvalue() assert "\nHPIFIT" in log assert "Computing 4385 HPI location guesses" in log _assert_quats(py_quats, mf_quats, dist_tol=0.001, angle_tol=0.7) # degenerate conditions raw_no_chpi = read_raw_fif(sample_fname) with pytest.raises(ValueError, match="No appropriate cHPI information"): _calculate_chpi_positions(raw_no_chpi) raw_bad = raw.copy() del raw_bad.info["hpi_meas"][0]["hpi_coils"][0]["coil_freq"] with pytest.raises(ValueError, match="No appropriate cHPI information"): _calculate_chpi_positions(raw_bad) raw_bad = raw.copy() for d in raw_bad.info["dig"]: if d["kind"] == FIFF.FIFFV_POINT_HPI: d["coord_frame"] = FIFF.FIFFV_COORD_UNKNOWN break with pytest.raises(RuntimeError, match="coordinate frame incorrect"): _calculate_chpi_positions(raw_bad) for d in raw_bad.info["dig"]: if d["kind"] == FIFF.FIFFV_POINT_HPI: d["coord_frame"] = FIFF.FIFFV_COORD_HEAD d["r"] = np.ones(3) raw_bad.crop(0, 1.0) picks = np.concatenate( [ np.arange(306, len(raw_bad.ch_names)), pick_types(raw_bad.info, meg=True)[::16], ] ) raw_bad.pick([raw_bad.ch_names[pick] for pick in picks]) with _record_warnings(), pytest.warns(RuntimeWarning, match="Discrepancy"): with catch_logging() as log_file: _calculate_chpi_positions(raw_bad, t_step_min=1.0, verbose=True) # ignore HPI info header and [done] footer assert "0/5 good HPI fits" in log_file.getvalue() # half the rate cuts off cHPI coils with raw.info._unlock(): raw.info["lowpass"] /= 2.0 with pytest.raises(RuntimeError, match="above the"): _calculate_chpi_positions(raw) @testing.requires_testing_data @pytest.mark.slowtest def test_calculate_chpi_snr(): """Test cHPI SNR calculation.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes") # include handling of NaN (when cHPI was off at the beginning) raw.load_data() raw.info["bads"] = ["MEG0342", "MEG1443"] stop = int(round(raw.info["sfreq"])) * 2 raw._data[raw.ch_names.index("STI201"), :stop] = 0 result = compute_chpi_snr(raw) # make sure all the entries are there keys = { f"{ch_type}_{key}" for ch_type in ("mag", "grad") for key in ("snr", "power", "resid") } assert set(result) == keys.union({"times", "freqs"}) # make sure the values are plausible, given the sample data file n_pts = len(raw.times) // int(round(raw.info["sfreq"] * 0.01)) # our logic in this test for this length is not perfect assert_allclose(result["mag_snr"].shape[0], n_pts, atol=5) n_nan = np.where(result["times"] <= raw.first_time + 2)[0][-1] assert_allclose(result["mag_snr"][:n_nan], np.nan) assert result["mag_snr"][n_nan:].min() > 1 assert result["mag_snr"][n_nan:].max() < 40 assert result["grad_snr"][n_nan:].min() > 1 assert result["grad_snr"][n_nan:].max() < 40 @testing.requires_testing_data @pytest.mark.slowtest def test_calculate_chpi_positions_artemis(): """Test on 5k artemis data.""" raw = read_raw_artemis123(art_fname, preload=True) mf_quats = read_head_pos(art_mc_fname) mf_quats[:, 8:] /= 100 # old code errantly had this factor py_quats = _calculate_chpi_positions(raw, t_step_min=2.0, verbose="debug") _assert_quats( py_quats, mf_quats, dist_tol=0.001, angle_tol=1.0, err_rtol=0.7, vel_atol=1e-2 ) @testing.requires_testing_data def test_warn_maxwell_filtered(): """Test that trying to compute locations on Maxwell filtered data warns.""" raw = read_raw_fif(sss_fif_fname).crop(0, 1) with pytest.warns(RuntimeWarning, match="Maxwell filter"): amps = compute_chpi_amplitudes(raw) assert len(amps["times"]) > 0 # but for this file, it does work! @testing.requires_testing_data def test_initial_fit_redo(): """Test that initial fits can be redone based on moments.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes") slopes = np.array([[c["slopes"] for c in raw.info["hpi_meas"][0]["hpi_coils"]]]) amps = np.linalg.norm(slopes, axis=-1) amps /= slopes.shape[-1] assert_array_less(amps, 5e-11) assert_array_less(1e-12, amps) proj, _, _ = _setup_ext_proj(raw.info, ext_order=1) chpi_amplitudes = dict(times=np.zeros(1), slopes=slopes, proj=proj) chpi_locs = compute_chpi_locs(raw.info, chpi_amplitudes) # check GOF coil_gof = raw.info["hpi_results"][0]["goodness"] assert_allclose(chpi_locs["gofs"][0], coil_gof, atol=0.3) # XXX not good # check moment # XXX our forward and theirs differ by an extra mult by _MAG_FACTOR coil_moment = raw.info["hpi_results"][0]["moments"] / _MAG_FACTOR py_moment = chpi_locs["moments"][0] coil_amp = np.linalg.norm(coil_moment, axis=-1, keepdims=True) py_amp = np.linalg.norm(py_moment, axis=-1, keepdims=True) assert_allclose(coil_amp, py_amp, rtol=0.2) coil_ori = coil_moment / coil_amp py_ori = py_moment / py_amp angles = np.rad2deg(np.arccos(np.abs(np.sum(coil_ori * py_ori, axis=1)))) assert_array_less(angles, 20) # check resulting dev_head_t head_pos = compute_head_pos(raw.info, chpi_locs) assert head_pos.shape == (1, 10) nm_pos = raw.info["dev_head_t"]["trans"] dist = 1000 * np.linalg.norm(nm_pos[:3, 3] - head_pos[0, 4:7]) assert 0.1 < dist < 2 angle = np.rad2deg( _angle_between_quats(rot_to_quat(nm_pos[:3, :3]), head_pos[0, 1:4]) ) assert 0.1 < angle < 2 gof = head_pos[0, 7] assert_allclose(gof, 0.9999, atol=1e-4) @testing.requires_testing_data def test_calculate_head_pos_chpi_on_chpi5_in_one_second_steps(): """Comparing estimated cHPI positions with MF results (one second).""" # Check to make sure our fits match MF decently mf_quats = read_head_pos(chpi5_pos_fname) raw = read_raw_fif(chpi5_fif_fname, allow_maxshield="yes") # the last two seconds contain a maxfilter problem! # fiff file timing: 26. to 43. s # maxfilter estimates a wrong head position for interval 16: 41.-42. s raw = _decimate_chpi(raw.crop(0.0, 10.0).load_data(), decim=8) # needs no interpolation, because maxfilter pos files comes with 1 s steps py_quats = _calculate_chpi_positions( raw, t_step_min=1.0, t_step_max=1.0, t_window=1.0, verbose="debug" ) _assert_quats( py_quats, mf_quats, dist_tol=0.002, angle_tol=1.2, vel_atol=3e-3 ) # 3 mm/s @pytest.mark.slowtest @testing.requires_testing_data def test_calculate_head_pos_chpi_on_chpi5_in_shorter_steps(): """Comparing estimated cHPI positions with MF results (smaller steps).""" # Check to make sure our fits match MF decently mf_quats = read_head_pos(chpi5_pos_fname) raw = read_raw_fif(chpi5_fif_fname, allow_maxshield="yes") raw = _decimate_chpi(raw.crop(0.0, 5.0).load_data(), decim=8) with pytest.warns(RuntimeWarning, match="cannot determine"): py_quats = _calculate_chpi_positions( raw, t_step_min=0.1, t_step_max=0.1, t_window=0.1, verbose="debug" ) # needs interpolation, tolerance must be increased _assert_quats( py_quats, mf_quats, dist_tol=0.002, angle_tol=1.2, vel_atol=0.02 ) # 2 cm/s is not great but probably fine def test_simulate_calculate_head_pos_chpi(): """Test calculation of cHPI positions with simulated data.""" # Read info dict from raw FIF file info = read_info(raw_fname) # Tune the info structure chpi_channel = "STI201" ncoil = len(info["hpi_results"][0]["order"]) coil_freq = 10 + np.arange(ncoil) * 5 hpi_subsystem = { "event_channel": chpi_channel, "hpi_coils": [ {"event_bits": np.array([256, 0, 256, 256], dtype=np.int32)}, {"event_bits": np.array([512, 0, 512, 512], dtype=np.int32)}, {"event_bits": np.array([1024, 0, 1024, 1024], dtype=np.int32)}, {"event_bits": np.array([2048, 0, 2048, 2048], dtype=np.int32)}, ], "ncoil": ncoil, } with info._unlock(): info["hpi_subsystem"] = hpi_subsystem for fi, freq in enumerate(coil_freq): info["hpi_meas"][0]["hpi_coils"][fi]["coil_freq"] = freq info["sfreq"] = 100.0 # this will speed it up a lot picks = pick_types(info, meg=True, stim=True, eeg=False, exclude=[]) info = pick_info(info, picks) info["chs"][info["ch_names"].index("STI 001")]["ch_name"] = "STI201" info._update_redundant() with info._unlock(): info["projs"] = [] info_trans = info["dev_head_t"]["trans"].copy() dev_head_pos_ini = np.concatenate( [rot_to_quat(info_trans[:3, :3]), info_trans[:3, 3]] ) ez = np.array([0, 0, 1]) # Unit vector in z-direction of head coordinates # Define some constants duration = 10 # Time / s # Quotient of head position sampling frequency # and raw sampling frequency head_pos_sfreq_quotient = 0.01 # Round number of head positions to the next integer S = int(duration * info["sfreq"] * head_pos_sfreq_quotient) assert S == 10 dz = 0.001 # Shift in z-direction is 0.1mm for each step dev_head_pos = np.zeros((S, 10)) dev_head_pos[:, 0] = np.arange(S) * info["sfreq"] * head_pos_sfreq_quotient dev_head_pos[:, 1:4] = dev_head_pos_ini[:3] dev_head_pos[:, 4:7] = dev_head_pos_ini[3:] + np.outer(np.arange(S) * dz, ez) dev_head_pos[:, 7] = 1.0 # m/s dev_head_pos[:, 9] = dz / (info["sfreq"] * head_pos_sfreq_quotient) # Round number of samples to the next integer raw_data = np.zeros((len(picks), int(duration * info["sfreq"] + 0.5))) raw = RawArray(raw_data, info) add_chpi(raw, dev_head_pos) quats = _calculate_chpi_positions( raw, t_step_min=raw.info["sfreq"] * head_pos_sfreq_quotient, t_step_max=raw.info["sfreq"] * head_pos_sfreq_quotient, t_window=1.0, ) _assert_quats( quats, dev_head_pos, dist_tol=0.001, angle_tol=1.0, vel_atol=4e-3 ) # 4 mm/s def _calculate_chpi_coil_locs(raw, verbose): """Wrap to facilitate change diff.""" chpi_amplitudes = compute_chpi_amplitudes(raw, verbose=verbose) chpi_locs = compute_chpi_locs(raw.info, chpi_amplitudes, verbose=verbose) return _chpi_locs_to_times_dig(chpi_locs) def _check_dists(info, cHPI_digs, n_bad=0, bad_low=0.02, bad_high=0.04): __tracebackhide__ = True orig = _get_hpi_initial_fit(info) hpi_coil_distances = cdist(orig, orig) new_pos = np.array([d["r"] for d in cHPI_digs]) mask, distances = _compute_good_distances(hpi_coil_distances, new_pos) good_idx = np.where(mask)[0] assert len(good_idx) >= 3 meds = np.empty(len(orig)) for ii in range(len(orig)): idx = np.setdiff1d(good_idx, ii) meds[ii] = np.median(distances[ii][idx]) meds = np.array(meds) assert_array_less(meds[good_idx], 0.003) bad_idx = np.where(~mask)[0] if len(bad_idx): bads = meds[bad_idx] assert_array_less(bad_low, bads) assert_array_less(bads, bad_high) @pytest.mark.slowtest @testing.requires_testing_data def test_calculate_chpi_coil_locs_artemis(): """Test computing just cHPI locations.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes", preload=True) # This is a little hack (aliasing while decimating) to make it much faster # for testing purposes only. We can relax this later if we find it breaks # something. raw_dec = _decimate_chpi(raw, 15) times, cHPI_digs = _calculate_chpi_coil_locs(raw_dec, verbose="debug") # spot check assert_allclose(times[0], 9.0, atol=1e-2) assert_allclose( cHPI_digs[0][2]["r"], [-0.01937833, 0.00346804, 0.06331209], atol=1e-3 ) assert_allclose(cHPI_digs[0][2]["gof"], 0.9957, atol=1e-3) assert_allclose(cHPI_digs[0][4]["r"], [-0.0655, 0.0755, 0.0004], atol=3e-3) assert_allclose(cHPI_digs[0][4]["gof"], 0.9323, atol=1e-3) _check_dists(raw.info, cHPI_digs[0], n_bad=1) # test on 5k artemis data raw = read_raw_artemis123(art_fname, preload=True) times, cHPI_digs = _calculate_chpi_coil_locs(raw, verbose="debug") assert len(np.setdiff1d(times, raw.times + raw.first_time)) == 0 # Should be somewhere around 1.5 s, depending on coil GOF values # around 0.98 it can change assert_allclose(times[5], 1.5, atol=2e-1) assert_allclose(cHPI_digs[5][0]["gof"], 0.995, atol=5e-3) assert_allclose(cHPI_digs[5][0]["r"], [-0.0157, 0.0655, 0.0018], atol=1e-3) _check_dists(raw.info, cHPI_digs[5]) coil_amplitudes = compute_chpi_amplitudes(raw) with pytest.raises(ValueError, match="too_close"): compute_chpi_locs(raw.info, coil_amplitudes, too_close="foo") # ensure values are in a reasonable range amps = np.linalg.norm(coil_amplitudes["slopes"], axis=-1) amps /= coil_amplitudes["slopes"].shape[-1] assert amps.shape == (len(coil_amplitudes["times"]), 3) assert_array_less(amps, 1e-11) assert_array_less(1e-13, amps) # with nan amplitudes (i.e., cHPI off) it should return an empty array, # but still one that is 3D coil_amplitudes["slopes"].fill(np.nan) chpi_locs = compute_chpi_locs(raw.info, coil_amplitudes) assert chpi_locs["rrs"].shape == (0, 3, 3) pos = compute_head_pos(raw.info, chpi_locs) assert pos.shape == (0, 10) def assert_suppressed(new, old, suppressed, retained): """Assert that some frequencies are suppressed and others aren't.""" __tracebackhide__ = True from scipy.signal import welch picks = pick_types(new.info, meg="grad") sfreq = new.info["sfreq"] new = new.get_data(picks) old = old.get_data(picks) f, new = welch(new, sfreq, "hann", nperseg=1024) _, old = welch(old, sfreq, "hann", nperseg=1024) new = np.median(new, axis=0) old = np.median(old, axis=0) for freqs, lim in ((suppressed, (10, 60)), (retained, (-3, 3))): for freq in freqs: fidx = np.argmin(np.abs(f - freq)) this_new = np.median(new[fidx]) this_old = np.median(old[fidx]) suppression = -10 * np.log10(this_new / this_old) assert lim[0] < suppression < lim[1], freq @testing.requires_testing_data def test_chpi_subtraction_filter_chpi(): """Test subtraction of cHPI signals.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes", preload=True) raw.info["bads"] = ["MEG0111"] raw.del_proj() raw_orig = raw.copy().crop(0, 16) with catch_logging() as log: filter_chpi(raw, include_line=False, t_window=0.2, verbose=True) log = log.getvalue() assert "No average EEG" not in log assert "5 cHPI" in log # MaxFilter doesn't do quite as well as our algorithm with the last bit raw.crop(0, 16) # remove cHPI status chans raw_c = read_raw_fif(sss_hpisubt_fname).crop(0, 16).load_data() raw_c.pick(["meg", "eeg", "eog", "ecg", "stim", "misc"]) assert_meg_snr(raw, raw_c, 143, 624) # cHPI suppressed but not line freqs (or others) assert_suppressed(raw, raw_orig, np.arange(83, 324, 60), [30, 60, 150]) raw = raw_orig.copy() with catch_logging() as log: filter_chpi(raw, include_line=True, t_window=0.2, verbose=True) log = log.getvalue() assert "5 cHPI" in log assert "6 line" in log # cHPI and line freqs suppressed suppressed = np.sort( np.concatenate( [ np.arange(83, 324, 60), np.arange(60, 301, 60), ] ) ) assert_suppressed(raw, raw_orig, suppressed, [30, 150]) # No HPI information raw = read_raw_fif(sample_fname, preload=True) raw_orig = raw.copy() assert raw.info["line_freq"] is None with pytest.raises(RuntimeError, match="line_freq.*consider setting it"): filter_chpi(raw, t_window=0.2) with raw.info._unlock(): raw.info["line_freq"] = 60.0 with pytest.raises(ValueError, match="No appropriate cHPI information"): filter_chpi(raw, t_window=0.2) # but this is allowed with catch_logging() as log: filter_chpi(raw, t_window="auto", allow_line_only=True, verbose=True) log = log.getvalue() assert "0 cHPI" in log assert "1 line" in log # Our one line freq suppressed but not others assert_suppressed(raw, raw_orig, [60], [30, 45, 75]) # When MaxFliter downsamples, like:: # $ maxfilter -nosss -ds 2 -f test_move_anon_raw.fif \ # -o test_move_anon_ds2_raw.fif # it can strip out some values of info, which we emulate here: raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes") raw = raw.crop(0, 1).load_data().resample(600.0, npad="auto") with raw.info._unlock(): raw.info["lowpass"] = 200.0 del raw.info["maxshield"] del raw.info["hpi_results"][0]["moments"] del raw.info["hpi_subsystem"]["event_channel"] with catch_logging() as log: filter_chpi(raw, t_window="auto", verbose=True) with pytest.raises(ValueError, match="must be > 0"): filter_chpi(raw, t_window=-1) assert "2 cHPI" in log.getvalue() @testing.requires_testing_data def test_calculate_head_pos_ctf(tmp_path): """Test extracting of cHPI positions from CTF data.""" raw = read_raw_ctf(ctf_chpi_fname) chpi_locs = extract_chpi_locs_ctf(raw) quats = compute_head_pos(raw.info, chpi_locs) mc_quats = read_head_pos(ctf_chpi_pos_fname) mc_quats[:, 9] /= 10000 # had old factor in there twice somehow... _assert_quats( quats, mc_quats, dist_tol=0.004, angle_tol=2.5, err_rtol=1.0, vel_atol=7e-3 ) # 7 mm/s plot_head_positions(quats, info=raw.info) with pytest.raises(RuntimeError, match="Could not find"): extract_chpi_locs_ctf(read_raw_fif(ctf_fname)) # save-load should not affect result fname_temp = tmp_path / "test_ctf_raw.fif" raw.save(fname_temp) raw_read = read_raw_fif(fname_temp) # the two attributes used by compute_head_pos assert_allclose( raw.info["dev_head_t"]["trans"], raw_read.info["dev_head_t"]["trans"] ) with pytest.warns(RuntimeWarning, match="is poor"): head_rrs = _get_hpi_initial_fit(raw.info, verbose="debug") with pytest.warns(RuntimeWarning, match="is poor"): head_rrs_2 = _get_hpi_initial_fit(raw_read.info, verbose="debug") assert_allclose(head_rrs, head_rrs_2, atol=1e-5) quats_2 = compute_head_pos(raw_read.info, chpi_locs) _assert_quats(quats, quats_2, dist_tol=1e-5, angle_tol=0.1) chpi_locs_2 = extract_chpi_locs_ctf(raw_read) assert_allclose(chpi_locs["rrs"], chpi_locs_2["rrs"], atol=1e-5) quats_3 = compute_head_pos(raw_read.info, chpi_locs_2) _assert_quats(quats, quats_3, dist_tol=1e-5, angle_tol=0.1) @testing.requires_testing_data def test_calculate_head_pos_kit(): """Test calculation of head position using KIT data.""" raw = read_raw_kit(con_fname, mrk_fname, elp_fname, hsp_fname) assert len(raw.info["hpi_results"]) == 1 chpi_locs = extract_chpi_locs_kit(raw) assert chpi_locs["rrs"].shape == (2, 5, 3) assert_array_less(chpi_locs["gofs"], 1.0) assert_array_less(0.98, chpi_locs["gofs"]) quats = compute_head_pos(raw.info, chpi_locs) assert quats.shape == (2, 10) # plotting works plot_head_positions(quats, info=raw.info) raw_berlin = read_raw_kit(berlin_fname) assert_allclose(raw_berlin.info["dev_head_t"]["trans"], np.eye(4)) assert len(raw_berlin.info["hpi_results"]) == 0 with pytest.raises(ValueError, match="Invalid value"): extract_chpi_locs_kit(raw_berlin) with pytest.raises(RuntimeError, match="not find appropriate"): extract_chpi_locs_kit(raw_berlin, "STI 014") with pytest.raises(RuntimeError, match="no initial cHPI"): compute_head_pos(raw_berlin.info, chpi_locs) @testing.requires_testing_data def test_get_active_chpi_ctf(): """Test extracting of cHPI positions from CTF data.""" raw = read_raw_ctf(ctf_chpi_fname) with pytest.raises(NotImplementedError, match="not implemented for other systems"): get_active_chpi(raw) @testing.requires_testing_data def test_get_active_chpi_neuromag(): """Test extracting of cHPI positions from neuromag data.""" raw = read_raw_fif(chpi_fif_fname, allow_maxshield="yes", preload=True) status_ch = raw.ch_names.index("STI201") # make artificial chpi signal first_three_on = 256 + 512 + 1024 all_on = 256 + 512 + 1024 + 2048 + 4096 raw._data[status_ch][:1000] = 0 raw._data[status_ch][1000:2000] = first_three_on raw._data[status_ch][2000:] = all_on # build target signal target_signal = 5 * np.ones_like(raw.times) target_signal[:1000] = 0 target_signal[1000:2000] = 3 assert_allclose(get_active_chpi(raw, on_missing="ignore"), target_signal) raw_no_chpi = read_raw_fif(sample_fname) assert_allclose( get_active_chpi(raw_no_chpi, on_missing="ignore"), np.zeros_like(raw_no_chpi.times), )