--- title: "Blink Detection" author: "Jason Geller" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Blink Detection} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, fig.width = 9, fig.height = 4.5, dpi = 110, out.width = "100%" ) ``` While it is generally agreed that blinks need to be taken care of during preprocessing, most eye tracking software does not include this event information. This has lead to the development of several great blink algorithms. gazeR includes two: 1. A **velocity-based** detector (`blink_velocity()`) — simple and intuitive; finds blinks by looking at how fast the pupil signal is changing. 2. A **noise-based** detector (`detect_blinks_noise()`) — implements Hershman, Henik & Cohen (2018); starts from the tracker's missing-data sentinel and refines each blink's onset and offset by walking outward on the smoothed first-difference signal. This vignette walks through both. The velocity-based steps below illustrate what is happening inside `blink_velocity()`; the final section shows how to call the noise-based detector directly. First load gazeR and some sample data. ```{r, warning=FALSE, message=FALSE} library(gazer) library(cowplot) library(data.table) library(patchwork) library(tidyverse) ``` ```{r,warning=FALSE, message=FALSE} pupil_path <- system.file("extdata", "pupil_sample_files_edf.xls.gz", package = "gazer") pupil_files1<-fread(pupil_path) pupil_files1 <- as_tibble(pupil_files1) summary(pupil_files1) ``` # Blink Detection First the uncorrected pupil signal for a single trial from a single subject is plotted as a continuous signal.The signal is measured in arbitrary units as outputted by Eyelink’s software (Please check your own data before using this). In the below figure, it is clear to see that the participant blinked three times during the trial. Eye blinks are characterized by a pronounced drop in the pupillary signal, followed by a full loss of signal, and then usually a recovery artifact when the signal comes back online. The data I am working with is a bit odd (positive deflections) which suggests some type of artifact. ```{r, warning=FALSE, message=FALSE} interp_graph <- pupil_files1 %>% dplyr::filter(subject=="11c.edf", trial=="1") pup_g<- ggplot(interp_graph, aes(x= time, y= pupil)) + geom_point()+ geom_line(colour="black") + ggtitle("Raw Pupil Signal") print(pup_g) ``` In order to detect these events systematically, I calculated a velocity profile for each trace. However, the original signal is too noisy to create a reliable profile. In order to detect blinks efficiently, I began by smoothing the signal using a weighted moving window average of 10 samples (i.e., 20 ms; 250 Hz tracker). ```{r, warning=FALSE, message=FALSE} pupil_blink_algo <- interp_graph %>% mutate(smooth_pupil=moving_average_pupil(pupil, n=10)) pup_g1<- ggplot(pupil_blink_algo, aes(x= time, y= smooth_pupil)) + geom_point()+ geom_line(colour="black") + ggtitle("Smoothed Pupil Signal") print(pup_g1) ``` Next I created a velocity profile by subtracting each sample from the immediately preceding sample in the signal. ```{r} pupil_blink_algo1 <- pupil_blink_algo %>% mutate(velocity_pupil=c(diff(smooth_pupil), NA)) pup_g2<- ggplot(pupil_blink_algo1, aes(x= time, y= velocity_pupil)) + geom_point()+ geom_line(colour="black") + ggtitle("Pupil Velocity") pup_g2 ``` Blink onsets were subsequently identified as occurring when the velocity crossed a predetermined negative threshold (I selected -5 based on Mathot's (2013) recommendation. This rapid decrease in pupil diameter corresponds to the apparent decrease in size of the pupil as the eyelid closes. Likewise, when the eyelid reopens there is a recovery artifact wherein pupil size rapidly gets larger. Thus, the algorithm detected the recovery period by indexing the time since onset that velocity exceed some positive threshold (I selected -5, again based off Matĥot (2013) recommendation). Finally, the offset was detected as the time at which velocity fell back down to 0. In this way, a blink corresponds to an onset, recovery, and offset index. According to Mathôt (2013), this detection algorithm underestimates the blink period by several milliseconds, thus I selected a margin value (10 ms) which was subtracted from the onset and added to the offset using the 'extend_blinks` function in gazeR. Finally, the pupil signal was then linearly interpolated. ```{r, warning=FALSE} pupil_blink_algo2 <- pupil_blink_algo1 %>% # Thrid, set neg and pos threshold here we are using -5 and 5 (probably good idea to look at data) mutate(blinks_onset_offset=ifelse(velocity_pupil <= -5 | velocity_pupil >= 5, 1, 0)) %>% # Four turn pupil values 0 if blink detection coded as 1 mutate(blinks_pupil=ifelse(blinks_onset_offset==1, pupil==NA, pupil)) %>% mutate(extendpupil=extend_blinks(blinks_pupil, fillback=10, fillforward=10, hz=250)) %>% #interpolate dplyr::mutate(interp = zoo::na.approx(extendpupil, na.rm = FALSE, rule=2)) ``` ```{r pressure, echo=FALSE, warning=FALSE} pup_g3 <- ggplot(pupil_blink_algo2, aes(x= time, y= pupil)) + geom_point()+ geom_line(colour="black") + geom_line(aes(x=time, y=interp), colour="purple") + xlab("Time (ms)") + ylab("Pupil Size (arbitrary units)") + ggtitle("Interpolated Pupil Signal") print(pup_g3) ``` By looking at interpolated data plotted over the raw data we can see that the algo did a pretty good job at finding the blinks and recovering the pupil signal. These steps are all contained in the `blink_mathot` function. ```{r, warning=FALSE, message=FALSE} pupil= plot_grid(pup_g + pup_g1 + pup_g2 + pup_g3) plot(pupil) ``` # Noise-based Blink Detection (Hershman, Henik & Cohen, 2018) The velocity-based detector above works well when blinks produce a clean drop / recovery pattern in the pupil signal. A complementary approach, due to Hershman, Henik & Cohen (2018), starts from the *missing-data* samples the tracker already emits (zeros on Tobii, `NA` on EyeLink) and refines each missing run outward — finding the last "calm" sample before the eyelid descent and the first calm sample after the recovery rise. Adjacent blinks that fall within a small gap are merged. `detect_blinks_noise()` is tracker-agnostic: by default it treats both `0` and `NA` as missing, so the same call works for Tobii and EyeLink-style exports without configuration. Before running detection, NA out samples whose pupil velocity exceeds what a real pupil can physically produce. `speed_pupil()` computes a per-sample dilation speed and `calc_mad()` derives a robust threshold (Kret & Sjak-Shie, 2019). Marking these high-velocity samples as missing *before* detection makes the eyelid descent / recovery part of each missing-data run, so the detector's onset/offset refinement has clean boundaries to anchor on. ```{r, warning=FALSE, message=FALSE} interp_graph <- interp_graph %>% dplyr::group_by(subject, trial) %>% dplyr::mutate(speed = speed_pupil(pupil, time)) %>% dplyr::ungroup() thr <- calc_mad(interp_graph$speed[!is.na(interp_graph$speed) & is.finite(interp_graph$speed)], n = 16) interp_graph <- interp_graph %>% dplyr::mutate(pupil = ifelse(!is.na(speed) & is.finite(speed) & speed > thr, NA_real_, pupil)) ``` Now run `detect_blinks_noise()`. It composes with `dplyr` the same way as `extend_blinks()` — call it inside `mutate()` under a `group_by()` so blink IDs restart per trial: ```{r, warning=FALSE, message=FALSE} blinks_noise <- interp_graph %>% dplyr::group_by(subject, trial) %>% dplyr::mutate(blink = detect_blinks_noise(pupil, hz = 250, smooth_ms = 20)) %>% dplyr::ungroup() ``` We're using `smooth_ms = 20` here (rather than the 10 ms default) because, on this recording, the default windows clip a bit of the eyelid descent and recovery — see the parameter sweep further down for the comparison. The result is an integer column `blink`: `0` outside detected blinks, `1, 2, …` inside each blink. ## What did the detector actually flag? Before looking at interpolated output, it's worth visualizing exactly which samples the algorithm marked as blink. The figure below shades each detected blink region in light orange on top of the raw trace. The dotted vertical lines mark where the tracker reported missing samples (`NA`s) — note how `detect_blinks_noise()` widens each missing-data run outward to capture the eyelid descent and recovery on either side. ```{r, warning=FALSE, message=FALSE, fig.height=5} library(dplyr) # One-row-per-blink summary: start time, end time, blink id. blink_windows <- blinks_noise %>% dplyr::filter(blink > 0) %>% dplyr::group_by(blink) %>% dplyr::summarise(start = min(time), end = max(time), .groups = "drop") # Times where the tracker itself reported a missing sample. na_times <- blinks_noise %>% dplyr::filter(is.na(pupil)) %>% dplyr::pull(time) pup_noise_regions <- ggplot(blinks_noise, aes(x = time, y = pupil)) + geom_rect(data = blink_windows, aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf), inherit.aes = FALSE, fill = "darkorange", alpha = 0.18) + geom_vline(xintercept = na_times, colour = "grey50", linetype = "dotted", alpha = 0.4) + geom_line(colour = "black") + geom_point(size = 0.6) + xlab("Time (ms)") + ylab("Pupil size (arbitrary units)") + ggtitle(sprintf("Detected blink regions (n = %d)", nrow(blink_windows))) + theme_bw(base_size = 13) print(pup_noise_regions) ``` The orange bands should bracket each blink: a stretch of pupil decrease, the missing-data run, and the recovery rise. ## Zooming in on one blink To see what the onset/offset refinement is doing, zoom into the first detected blink window with ±200 ms of context: ```{r, warning=FALSE, message=FALSE, fig.height=4} focus <- blink_windows[1, ] window <- blinks_noise %>% dplyr::filter(time >= focus$start - 200, time <= focus$end + 200) pup_zoom <- ggplot(window, aes(x = time, y = pupil)) + geom_rect(data = focus, aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf), inherit.aes = FALSE, fill = "darkorange", alpha = 0.22) + geom_vline(xintercept = window$time[is.na(window$pupil)], colour = "grey40", linetype = "dotted") + geom_line(colour = "black") + geom_point(size = 1.4) + xlab("Time (ms)") + ylab("Pupil size (arbitrary units)") + ggtitle("Blink #1: refined window (orange) vs. missing samples (dotted)") + theme_bw(base_size = 13) print(pup_zoom) ``` The dotted vertical lines mark the raw missing samples; the orange band shows the *refined* blink — the algorithm walks outward on the smoothed first-difference signal until it finds the last calm sample before the descent and the first calm sample after the recovery. ## Interpolated output To interpolate across the detected blinks, route the data through gazeR's `interpolate_pupil()`. The function expects a binary `blink` column (`blink == 1` inside blinks), so binarize the integer IDs from `detect_blinks_noise()` first. The interpolated values land in a new `interp` column. ```{r, warning=FALSE, message=FALSE, fig.height=5} blinks_noise <- blinks_noise %>% dplyr::mutate(blink = as.integer(blink > 0)) %>% interpolate_pupil(pupil = "pupil", subject = "subject", trial = "trial", extendblinks = FALSE, type = "linear", hz = 250) ``` Now plot the interpolated signal over the raw trace: ```{r, warning=FALSE, message=FALSE, fig.height=5} pup_noise_interp <- ggplot(blinks_noise, aes(x = time, y = pupil)) + geom_rect(data = blink_windows, aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf), inherit.aes = FALSE, fill = "darkorange", alpha = 0.12) + geom_line(colour = "grey60") + geom_line(aes(y = interp), colour = "darkorange", linewidth = 0.7) + xlab("Time (ms)") + ylab("Pupil size (arbitrary units)") + ggtitle("Speed-filtered + noise-based detection + interpolation") + theme_bw(base_size = 13) print(pup_noise_interp) ``` ## Velocity vs. noise side by side ```{r, warning=FALSE, message=FALSE, fig.height=7, fig.width=10} pup_velocity_panel <- pup_g3 + ggtitle("Velocity-based") pup_noise_panel <- pup_noise_interp + ggtitle("Noise-based") plot_grid(pup_velocity_panel, pup_noise_panel, ncol = 1) ``` ## Tuning Two parameters are worth knowing about: - `smooth_ms` (default 10 ms) — the centered moving-average window applied before the algorithm searches for onset / offset. Wider windows are more robust to single-sample noise but can over-smooth short blinks. - `gap_ms` (default 100 ms) — adjacent blinks separated by less than this are merged into a single blink. Useful for double-blinks where the eyelid only briefly re-opens. If your tracker uses a different missing-data sentinel, pass it explicitly via `missing` (e.g. `missing = -1` for a tracker that flags missing samples with `-1`). ### Effect of `smooth_ms` on this trial The algorithm anchors on missing-data runs, so widening the smoother doesn't *find* additional blinks — it widens the window around each one. Below, the same trial is shown at three settings: ```{r, warning=FALSE, message=FALSE, fig.height=8, fig.width=10} sweep_one <- function(sm) { d <- interp_graph %>% dplyr::mutate(blink = detect_blinks_noise(pupil, hz = 250, smooth_ms = sm)) windows <- d %>% dplyr::filter(blink > 0) %>% dplyr::group_by(blink) %>% dplyr::summarise(start = min(time), end = max(time), .groups = "drop") list( data = d, windows = windows, label = sprintf("smooth_ms = %d (%d blinks, mean width = %.0f ms)", sm, nrow(windows), if (nrow(windows)) mean(windows$end - windows$start) else 0) ) } panels <- lapply(c(10, 20, 40), sweep_one) make_panel <- function(p) { ggplot(p$data, aes(x = time, y = pupil)) + geom_rect(data = p$windows, aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf), inherit.aes = FALSE, fill = "darkorange", alpha = 0.20) + geom_line(colour = "black") + geom_point(size = 0.5) + xlab("Time (ms)") + ylab("Pupil size") + ggtitle(p$label) + theme_bw(base_size = 12) } plot_grid(plotlist = lapply(panels, make_panel), ncol = 1) ``` At the default 10 ms the orange bands hug the missing-data runs tightly — you can see the eyelid descent and recovery just outside the shaded region. At 20 ms the windows widen slightly to absorb more of the descent. At 40 ms the windows roughly double in width, capturing the full eyelid trajectory before and after each missing run. Pushing the smoother to 80 ms (not shown) starts producing ~280 ms blink windows that bleed into surrounding pupil activity, so somewhere in the 20–40 ms range is the sweet spot for this recording. Pick the smallest value that produces clean interpolation downstream — wider isn't automatically better. ## When detection looks "off" The noise-based detector anchors on the tracker's missing-data flag, so a blink without a corresponding NA / 0 run in the raw signal will *not* be flagged. If you eyeball the trace and see what looks like a blink but the orange band doesn't cover it, check: 1. **Does the raw pupil column actually contain `NA`s or `0`s during that event?** If your tracker has already been "cleaned" (e.g. blinks filled in upstream), the noise-based method has nothing to anchor on — use `blink_velocity()` instead. 2. **Sentinel mismatch.** Some pipelines flag missing data with `-1` or `999`. Pass `missing = -1` (etc.) so the algorithm knows what to look for. 3. **Too-narrow `gap_ms`.** A double-blink can look like two separate detections; raise `gap_ms` to merge them. ## Comparison The two detectors make different tradeoffs: - `blink_velocity()` finds blink onsets *before* the signal goes missing, useful when the tracker is slow to report zeros / `NA`s. It needs a usable threshold (sensitive to recording gain). - `detect_blinks_noise()` anchors on the tracker's own missing-data flag, so it works without per-recording threshold tuning, but it can miss blinks that don't produce a missing-data run. In practice it's reasonable to try both on a representative trial and pick whichever recovers the cleaner signal for your data. # Worked example: Tobii Titta data + noise-based detection The walkthrough above used a 250 Hz EyeLink recording. To show the noise-based detector end-to-end on a different tracker, the package also ships one trimmed trial recorded with a Tobii spectrum at 120 Hz, captured via the Python `Titta` package. The gaze and message CSVs live in `inst/extdata/titta_gaze/` and `inst/extdata/titta_msg/` (paired by basename). `merge_tobii_titta()` reads the two folders, attaches each message to the nearest gaze sample, assigns trials from the `onset_` / `offset_` markers, and rebases time to be trial-relative in ms. ```{r, warning=FALSE, message=FALSE} gaze_dir <- system.file("extdata", "titta_gaze", package = "gazer") msg_dir <- system.file("extdata", "titta_msg", package = "gazer") titta_raw <- merge_tobii_titta(gaze = gaze_dir, msg = msg_dir) ``` Tobii reports left and right eyes separately. Average them with `compute_monocular_mean()` (which returns `NA` only when both eyes are missing), then run the result through `make_gazer()` so the standard `subject` / `trial` / `time` / `x` / `y` / `pupil` columns end up at the front of the data frame: ```{r, warning=FALSE, message=FALSE} titta <- titta_raw %>% mutate( pupil = compute_monocular_mean(left_pupil_diameter, right_pupil_diameter), x = compute_monocular_mean(left_gaze_point_on_display_area_x, right_gaze_point_on_display_area_x), y = compute_monocular_mean(left_gaze_point_on_display_area_y, right_gaze_point_on_display_area_y) ) %>% make_gazer(subject = "subject", trial = "trial", time = "time", x = "x", y = "y", pupil = "pupil") head(titta[, c("subject", "trial", "time", "x", "y", "pupil")]) ``` From here, the rest of the pipeline is identical to the EyeLink walkthrough above — only `hz = 120` differs. Apply the speed filter, run `detect_blinks_noise()`, then `interpolate_pupil()`: ```{r, warning=FALSE, message=FALSE} # 1. speed_pupil + MAD filter titta <- titta %>% dplyr::group_by(subject, trial) %>% dplyr::mutate(speed = speed_pupil(pupil, time)) %>% dplyr::ungroup() thr <- calc_mad(titta$speed[!is.na(titta$speed) & is.finite(titta$speed)], n = 16) titta <- titta %>% dplyr::mutate(pupil = ifelse(!is.na(speed) & is.finite(speed) & speed > thr, NA_real_, pupil)) # 2. detect blinks + 3. interpolate titta <- titta %>% dplyr::group_by(subject, trial) %>% dplyr::mutate(blink = as.integer(detect_blinks_noise(pupil, hz = 120) > 0)) %>% dplyr::ungroup() %>% interpolate_pupil(pupil = "pupil", subject = "subject", trial = "trial", extendblinks = FALSE, type = "linear", hz = 120) ``` ```{r, warning=FALSE, message=FALSE, fig.height=5} titta_windows <- titta %>% dplyr::filter(blink > 0) %>% dplyr::mutate(grp = cumsum(c(1, diff(time) > 1000 / 120 + 1))) %>% dplyr::group_by(grp) %>% dplyr::summarise(start = min(time), end = max(time), .groups = "drop") ggplot(titta, aes(x = time, y = pupil)) + geom_rect(data = titta_windows, aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf), inherit.aes = FALSE, fill = "darkorange", alpha = 0.12) + geom_line(colour = "grey60") + geom_line(aes(y = interp), colour = "darkorange", linewidth = 0.7) + xlab("Time (ms)") + ylab("Pupil size (mm)") + ggtitle("Titta trial: speed filter + noise-based detection + interpolation") + theme_bw(base_size = 13) ``` The same `detect_blinks_noise()` call worked on both trackers without changing the algorithm — only `hz` differed (250 for EyeLink, 120 for Tobii).