Move the gaze calibration code into gazecal.py

The correction models, filters, blink filter, and SteamVR log reader move
out of the probe so the gaze service can use them too. Two additions on the
way: the log reader follows the headset going on and off (SteamVR starts its
eye model over each time it goes on), and LiveCorrection counts lessons from
before the last time it went on less (OLD_WEAR), so the first few after
relearn the offset while the shape is kept.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
DeeJanuzandClaude Opus 5.5 committed 2026-09-27 22:51:24 -06:00
1 parent 526dfa3511
commit 1d0211dc5e
2 files changed
+524 -477

No files matched your search

+520
View File
@@ -0,0 +1,520 @@
"""gazecal: gaze calibration shared by ft-gazeprobe and ft-gazed.
The correction models (Correction: the calibration fitted from calibration dots;
LiveCorrection: what clicks teach on the fly, on top of it), the smoothing filters, the
blink and dropout filter for one look at a spot, and SteamEyeLog, which follows SteamVR's
eye tracking log. Angles are head-relative degrees (yaw +left, pitch +up), as ft-gaze
reports them.
"""
import math
import os
import statistics
import time
from pathlib import Path
STATE = Path.home() / ".local" / "state" / "frametop" / "gaze"
# --- Small math ---------------------------------------------------------------------
def px_from_deg(j, dy, dp):
"""Pixels for a head-relative change of (yaw, pitch) degrees, from ft-gaze's Jacobian."""
return j[0] * dy + j[2] * dp, j[1] * dy + j[3] * dp
def deg_from_px(j, dx, dy):
"""Head-relative (yaw, pitch) degrees for a pixel offset: the Jacobian's inverse."""
det = j[0] * j[3] - j[2] * j[1]
if abs(det) < 1e-9:
return 0.0, 0.0
return (j[3] * dx - j[2] * dy) / det, (-j[1] * dx + j[0] * dy) / det
class OneEuro:
"""One Euro filter (Casiez et al. 2012): smooth when still, quick when moving.
`scale` turns the input's units into degrees, so beta is per degree a second."""
def __init__(self, min_cutoff=1.0, beta=0.01, d_cutoff=1.0):
self.min_cutoff, self.beta, self.d_cutoff = min_cutoff, beta, d_cutoff
self.x = self.dx = self.t = None
@staticmethod
def alpha(cutoff, dt):
tau = 1.0 / (2 * math.pi * cutoff)
return 1.0 / (1.0 + tau / dt)
def __call__(self, x, t, scale=1.0):
if self.t is None or t <= self.t or t - self.t > 0.5:
self.x, self.dx, self.t = x, 0.0, t
return x
dt = t - self.t
dx = (x - self.x) / dt
a_d = self.alpha(self.d_cutoff, dt)
self.dx = a_d * dx + (1 - a_d) * self.dx
cutoff = self.min_cutoff + self.beta * abs(self.dx) * scale
a = self.alpha(cutoff, dt)
self.x = a * x + (1 - a) * self.x
self.t = t
return self.x
class Fixation:
"""Dispersion-based fixations: while gaze stays within `radius` degrees of the current
fixation's mean, the output is that mean, so the dot sits still; two samples in a row
outside it start a new fixation there, so a glance elsewhere moves the dot at once."""
def __init__(self, radius=1.0):
self.radius = radius
self.reset()
def reset(self):
self.sum = [0.0, 0.0]
self.count = 0
self.outside = []
self.last_t = None
def __call__(self, x, y, t, dpp):
if self.last_t is not None and (t <= self.last_t or t - self.last_t > 0.5):
self.reset()
self.last_t = t
if self.count:
mx, my = self.sum[0] / self.count, self.sum[1] / self.count
if math.hypot(x - mx, y - my) * dpp > self.radius:
self.outside.append((x, y))
if len(self.outside) < 2:
return mx, my # one stray sample: probably noise
self.sum = [sum(p[0] for p in self.outside), sum(p[1] for p in self.outside)]
self.count = len(self.outside)
self.outside = []
return self.sum[0] / self.count, self.sum[1] / self.count
self.outside = []
if self.count >= 90: # the last second or so: a slow drift still gets followed
self.sum = [self.sum[0] * 89 / 90, self.sum[1] * 89 / 90]
self.count = 89
self.sum[0] += x
self.sum[1] += y
self.count += 1
return self.sum[0] / self.count, self.sum[1] / self.count
MODELS = ["none", "offset", "affine", "affine+grid", "quadratic", "quadratic+grid"]
DEFAULT_MODEL = "quadratic"
class Correction:
"""Gaze correction in degrees, looked up by where in your view you're looking
(head-relative yaw and pitch, hy and hp):
offset one (yaw, pitch) offset everywhere
affine plus a straight-line change across the view: a gain and a tilt
quadratic plus curvature (hy*hp, hy^2, hp^2): the second-order polynomial video
eye trackers usually calibrate with. The tracker's error grows as
the eye turns away from the centre (on the Frame it overstates
vertical movement, more the further up or down you look, and looking
up adds a sideways error), and a straight line can only follow part
of that
...+grid plus a bilinear grid of what's left every 10 degrees
Coefficients: C @ f, f = [1, x, y, x*y, x^2, y^2] with x = hy/30, y = hp/30; the
terms a model doesn't use are 0. A polynomial runs away outside the spots it was fitted
on, so its input is clamped to the range of view it has seen, plus a margin.
"""
YAWS = list(range(-40, 41, 10))
PITCHES = list(range(-30, 31, 10))
NF = 6
MARGIN = 3.0 # degrees past the fitted range that the polynomial still follows
def __init__(self):
self.reset()
def reset(self):
self.a = [[0.0] * self.NF, [0.0] * self.NF]
self.grid = [[[0.0, 0.0] for _ in self.PITCHES] for _ in self.YAWS]
self.samples = 0
self.range = None # [hy min, hy max, hp min, hp max] of the samples so far
@staticmethod
def base(mode):
return mode.split("+")[0]
def clamp(self, hy, hp):
if not self.range:
return hy, hp
y0, y1, p0, p1 = self.range
m = self.MARGIN
return min(max(hy, y0 - m), y1 + m), min(max(hp, p0 - m), p1 + m)
def features(self, hy, hp, mode):
kind = self.base(mode)
if kind == "offset":
return [1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
hy, hp = self.clamp(hy, hp)
x, y = hy / 30.0, hp / 30.0
if kind == "affine":
return [1.0, x, y, 0.0, 0.0, 0.0]
return [1.0, x, y, x * y, x * x, y * y]
def extend(self, hy, hp):
if self.range is None:
self.range = [hy, hy, hp, hp]
else:
r = self.range
self.range = [min(r[0], hy), max(r[1], hy), min(r[2], hp), max(r[3], hp)]
def weights(self, hy, hp):
def cell(v, axis):
v = min(max(v, axis[0]), axis[-1])
i = min(int((v - axis[0]) // 10), len(axis) - 2)
return i, (v - axis[i]) / 10.0
i, fy = cell(hy, self.YAWS)
k, fp = cell(hp, self.PITCHES)
return [((i, k), (1 - fy) * (1 - fp)), ((i + 1, k), fy * (1 - fp)),
((i, k + 1), (1 - fy) * fp), ((i + 1, k + 1), fy * fp)]
def get(self, hy, hp, mode):
if mode == "none":
return 0.0, 0.0
f = self.features(hy, hp, mode)
cy = sum(a * b for a, b in zip(self.a[0], f))
cp = sum(a * b for a, b in zip(self.a[1], f))
if mode.endswith("+grid"):
for (i, k), w in self.weights(hy, hp):
cy += w * self.grid[i][k][0]
cp += w * self.grid[i][k][1]
return cy, cp
def learn(self, hy, hp, dy, dp, mode, rate):
"""One sample: the correction here should have been (dy, dp) degrees more.
Normalized LMS for the polynomial; the grid takes half when it's on."""
if mode == "none":
return
self.samples += 1
self.extend(hy, hp)
grid = mode.endswith("+grid")
share = rate * 0.5 if grid else rate
f = self.features(hy, hp, mode)
norm = sum(v * v for v in f)
for row, d in ((self.a[0], dy), (self.a[1], dp)):
for n in range(self.NF):
row[n] += share * d * f[n] / norm
if grid:
rest = rate - share
for (i, k), w in self.weights(hy, hp):
self.grid[i][k][0] += rest * w * dy
self.grid[i][k][1] += rest * w * dp
def fit(self, points, mode, ridge=0.05, smooth=0.3):
"""Batch fit from (hy, hp, dy, dp) points, each the whole error there (degrees)."""
self.reset()
if mode == "none" or not points:
return
self.samples = len(points)
for p in points:
self.extend(p[0], p[1])
kind = self.base(mode)
used = {"offset": 1, "affine": 3, "quadratic": 6}[kind]
if used > 1 and len(points) < used + 2: # too few spots for this many terms
kind, used = ("affine", 3) if len(points) >= 5 else ("offset", 1)
if kind == "offset":
self.a[0][0] = statistics.fmean(p[2] for p in points)
self.a[1][0] = statistics.fmean(p[3] for p in points)
else:
# Least squares, with a little ridge on everything but the offset, so a
# lopsided set of spots can't bend it far.
X = [self.features(p[0], p[1], kind)[:used] for p in points]
M = [[sum(x[r] * x[c] for x in X) + (ridge * len(X) if r == c and r else 0.0) for c in range(used)]
for r in range(used)]
for out, col in ((self.a[0], 2), (self.a[1], 3)):
b = [sum(x[r] * p[col] for x, p in zip(X, points)) for r in range(used)]
out[:used] = solve(M, b)
if not mode.endswith("+grid"):
return
# Each node: the weighted mean of what the polynomial left over near it, shrunk toward 0.
acc = [[[0.0, 0.0, 0.0] for _ in self.PITCHES] for _ in self.YAWS]
for hy, hp, dy, dp in points:
ly, lp = self.get(hy, hp, kind)
for (i, k), w in self.weights(hy, hp):
acc[i][k][0] += w * (dy - ly)
acc[i][k][1] += w * (dp - lp)
acc[i][k][2] += w
for i in range(len(self.YAWS)):
for k in range(len(self.PITCHES)):
sy, sp, sw = acc[i][k]
self.grid[i][k] = [sy / (sw + smooth), sp / (sw + smooth)]
def offset(self):
return self.a[0][0], self.a[1][0]
def to_json(self):
return {"coef": self.a, "range": self.range, "grid": self.grid, "samples": self.samples}
def from_json(self, d):
self.reset()
a = d.get("coef") or d.get("affine") # "affine": the 3-term version of this file
if a and len(a) == 2 and all(len(r) in (3, self.NF) for r in a):
self.a = [list(map(float, r)) + [0.0] * (self.NF - len(r)) for r in a]
grid = d.get("grid")
if grid and len(grid) == len(self.YAWS) and all(len(r) == len(self.PITCHES) for r in grid):
self.grid = [[list(c) for c in row] for row in grid]
r = d.get("range")
self.range = list(map(float, r)) if r and len(r) == 4 else None
self.samples = d.get("samples", 0)
def solve(M, b):
"""Solve a small linear system (Gaussian elimination with pivoting); zeros if singular."""
n = len(b)
A = [row[:] + [b[i]] for i, row in enumerate(M)]
for c in range(n):
piv = max(range(c, n), key=lambda r: abs(A[r][c]))
if abs(A[piv][c]) < 1e-12:
return [0.0] * n
A[c], A[piv] = A[piv], A[c]
for r in range(n):
if r != c:
f = A[r][c] / A[c][c]
for k in range(c, n + 1):
A[r][k] -= f * A[c][k]
return [A[i][n] / A[i][i] for i in range(n)]
class LiveCorrection:
"""Corrections learned on the fly from snapped clicks, on top of the calibration.
Each click on an element is a measurement: you were looking at that element when you
pressed, and the tracker put your gaze at the raw point, so the gap between them is the
whole error there. Whatever the calibration doesn't already explain (the residual) is
fitted with the same quadratic terms. Every term but the offset is held close to zero
(ridge), so one click shifts the whole correction and more clicks bend it. Whatever is
still left near a click is added within a few degrees of it: on the Frame, errors less
than 3 degrees apart are alike, and ones further apart are unrelated. Recent clicks count more (a
half-life counted in clicks), so it follows SteamVR's gaze as that drifts or relearns.
Replayed on logged points: a calibration from an earlier session was 4.95 degrees off;
one click brought that to 2.3, five to 1.6, twenty to 1.2.
A big element says little about where on it you looked, so each axis is weighted by the
element's size along it: a list row 12 degrees wide barely counts sideways.
Putting the headset back on moves the error (SteamVR starts its eye model over each
time, and the headset sits a little differently), so clicks from before the last time
it went on (`wear_time`, when set) count OLD_WEAR as much: the offset is relearned from
the first few clicks after, and the shape is kept meanwhile."""
RIDGE = [0.01, 0.5, 0.5, 0.5, 0.5, 0.5]
KERNEL = 2.5 # degrees: how far a click's leftover reaches
SHRINK = 0.5 # near one click, half its leftover; near several, nearly all
HALF_LIFE = 40 # clicks
KEEP = 150
MARGIN = 3.0
OLD_WEAR = 0.3
def __init__(self):
self.samples = [] # dicts: time, hy, hp, dy, dp (the whole error), wy, wp
self.wear_time = None
self.reset_fit()
def reset_fit(self):
self.cy = [0.0] * 6
self.cp = [0.0] * 6
self.left = [] # (hy, hp, leftover yaw, leftover pitch, wy, wp, decay)
self.range = None
def features(self, hy, hp):
if self.range:
y0, y1, p0, p1 = self.range
hy = min(max(hy, y0 - self.MARGIN), y1 + self.MARGIN)
hp = min(max(hp, p0 - self.MARGIN), p1 + self.MARGIN)
x, y = hy / 30.0, hp / 30.0
return [1.0, x, y, x * y, x * x, y * y]
def add(self, sample, base, mode):
self.samples = (self.samples + [sample])[-self.KEEP:]
self.refit(base, mode)
def undo(self, base, mode):
if self.samples:
self.samples.pop()
self.refit(base, mode)
def refit(self, base, mode):
"""Refit from the samples against the calibration as it is now."""
self.reset_fit()
n = len(self.samples)
if not n:
return
self.range = [min(s["hy"] for s in self.samples), max(s["hy"] for s in self.samples),
min(s["hp"] for s in self.samples), max(s["hp"] for s in self.samples)]
rows = []
for k, s in enumerate(self.samples):
decay = 0.5 ** ((n - 1 - k) / self.HALF_LIFE)
if self.wear_time and s.get("time", 0) < self.wear_time:
decay *= self.OLD_WEAR
by, bp = base.get(s["hy"], s["hp"], mode)
rows.append((s, self.features(s["hy"], s["hp"]), s["dy"] - by, s["dp"] - bp, decay))
for out, ri, wi in ((self.cy, 2, "wy"), (self.cp, 3, "wp")):
M = [[self.RIDGE[r] if r == c else 0.0 for c in range(6)] for r in range(6)]
b = [0.0] * 6
for row in rows:
w = row[4] * row[0][wi]
f = row[1]
for r in range(6):
b[r] += w * f[r] * row[ri]
for c in range(6):
M[r][c] += w * f[r] * f[c]
out[:] = solve(M, b)
for s, f, ry, rp, decay in rows:
ly = ry - sum(a * v for a, v in zip(self.cy, f))
lp = rp - sum(a * v for a, v in zip(self.cp, f))
self.left.append((s["hy"], s["hp"], ly, lp, s["wy"] * decay, s["wp"] * decay))
def get(self, hy, hp):
if not self.samples:
return 0.0, 0.0
f = self.features(hy, hp)
cy = sum(a * v for a, v in zip(self.cy, f))
cp = sum(a * v for a, v in zip(self.cp, f))
k2 = 2 * self.KERNEL ** 2
sy = sp = wy = wp = 0.0
for y, p, ly, lp, ay, ap in self.left:
d2 = (y - hy) ** 2 + (p - hp) ** 2
if d2 > 9 * k2:
continue
g = math.exp(-d2 / k2)
sy += g * ay * ly
wy += g * ay
sp += g * ap * lp
wp += g * ap
return cy + sy / (wy + self.SHRINK), cp + sp / (wp + self.SHRINK)
def offset(self):
return self.cy[0], self.cp[0]
class SteamEyeLog:
"""Follows SteamVR's eye tracking log (read only) for what moves the raw gaze under a
calibration.
SteamVR's eye tracker calibrates itself from clicks: a quick mouse-button down and up
(the laser or the Frametop pointer), with the gaze within 5 degrees of the click and
held still, is taken as "you were looking there" ("Accept usercal"). Accepted clicks
were all under 0.14 s; 0.38 s was "too slow", and one that moved was refused. It keeps that in
the running `eyetracking` process and saves nothing, so when the process starts again
(SteamVR restarting), its calibration starts over. Both events are counted here, and
each time the headset goes on ("HMD on"): the eye model starts over then too."""
PATH = Path.home() / ".local" / "share" / "Steam" / "logs" / "eyetracking.txt"
def __init__(self):
self.pos = 0
self.inode = None
self.partial = ""
self.starts = [] # when the eyetracking process started
self.accepts = [] # when it learned from a click
self.rejects = []
self.wears = [] # when the headset went on
self.offs = [] # ... and off
@staticmethod
def stamp(line):
head = line.split(" [", 1)[0]
main, _, frac = head.partition(".")
try:
return time.mktime(time.strptime(main.strip(), "%a %b %d %Y %H:%M:%S")) + float("0." + (frac or "0"))
except ValueError:
return None
def poll(self):
"""Read what's new. True if the eye tracker started again since the last poll."""
try:
st = os.stat(self.PATH)
except OSError:
return False
if st.st_ino != self.inode or st.st_size < self.pos:
self.inode, self.pos, self.partial = st.st_ino, 0, ""
if st.st_size == self.pos:
return False
first = self.pos == 0 and not self.starts
try:
with open(self.PATH, "rb") as f:
f.seek(self.pos)
data = f.read()
except OSError:
return False
self.pos += len(data)
lines = (self.partial + data.decode("utf-8", "replace")).split("\n")
self.partial = lines.pop()
restarted = False
for line in lines:
if "usercal" not in line and "startup with PID" not in line and "HMD o" not in line:
continue
t = self.stamp(line)
if t is None:
continue
if "startup with PID" in line:
self.starts.append(t)
restarted = not first
elif "HMD on" in line:
self.wears.append(t)
elif "HMD off" in line:
self.offs.append(t)
elif "Accept usercal" in line:
self.accepts.append(t)
elif "Reject usercal" in line:
self.rejects.append(t)
return restarted
def worn(self):
"""When the headset last went on (None if not in this log)."""
return self.wears[-1] if self.wears else None
def wearing(self):
"""Whether the headset is on, as far as the log says (None: it doesn't say)."""
if not self.wears and not self.offs:
return None
return bool(self.wears) and (not self.offs or self.wears[-1] > self.offs[-1])
def started(self):
return self.starts[-1] if self.starts else None
def accepted_since(self, t):
return sum(1 for a in self.accepts if a >= t)
def cross_validate(points, mode):
"""Leave-one-out: each point's error under a model fitted on all the others (degrees)."""
errs = []
for i in range(len(points)):
c = Correction()
c.fit(points[:i] + points[i + 1:], mode)
cy, cp = c.get(points[i][0], points[i][1], mode)
errs.append(math.hypot(points[i][2] - cy, points[i][3] - cp))
return errs
def steady_samples(samples, vergence_jump=1.5):
"""The samples of one look at one spot where the tracker had both eyes: none in a blink
(openness under half its median over the samples), and none where the angle between the eyes' directions (`lr`, the
vergence) is more than `vergence_jump` degrees from its median over the samples. The
vergence itself depends on distance (about 2.8 degrees for a screen 1.3 m away, a
fraction of one far off), so only a jump away from what it was during this look means
the tracker lost an eye. Without the mmap there's nothing to judge by: all are kept."""
# Openness: a blink is a sharp drop from what it was during this look. Not a fixed
# level: looking down, the upper lids come down with the eyes, and in bright light you
# squint, so the reading can stay under 0.5 for the whole look while the tracker follows
# the eyes fine (a calibration dot at the bottom of the bright round failed that way).
opens = [min(o) for o in ((smp["src"].get("mmap1") or {}).get("open") for smp in samples) if o]
floor = max(0.12, 0.5 * statistics.median(opens)) if len(opens) >= 5 else 0.12
opened = []
for smp in samples:
o = (smp["src"].get("mmap1") or {}).get("open")
if not o or min(o) >= floor:
opened.append(smp)
def vergence(smp):
return (smp["src"].get("mmap1") or {}).get("lr", (smp["src"].get("mmap2") or {}).get("lr"))
have = [v for v in map(vergence, opened) if v is not None]
if len(have) < 5:
return opened
med = statistics.median(have)
return [smp for smp in opened if vergence(smp) is None or abs(vergence(smp) - med) <= vergence_jump]
+4 -477
View File
@@ -78,444 +78,10 @@ SOURCE_COLORS = {"action": (0.2, 0.8, 1.0), "mmap1": (1.0, 0.6, 0.1), "mmap2": (
TRIGGER_KEYS = {Gdk.KEY_Return, Gdk.KEY_KP_Enter, Gdk.KEY_space}
# --- Small math ---------------------------------------------------------------------
def px_from_deg(j, dy, dp):
"""Pixels for a head-relative change of (yaw, pitch) degrees, from ft-gaze's Jacobian."""
return j[0] * dy + j[2] * dp, j[1] * dy + j[3] * dp
def deg_from_px(j, dx, dy):
"""Head-relative (yaw, pitch) degrees for a pixel offset: the Jacobian's inverse."""
det = j[0] * j[3] - j[2] * j[1]
if abs(det) < 1e-9:
return 0.0, 0.0
return (j[3] * dx - j[2] * dy) / det, (-j[1] * dx + j[0] * dy) / det
class OneEuro:
"""One Euro filter (Casiez et al. 2012): smooth when still, quick when moving.
`scale` turns the input's units into degrees, so beta is per degree a second."""
def __init__(self, min_cutoff=1.0, beta=0.01, d_cutoff=1.0):
self.min_cutoff, self.beta, self.d_cutoff = min_cutoff, beta, d_cutoff
self.x = self.dx = self.t = None
@staticmethod
def alpha(cutoff, dt):
tau = 1.0 / (2 * math.pi * cutoff)
return 1.0 / (1.0 + tau / dt)
def __call__(self, x, t, scale=1.0):
if self.t is None or t <= self.t or t - self.t > 0.5:
self.x, self.dx, self.t = x, 0.0, t
return x
dt = t - self.t
dx = (x - self.x) / dt
a_d = self.alpha(self.d_cutoff, dt)
self.dx = a_d * dx + (1 - a_d) * self.dx
cutoff = self.min_cutoff + self.beta * abs(self.dx) * scale
a = self.alpha(cutoff, dt)
self.x = a * x + (1 - a) * self.x
self.t = t
return self.x
class Fixation:
"""Dispersion-based fixations: while gaze stays within `radius` degrees of the current
fixation's mean, the output is that mean, so the dot sits still; two samples in a row
outside it start a new fixation there, so a glance elsewhere moves the dot at once."""
def __init__(self, radius=1.0):
self.radius = radius
self.reset()
def reset(self):
self.sum = [0.0, 0.0]
self.count = 0
self.outside = []
self.last_t = None
def __call__(self, x, y, t, dpp):
if self.last_t is not None and (t <= self.last_t or t - self.last_t > 0.5):
self.reset()
self.last_t = t
if self.count:
mx, my = self.sum[0] / self.count, self.sum[1] / self.count
if math.hypot(x - mx, y - my) * dpp > self.radius:
self.outside.append((x, y))
if len(self.outside) < 2:
return mx, my # one stray sample: probably noise
self.sum = [sum(p[0] for p in self.outside), sum(p[1] for p in self.outside)]
self.count = len(self.outside)
self.outside = []
return self.sum[0] / self.count, self.sum[1] / self.count
self.outside = []
if self.count >= 90: # the last second or so: a slow drift still gets followed
self.sum = [self.sum[0] * 89 / 90, self.sum[1] * 89 / 90]
self.count = 89
self.sum[0] += x
self.sum[1] += y
self.count += 1
return self.sum[0] / self.count, self.sum[1] / self.count
MODELS = ["none", "offset", "affine", "affine+grid", "quadratic", "quadratic+grid"]
DEFAULT_MODEL = "quadratic"
class Correction:
"""Gaze correction in degrees, looked up by where in your view you're looking
(head-relative yaw and pitch, hy and hp):
offset one (yaw, pitch) offset everywhere
affine plus a straight-line change across the view: a gain and a tilt
quadratic plus curvature (hy*hp, hy^2, hp^2): the second-order polynomial video
eye trackers usually calibrate with. The tracker's error grows as
the eye turns away from the centre (on the Frame it overstates
vertical movement, more the further up or down you look, and looking
up adds a sideways error), and a straight line can only follow part
of that
...+grid plus a bilinear grid of what's left every 10 degrees
Coefficients: C @ f, f = [1, x, y, x*y, x^2, y^2] with x = hy/30, y = hp/30; the
terms a model doesn't use are 0. A polynomial runs away outside the spots it was fitted
on, so its input is clamped to the range of view it has seen, plus a margin.
"""
YAWS = list(range(-40, 41, 10))
PITCHES = list(range(-30, 31, 10))
NF = 6
MARGIN = 3.0 # degrees past the fitted range that the polynomial still follows
def __init__(self):
self.reset()
def reset(self):
self.a = [[0.0] * self.NF, [0.0] * self.NF]
self.grid = [[[0.0, 0.0] for _ in self.PITCHES] for _ in self.YAWS]
self.samples = 0
self.range = None # [hy min, hy max, hp min, hp max] of the samples so far
@staticmethod
def base(mode):
return mode.split("+")[0]
def clamp(self, hy, hp):
if not self.range:
return hy, hp
y0, y1, p0, p1 = self.range
m = self.MARGIN
return min(max(hy, y0 - m), y1 + m), min(max(hp, p0 - m), p1 + m)
def features(self, hy, hp, mode):
kind = self.base(mode)
if kind == "offset":
return [1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
hy, hp = self.clamp(hy, hp)
x, y = hy / 30.0, hp / 30.0
if kind == "affine":
return [1.0, x, y, 0.0, 0.0, 0.0]
return [1.0, x, y, x * y, x * x, y * y]
def extend(self, hy, hp):
if self.range is None:
self.range = [hy, hy, hp, hp]
else:
r = self.range
self.range = [min(r[0], hy), max(r[1], hy), min(r[2], hp), max(r[3], hp)]
def weights(self, hy, hp):
def cell(v, axis):
v = min(max(v, axis[0]), axis[-1])
i = min(int((v - axis[0]) // 10), len(axis) - 2)
return i, (v - axis[i]) / 10.0
i, fy = cell(hy, self.YAWS)
k, fp = cell(hp, self.PITCHES)
return [((i, k), (1 - fy) * (1 - fp)), ((i + 1, k), fy * (1 - fp)),
((i, k + 1), (1 - fy) * fp), ((i + 1, k + 1), fy * fp)]
def get(self, hy, hp, mode):
if mode == "none":
return 0.0, 0.0
f = self.features(hy, hp, mode)
cy = sum(a * b for a, b in zip(self.a[0], f))
cp = sum(a * b for a, b in zip(self.a[1], f))
if mode.endswith("+grid"):
for (i, k), w in self.weights(hy, hp):
cy += w * self.grid[i][k][0]
cp += w * self.grid[i][k][1]
return cy, cp
def learn(self, hy, hp, dy, dp, mode, rate):
"""One sample: the correction here should have been (dy, dp) degrees more.
Normalized LMS for the polynomial; the grid takes half when it's on."""
if mode == "none":
return
self.samples += 1
self.extend(hy, hp)
grid = mode.endswith("+grid")
share = rate * 0.5 if grid else rate
f = self.features(hy, hp, mode)
norm = sum(v * v for v in f)
for row, d in ((self.a[0], dy), (self.a[1], dp)):
for n in range(self.NF):
row[n] += share * d * f[n] / norm
if grid:
rest = rate - share
for (i, k), w in self.weights(hy, hp):
self.grid[i][k][0] += rest * w * dy
self.grid[i][k][1] += rest * w * dp
def fit(self, points, mode, ridge=0.05, smooth=0.3):
"""Batch fit from (hy, hp, dy, dp) points, each the whole error there (degrees)."""
self.reset()
if mode == "none" or not points:
return
self.samples = len(points)
for p in points:
self.extend(p[0], p[1])
kind = self.base(mode)
used = {"offset": 1, "affine": 3, "quadratic": 6}[kind]
if used > 1 and len(points) < used + 2: # too few spots for this many terms
kind, used = ("affine", 3) if len(points) >= 5 else ("offset", 1)
if kind == "offset":
self.a[0][0] = statistics.fmean(p[2] for p in points)
self.a[1][0] = statistics.fmean(p[3] for p in points)
else:
# Least squares, with a little ridge on everything but the offset, so a
# lopsided set of spots can't bend it far.
X = [self.features(p[0], p[1], kind)[:used] for p in points]
M = [[sum(x[r] * x[c] for x in X) + (ridge * len(X) if r == c and r else 0.0) for c in range(used)]
for r in range(used)]
for out, col in ((self.a[0], 2), (self.a[1], 3)):
b = [sum(x[r] * p[col] for x, p in zip(X, points)) for r in range(used)]
out[:used] = solve(M, b)
if not mode.endswith("+grid"):
return
# Each node: the weighted mean of what the polynomial left over near it, shrunk toward 0.
acc = [[[0.0, 0.0, 0.0] for _ in self.PITCHES] for _ in self.YAWS]
for hy, hp, dy, dp in points:
ly, lp = self.get(hy, hp, kind)
for (i, k), w in self.weights(hy, hp):
acc[i][k][0] += w * (dy - ly)
acc[i][k][1] += w * (dp - lp)
acc[i][k][2] += w
for i in range(len(self.YAWS)):
for k in range(len(self.PITCHES)):
sy, sp, sw = acc[i][k]
self.grid[i][k] = [sy / (sw + smooth), sp / (sw + smooth)]
def offset(self):
return self.a[0][0], self.a[1][0]
def to_json(self):
return {"coef": self.a, "range": self.range, "grid": self.grid, "samples": self.samples}
def from_json(self, d):
self.reset()
a = d.get("coef") or d.get("affine") # "affine": the 3-term version of this file
if a and len(a) == 2 and all(len(r) in (3, self.NF) for r in a):
self.a = [list(map(float, r)) + [0.0] * (self.NF - len(r)) for r in a]
grid = d.get("grid")
if grid and len(grid) == len(self.YAWS) and all(len(r) == len(self.PITCHES) for r in grid):
self.grid = [[list(c) for c in row] for row in grid]
r = d.get("range")
self.range = list(map(float, r)) if r and len(r) == 4 else None
self.samples = d.get("samples", 0)
def solve(M, b):
"""Solve a small linear system (Gaussian elimination with pivoting); zeros if singular."""
n = len(b)
A = [row[:] + [b[i]] for i, row in enumerate(M)]
for c in range(n):
piv = max(range(c, n), key=lambda r: abs(A[r][c]))
if abs(A[piv][c]) < 1e-12:
return [0.0] * n
A[c], A[piv] = A[piv], A[c]
for r in range(n):
if r != c:
f = A[r][c] / A[c][c]
for k in range(c, n + 1):
A[r][k] -= f * A[c][k]
return [A[i][n] / A[i][i] for i in range(n)]
class LiveCorrection:
"""Corrections learned on the fly from snapped clicks, on top of the calibration.
Each click on an element is a measurement: you were looking at that element when you
pressed, and the tracker put your gaze at the raw point, so the gap between them is the
whole error there. Whatever the calibration doesn't already explain (the residual) is
fitted with the same quadratic terms. Every term but the offset is held close to zero
(ridge), so one click shifts the whole correction and more clicks bend it. Whatever is
still left near a click is added within a few degrees of it: on the Frame, errors less
than 3 degrees apart are alike, and ones further apart are unrelated. Recent clicks count more (a
half-life counted in clicks), so it follows SteamVR's gaze as that drifts or relearns.
Replayed on logged points: a calibration from an earlier session was 4.95 degrees off;
one click brought that to 2.3, five to 1.6, twenty to 1.2.
A big element says little about where on it you looked, so each axis is weighted by the
element's size along it: a list row 12 degrees wide barely counts sideways."""
RIDGE = [0.01, 0.5, 0.5, 0.5, 0.5, 0.5]
KERNEL = 2.5 # degrees: how far a click's leftover reaches
SHRINK = 0.5 # near one click, half its leftover; near several, nearly all
HALF_LIFE = 40 # clicks
KEEP = 150
MARGIN = 3.0
def __init__(self):
self.samples = [] # dicts: time, hy, hp, dy, dp (the whole error), wy, wp
self.reset_fit()
def reset_fit(self):
self.cy = [0.0] * 6
self.cp = [0.0] * 6
self.left = [] # (hy, hp, leftover yaw, leftover pitch, wy, wp, decay)
self.range = None
def features(self, hy, hp):
if self.range:
y0, y1, p0, p1 = self.range
hy = min(max(hy, y0 - self.MARGIN), y1 + self.MARGIN)
hp = min(max(hp, p0 - self.MARGIN), p1 + self.MARGIN)
x, y = hy / 30.0, hp / 30.0
return [1.0, x, y, x * y, x * x, y * y]
def add(self, sample, base, mode):
self.samples = (self.samples + [sample])[-self.KEEP:]
self.refit(base, mode)
def undo(self, base, mode):
if self.samples:
self.samples.pop()
self.refit(base, mode)
def refit(self, base, mode):
"""Refit from the samples against the calibration as it is now."""
self.reset_fit()
n = len(self.samples)
if not n:
return
self.range = [min(s["hy"] for s in self.samples), max(s["hy"] for s in self.samples),
min(s["hp"] for s in self.samples), max(s["hp"] for s in self.samples)]
rows = []
for k, s in enumerate(self.samples):
decay = 0.5 ** ((n - 1 - k) / self.HALF_LIFE)
by, bp = base.get(s["hy"], s["hp"], mode)
rows.append((s, self.features(s["hy"], s["hp"]), s["dy"] - by, s["dp"] - bp, decay))
for out, ri, wi in ((self.cy, 2, "wy"), (self.cp, 3, "wp")):
M = [[self.RIDGE[r] if r == c else 0.0 for c in range(6)] for r in range(6)]
b = [0.0] * 6
for row in rows:
w = row[4] * row[0][wi]
f = row[1]
for r in range(6):
b[r] += w * f[r] * row[ri]
for c in range(6):
M[r][c] += w * f[r] * f[c]
out[:] = solve(M, b)
for s, f, ry, rp, decay in rows:
ly = ry - sum(a * v for a, v in zip(self.cy, f))
lp = rp - sum(a * v for a, v in zip(self.cp, f))
self.left.append((s["hy"], s["hp"], ly, lp, s["wy"] * decay, s["wp"] * decay))
def get(self, hy, hp):
if not self.samples:
return 0.0, 0.0
f = self.features(hy, hp)
cy = sum(a * v for a, v in zip(self.cy, f))
cp = sum(a * v for a, v in zip(self.cp, f))
k2 = 2 * self.KERNEL ** 2
sy = sp = wy = wp = 0.0
for y, p, ly, lp, ay, ap in self.left:
d2 = (y - hy) ** 2 + (p - hp) ** 2
if d2 > 9 * k2:
continue
g = math.exp(-d2 / k2)
sy += g * ay * ly
wy += g * ay
sp += g * ap * lp
wp += g * ap
return cy + sy / (wy + self.SHRINK), cp + sp / (wp + self.SHRINK)
def offset(self):
return self.cy[0], self.cp[0]
class SteamEyeLog:
"""Follows SteamVR's eye tracking log (read only) for what moves the raw gaze under a
calibration.
SteamVR's eye tracker calibrates itself from clicks: a quick mouse-button down and up
(the laser or the Frametop pointer), with the gaze within 5 degrees of the click and
held still, is taken as "you were looking there" ("Accept usercal"). Accepted clicks
were all under 0.14 s; 0.38 s was "too slow", and one that moved was refused. It keeps that in
the running `eyetracking` process and saves nothing, so when the process starts again
(SteamVR restarting), its calibration starts over. Both events are counted here."""
PATH = Path.home() / ".local" / "share" / "Steam" / "logs" / "eyetracking.txt"
def __init__(self):
self.pos = 0
self.inode = None
self.partial = ""
self.starts = [] # when the eyetracking process started
self.accepts = [] # when it learned from a click
self.rejects = []
@staticmethod
def stamp(line):
head = line.split(" [", 1)[0]
main, _, frac = head.partition(".")
try:
return time.mktime(time.strptime(main.strip(), "%a %b %d %Y %H:%M:%S")) + float("0." + (frac or "0"))
except ValueError:
return None
def poll(self):
"""Read what's new. True if the eye tracker started again since the last poll."""
try:
st = os.stat(self.PATH)
except OSError:
return False
if st.st_ino != self.inode or st.st_size < self.pos:
self.inode, self.pos, self.partial = st.st_ino, 0, ""
if st.st_size == self.pos:
return False
first = self.pos == 0 and not self.starts
try:
with open(self.PATH, "rb") as f:
f.seek(self.pos)
data = f.read()
except OSError:
return False
self.pos += len(data)
lines = (self.partial + data.decode("utf-8", "replace")).split("\n")
self.partial = lines.pop()
restarted = False
for line in lines:
if "usercal" not in line and "startup with PID" not in line:
continue
t = self.stamp(line)
if t is None:
continue
if "startup with PID" in line:
self.starts.append(t)
restarted = not first
elif "Accept usercal" in line:
self.accepts.append(t)
elif "Reject usercal" in line:
self.rejects.append(t)
return restarted
def started(self):
return self.starts[-1] if self.starts else None
def accepted_since(self, t):
return sum(1 for a in self.accepts if a >= t)
# The math, filters, correction models, and SteamVR log reader are shared with ft-gazed.
sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
from gazecal import (DEFAULT_MODEL, MODELS, Correction, Fixation, LiveCorrection, OneEuro, # noqa: E402
SteamEyeLog, cross_validate, deg_from_px, px_from_deg, steady_samples)
# --- Gaze from ft-gaze ----------------------------------------------------------------
@@ -2965,45 +2531,6 @@ def region_errors(targets, source, key="cerr_deg"):
return sorted(((r, statistics.fmean(v), len(v)) for r, v in by.items()), key=lambda x: -x[1])
def cross_validate(points, mode):
"""Leave-one-out: each point's error under a model fitted on all the others (degrees)."""
errs = []
for i in range(len(points)):
c = Correction()
c.fit(points[:i] + points[i + 1:], mode)
cy, cp = c.get(points[i][0], points[i][1], mode)
errs.append(math.hypot(points[i][2] - cy, points[i][3] - cp))
return errs
def steady_samples(samples, vergence_jump=1.5):
"""The samples of one look at one spot where the tracker had both eyes: none in a blink
(openness under half its median over the samples), and none where the angle between the eyes' directions (`lr`, the
vergence) is more than `vergence_jump` degrees from its median over the samples. The
vergence itself depends on distance (about 2.8 degrees for a screen 1.3 m away, a
fraction of one far off), so only a jump away from what it was during this look means
the tracker lost an eye. Without the mmap there's nothing to judge by: all are kept."""
# Openness: a blink is a sharp drop from what it was during this look. Not a fixed
# level: looking down, the upper lids come down with the eyes, and in bright light you
# squint, so the reading can stay under 0.5 for the whole look while the tracker follows
# the eyes fine (a calibration dot at the bottom of the bright round failed that way).
opens = [min(o) for o in ((smp["src"].get("mmap1") or {}).get("open") for smp in samples) if o]
floor = max(0.12, 0.5 * statistics.median(opens)) if len(opens) >= 5 else 0.12
opened = []
for smp in samples:
o = (smp["src"].get("mmap1") or {}).get("open")
if not o or min(o) >= floor:
opened.append(smp)
def vergence(smp):
return (smp["src"].get("mmap1") or {}).get("lr", (smp["src"].get("mmap2") or {}).get("lr"))
have = [v for v in map(vergence, opened) if v is not None]
if len(have) < 5:
return opened
med = statistics.median(have)
return [smp for smp in opened if vergence(smp) is None or abs(vergence(smp) - med) <= vergence_jump]
def main():
ap = argparse.ArgumentParser(description="Eye tracking playground for the Frametop desktop")
ap.add_argument("--screen", type=int, default=0, help="Frametop screen to open on (default: where it opens)")