{ "nbformat": 4, "nbformat_minor": 5, "metadata": {}, "cells": [ { "metadata": {}, "source": [ "# Detecting Blinks from the Pupil Signal\n", "\n", "Many eye trackers (e.g., EyeLink) provide hardware-detected blink markers, but some\n", "do not. Even when they do, you may want an independent detection method for validation\n", "or for trackers that only provide a raw pupil signal.\n", "\n", "pymovements includes a `blink` detection function that identifies blinks directly from\n", "the pupil size time series. The algorithm adapts the differential detection approach of\n", "[Hershman et al. (2018)](https://doi.org/10.3758/s13428-017-1008-1)\n", "with a two-stage pipeline inspired by\n", "[PupilPre](https://CRAN.R-project.org/package=PupilPre)\n", "(Kyrolainen et al., 2019):\n", "\n", "1. **Flagging** -- samples where pupil is NaN/zero or shows rapid changes (large absolute difference)\n", "2. **Island absorption** -- short unflagged gaps between flagged regions are absorbed\n", "\n", "Duration defaults (50--500 ms) follow\n", "[Nystrom et al. (2024)](https://doi.org/10.3758/s13428-023-02333-9).\n", "\n", "This tutorial demonstrates:\n", "- Loading EyeLink data with `from_asc()`\n", "- Inspecting the recording structure (participant, session, duration)\n", "- Running algorithmic blink detection with `gaze.detect('blink')`\n", "- Listing **all** detected blink instances with timing and duration\n", "- Visualizing each blink in context\n", "- Comparing with EyeLink hardware blink events\n", "- Tuning detection parameters" ], "cell_type": "markdown", "id": "a0" }, { "metadata": {}, "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import polars as pl\n", "\n", "import pymovements as pm\n", "from pymovements.events import blink as blink_fn\n", "from pymovements.gaze.io import from_asc" ], "cell_type": "code", "outputs": [], "id": "a1", "execution_count": null }, { "metadata": {}, "source": [ "## 1. Load EyeLink Data\n", "\n", "We use the `ToyDatasetEyeLink` dataset. The file name encodes the\n", "**participant** and **session**: `subject_1_session_1.asc`.\n", "\n", "We load with `events=True` so that hardware blink events (`blink_eyelink`)\n", "are also available for comparison later." ], "cell_type": "markdown", "id": "b0" }, { "metadata": {}, "source": [ "# Download the dataset\n", "dataset = pm.Dataset('ToyDatasetEyeLink', path='data/ToyDataset')\n", "dataset.download()\n", "\n", "# Load the first ASC file\n", "raw_dir = dataset.paths.raw / 'pymovements-toy-dataset-eyelink-main'\n", "asc_file = raw_dir / 'raw' / 'subject_1_session_1.asc'\n", "\n", "gaze = from_asc(\n", " asc_file,\n", " patterns='eyelink',\n", " encoding='ascii',\n", " events=True,\n", ")\n", "\n", "print(f'File: {asc_file.name}')\n", "print('Subject: 1, Session: 1')\n", "print(f'Samples: {gaze.samples.shape}')\n", "print(f'Columns: {gaze.samples.columns}')\n", "gaze.samples.head()" ], "cell_type": "code", "outputs": [], "id": "b1", "execution_count": null }, { "metadata": {}, "source": [ "## 2. Understand the Recording\n", "\n", "This ASC file contains a **single continuous recording** of one participant in one session.\n", "\n", "The `time` column is in **milliseconds** (EyeLink native), but the values are large\n", "absolute timestamps (ms since the tracker started). We convert to **seconds relative\n", "to the recording start** for readable plots." ], "cell_type": "markdown", "id": "b2" }, { "metadata": {}, "source": [ "time_arr = gaze.samples['time'].to_numpy()\n", "pupil_arr = gaze.samples['pupil'].to_numpy()\n", "\n", "# Recording timing\n", "t0 = time_arr[0]\n", "duration_s = (time_arr[-1] - t0) / 1000\n", "\n", "print(f'EyeLink time range: {t0:.0f} \u2013 {time_arr[-1]:.0f} ms')\n", "print(f'Recording duration: {duration_s:.1f} s ({duration_s / 60:.1f} min)')\n", "print(f'Sampling rate: ~{len(time_arr) / duration_s:.0f} Hz')\n", "print(f'Total samples: {len(time_arr)}')" ], "cell_type": "code", "outputs": [], "id": "b3", "execution_count": null }, { "metadata": {}, "source": [ "## 3. Inspect the Raw Pupil Signal\n", "\n", "The pupil signal typically drops to zero or NaN during blinks. Let's plot the\n", "full recording using time relative to the start (in seconds) so the x-axis\n", "is easy to interpret." ], "cell_type": "markdown", "id": "c0" }, { "metadata": {}, "source": [ "# Relative time in seconds for all plots\n", "time_s = (time_arr - t0) / 1000\n", "\n", "fig, ax = plt.subplots(figsize=(14, 3))\n", "ax.plot(time_s, pupil_arr, color='mediumpurple', linewidth=0.5)\n", "ax.set_xlabel('Time since recording start (s)')\n", "ax.set_ylabel('Pupil Size')\n", "ax.set_title(\n", " f'Raw Pupil Signal \u2014 Subject 1, Session 1 '\n", " f'(full recording, {duration_s:.0f} s)',\n", ")\n", "ax.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", "plt.show()" ], "cell_type": "code", "outputs": [], "id": "c1", "execution_count": null }, { "metadata": {}, "source": [ "## 4. Detect Blinks\n", "\n", "We use `gaze.detect('blink')` which calls the blink detection algorithm. The function\n", "automatically reads the `pupil` column from the gaze dataframe.\n", "\n", "The detection algorithm adapts the differential approach of\n", "[Hershman et al. (2018)](https://doi.org/10.3758/s13428-017-1008-1)\n", "as implemented in [PupilPre](https://CRAN.R-project.org/package=PupilPre)\n", "(Kyrolainen et al., 2019).\n", "\n", "Default parameters:\n", "- `delta=None` -- auto-estimated from the 95th percentile of absolute pupil differences\n", "- `minimum_gap=3` -- short gaps up to 3 samples are absorbed\n", "- `minimum_candidates_around_gap=2` -- at least 2 flagged samples on each side to absorb a gap\n", "- `minimum_duration=50` -- events shorter than 50 ms are discarded\n", "- `maximum_duration=500` -- events longer than 500 ms are discarded\n", "\n", "The 50--500 ms duration defaults follow\n", "[Nystrom et al. (2024)](https://doi.org/10.3758/s13428-023-02333-9),\n", "who established that typical blinks fall within this range. Events outside\n", "this window are likely artifacts (too short) or extended eye closures (too long)." ], "cell_type": "markdown", "id": "d0" }, { "metadata": {}, "source": [ "gaze.detect('blink')\n", "\n", "detected_blinks = gaze.events.frame.filter(pl.col('name') == 'blink')\n", "print(f'Detected {len(detected_blinks)} blink events')" ], "cell_type": "code", "outputs": [], "id": "d1", "execution_count": null }, { "metadata": {}, "source": [ "### Complete list of all detected blinks\n", "\n", "Each row is one blink event. We show onset/offset (original EyeLink ms),\n", "duration, and time relative to recording start for easy reference." ], "cell_type": "markdown", "id": "d2" }, { "metadata": {}, "source": [ "blink_table = detected_blinks.select(\n", " pl.col('onset'),\n", " pl.col('offset'),\n", ").with_columns(\n", " (pl.col('offset') - pl.col('onset')).alias('duration_ms'),\n", " ((pl.col('onset') - t0) / 1000).round(2).alias('onset_s'),\n", ").with_row_index('blink_nr', offset=1)\n", "\n", "# Show ALL rows \u2014 no truncation\n", "with pl.Config(tbl_rows=-1):\n", " print(blink_table)" ], "cell_type": "code", "outputs": [], "id": "d3", "execution_count": null }, { "metadata": {}, "source": [ "## 5. Blink-by-Blink Visualization\n", "\n", "We plot each detected blink with ~500 ms of context on each side, so you\n", "can see the pupil drop in detail. Time on the x-axis is relative to the\n", "blink onset (in ms)." ], "cell_type": "markdown", "id": "e0" }, { "metadata": {}, "source": [ "blink_onsets = detected_blinks['onset'].to_list()\n", "blink_offsets = detected_blinks['offset'].to_list()\n", "\n", "n_blinks = len(blink_onsets)\n", "ncols = min(5, n_blinks)\n", "nrows = max(1, int(np.ceil(n_blinks / ncols)))\n", "context_ms = 500 # ms of context before/after blink\n", "\n", "fig, axes = plt.subplots(\n", " nrows, ncols, figsize=(ncols * 3, nrows * 2.2), squeeze=False,\n", ")\n", "\n", "for idx in range(n_blinks):\n", " r, c = divmod(idx, ncols)\n", " ax = axes[r, c]\n", "\n", " onset = blink_onsets[idx]\n", " offset = blink_offsets[idx]\n", " dur = offset - onset\n", " win_start = onset - context_ms\n", " win_end = offset + context_ms\n", "\n", " mask = (time_arr >= win_start) & (time_arr <= win_end)\n", " # Time relative to blink onset (ms)\n", " t_rel = time_arr[mask] - onset\n", "\n", " ax.plot(t_rel, pupil_arr[mask], color='mediumpurple', linewidth=0.8)\n", " ax.axvspan(0, dur, alpha=0.3, color='coral')\n", "\n", " ax.set_title(f'#{idx + 1} ({dur} ms)', fontsize=8)\n", " ax.tick_params(labelsize=6)\n", " ax.set_yticks([])\n", "\n", "# Hide unused subplots\n", "for idx in range(n_blinks, nrows * ncols):\n", " r, c = divmod(idx, ncols)\n", " axes[r, c].set_visible(False)\n", "\n", "fig.suptitle(\n", " f'All {n_blinks} Detected Blinks \u2014 Pupil Signal '\n", " f'(coral = blink, x = ms from onset)',\n", " fontsize=11, fontweight='bold',\n", ")\n", "plt.tight_layout()\n", "plt.show()" ], "cell_type": "code", "outputs": [], "id": "e1", "execution_count": null }, { "metadata": {}, "source": [ "## 6. Compare with EyeLink Hardware Blinks\n", "\n", "Since we loaded with `events=True`, we also have `blink_eyelink` events from the\n", "EyeLink hardware detector. Let's compare both." ], "cell_type": "markdown", "id": "f0" }, { "metadata": {}, "source": [ "hw_blinks = gaze.events.frame.filter(pl.col('name') == 'blink_eyelink')\n", "print(f'Hardware (EyeLink) blinks: {len(hw_blinks)}')\n", "print(f'Algorithmic blinks: {len(detected_blinks)}')\n", "\n", "hw_onsets = hw_blinks['onset'].to_list()\n", "hw_offsets = hw_blinks['offset'].to_list()\n", "\n", "if len(hw_onsets) > 0:\n", " # Window around the first few hardware blinks\n", " focus_start = hw_onsets[0] - 500\n", " focus_end = hw_offsets[min(2, len(hw_offsets) - 1)] + 500\n", " mask = (time_arr >= focus_start) & (time_arr <= focus_end)\n", " focus_t_s = (time_arr[mask] - focus_start) / 1000\n", "\n", " fig, ax = plt.subplots(figsize=(14, 3.5))\n", " ax.plot(focus_t_s, pupil_arr[mask], color='mediumpurple', linewidth=0.8)\n", "\n", " for i, (on, off) in enumerate(zip(hw_onsets, hw_offsets)):\n", " if off >= focus_start and on <= focus_end:\n", " ax.axvspan(\n", " (on - focus_start) / 1000, (off - focus_start) / 1000,\n", " alpha=0.2, color='steelblue',\n", " label='EyeLink' if i == 0 else None,\n", " )\n", "\n", " for i, (on, off) in enumerate(zip(blink_onsets, blink_offsets)):\n", " if off >= focus_start and on <= focus_end:\n", " ax.axvspan(\n", " (on - focus_start) / 1000, (off - focus_start) / 1000,\n", " alpha=0.2, color='coral',\n", " label='Algorithmic' if i == 0 else None,\n", " )\n", "\n", " ax.legend(loc='upper right')\n", " ax.set_xlabel('Time (s)')\n", " ax.set_ylabel('Pupil Size')\n", " ax.set_title('Hardware vs. Algorithmic Blink Detection')\n", " ax.grid(True, alpha=0.3)\n", " plt.tight_layout()\n", " plt.show()\n", "else:\n", " print('No hardware blink events found for comparison.')" ], "cell_type": "code", "outputs": [], "id": "f1", "execution_count": null }, { "metadata": {}, "source": [ "## 7. Parameter Tuning\n", "\n", "### What is `delta`?\n", "\n", "`delta` is a threshold on **sample-to-sample pupil change**: if `|pupil[i+1] - pupil[i]| > delta`,\n", "both samples get flagged as part of a blink.\n", "\n", "**Why it matters:** The algorithm first flags samples where pupil is 0 or NaN (the *middle* of a\n", "blink where the tracker lost the pupil entirely). But blinks also have *ramps* \u2014 the pupil signal\n", "drops rapidly as the eyelid closes and recovers rapidly as it opens. `delta` catches these\n", "transitional samples at the blink edges.\n", "\n", "**When `delta=None` (default):** It is auto-estimated as `5 x 95th percentile` of all valid\n", "`|diff(pupil)|`. This adapts to each recording \u2014 noisy data gets a higher threshold, clean data\n", "gets a lower one.\n", "\n", "**When to change it:**\n", "\n", "| Situation | What you see | What to do |\n", "|-----------|-------------|------------|\n", "| False positives (non-blink artifacts flagged) | Pupil constriction or noise marked as blinks | **Increase** delta (e.g., 200-500) |\n", "| Missed blinks (edges not captured) | Blink regions too narrow, missing the onset/offset ramps | **Decrease** delta (e.g., 20-50) |\n", "| Very noisy recording | Auto-delta is too high due to overall noise | Set delta **explicitly** to a fixed value |\n", "| Clean recording | Auto-delta works fine | Leave as `None` |" ], "cell_type": "markdown", "id": "g0" }, { "metadata": {}, "source": [ "deltas = [None, 50.0, 200.0, 500.0]\n", "\n", "# Focus window around first few blinks\n", "if len(blink_onsets) > 0:\n", " focus_start = blink_onsets[0] - 500\n", " focus_end = blink_offsets[min(2, len(blink_offsets) - 1)] + 500\n", "else:\n", " focus_start = t0\n", " focus_end = t0 + 5000\n", "\n", "mask = (time_arr >= focus_start) & (time_arr <= focus_end)\n", "focus_t_s = (time_arr[mask] - focus_start) / 1000\n", "\n", "fig, axes = plt.subplots(len(deltas), 1, figsize=(14, 2.5 * len(deltas)), sharex=True)\n", "\n", "for ax, d in zip(axes, deltas):\n", " events = blink_fn(\n", " pupil=pupil_arr,\n", " timesteps=time_arr.astype(int),\n", " delta=d,\n", " )\n", "\n", " ev_df = events.frame\n", " ax.plot(focus_t_s, pupil_arr[mask], color='mediumpurple', linewidth=0.8)\n", "\n", " if len(ev_df) > 0:\n", " for row in ev_df.to_dicts():\n", " on, off = row['onset'], row['offset']\n", " if off >= focus_start and on <= focus_end:\n", " ax.axvspan(\n", " (on - focus_start) / 1000, (off - focus_start) / 1000,\n", " alpha=0.3, color='coral',\n", " )\n", "\n", " label = f'delta={d}' if d is not None else 'delta=auto'\n", " ax.set_ylabel('Pupil', fontsize=9)\n", " ax.set_title(f'{label} \u2014 {len(ev_df)} events total', fontsize=10)\n", " ax.grid(True, alpha=0.3)\n", "\n", "axes[-1].set_xlabel('Time (s)')\n", "fig.suptitle('Effect of delta on Blink Detection', fontsize=13, fontweight='bold')\n", "plt.tight_layout()\n", "plt.show()" ], "cell_type": "code", "outputs": [], "id": "g1", "execution_count": null }, { "metadata": {}, "source": [ "### What is `minimum_gap`?\n", "\n", "`minimum_gap` controls **island absorption** \u2014 merging short unflagged gaps between flagged\n", "regions. During a blink, the tracker sometimes recovers valid pupil readings for 1-3 samples\n", "in the middle of an otherwise complete blink. Without absorption, this splits one real blink\n", "into two or more events.\n", "\n", "- `minimum_gap=0` \u2014 disables absorption entirely (no merging)\n", "- `minimum_gap=3` (default) \u2014 gaps of up to 3 unflagged samples are absorbed if surrounded\n", " by at least `minimum_candidates_around_gap=2` flagged samples on each side\n", "- `minimum_gap=10` \u2014 aggressive merging, useful for very noisy data\n", "\n", "**When to change it:** If you see single blinks split into multiple events, increase it.\n", "If distinct close blinks are being incorrectly merged, decrease it." ], "cell_type": "markdown", "id": "h0" }, { "metadata": {}, "source": [ "runs = [0, 1, 3, 10]\n", "\n", "fig, axes = plt.subplots(len(runs), 1, figsize=(14, 2.5 * len(runs)), sharex=True)\n", "\n", "for ax, mvr in zip(axes, runs):\n", " events = blink_fn(\n", " pupil=pupil_arr,\n", " timesteps=time_arr.astype(int),\n", " minimum_gap=mvr,\n", " )\n", "\n", " ev_df = events.frame\n", " ax.plot(focus_t_s, pupil_arr[mask], color='mediumpurple', linewidth=0.8)\n", "\n", " if len(ev_df) > 0:\n", " for row in ev_df.to_dicts():\n", " on, off = row['onset'], row['offset']\n", " if off >= focus_start and on <= focus_end:\n", " ax.axvspan(\n", " (on - focus_start) / 1000, (off - focus_start) / 1000,\n", " alpha=0.3, color='coral',\n", " )\n", "\n", " ax.set_ylabel('Pupil', fontsize=9)\n", " ax.set_title(f'minimum_gap={mvr} \u2014 {len(ev_df)} events total', fontsize=10)\n", " ax.grid(True, alpha=0.3)\n", "\n", "axes[-1].set_xlabel('Time (s)')\n", "fig.suptitle(\n", " 'Effect of minimum_gap on Blink Detection', fontsize=13, fontweight='bold',\n", ")\n", "plt.tight_layout()\n", "plt.show()" ], "cell_type": "code", "outputs": [], "id": "h1", "execution_count": null }, { "metadata": {}, "source": [ "### Duration filtering: `minimum_duration` and `maximum_duration`\n", "\n", "Following [Nystrom et al. (2024)](https://doi.org/10.3758/s13428-023-02333-9),\n", "typical blinks last **50--500 ms**. The defaults `minimum_duration=50` and\n", "`maximum_duration=500` filter out events outside this range:\n", "\n", "- Events **< 50 ms** are likely partial tracking losses or noise, not true blinks\n", "- Events **> 500 ms** are likely extended eye closures or prolonged tracking loss\n", "\n", "Set `maximum_duration=None` to disable the upper bound (e.g., if you want to\n", "capture extended eye closures too)." ], "cell_type": "markdown", "id": "i0" }, { "metadata": {}, "source": [ "# Compare: default duration filter vs. no upper bound\n", "events_default = blink_fn(\n", " pupil=pupil_arr,\n", " timesteps=time_arr.astype(int),\n", ")\n", "events_no_max = blink_fn(\n", " pupil=pupil_arr,\n", " timesteps=time_arr.astype(int),\n", " maximum_duration=None,\n", ")\n", "events_strict = blink_fn(\n", " pupil=pupil_arr,\n", " timesteps=time_arr.astype(int),\n", " minimum_duration=80,\n", " maximum_duration=300,\n", ")\n", "\n", "print(f'Default (50-500 ms): {len(events_default.frame)} blinks')\n", "print(f'No upper bound (50+ ms): {len(events_no_max.frame)} blinks')\n", "print(f'Strict (80-300 ms): {len(events_strict.frame)} blinks')\n", "\n", "# Show the duration distribution\n", "if len(events_no_max.frame) > 0:\n", " durations = events_no_max.frame['duration'].to_numpy()\n", " fig, ax = plt.subplots(figsize=(10, 3))\n", " ax.hist(durations, bins=30, color='mediumpurple', edgecolor='white', alpha=0.8)\n", " ax.axvline(50, color='coral', linestyle='--', linewidth=2, label='min=50 ms')\n", " ax.axvline(500, color='steelblue', linestyle='--', linewidth=2, label='max=500 ms')\n", " ax.set_xlabel('Blink Duration (ms)')\n", " ax.set_ylabel('Count')\n", " ax.set_title('Blink Duration Distribution (no upper bound)')\n", " ax.legend()\n", " ax.grid(True, alpha=0.3)\n", " plt.tight_layout()\n", " plt.show()" ], "cell_type": "code", "outputs": [], "id": "wo8lim1otpr", "execution_count": null }, { "metadata": {}, "source": [ "## Summary\n", "\n", "- `blink()` detects blinks from the pupil signal without relying on eye-tracker markers\n", "- The algorithm adapts [Hershman et al. (2018)](https://doi.org/10.3758/s13428-017-1008-1) as implemented in [PupilPre](https://CRAN.R-project.org/package=PupilPre) (Kyrolainen et al., 2019)\n", "- Duration defaults follow [Nystrom et al. (2024)](https://doi.org/10.3758/s13428-023-02333-9): **50--500 ms**\n", "- Key parameters: `delta` (sensitivity), `minimum_gap` (gap merging), `minimum_duration` / `maximum_duration` (duration filtering)\n", "- Use `gaze.detect('blink')` for seamless integration with the pymovements pipeline\n", "- For cleaning blink artifacts from gaze data, see the [Blink Cleaning](blink-cleaning.ipynb) tutorial" ], "cell_type": "markdown", "id": "cneiff3c14" } ] }