code wiki / _hdl_build / nx_svmetrology.nx
nx_svmetrology.nx source
↩ module page · 371 lines · 18392 B
1// nx_svmetrology.nx -- THE MULTI-CUE SCALE SOLVER (debts 1785958292 + 1785957867; operator 2026-08-05:
2// "be clever but stay reproducible"). Single-view metrology done as an OVER-DETERMINED system.
3//
4// THE DOCTRINE, IN CODE: never choose a reference. Every scale cue in a scene is an independent estimate
5// of ONE unknown (camera height, in micrometres), so N cues leave N-1 degrees of freedom of VALIDATION.
6// Fuse by inverse-variance weighting, then TEST the agreement with chi-square. Below threshold: report
7// the interval. Above: REFUSE and NAME THE OUTLIER -- averaging away a disagreement hides the broken
8// assumption that caused it (a partition is a claim; check the parts sum).
9// A SINGLE CUE IS NEVER VALIDATED: 0 dof means there is nothing to check it against, and the organ says
10// so out loud rather than reporting a confident number.
11//
12// GEOMETRY (Criminisi/Reid/Zisserman single-view metrology): for a flat ground plane the horizon sits at
13// camera-optical-axis height, so for any object standing ON that plane, with no focal length needed:
14// h_obj / h_cam = (y_base - y_top) / (y_base - y_horizon)
15// LEVERAGE LAW (measured 2026-08-05): the error scales INVERSELY with (y_base - y_horizon), so a NEAR
16// reference beats a distant one however crisp the distant one looks. The organ reports that denominator
17// as `leverage_px` on every cue it derives, so a caller can see which reference is worth having.
18//
19// nx_svmetrology camheight <y_base> <y_top> <y_horizon> <h_obj_um> <sigma_permil> -> one cue row
20// nx_svmetrology height <y_base> <y_top> <y_horizon> <h_cam_um> <sigma_permil> -> a height, given scale
21// nx_svmetrology fuse <cuefile> [max_reduced_chi2_x100] -> posterior + chi2 + residuals
22// nx_svmetrology selftest
23// CUE ROWS (tab-separated; the cue registry plane uses the same shape):
24// cue <TAB> <name> <TAB> <value_um> <TAB> <sigma_um> [<TAB> provenance [<TAB> source]]
25// license_tier: ORIGINAL expect_exit: 0
26import "nx_syscalls.nx"
27const SV_MAGIC_2000: i64 = 2000
28const SV_MAGIC_1700000: i64 = 1700000
29
30const SV_CAP: i64 = 262144
31const SV_MAXCUE: i64 = 256
32const SV_WK: i64 = 1000000000000
33const SV_DEFCHI: i64 = 400
34const SV_TAB: i64 = 9
35const SV_NL: i64 = 10
36
37func hw(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
38func pn(v: i64) -> i64 { let b: *u8=sys_mmap(32) as *u8; var x: i64=v; var ng: i64=0; if x<0{ng=1;x=0-x} var i: i64=31; if x==0{b[i]=48 as u8;i=i-1} while x>0{b[i]=(48+x%10) as u8;x=x/10;i=i-1} if ng==1{b[i]=45 as u8;i=i-1} sys_write(1,(b as i64+i+1) as *u8,31-i); return 0 }
39func sv_pr(b: *u8, a: i64, e: i64) -> i64 { if e > a { sys_write(1, ((b as i64)+a) as *u8, e-a) } return 0 }
40func sv_isqrt(x: i64) -> i64 {
41 if x <= 0 { return 0 }
42 var r: i64 = x
43 var q: i64 = (x/2) + 1
44 while q < r { r = q; q = (r + x/r)/2 }
45 return r
46}
47func sv_refuse(reason: *u8) -> i64 { hw("SVMETROLOGY REFUSED: " as *u8); hw(reason); hw("\n" as *u8); return 0 }
48
49// ---- text helpers (flag-terminated scans: a break-by-sentinel that overwrites the cursor destroys
50// the position it exists to report -- banked law, re-broken twice in this estate) ----
51func sv_eol(b: *u8, n: i64, i: i64) -> i64 {
52 var e: i64 = i
53 var f: i64 = 0
54 while f == 0 { if e >= n { f = 1 } else { if b[e] == (SV_NL as u8) { f = 1 } else { e = e + 1 } } }
55 return e
56}
57func sv_cols(b: *u8, ls: i64, le: i64, col: *i64, maxc: i64) -> i64 {
58 var c: i64 = 0
59 var p: i64 = ls
60 while c < maxc {
61 var q: i64 = p
62 var f: i64 = 0
63 while f == 0 { if q >= le { f = 1 } else { if b[q] == (SV_TAB as u8) { f = 1 } else { q = q + 1 } } }
64 col[c*2] = p
65 col[c*2+1] = q
66 if q >= le { c = c + 1; while c < maxc { col[c*2] = le; col[c*2+1] = le; c = c + 1 } return c }
67 p = q + 1
68 c = c + 1
69 }
70 return c
71}
72func sv_int(b: *u8, a: i64, e: i64) -> i64 {
73 var v: i64 = 0
74 var neg: i64 = 0
75 var i: i64 = a
76 if i < e { if b[i] == (45 as u8) { neg = 1; i = i + 1 } }
77 while i < e {
78 let ch: i64 = b[i] as i64
79 if ch >= 48 { if ch <= 57 { v = v*10 + (ch-48) } }
80 i = i + 1
81 }
82 if neg == 1 { return 0 - v }
83 return v
84}
85func sv_same(b: *u8, a: i64, e: i64, lit: *u8) -> i64 {
86 var l: i64 = 0
87 while lit[l] != (0 as u8) { l = l + 1 }
88 if e - a != l { return 0 }
89 var i: i64 = 0
90 while i < l { if b[a+i] != lit[i] { return 0 } i = i + 1 }
91 return 1
92}
93func sv_argint(s: *u8) -> i64 {
94 var v: i64 = 0
95 var neg: i64 = 0
96 var i: i64 = 0
97 if s[0] == (45 as u8) { neg = 1; i = 1 }
98 while s[i] != (0 as u8) {
99 let ch: i64 = s[i] as i64
100 if ch >= 48 { if ch <= 57 { v = v*10 + (ch-48) } }
101 i = i + 1
102 }
103 if neg == 1 { return 0 - v }
104 return v
105}
106func sv_read(path: *u8, buf: *u8, cap: i64) -> i64 {
107 let fd: i64 = sys_openat_rd(path)
108 if fd < 0 { return 0 - 1 }
109 var n: i64 = 0
110 var go: i64 = 1
111 while go == 1 {
112 let r: i64 = sys_read(fd, ((buf as i64) + n) as *u8, cap - n)
113 if r <= 0 { go = 0 } else { n = n + r }
114 if n >= cap { go = 0 }
115 }
116 sys_close(fd)
117 return n
118}
119
120// ---- the cross-ratio: camera height from a reference of known height ----
121// h_obj/h_cam = (y_base - y_top)/(y_base - y_horizon) => h_cam = h_obj * (y_base - y_hor)/(y_base - y_top)
122func sv_camheight(yb: i64, yt: i64, yh: i64, hobj: i64, sigpm: i64) -> i64 {
123 let span: i64 = yb - yt
124 let lev: i64 = yb - yh
125 if span <= 0 { sv_refuse("y_top must be above y_base in image coordinates (span <= 0)" as *u8); return 3 }
126 if lev <= 0 { sv_refuse("reference base must be BELOW the horizon (leverage <= 0) -- an object at or above the horizon line is at/over camera height and carries no scale" as *u8); return 3 }
127 let hcam: i64 = hobj * lev / span
128 // uncertainty: the prior's relative sigma, plus a 2px read error amplified by 1/leverage
129 let pxrel: i64 = SV_MAGIC_2000 / lev
130 var rel: i64 = sigpm + pxrel
131 let sig: i64 = hcam * rel / 1000
132 hw("cue\tcamheight_from_ref\t" as *u8); pn(hcam)
133 hw("\t" as *u8); pn(sig)
134 hw("\tMEASURED\tcross-ratio horizon=" as *u8); pn(yh)
135 hw(" base=" as *u8); pn(yb)
136 hw(" top=" as *u8); pn(yt)
137 hw(" leverage_px=" as *u8); pn(lev)
138 hw(" prior_permil=" as *u8); pn(sigpm)
139 hw(" px_term_permil=" as *u8); pn(pxrel)
140 hw("\n" as *u8)
141 return 0
142}
143// object height given a known camera height
144func sv_height(yb: i64, yt: i64, yh: i64, hcam: i64, sigpm: i64) -> i64 {
145 let span: i64 = yb - yt
146 let lev: i64 = yb - yh
147 if span <= 0 { sv_refuse("span <= 0" as *u8); return 3 }
148 if lev <= 0 { sv_refuse("base at or above the horizon -- no scale" as *u8); return 3 }
149 let hobj: i64 = hcam * span / lev
150 let pxrel: i64 = SV_MAGIC_2000 / lev
151 let sig: i64 = hobj * (sigpm + pxrel) / 1000
152 hw("{\x22organ\x22:\x22nx_svmetrology\x22,\x22verb\x22:\x22height\x22,\x22height_um\x22:" as *u8); pn(hobj)
153 hw(",\x22sigma_um\x22:" as *u8); pn(sig)
154 hw(",\x22leverage_px\x22:" as *u8); pn(lev)
155 hw(",\x22note\x22:\x22error scales INVERSELY with leverage_px: prefer the NEAREST reference, not the clearest\x22}\n" as *u8)
156 return 0
157}
158
159// ---- fusion ----
160func sv_fuse(buf: *u8, n: i64, maxchi: i64) -> i64 {
161 let c0: *i64 = sys_mmap(SV_MAXCUE*8) as *i64
162 let c1: *i64 = sys_mmap(SV_MAXCUE*8) as *i64
163 let val: *i64 = sys_mmap(SV_MAXCUE*8) as *i64
164 let sig: *i64 = sys_mmap(SV_MAXCUE*8) as *i64
165 let col: *i64 = sys_mmap(12*2*8) as *i64
166 var nc: i64 = 0
167 var i: i64 = 0
168 while i < n {
169 let le: i64 = sv_eol(buf, n, i)
170 if le > i {
171 sv_cols(buf, i, le, col, 6)
172 if sv_same(buf, col[0], col[1], "cue" as *u8) == 1 {
173 if nc < SV_MAXCUE {
174 c0[nc] = col[2]
175 c1[nc] = col[3]
176 val[nc] = sv_int(buf, col[4], col[5])
177 var s: i64 = sv_int(buf, col[6], col[7])
178 if s <= 0 { s = 1 }
179 sig[nc] = s
180 nc = nc + 1
181 }
182 }
183 }
184 i = le + 1
185 }
186 if nc == 0 { sv_refuse("no cue rows -- refusing to invent a scale" as *u8); return 4 }
187 // inverse-variance fusion, weights scaled by SV_WK to stay integral
188 var wsum: i64 = 0
189 var wxsum: i64 = 0
190 var k: i64 = 0
191 while k < nc {
192 let w: i64 = SV_WK / (sig[k] * sig[k])
193 wsum = wsum + w
194 wxsum = wxsum + w * (val[k]/1000)
195 k = k + 1
196 }
197 if wsum <= 0 { sv_refuse("all cue sigmas overflow the weight scale (sigma too large to be informative)" as *u8); return 5 }
198 let post: i64 = (wxsum / wsum) * 1000
199 let postsig: i64 = sv_isqrt(SV_WK / wsum)
200 // chi-square over the residuals
201 var chi100: i64 = 0
202 var worst: i64 = 0 - 1
203 var worstr: i64 = 0
204 k = 0
205 while k < nc {
206 var d: i64 = val[k] - post
207 if d < 0 { d = 0 - d }
208 let r100: i64 = d * 100 / sig[k]
209 chi100 = chi100 + (r100 * r100) / 100
210 if r100 > worstr { worstr = r100; worst = k }
211 k = k + 1
212 }
213 let dof: i64 = nc - 1
214 hw("{\x22organ\x22:\x22nx_svmetrology\x22,\x22verb\x22:\x22fuse\x22,\x22cues\x22:" as *u8); pn(nc)
215 hw(",\x22posterior_um\x22:" as *u8); pn(post)
216 hw(",\x22sigma_um\x22:" as *u8); pn(postsig)
217 hw(",\x22ci95_lo_um\x22:" as *u8); pn(post - 2*postsig)
218 hw(",\x22ci95_hi_um\x22:" as *u8); pn(post + 2*postsig)
219 hw(",\x22dof\x22:" as *u8); pn(dof)
220 hw(",\x22chi2_x100\x22:" as *u8); pn(chi100)
221 hw(",\x22residuals\x22:[" as *u8)
222 k = 0
223 while k < nc {
224 if k > 0 { hw("," as *u8) }
225 hw("{\x22cue\x22:\x22" as *u8); sv_pr(buf, c0[k], c1[k])
226 hw("\x22,\x22value_um\x22:" as *u8); pn(val[k])
227 hw(",\x22sigma_um\x22:" as *u8); pn(sig[k])
228 var d2: i64 = val[k] - post
229 if d2 < 0 { d2 = 0 - d2 }
230 hw(",\x22residual_sigma_x100\x22:" as *u8); pn(d2 * 100 / sig[k])
231 hw("}" as *u8)
232 k = k + 1
233 }
234 hw("]" as *u8)
235 if dof == 0 {
236 hw(",\x22verdict\x22:\x22UNVALIDATED-SINGLE-CUE\x22,\x22why\x22:\x22one cue leaves ZERO degrees of freedom: nothing can check it. The interval is the prior's, not evidence of agreement -- acquire a second independent cue (debt 1785958362)\x22}\n" as *u8)
237 return 2
238 }
239 let red100: i64 = chi100 / dof
240 hw(",\x22reduced_chi2_x100\x22:" as *u8); pn(red100)
241 hw(",\x22threshold_x100\x22:" as *u8); pn(maxchi)
242 if red100 > maxchi {
243 hw(",\x22verdict\x22:\x22INCONSISTENT-REFUSED\x22,\x22outlier_cue\x22:\x22" as *u8); sv_pr(buf, c0[worst], c1[worst])
244 hw("\x22,\x22outlier_residual_sigma_x100\x22:" as *u8); pn(worstr)
245 // ★REFUSING IS NECESSARY BUT NOT SUFFICIENT: A DIAGNOSIS BEATS A VERDICT. Leave-one-out the named
246 // outlier and report what the evidence says WITHOUT it, so the caller learns whether ONE broken
247 // assumption explains the disagreement (repairable) or several do (the model itself is wrong).
248 if nc >= 3 {
249 var w2: i64 = 0
250 var wx2: i64 = 0
251 var j: i64 = 0
252 while j < nc {
253 if j != worst {
254 let ww: i64 = SV_WK / (sig[j] * sig[j])
255 w2 = w2 + ww
256 wx2 = wx2 + ww * (val[j]/1000)
257 }
258 j = j + 1
259 }
260 if w2 > 0 {
261 let post2: i64 = (wx2 / w2) * 1000
262 let psig2: i64 = sv_isqrt(SV_WK / w2)
263 var chi2b: i64 = 0
264 j = 0
265 while j < nc {
266 if j != worst {
267 var dd: i64 = val[j] - post2
268 if dd < 0 { dd = 0 - dd }
269 let rr: i64 = dd * 100 / sig[j]
270 chi2b = chi2b + (rr * rr) / 100
271 }
272 j = j + 1
273 }
274 let red2: i64 = chi2b / (nc - 2)
275 hw(",\x22repair\x22:{\x22dropped\x22:\x22" as *u8); sv_pr(buf, c0[worst], c1[worst])
276 hw("\x22,\x22posterior_um\x22:" as *u8); pn(post2)
277 hw(",\x22sigma_um\x22:" as *u8); pn(psig2)
278 hw(",\x22reduced_chi2_x100\x22:" as *u8); pn(red2)
279 if red2 <= maxchi { hw(",\x22diagnosis\x22:\x22ONE-BROKEN-ASSUMPTION: the remaining cues agree, so the named cue alone is at fault -- re-examine ITS assumption (posture, identity, flatness), then re-fuse\x22}" as *u8) } else { hw(",\x22diagnosis\x22:\x22MULTIPLE-BROKEN-ASSUMPTIONS: the cues still disagree after dropping the worst, so the MODEL is suspect (ground plane, horizon estimate or lens), not one reference\x22}" as *u8) }
280 }
281 }
282 hw(",\x22why\x22:\x22the cues disagree beyond their stated uncertainties, so an assumption is broken (misidentified reference, non-flat ground, lens distortion). Averaging would HIDE it.\x22}\n" as *u8)
283 return 1
284 }
285 hw(",\x22verdict\x22:\x22CONSISTENT\x22,\x22why\x22:\x22" as *u8); pn(dof)
286 hw(" degrees of freedom of validation passed: independent cues agree within their uncertainties\x22}\n" as *u8)
287 return 0
288}
289
290// ---- teeth ----
291func sv_wr(path: *u8, s: *u8) -> i64 {
292 var n: i64 = 0
293 while s[n] != (0 as u8) { n = n + 1 }
294 let fd: i64 = sys_openat_wr(path, 420)
295 if fd < 0 { return 0 - 1 }
296 sys_write(fd, s, n)
297 sys_close(fd)
298 return 0
299}
300func sv_runfile(path: *u8, maxchi: i64) -> i64 {
301 let b: *u8 = sys_mmap(SV_CAP)
302 let n: i64 = sv_read(path, b, SV_CAP)
303 if n <= 0 { sv_refuse("cue file unreadable or empty" as *u8); return 4 }
304 return sv_fuse(b, n, maxchi)
305}
306func sv_selftest() -> i64 {
307 var fails: i64 = 0
308 // T0 three agreeing cues around 850mm
309 sv_wr("/tmp/sv_ok.txt" as *u8, "cue\tfoot_ratio\t850000\t40000\tMEASURED\tfoot=0.152*stature\ncue\tfigure_population\t862000\t70000\tPRIOR\tadult height distribution\ncue\tfootprint_plane\t845000\t35000\tMEASURED\tground rectification\n" as *u8)
310 hw("T0 three agreeing cues -> CONSISTENT with a tightened interval:\n" as *u8)
311 if sv_runfile("/tmp/sv_ok.txt" as *u8, SV_DEFCHI) != 0 { fails = fails + 1; hw("T0 FAIL\n" as *u8) } else { hw("T0 PASS\n" as *u8) }
312 // T1 same three plus a gross outlier (the seated-vs-standing misread: 2x)
313 sv_wr("/tmp/sv_bad.txt" as *u8, "cue\tfoot_ratio\t850000\t40000\tMEASURED\tfoot=0.152*stature\ncue\tfigure_population\t862000\t70000\tPRIOR\tadult height distribution\ncue\tfootprint_plane\t845000\t35000\tMEASURED\tground rectification\ncue\tdistant_figure_STANDING\t1700000\t90000\tINFERRED\tposture assumed standing\n" as *u8)
314 hw("T1 add a 2x outlier (posture misread) -> must REFUSE and NAME it:\n" as *u8)
315 if sv_runfile("/tmp/sv_bad.txt" as *u8, SV_DEFCHI) != 1 { fails = fails + 1; hw("T1 FAIL outlier not refused\n" as *u8) } else { hw("T1 PASS refused + named\n" as *u8) }
316 // T2 single cue -> unvalidated
317 sv_wr("/tmp/sv_one.txt" as *u8, "cue\tfoot_ratio\t850000\t40000\tMEASURED\tsole reference\n" as *u8)
318 hw("T2 single cue -> UNVALIDATED (0 dof, nothing can check it):\n" as *u8)
319 if sv_runfile("/tmp/sv_one.txt" as *u8, SV_DEFCHI) != 2 { fails = fails + 1; hw("T2 FAIL single cue not flagged\n" as *u8) } else { hw("T2 PASS\n" as *u8) }
320 // T3 empty -> refuse
321 sv_wr("/tmp/sv_none.txt" as *u8, "# no cues here\n" as *u8)
322 hw("T3 no cues -> REFUSE (never invent a scale):\n" as *u8)
323 if sv_runfile("/tmp/sv_none.txt" as *u8, SV_DEFCHI) == 0 { fails = fails + 1; hw("T3 FAIL\n" as *u8) } else { hw("T3 PASS\n" as *u8) }
324 // T4 the leverage law: a NEAR reference must beat a distant one
325 hw("T4 leverage law -- same 2px read error, near vs distant reference:\n" as *u8)
326 hw(" near (base 900, horizon 655, leverage 245): " as *u8)
327 sv_camheight(900, 300, 655, SV_MAGIC_1700000, 40)
328 hw(" distant(base 680, horizon 655, leverage 25): " as *u8)
329 sv_camheight(680, 635, 655, SV_MAGIC_1700000, 40)
330 hw(" (the distant cue's px term is ~10x the near one -- that is the law, printed)\n" as *u8)
331 // T5 a reference at/above the horizon carries no scale
332 hw("T5 base above the horizon -> REFUSE:\n" as *u8)
333 if sv_camheight(600, 300, 655, SV_MAGIC_1700000, 40) == 0 { fails = fails + 1; hw("T5 FAIL\n" as *u8) } else { hw("T5 PASS\n" as *u8) }
334 if fails == 0 { hw("SVMETROLOGY-SELFTEST GREEN 5/5 (+ leverage demonstration)\n" as *u8); return 0 }
335 hw("SVMETROLOGY-SELFTEST RED fails=" as *u8); pn(fails); hw("\n" as *u8)
336 return 1
337}
338
339func main(argc: i64, argv: *i64) -> i64 {
340 if argc < 2 {
341 hw("usage: nx_svmetrology camheight <y_base> <y_top> <y_horizon> <h_obj_um> <sigma_permil>\n nx_svmetrology height <y_base> <y_top> <y_horizon> <h_cam_um> <sigma_permil>\n nx_svmetrology fuse <cuefile> [max_reduced_chi2_x100]\n nx_svmetrology selftest\n" as *u8)
342 sys_exit(2)
343 return 2
344 }
345 let a1: *u8 = argv[1] as *u8
346 if sv_same(a1, 0, 8, "selftest" as *u8) == 1 { let rc: i64 = sv_selftest(); sys_exit(rc); return rc }
347 if sv_same(a1, 0, 10, "camheight" as *u8) == 1 {
348 if argc < 7 { hw("usage: nx_svmetrology camheight <y_base> <y_top> <y_horizon> <h_obj_um> <sigma_permil>\n" as *u8); sys_exit(2); return 2 }
349 let rc2: i64 = sv_camheight(sv_argint(argv[2] as *u8), sv_argint(argv[3] as *u8), sv_argint(argv[4] as *u8), sv_argint(argv[5] as *u8), sv_argint(argv[6] as *u8))
350 sys_exit(rc2)
351 return rc2
352 }
353 if sv_same(a1, 0, 7, "height" as *u8) == 1 {
354 if argc < 7 { hw("usage: nx_svmetrology height <y_base> <y_top> <y_horizon> <h_cam_um> <sigma_permil>\n" as *u8); sys_exit(2); return 2 }
355 let rc3: i64 = sv_height(sv_argint(argv[2] as *u8), sv_argint(argv[3] as *u8), sv_argint(argv[4] as *u8), sv_argint(argv[5] as *u8), sv_argint(argv[6] as *u8))
356 sys_exit(rc3)
357 return rc3
358 }
359 if sv_same(a1, 0, 5, "fuse" as *u8) == 1 {
360 if argc < 3 { hw("usage: nx_svmetrology fuse <cuefile> [max_reduced_chi2_x100]\n" as *u8); sys_exit(2); return 2 }
361 var mc: i64 = SV_DEFCHI
362 if argc >= 4 { mc = sv_argint(argv[3] as *u8) }
363 if mc <= 0 { mc = SV_DEFCHI }
364 let rc4: i64 = sv_runfile(argv[2] as *u8, mc)
365 sys_exit(rc4)
366 return rc4
367 }
368 hw("unknown verb\n" as *u8)
369 sys_exit(2)
370 return 2
371}