The current contact detector in kinematic_regimes.py works by scanning per-axis velocity sign reversals. It is a cheap proxy for the floor-bounce case and fails in geometries where the contact normal is not aligned with a coordinate axis. This note works out the principled replacement: a frame-invariant, impulse-magnitude detector with a closed-form normal estimator, plus a paired variant for two-body contacts. It also sketches a wavelet-based extension when SNR is poor.
The current scanner looks for $v_y$ flipping sign (floor bounce) or $v_x$ flipping sign (wall on the side). Cases it cannot see, in increasing severity:
| scenario | what the velocity does | v1 verdict |
|---|---|---|
| Vertical floor bounce | $v_y$ flips, $v_x$ ≈ const | detected |
| Angled wall, ball continues down | $v_x$ flips, $v_y$ continues (gravity) | only if axis="x" is scanned; missed by axis="y" default |
| Glancing blow on a slope | $v_x$ and $v_y$ both change direction but neither reverses sign | missed entirely |
| Galilean cannon (ball A hits ball B mid-air) | $v_A$ drops, $v_B$ jumps; no per-axis sign flip on either | missed |
| Pendulum string going taut | radial $v$ jumps to zero; tangential $v$ continues | missed |
Two issues compound: the detector is coordinate-system dependent (rotate the frame and the answer changes), and it confuses direction reversal with sign reversal of one component — those are different things.
Mechanically, a contact event is the application of an impulsive force over a vanishingly small interval. By the impulse-momentum theorem,
$$ \Delta \mathbf{p} = \int_{t-\epsilon}^{t+\epsilon} \mathbf{F}(\tau)\,d\tau = \mathbf{J}, \qquad \epsilon \to 0. $$Position is continuous through the event; velocity has a jump discontinuity; acceleration is approximately a Dirac delta. At the sample scale this manifests as a brief, anomalously large spike in $\|\mathbf{a}\|$ — and crucially, the spike is invariant under rotation of the coordinate frame because $\|\mathbf{a}\|$ is a scalar.
This gives the universal detector statement:
From samples $\{(t_i, \mathbf{r}_i)\}$, estimate velocity by smoothed central difference:
$$ \mathbf{v}_i = \frac{\mathbf{r}_{i+1} - \mathbf{r}_{i-1}}{t_{i+1} - t_{i-1}}, \qquad \mathbf{a}_i = \frac{\mathbf{v}_{i+1} - \mathbf{v}_{i-1}}{t_{i+1} - t_{i-1}}. $$Define the impulse-norm time series
$$ \alpha_i \;=\; \|\mathbf{a}_i\| \;=\; \sqrt{a_{x,i}^2 + a_{y,i}^2}. $$This is the rotation-invariant scalar we threshold.
Standard deviation is poisoned by the very impacts we are trying to detect (a single 50× spike inflates $\sigma$ enormously). Use median absolute deviation:
$$ \tilde\alpha = \mathrm{median}(\alpha_i), \qquad \mathrm{MAD} = \mathrm{median}\big(|\alpha_i - \tilde\alpha|\big), \qquad \hat\sigma_\alpha = 1.4826 \cdot \mathrm{MAD}. $$The factor 1.4826 makes $\hat\sigma_\alpha$ a consistent estimator of $\sigma$ for Gaussian noise. The MAD is robust to up to 50% of the samples being outliers — far more than any realistic impact rate.
Detect a contact when
$$ \alpha_i \;>\; \tilde\alpha + k\,\hat\sigma_\alpha \quad\text{and}\quad \alpha_i \;=\; \max_{|j-i|\le w}\alpha_j, $$with $k \in [6, 10]$ and a non-maximum-suppression window $w$ of a few frames.
At a detected peak, sample velocity in two clean windows skipping the central straddling frame:
$$ \mathbf{v}^- = \mathrm{mean}(\mathbf{v}_{i-s-2 \,..\, i-s}), \qquad \mathbf{v}^+ = \mathrm{mean}(\mathbf{v}_{i+s \,..\, i+s+2}). $$The contact normal is the unit vector of the velocity jump:
$$ \hat{\mathbf{n}} \;=\; \frac{\mathbf{v}^+ - \mathbf{v}^-}{\|\mathbf{v}^+ - \mathbf{v}^-\|}. $$The general (oblique) coefficient of restitution is the ratio of normal components, not full magnitudes:
$$ e \;=\; -\frac{(\mathbf{v}^+ - \mathbf{v}_\text{wall}) \cdot \hat{\mathbf{n}}}{(\mathbf{v}^- - \mathbf{v}_\text{wall}) \cdot \hat{\mathbf{n}}}. $$For a fixed wall ($\mathbf{v}_\text{wall} = \mathbf{0}$) this collapses to $e = -\mathbf{v}^+\cdot\hat{\mathbf{n}} \,/\, \mathbf{v}^-\cdot\hat{\mathbf{n}}$, which reduces to the familiar 1D form when $\hat{\mathbf{n}} = \hat{\mathbf{y}}$.
Decompose $\Delta\mathbf{v}$ into normal and tangential parts:
$$ \Delta\mathbf{v}_n = (\Delta\mathbf{v}\cdot\hat{\mathbf{n}})\hat{\mathbf{n}}, \qquad \Delta\mathbf{v}_t = \Delta\mathbf{v} - \Delta\mathbf{v}_n. $$If the impact slipped, the kinetic friction coefficient is
$$ \mu_k \;\approx\; \frac{\|\Delta\mathbf{v}_t\|}{(1+e)\,|\Delta\mathbf{v}\cdot\hat{\mathbf{n}}|}. $$If $\|\Delta\mathbf{v}_t\|$ is small relative to that bound, the contact stuck — useful information the v1 detector simply cannot extract.
For two trajectories $\mathbf{r}_A(t), \mathbf{r}_B(t)$ with radii $R_A, R_B$ and (optionally) masses $m_A, m_B$, declare a pair contact at $t^*$ when all of the following hold simultaneously:
Without masses, criteria 1–3 detect the event. With masses (or with a known $e$), the system has enough constraints to identify the missing parameter:
$$ \frac{m_A}{m_B} \;=\; \frac{\|\Delta\mathbf{v}_B\|}{\|\Delta\mathbf{v}_A\|} \quad\text{(1D collinear)}, \qquad e \;=\; -\frac{(\mathbf{v}_A^+ - \mathbf{v}_B^+)\cdot\hat{\mathbf{n}}}{(\mathbf{v}_A^- - \mathbf{v}_B^-)\cdot\hat{\mathbf{n}}}. $$This is exactly the situation in the Galilean cannon: the v1 detector misses the contact because no axis flips, but the proximity-+-coincident-spike test is unambiguous.
A constraint event is not an instantaneous collision but the sudden activation of a holonomic constraint $g(\mathbf{r}) = 0$. The signature is selective: the velocity component along $\nabla g$ jumps to zero (or rebounds), while the orthogonal components continue smoothly.
For a pendulum bob with pivot $\mathbf{p}$ and string length $L$, the constraint is $\|\mathbf{r} - \mathbf{p}\| = L$. Define
$$ \hat{\mathbf{r}}(t) \;=\; \frac{\mathbf{r}(t) - \mathbf{p}}{\|\mathbf{r}(t) - \mathbf{p}\|}, \qquad v_\parallel(t) \;=\; \mathbf{v}(t)\cdot \hat{\mathbf{r}}(t). $$String-taut events are local extrema of $\|\mathbf{r} - \mathbf{p}\|$ where $v_\parallel$ jumps from outward-positive to non-positive. The pivot $\mathbf{p}$ is recovered from the trajectory itself by least-squares circle fitting (the pipeline already does this in derived.json's pivot_guess_px).
For tracker traces that are noisy enough that $\hat\sigma_\alpha$ is dominated by jitter rather than smooth motion, the threshold-on-MAD scheme degrades. The mathematically clean fall-back is wavelet-based singularity detection (Mallat & Hwang 1992):
This is overkill for clean tracker output but it gives a principled answer when the data is bad. It also handles the case where the tracker emits a brief jitter spike that mimics an impulse — the jitter has $\alpha < 0$ (worse than a step) and its modulus maxima grow at fine scales rather than staying constant.
| case | v1 (per-axis sign flip) | v2 (Δv-magnitude + MAD) |
|---|---|---|
| Vertical floor bounce | ✓ | ✓ |
| Wall hit, axis-aligned | ✓ if scanning that axis | ✓ |
| Angled wall, partial reversal | missed | ✓ + correct normal |
| Glancing blow | missed | ✓ |
| Galilean cannon (pair contact) | missed | ✓ via paired test |
| String going taut | missed | ✓ via projection variant |
| Sliding → rolling onset (smooth) | false-positive risk on jitter | correctly rejects |
| Frame-invariant | no | yes |
| Returns contact normal | 2D heuristic merge | closed-form $\Delta\mathbf{v}/\|\Delta\mathbf{v}\|$ |
| Returns friction estimate | no | $\mu_k$ from tangential Δv |
def contacts_v2_from_trajectory(traj, k=8.0, skip=1, w=2):
ts, rs = load(traj) # rs is shape (N, 2)
vs = central_diff(rs, ts) # (N, 2)
accs = central_diff(vs, ts) # (N, 2)
alpha = np.linalg.norm(accs, axis=1)
a_med = np.median(alpha)
mad = np.median(np.abs(alpha - a_med))
sigma = 1.4826 * mad
thr = a_med + k * sigma
# local-maxima with non-max suppression of width w
peaks = nms_peaks(alpha, threshold=thr, half_width=w)
events = []
for i in peaks:
v_pre = vs[i - skip - 2 : i - skip].mean(axis=0)
v_post = vs[i + skip : i + skip + 2].mean(axis=0)
dv = v_post - v_pre
nhat = dv / (np.linalg.norm(dv) + 1e-12)
v_norm_pre = float(v_pre @ nhat)
v_norm_post = float(v_post @ nhat)
e = clip(-v_norm_post / v_norm_pre, 0.0, 1.0) if v_norm_pre != 0 else None
dv_t = dv - (dv @ nhat) * nhat
events.append({
"t": float(ts[i]),
"alpha_peak": float(alpha[i]),
"alpha_threshold": float(thr),
"normal": (float(nhat[0]), float(nhat[1])),
"v_before": tuple(map(float, v_pre)),
"v_after": tuple(map(float, v_post)),
"restitution_normal": e,
"tangential_dv_px_per_s": float(np.linalg.norm(dv_t)),
"friction_lower_bound": (
float(np.linalg.norm(dv_t) / ((1.0 + e) * abs(v_norm_pre)))
if e is not None and v_norm_pre != 0 else None
),
})
return events
contacts_v2_from_trajectory next to the existing function. Do not delete v1.compute_contacts_v2 orchestrator tool with the same interface plus an include_friction flag.compute_pair_contact(traj_a, traj_b); validate on dual-ball-drop.infer_regimes / regimes_from_trajectories to call v2. Keep v1 reachable behind a flag for one cycle for parity diffing.Written 2026-04-27 as a design doc for the cradle pipeline's contact detector. The implementation is not yet landed; this note exists to make the physics explicit before code is written.