diff --git a/mne_bids_pipeline/_config.py b/mne_bids_pipeline/_config.py index 17b9052b2..ef3a0fff2 100644 --- a/mne_bids_pipeline/_config.py +++ b/mne_bids_pipeline/_config.py @@ -1136,11 +1136,85 @@ ``` """ +# ## Sync with eyelink +# +# These parameters determine how eyelink syncing is performed and reported + sync_eyelink: bool = False +""" +Carry out eyelink syncing if True +""" remove_blink_saccades: bool = True -sync_eventtype_regex: str = ".*" +""" +Currently not implemented +""" + +sync_eventtype_regex: str = "" +""" +Regular expression which will be used to select events in the EEG file for synchronisation + +???+ example "Example" + ```python + sync_eventtype_regex = ".*_grab_this" # select all events with some text, immediately followed by "_grab_this" + ``` +""" + sync_eventtype_regex_et: str = "" +""" +Regular expression which will be used to select events in the eye-tracking file for synchronisation + +???+ example "Example" + ```python + sync_eventtype_regex_et = ".*_grab_that" # select all events with some text, immediately followed by "_grab_that" + ``` +""" +sync_heog_ch: tuple[str, str] | str | None = None +""" +The HEOG channel(s) to use for assessing synchronisation by means of cross correlation + +???+ example "Example" + ```python + sync_heog_ch = ["HEOGL", "HEOGR"] # subtract the second channel from the first to make a bipolar HEOG channel + sync_heog_ch = "HEOGL" # use the HEOGL channel + sync_heog_ch = None # skip HEOG cross correlation + ``` +""" + +sync_et_ch: tuple[str, str] | str | None = ("xpos_left", "xpos_right") +""" +The ET channel to use for assessing synchronisation by means of cross correlation + +???+ example "Example" + ```python + + sync_et_ch = "xpos_left" # use the xpos_left channel + sync_et_ch = None # skip HEOG cross correlation + ``` +""" + +sync_heog_highpass: float = 0.1 +""" +High-pass filter to apply to the HEOG channel. Filtered signal will not be saved. +""" + +sync_heog_lowpass: float = 10 +""" +Low-pass filter to apply to the HEOG channel. Filtered signal will not be saved. +""" + +sync_plot_samps: int | tuple[int, int] | None = None +""" +When displaying HEOG-ET cross-correlation, constrict plotting to sync_plot_samps[0] to sync_plot_samps[1] around the midpoint + +???+ example "Example" + ```python + + sync_plot_samps = 500 # plot -500 to 500 around the midpoint + sync_plot_samps = (-300, 500) # plot -300 to 500 around the midpoint + sync_plot_samps = None # plot the entire cross correlation + ``` +""" # ### SSP, ICA, and artifact regression diff --git a/mne_bids_pipeline/steps/preprocessing/_05b_sync_eyelink.py b/mne_bids_pipeline/steps/preprocessing/_05b_sync_eyelink.py index 136aecf2a..d8ec236f9 100644 --- a/mne_bids_pipeline/steps/preprocessing/_05b_sync_eyelink.py +++ b/mne_bids_pipeline/steps/preprocessing/_05b_sync_eyelink.py @@ -19,6 +19,43 @@ from ..._report import _open_report from ..._run import _prep_out_files, _update_for_splits, failsafe_run, save_logs +def _check_HEOG_ET_vars(cfg): + # helper function for sorting out heog and et channels + bipolar = False + if isinstance(cfg.sync_heog_ch, tuple): + heog_ch = "bi_HEOG" + bipolar = True + else: + heog_ch = cfg.sync_heog_ch + + if isinstance(cfg.sync_et_ch, tuple): + et_ch = list(cfg.sync_et_ch) + else: + et_ch = [cfg.sync_et_ch] + + return heog_ch, et_ch, bipolar + +def _mark_calibration_as_bad(raw): + # marks recalibration beginnings and ends as one bad segment + cur_idx = None + cur_start_time = 0. + for annot in raw.annotations: + calib_match = re.match(".* Recalibration (start|end) \\| (\\d*)", annot["description"]) + if not calib_match: continue + + calib_status, calib_idx = calib_match.group(1), calib_match.group(2) + if calib_idx == cur_idx and calib_status == "end": + duration = annot["onset"] - cur_start_time + raw.annotations.append(cur_start_time, duration, f"BAD_Recalibrate {calib_idx}") + cur_idx, cur_start_time = None, 0. + elif calib_status == "start" and cur_idx is None: + cur_idx = calib_idx + cur_start_time = annot["onset"] + elif calib_status == "start" and cur_idx is not None: + raise ValueError(f"Annotation {annot["description"]} could not be assigned membership") + + return raw + def get_input_fnames_sync_eyelink( *, @@ -104,6 +141,7 @@ def sync_eyelink( ) -> dict: """Run Sync for Eyelink.""" import matplotlib.pyplot as plt + from scipy.signal import correlate raw_fnames = [in_files.pop(f"raw_run-{run}") for run in cfg.runs] et_fnames = [in_files.pop(f"et_run-{run}") for run in cfg.runs] @@ -114,8 +152,6 @@ def sync_eyelink( bids_basename = raw_fnames[0].copy().update(processing=None, split=None, run=None) out_files["eyelink"] = bids_basename.copy().update(processing="eyelink", suffix="raw") del bids_basename - - for idx, (run, raw_fname,et_fname,et_edf_fname) in enumerate(zip(cfg.runs, raw_fnames,et_fnames,et_edf_fnames)): msg = f"Syncing eyelink data (fake for now) {raw_fname.basename}" @@ -128,7 +164,7 @@ def sync_eyelink( import subprocess subprocess.run(["edf2asc", et_edf_fname]) # TODO: Still needs to be tested - raw_et = mne.io.read_raw_eyelink(et_fname,find_overlaps=True) + raw_et = mne.io.read_raw_eyelink(et_fname, find_overlaps=True) # If the user did not specify a regular expression for the eye-tracking sync events, it is assumed that it's # identical to the regex for the EEG sync events @@ -137,7 +173,6 @@ def sync_eyelink( et_sync_times = [annotation["onset"] for annotation in raw_et.annotations if re.search(cfg.sync_eventtype_regex_et,annotation["description"])] sync_times = [annotation["onset"] for annotation in raw.annotations if re.search(cfg.sync_eventtype_regex, annotation["description"])] - assert len(et_sync_times) == len(sync_times),f"Detected eyetracking and EEG sync events were not of equal size ({len(et_sync_times)} vs {len(sync_times)}). Adjust your regular expressions via 'sync_eventtype_regex_et' and 'sync_eventtype_regex' accordingly" #logger.info(**gen_log_kwargs(message=f"{et_sync_times}")) #logger.info(**gen_log_kwargs(message=f"{sync_times}")) @@ -154,18 +189,24 @@ def sync_eyelink( #mne.preprocessing.eyetracking.interpolate_blinks(raw_et, buffer=(0.05, 0.05), interpolate_gaze=True) - # Align the data mne.preprocessing.realign_raw(raw, raw_et, sync_times, et_sync_times) - # Add ET data to EEG raw.add_channels([raw_et], force_update_info=True) raw._raw_extras.append(raw_et._raw_extras) # Also add ET annotations to EEG - raw.set_annotations(mne.annotations._combine_annotations(raw.annotations,raw_et.annotations,0,raw.first_samp,raw_et.first_samp,raw.info["sfreq"])) - + # first mark et sync event descriptions so we can differentiate them later + for idx, desc in enumerate(raw_et.annotations.description): + if re.search(cfg.sync_eventtype_regex_et, desc): + raw_et.annotations.description[idx] = "ET_" + desc + raw.set_annotations(mne.annotations._combine_annotations(raw.annotations, + raw_et.annotations, + 0, + raw.first_samp, + raw_et.first_samp, + raw.info["sfreq"])) msg = f"Saving synced data to disk." logger.info(**gen_log_kwargs(message=msg)) @@ -178,11 +219,91 @@ def sync_eyelink( # no idea what the split stuff is... _update_for_splits(out_files, "eyelink") # TODO: Find out if we need to add this or not - # Add to report + fig, axes = plt.subplots(2, 2, figsize=(19.2, 19.2)) + msg = f"Adding figure to report." + logger.info(**gen_log_kwargs(message=msg)) tags = ("sync", "eyelink") title = "Synchronize Eyelink" + caption = ( + f"The `realign_raw` function from MNE was used to align an Eyelink `asc` file to the M/EEG file." + f"The Eyelink-data was added as annotations and appended as new channels." + ) + if cfg.sync_heog_ch is None or cfg.sync_et_ch is None: + # we need both an HEOG channel and ET channel specified to do cross-correlation + msg = f"HEOG and/or ET channel not specified; cannot produce cross-correlation for report." + logger.info(**gen_log_kwargs(message=msg)) + caption += "\nHEOG and/or eye tracking channels were not specified and no cross-correlation was performed." + axes[0,0].text(0.5, 0.5, 'HEOG/ET cross-correlation unavailable', fontsize=34, + horizontalalignment='center', verticalalignment='center', transform=axes[0,0].transAxes) + axes[0,0].axis("off") + else: + # return _prep_out_files(exec_params=exec_params, out_files=out_files) + # calculate cross correlation of HEOG with ET + heog_ch, et_ch, bipolar = _check_HEOG_ET_vars(cfg) + if bipolar: + # create bipolar HEOG + raw = mne.set_bipolar_reference(raw, *cfg.sync_heog_ch, ch_name=heog_ch, drop_refs=False) + raw.filter(l_freq=cfg.sync_heog_highpass, h_freq=cfg.sync_heog_lowpass, picks=heog_ch) # get rid of drift and high freq noise + _mark_calibration_as_bad(raw) + # extract HEOG and ET as arrays + heog_array = raw.get_data(picks=[heog_ch], reject_by_annotation="omit") + et_array = raw.get_data(picks=et_ch, reject_by_annotation="omit") + if len(et_array) > 1: + et_array = et_array.mean(axis=0, keepdims=True) + # cross correlate them + corr = correlate(heog_array[0], et_array[0], mode="same") / heog_array.shape[1] + # plot cross correlation + # figure out how much we plot + midpoint = len(corr) // 2 + plot_samps = (-cfg.sync_plot_samps, cfg.sync_plot_samps) if isinstance(cfg.sync_plot_samps, int) else cfg.sync_plot_samps + if isinstance(plot_samps, tuple): + x_range = np.arange(plot_samps[0], plot_samps[1]) + y_range = np.arange(midpoint+plot_samps[0], midpoint+plot_samps[1]) + else: # None + y_range = np.arange(len(corr)) + x_range = y_range - midpoint + # plot + axes[0,0].plot(x_range, corr[y_range], color="black") + axes[0,0].axvline(linestyle="--", alpha=0.3) + axes[0,0].set_title("Cross correlation HEOG and ET") + axes[0,0].set_xlabel("Samples") + axes[0,0].set_ylabel("X correlation") + # calculate delay + delay_idx = abs(corr).argmax() - midpoint + delay_time = delay_idx * (raw.times[1] - raw.times[0]) + caption += f"\nThere was an estimated synchronisation delay of {delay_idx} samples ({delay_time:.3f} seconds.)" + + # regression between synced events + # we assume here that these annotations are sequential pairs of the same event in raw and et. otherwise this will break + raw_onsets = [annot["onset"] for annot in raw.annotations if re.match(cfg.sync_eventtype_regex, annot["description"])] + et_onsets = [annot["onset"] for annot in raw.annotations if re.match("ET_"+cfg.sync_eventtype_regex_et, annot["description"])] + + if len(raw_onsets) != len(et_onsets): + raise ValueError(f"Lengths of raw {len(raw_onsets)} and ET {len(et_onsets)} onsets do not match.") + # regress and plot + coef = np.polyfit(raw_onsets, et_onsets, 1) + preds = np.poly1d(coef)(raw_onsets) + resids = et_onsets - preds + axes[0,1].plot(raw_onsets, et_onsets, "o", alpha=0.3, color="black") + axes[0,1].plot(raw_onsets, preds, "--k") + axes[0,1].set_title("Regression") + axes[0,1].set_xlabel("Raw onsets (seconds)") + axes[0,1].set_ylabel("ET onsets (seconds)") + # residuals + axes[1,0].plot(np.arange(len(resids)), resids, "o", alpha=0.3, color="black") + axes[1,0].axhline(linestyle="--") + axes[1,0].set_title("Residuals") + axes[1,0].set_ylabel("Residual (seconds)") + axes[1,0].set_xlabel("Samples") + # histogram of distances between events in time + axes[1,1].hist(np.array(raw_onsets) - np.array(et_onsets), bins=11, range=(-5,5), color="black") + axes[1,1].set_title("Raw - ET event onset distances histogram") + axes[1,1].set_xlabel("milliseconds") + # this doesn't seem to help, though it should... + fig.tight_layout() + with _open_report( cfg=cfg, exec_params=exec_params, @@ -190,13 +311,7 @@ def sync_eyelink( session=session, task=cfg.task, ) as report: - - - caption = ( - f"The `realign_raw` function from MNE was used to align an Eyelink `asc` file to the M/EEG file." - f"The Eyelink-data was added as annotations and appended as new channels." - ) - fig = raw_et.plot(scalings=dict(eyegaze=1e3)) + caption = caption report.add_figure( fig=fig, title="Eyelink data", @@ -211,9 +326,6 @@ def sync_eyelink( - - - def get_config( *, config: SimpleNamespace, @@ -227,6 +339,11 @@ def get_config( remove_blink_saccades = config.remove_blink_saccades, sync_eventtype_regex = config.sync_eventtype_regex, sync_eventtype_regex_et = config.sync_eventtype_regex_et, + sync_heog_ch = config.sync_heog_ch, + sync_et_ch = config.sync_et_ch, + sync_heog_highpass = config.sync_heog_highpass, + sync_heog_lowpass = config.sync_heog_lowpass, + sync_plot_samps = config.sync_plot_samps, processing= "filt" if config.regress_artifact is None else "regress", _raw_split_size=config._raw_split_size, @@ -255,7 +372,4 @@ def main(*, config: SimpleNamespace) -> None: for subject in get_subjects(config) for session in get_sessions(config) ) - save_logs(config=config, logs=logs) - - - + save_logs(config=config, logs=logs) \ No newline at end of file