nx_motion.nx source
↩ module page · 459 lines · 18817 B
1// nx_motion.nx -- frame-to-frame motion (Lucas-Kanade optical flow).
2//
3// Brightness-constancy + small-motion assumption:
4// I(x, y, t) = I(x + u, y + v, t + 1)
5// Taylor expand:
6// I_x * u + I_y * v + I_t = 0
7// Over an N-pixel window this is N equations in 2 unknowns; solve
8// via 2x2 least-squares (normal equations + Cramer's rule).
9//
10// What this unlocks:
11// + motion detection (frame diff)
12// + tracking (per-pixel flow vectors)
13// + camera-motion estimation (mean flow direction)
14// + scene-cut detection (flow magnitude spike)
15// + activity recognition substrate
16//
17// Returns flow in Q8 sub-pixel precision (1 unit = 1/256 pixel).
18// 5x5 integration window default.
19//
20// genealogy_id: lucas_kanade_1981 + horn_schunck_1981
21// lineage_id: gradient_based_optical_flow + 2x2_least_squares
22
23// nx_safety_envelope:
24// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
25// sil_target: SIL1
26// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
27// verdict: NOT_YET_EVALUATED
28
29import "syscalls.nx"
30import "nx_image.nx"
31
32const NX_MOTION_Q: i64 = 256 // Q8 sub-pixel
33const NX_MOTION_WIN: i64 = 5 // 5x5 integration
34// THE ESTIMATOR'S OWN DOMAIN OF VALIDITY, DERIVED FROM THE WINDOW AND NOTHING ELSE. A 5x5
35// integration window cannot observe a displacement larger than its own radius: beyond that the two
36// frames share no structure inside the window and the 2x2 solve is reading noise. So a per-pixel
37// result larger than the half-window is not a large measurement, it is the system being
38// rank-deficient -- the aperture problem -- and `det != 0` is the only guard against it, which
39// passes for a determinant of one.
40const NX_MOTION_LK_DOMAIN_Q8: i64 = (NX_MOTION_WIN / 2) * NX_MOTION_Q
41
42// ===== Frame difference =================================================
43//
44// Absolute difference between two same-size grayscale images.
45
46func nx_motion_frame_diff(a: *Image, b: *Image) -> *Image {
47 let w: i64 = a.width
48 let h: i64 = a.height
49 let dst: *Image = nx_image_alloc(w, h, 1)
50 var y: i64 = 0
51 while y < h {
52 var x: i64 = 0
53 while x < w {
54 let av: i64 = nx_image_get(a, x, y, 0)
55 let bv: i64 = nx_image_get(b, x, y, 0)
56 var d: i64 = av - bv
57 if d < 0 { d = -d }
58 if d > 255 { d = 255 }
59 nx_image_set(dst, x, y, 0, d)
60 x = x + 1
61 }
62 y = y + 1
63 }
64 return dst
65}
66
67// ===== Temporal derivative ============================================
68//
69// I_t = frame2 - frame1. Signed output in ImageS64.
70
71func nx_motion_temporal_grad_into(frame1: *Image, frame2: *Image, dst: *ImageS64) -> *ImageS64 {
72 let w: i64 = frame1.width
73 let h: i64 = frame1.height
74 // 'dst' is CALLER-OWNED scratch. The allocating wrapper sits at the end of this file.
75 var y: i64 = 0
76 while y < h {
77 var x: i64 = 0
78 while x < w {
79 let v: i64 = nx_image_get(frame2, x, y, 0) - nx_image_get(frame1, x, y, 0)
80 nx_image_s64_set(dst, x, y, v)
81 x = x + 1
82 }
83 y = y + 1
84 }
85 return dst
86}
87
88// ===== Lucas-Kanade per-pixel flow =====================================
89//
90// Computes flow (u, v) at each pixel by solving the 2x2 normal
91// equations over a 5x5 window of spatial+temporal gradients.
92//
93// Pixels with insufficient texture (det close to zero) get u=v=0.
94//
95// Returns two ImageS64 maps: flow_u (Q8) and flow_v (Q8).
96
97struct FlowField {
98 u: *ImageS64, // Q8 horizontal flow
99 v: *ImageS64, // Q8 vertical flow
100}
101
102// SCRATCH-REUSING VARIANT. Ix/Iy/It and field.u/field.v are ALL caller-owned. The allocating form
103// at the end of this file delegates here, so there is one Lucas-Kanade in the estate.
104// WHY THE CLEAR LOOP BELOW IS NOT OPTIONAL: LK writes a pixel only where the 2x2 determinant is
105// non-zero, and it skips a half-window border entirely. On a FRESH mmap those pixels are zero and
106// the result is correct by accident. On a REUSED buffer they hold the PREVIOUS frame's flow, so a
107// low-texture region would silently report the last frame it could see. That is a stale-value
108// defect that produces plausible motion instead of an error, so the buffers are cleared up front.
109func nx_motion_lucas_kanade_into(frame1: *Image, frame2: *Image,
110 Ix: *ImageS64, Iy: *ImageS64, It: *ImageS64,
111 field: *FlowField) -> *FlowField {
112 let w: i64 = frame1.width
113 let h: i64 = frame1.height
114
115 // Compute spatial + temporal gradients on frame1.
116 nx_image_sobel_x_into(frame1, Ix)
117 nx_image_sobel_y_into(frame1, Iy)
118 nx_motion_temporal_grad_into(frame1, frame2, It)
119
120 let flow_u: *ImageS64 = field.u
121 let flow_v: *ImageS64 = field.v
122 // Clear both flow buffers: see the header note. A reused buffer must not inherit last frame's flow.
123 var cy: i64 = 0
124 while cy < h {
125 var cx: i64 = 0
126 while cx < w {
127 nx_image_s64_set(flow_u, cx, cy, 0)
128 nx_image_s64_set(flow_v, cx, cy, 0)
129 cx = cx + 1
130 }
131 cy = cy + 1
132 }
133
134
135 let half: i64 = NX_MOTION_WIN / 2
136 var y: i64 = half
137 while y < h - half {
138 var x: i64 = half
139 while x < w - half {
140 // Accumulate over 5x5 window.
141 var sxx: i64 = 0
142 var sxy: i64 = 0
143 var syy: i64 = 0
144 var sxt: i64 = 0
145 var syt: i64 = 0
146 var dy: i64 = -half
147 while dy <= half {
148 var dx: i64 = -half
149 while dx <= half {
150 let ix: i64 = nx_image_s64_get(Ix, x + dx, y + dy)
151 let iy: i64 = nx_image_s64_get(Iy, x + dx, y + dy)
152 let it: i64 = nx_image_s64_get(It, x + dx, y + dy)
153 sxx = sxx + ix * ix
154 sxy = sxy + ix * iy
155 syy = syy + iy * iy
156 sxt = sxt + ix * it
157 syt = syt + iy * it
158 dx = dx + 1
159 }
160 dy = dy + 1
161 }
162 // 2x2 system: [sxx sxy] [u] [-sxt]
163 // [sxy syy] [v] = [-syt]
164 // det = sxx * syy - sxy^2
165 let det: i64 = sxx * syy - sxy * sxy
166 if det != 0 {
167 // u = (syy * (-sxt) - sxy * (-syt)) / det
168 // v = (sxx * (-syt) - sxy * (-sxt)) / det
169 // Scale up by Q for sub-pixel precision.
170 let u_num: i64 = (sxy * syt - syy * sxt) * NX_MOTION_Q
171 let v_num: i64 = (sxy * sxt - sxx * syt) * NX_MOTION_Q
172 nx_image_s64_set(flow_u, x, y, u_num / det)
173 nx_image_s64_set(flow_v, x, y, v_num / det)
174 }
175 x = x + 1
176 }
177 y = y + 1
178 }
179
180 return field
181}
182
183// ===== Mean flow ========================================================
184//
185// Mean (u, v) over high-confidence pixels (those where |flow| >= threshold).
186// Useful for global camera-motion estimate.
187
188// ROI VARIANT (added 2026-08-25 for nx_dynaoracle). Mean flow over a RECTANGLE, so a caller can
189// ask about one region -- a tissue patch versus a rigid torso patch -- instead of the whole frame.
190// nx_motion_mean_flow below now DELEGATES to this with the full-frame rect, so there is exactly
191// one flow-averaging implementation in the estate and every existing caller is bit-identical BY
192// CONSTRUCTION. The rect is clamped to the image, and the returned count is the caller's coverage
193// number: n==0 means "no pixel in this rect cleared the confidence threshold", which is an
194// UNOBSERVED region and must never be read as "this region did not move".
195// BOUNDED VARIANT (2026-08-25). max_mag_q8 <= 0 means NO BOUND, which is the incumbent contract
196// exactly; nx_motion_mean_flow_rect below delegates here with 0, so every existing caller is
197// bit-identical BY CONSTRUCTION and only the pyramid opts in. out_rejected may be 0.
198//
199// WHY THIS EXISTS, MEASURED. The pyramid DOUBLES each level's estimate on the way down, so a single
200// rank-deficient window at the coarsest level arrives at the finest level multiplied by 32. On a
201// synthetic clip textured only at 17, 13 and 29 pixels -- every one of which the binomial decimation
202// annihilates within two halvings -- the coarse levels see an almost one-dimensional field, the 2x2
203// determinant collapses toward zero, and the quotient explodes. Measured end-to-end on that clip:
204// a rigid 0.5 Hz reference drift was recovered as 2.0 Hz at ALL EIGHTEEN grid points, and the
205// integrated series carried 2 pixels per frame of manufactured drift against a true 0.42. Holding
206// the pyramid to ONE level -- the same code, the same clip, only the depth changed -- recovered
207// 0.535 Hz and cut the manufactured drift by a factor of 4,148. That is what indicted the coarse
208// levels rather than the estimator.
209//
210// A REJECTED PIXEL IS COUNTED, NEVER SILENTLY DROPPED: a population that shrinks without saying so
211// is a measurement nobody knows is partial.
212func nx_motion_mean_flow_rect_bounded(flow: *FlowField, mag_threshold: i64, max_mag_q8: i64,
213 rx0: i64, ry0: i64, rx1: i64, ry1: i64,
214 out_u: *i64, out_v: *i64, out_rejected: *i64) -> i64 {
215 let w: i64 = flow.u.width
216 let h: i64 = flow.u.height
217 var x0: i64 = rx0
218 var y0: i64 = ry0
219 var x1: i64 = rx1
220 var y1: i64 = ry1
221 if x0 < 0 { x0 = 0 }
222 if y0 < 0 { y0 = 0 }
223 if x1 > w { x1 = w }
224 if y1 > h { y1 = h }
225 var sum_u: i64 = 0
226 var sum_v: i64 = 0
227 var n: i64 = 0
228 var rej: i64 = 0
229 var y: i64 = y0
230 while y < y1 {
231 var x: i64 = x0
232 while x < x1 {
233 let uv: i64 = nx_image_s64_get(flow.u, x, y)
234 let vv: i64 = nx_image_s64_get(flow.v, x, y)
235 var au: i64 = uv
236 if au < 0 { au = -au }
237 var av: i64 = vv
238 if av < 0 { av = -av }
239 if au + av >= mag_threshold {
240 var admit: i64 = 1
241 if max_mag_q8 > 0 {
242 if au > max_mag_q8 { admit = 0 }
243 if av > max_mag_q8 { admit = 0 }
244 }
245 if admit == 1 {
246 sum_u = sum_u + uv
247 sum_v = sum_v + vv
248 n = n + 1
249 }
250 else { rej = rej + 1 }
251 }
252 x = x + 1
253 }
254 y = y + 1
255 }
256 if (out_rejected as i64) != 0 { out_rejected[0] = rej }
257 if n == 0 {
258 out_u[0] = 0
259 out_v[0] = 0
260 return 0
261 }
262 out_u[0] = sum_u / n
263 out_v[0] = sum_v / n
264 return n
265}
266
267// THE ORIGINAL CONTRACT, expressed over the bounded one with NO bound. Delegation is what makes
268// "every existing caller is unchanged" a fact rather than a hope.
269func nx_motion_mean_flow_rect(flow: *FlowField, mag_threshold: i64,
270 rx0: i64, ry0: i64, rx1: i64, ry1: i64,
271 out_u: *i64, out_v: *i64) -> i64 {
272 return nx_motion_mean_flow_rect_bounded(flow, mag_threshold, 0, rx0, ry0, rx1, ry1,
273 out_u, out_v, 0 as *i64)
274}
275
276// Whole-frame mean flow: the ORIGINAL contract, now expressed as the full-frame rect. Callers of
277// this function are unchanged -- the delegation is what makes that a fact rather than a hope.
278func nx_motion_mean_flow(flow: *FlowField, mag_threshold: i64,
279 out_u: *i64, out_v: *i64) -> i64 {
280 return nx_motion_mean_flow_rect(flow, mag_threshold,
281 0, 0, flow.u.width, flow.u.height, out_u, out_v)
282}
283
284// ===== ALLOCATING FORMS =================================================
285//
286// These are the ORIGINAL contracts, unchanged for every existing caller. Each allocates the
287// scratch the _into variant needs and delegates. A caller that runs per-frame should call the
288// _into form and hoist the buffers; a caller that runs once should keep using these.
289
290func nx_motion_temporal_grad(frame1: *Image, frame2: *Image) -> *ImageS64 {
291 let dst: *ImageS64 = nx_image_s64_alloc(frame1.width, frame1.height)
292 return nx_motion_temporal_grad_into(frame1, frame2, dst)
293}
294
295func nx_motion_lucas_kanade(frame1: *Image, frame2: *Image) -> *FlowField {
296 let w: i64 = frame1.width
297 let h: i64 = frame1.height
298 let Ix: *ImageS64 = nx_image_s64_alloc(w, h)
299 let Iy: *ImageS64 = nx_image_s64_alloc(w, h)
300 let It: *ImageS64 = nx_image_s64_alloc(w, h)
301 let field: *FlowField = (sys_mmap(16)) as *FlowField
302 field.u = nx_image_s64_alloc(w, h)
303 field.v = nx_image_s64_alloc(w, h)
304 return nx_motion_lucas_kanade_into(frame1, frame2, Ix, Iy, It, field)
305}
306
307// ===== PYRAMIDAL MEAN FLOW (2026-08-25) =================================
308//
309// WHY: the single-scale estimator above is correct and BOUNDED -- its own header declares the
310// small-motion assumption. MEASURED on a synthetic fixture with a 14-pixel oscillation at 2 Hz and
311// 30 fps (peak 5.9 px/frame), it recovered about 1 px of a 14 px signal and the consumer correctly
312// reported UNOBSERVABLE. That is the declared limit doing its job, not a defect -- and the remedy
313// is the standard one: estimate coarse-to-fine so every level sees a small displacement.
314//
315// This COMPOSES nx_motion_lucas_kanade_into and nx_motion_mean_flow_rect. There is still exactly
316// one optical-flow implementation in the estate; this only chooses what to feed it.
317//
318// The workspace is allocated ONCE per region size and reused for every frame pair, because a
319// pyramid rebuilt per frame is the hot-loop allocation this file was just refactored to remove.
320
321const NX_MOTION_PYR_MAX: i64 = 6
322
323struct MotionPyr {
324 levels: i64,
325 a: *i64,
326 b: *i64,
327 warp: *i64,
328 Ix: *i64,
329 Iy: *i64,
330 It: *i64,
331 field: *FlowField,
332 rejected: i64,
333}
334const NX_MOTION_PYR_BYTES: i64 = 72
335
336// DERIVED level count: keep halving while both axes stay at or above the integration window, and
337// never exceed NX_MOTION_PYR_MAX. A level smaller than the window has no interior pixels and would
338// contribute a structurally empty estimate.
339func nx_motion_pyr_levels(w: i64, h: i64) -> i64 {
340 var n: i64 = 1
341 var cw: i64 = w
342 var ch: i64 = h
343 while n < NX_MOTION_PYR_MAX {
344 let nw: i64 = cw / 2
345 let nh: i64 = ch / 2
346 if nw < NX_MOTION_WIN { return n }
347 if nh < NX_MOTION_WIN { return n }
348 cw = nw
349 ch = nh
350 n = n + 1
351 }
352 return n
353}
354
355func nx_motion_pyr_new(w: i64, h: i64) -> *MotionPyr {
356 let p: *MotionPyr = (sys_mmap(NX_MOTION_PYR_BYTES)) as *MotionPyr
357 let n: i64 = nx_motion_pyr_levels(w, h)
358 p.levels = n
359 p.a = (sys_mmap(n * 8 + 16)) as *i64
360 p.b = (sys_mmap(n * 8 + 16)) as *i64
361 p.warp = (sys_mmap(n * 8 + 16)) as *i64
362 p.Ix = (sys_mmap(n * 8 + 16)) as *i64
363 p.Iy = (sys_mmap(n * 8 + 16)) as *i64
364 p.It = (sys_mmap(n * 8 + 16)) as *i64
365 let f: *FlowField = (sys_mmap(16)) as *FlowField
366 f.u = nx_image_s64_alloc(w, h)
367 f.v = nx_image_s64_alloc(w, h)
368 p.field = f
369 p.rejected = 0
370 var i: i64 = 0
371 var cw: i64 = w
372 var ch: i64 = h
373 while i < n {
374 p.a[i] = (nx_image_alloc(cw, ch, 1)) as i64
375 p.b[i] = (nx_image_alloc(cw, ch, 1)) as i64
376 p.warp[i] = (nx_image_alloc(cw, ch, 1)) as i64
377 p.Ix[i] = (nx_image_s64_alloc(cw, ch)) as i64
378 p.Iy[i] = (nx_image_s64_alloc(cw, ch)) as i64
379 p.It[i] = (nx_image_s64_alloc(cw, ch)) as i64
380 cw = cw / 2
381 ch = ch / 2
382 i = i + 1
383 }
384 return p
385}
386
387func nx_motion_copy_into(src: *Image, dst: *Image) -> *Image {
388 var y: i64 = 0
389 while y < dst.height {
390 var x: i64 = 0
391 while x < dst.width {
392 nx_image_set(dst, x, y, 0, nx_image_get(src, x, y, 0))
393 x = x + 1
394 }
395 y = y + 1
396 }
397 return dst
398}
399
400// Mean flow over the whole pair, coarse-to-fine. Returns the pixel count that cleared the threshold
401// at the BEST-SUPPORTED level -- see the note inside. out_u/out_v are Q8.
402func nx_motion_mean_flow_pyr(p: *MotionPyr, f1: *Image, f2: *Image, mag_threshold: i64,
403 out_u: *i64, out_v: *i64) -> i64 {
404 let n: i64 = p.levels
405 nx_motion_copy_into(f1, (p.a[0]) as *Image)
406 nx_motion_copy_into(f2, (p.b[0]) as *Image)
407 var i: i64 = 1
408 while i < n {
409 nx_image_downsample2_into((p.a[i - 1]) as *Image, (p.a[i]) as *Image)
410 nx_image_downsample2_into((p.b[i - 1]) as *Image, (p.b[i]) as *Image)
411 i = i + 1
412 }
413 var gu: i64 = 0
414 var gv: i64 = 0
415 var cnt: i64 = 0
416 let lu: *i64 = (sys_mmap(16)) as *i64
417 let lv: *i64 = (sys_mmap(16)) as *i64
418 let lrej: *i64 = (sys_mmap(16)) as *i64
419 lrej[0] = 0
420 var acc_u_q8: i64 = 0
421 var acc_v_q8: i64 = 0
422 var l: i64 = n - 1
423 while l >= 0 {
424 let A: *Image = (p.a[l]) as *Image
425 let B: *Image = (p.b[l]) as *Image
426 let Wp: *Image = (p.warp[l]) as *Image
427 nx_image_shift_into(B, gu, gv, Wp)
428 let fld: *FlowField = p.field
429 fld.u.width = A.width
430 fld.u.height = A.height
431 fld.v.width = A.width
432 fld.v.height = A.height
433 nx_motion_lucas_kanade_into(A, Wp, (p.Ix[l]) as *ImageS64, (p.Iy[l]) as *ImageS64,
434 (p.It[l]) as *ImageS64, fld)
435 // COVERAGE IS THE BEST-SUPPORTED LEVEL, NOT THE FINEST. In coarse-to-fine, a near-zero
436 // residual at the finest level means the COARSER levels already captured the motion -- that
437 // is the estimator succeeding, and reporting the finest count as coverage would mark a good
438 // estimate as unobserved.
439 // HOLD EVERY LEVEL TO THE ESTIMATOR'S DECLARED DOMAIN. A residual larger than the window's
440 // own radius is the aperture problem, not motion, and at level l it would be multiplied by
441 // 2^l before it reached the answer. A level whose pixels are all refused contributes a mean
442 // of zero, which leaves the coarser guess untouched -- the level ABSTAINS instead of
443 // poisoning, which is the right behaviour for an axis that cannot see.
444 let lvl_cnt: i64 = nx_motion_mean_flow_rect_bounded(fld, mag_threshold,
445 NX_MOTION_LK_DOMAIN_Q8, 0, 0, A.width, A.height, lu, lv, lrej)
446 p.rejected = p.rejected + lrej[0]
447 if lvl_cnt > cnt { cnt = lvl_cnt }
448 acc_u_q8 = gu * NX_MOTION_Q + lu[0]
449 acc_v_q8 = gv * NX_MOTION_Q + lv[0]
450 if l > 0 {
451 gu = acc_u_q8 * 2 / NX_MOTION_Q
452 gv = acc_v_q8 * 2 / NX_MOTION_Q
453 }
454 l = l - 1
455 }
456 out_u[0] = acc_u_q8
457 out_v[0] = acc_v_q8
458 return cnt
459}