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.

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).

  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.

  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

## Warning: Removed 31 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

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.

 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))

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.

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.

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:

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.

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)

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)

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)

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)

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)

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)

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.

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:

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")])

## # A tibble: 6 × 6
##   subject   trial  time     x     y pupil
##   <chr>     <int> <dbl> <dbl> <dbl> <dbl>
## 1 DS21_2120     1  0    0.506 0.478  3.21
## 2 DS21_2120     1  8.33 0.507 0.475  3.22
## 3 DS21_2120     1 16.7  0.507 0.477  3.22
## 4 DS21_2120     1 25.0  0.509 0.477  3.22
## 5 DS21_2120     1 33.3  0.508 0.478  3.22
## 6 DS21_2120     1 41.7  0.506 0.477  3.23

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():

# 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)

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).