eyelink sync updates, first commit

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
 

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,25 +219,99 @@ 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,
         subject=subject,
         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)