nx_meshdist.nx source

1// nx_meshdist.nx -- THE METRIC GAP RULER (form ladder 1785935557; operator 2026-08-05: measure the gap 2// between two objects below the millimetre -- bodies today, screw-to-screw tomorrow). 3// Bidirectional point-to-surface distance between two NXMSH2 meshes, the METRO / DTU convention done 4// sovereignly: candidate->oracle = ACCURACY, oracle->candidate = COMPLETENESS (one direction alone is 5// gameable: a sphere inside the body has perfect accuracy and no completeness). Reported per direction: 6// mean / rms / p95 (the medical HD95) / max, in MICROMETRES, plus F-score-style fractions within 1mm and 7// 0.5mm (the Tanks-and-Temples F@tau shape). All integer math in 10um quanta -- products stay inside i64 8// by construction (coords <= ~2e5 q, near-cell evaluation bounds |p-a|), so the quantisation floor is 9// 10um, declared in the output. Screws and finer work pass a tighter unit via a future scale arg. 10// ALIGNMENT (declared, never silent): candidate is translated centroid-to-centroid and uniformly scaled 11// by the bbox-height ratio (permil printed). No rotation search in v1 -- both estates' bodies stand y-up. 12// Sampling v1: every k-th triangle CENTROID, k chosen to cap samples (bias declared; area-weighted 13// sampling is the named v2 rung). 14// nx_meshdist <cand.nxmesh> <oracle.nxmesh> [samples] | selftest 15// license_tier: ORIGINAL expect_exit: 0 16import "nx_syscalls.nx" 17 18const MD_CAP: i64 = 16777216 19const MD_MAXTRI: i64 = 400000 20const MD_GRID: i64 = 48 21const MD_MAXSAMP: i64 = 16384 22const MD_HISTMAX: i64 = 65536 23const MD_Q_PER_MM: i64 = 100 24const MD_UM_PER_Q: i64 = 10 25const MD_M8388607: i64 = 8388607 26const MD_M8388608: i64 = 8388608 27const MD_FARQ: i64 = 60000 28const MD_BIG: i64 = 4611686018427387903 29 30func hw(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 31func 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 } 32func md_u32(b: *u8, o: i64) -> i64 { return (b[o] as i64) + ((b[o+1] as i64)<<8) + ((b[o+2] as i64)<<16) + ((b[o+3] as i64)<<24) } 33// IEEE754 f32 (millimetres) -> integer 10um quanta (value*100) 34func md_f32q(w: i64) -> i64 { 35 let sign: i64 = (w >> 31) & 1 36 let expo: i64 = (w >> 23) & 255 37 if expo == 0 { return 0 } 38 var mant: i64 = (w & MD_M8388607) | MD_M8388608 39 let sh: i64 = expo - 127 40 var v: i64 = 0 41 if sh >= 23 { if sh - 23 > 30 { return 0 } } 42 if sh >= 23 { v = mant * MD_Q_PER_MM * (1 << (sh - 23)) } 43 if sh < 23 { if 23 - sh > 62 { return 0 } } 44 if sh < 23 { v = (mant * MD_Q_PER_MM) >> (23 - sh) } 45 if sign == 1 { return 0 - v } 46 return v 47} 48func md_isqrt(x: i64) -> i64 { 49 if x <= 0 { return 0 } 50 var r: i64 = x 51 var q: i64 = (x/2) + 1 52 while q < r { r = q; q = (r + x/r)/2 } 53 return r 54} 55func md_refuse(reason: *u8) -> i64 { hw("MESHDIST REFUSED: " as *u8); hw(reason); hw("\n" as *u8); return 0 } 56 57// load an NXMSH2 into vx (9 coords per tri, q10um). returns ntris or -1. 58func md_load(path: *u8, vx: *i64) -> i64 { 59 let fd: i64 = sys_openat_rd(path) 60 if fd < 0 { return 0 - 1 } 61 let b: *u8 = sys_mmap(MD_CAP + 64) 62 var n: i64 = 0 63 var go: i64 = 1 64 while go == 1 { 65 let r: i64 = sys_read(fd, ((b as i64) + n) as *u8, MD_CAP - n) 66 if r <= 0 { go = 0 } else { n = n + r } 67 if n >= MD_CAP { go = 0 } 68 } 69 sys_close(fd) 70 if n < 44 { return 0 - 1 } 71 if b[0] != (78 as u8) { return 0 - 1 } 72 if b[5] != (50 as u8) { return 0 - 1 } 73 let nlay: i64 = md_u32(b, 8) 74 let nt: i64 = md_u32(b, 12) 75 if nt <= 0 { return 0 - 1 } 76 if nt > MD_MAXTRI { return 0 - 1 } 77 let hdr: i64 = 16 + nlay*24 78 if hdr + nt*84 > n { return 0 - 1 } 79 var t: i64 = 0 80 while t < nt { 81 var c: i64 = 0 82 while c < 9 { vx[t*9+c] = md_f32q(md_u32(b, hdr + t*84 + c*4)); c = c + 1 } 83 t = t + 1 84 } 85 return nt 86} 87 88// squared distance point->triangle, all q10um (Ericson clamp; t in 1/16384) 89func md_pt(px: i64, py: i64, pz: i64, v: *i64, o: i64) -> i64 { 90 let ax: i64 = v[o] 91 let ay: i64 = v[o+1] 92 let az: i64 = v[o+2] 93 var abx: i64 = v[o+3]-ax 94 var aby: i64 = v[o+4]-ay 95 var abz: i64 = v[o+5]-az 96 var acx: i64 = v[o+6]-ax 97 var acy: i64 = v[o+7]-ay 98 var acz: i64 = v[o+8]-az 99 var apx: i64 = px-ax 100 var apy: i64 = py-ay 101 var apz: i64 = pz-az 102 // far guard: beyond ~6m component the exact solve could overflow -- return a conservative vertex bound 103 if apx > MD_FARQ { return MD_BIG/4 } 104 if apx < 0-MD_FARQ { return MD_BIG/4 } 105 if apy > MD_FARQ { return MD_BIG/4 } 106 if apy < 0-MD_FARQ { return MD_BIG/4 } 107 if apz > MD_FARQ { return MD_BIG/4 } 108 if apz < 0-MD_FARQ { return MD_BIG/4 } 109 let d1: i64 = abx*apx + aby*apy + abz*apz 110 let d2: i64 = acx*apx + acy*apy + acz*apz 111 var cx: i64 = 0 112 var cy: i64 = 0 113 var cz: i64 = 0 114 var done: i64 = 0 115 if d1 <= 0 { if d2 <= 0 { cx = ax; cy = ay; cz = az; done = 1 } } 116 let bpx: i64 = px-v[o+3] 117 let bpy: i64 = py-v[o+4] 118 let bpz: i64 = pz-v[o+5] 119 let d3: i64 = abx*bpx + aby*bpy + abz*bpz 120 let d4: i64 = acx*bpx + acy*bpy + acz*bpz 121 if done == 0 { if d3 >= 0 { if d4 <= d3 { cx = v[o+3]; cy = v[o+4]; cz = v[o+5]; done = 1 } } } 122 if done == 0 { 123 let vc: i64 = d1*d4 - d3*d2 124 if vc <= 0 { if d1 >= 0 { if d3 <= 0 { 125 var den: i64 = d1 - d3 126 if den == 0 { den = 1 } 127 let t1: i64 = d1*16384/den 128 cx = ax + abx*t1/16384; cy = ay + aby*t1/16384; cz = az + abz*t1/16384; done = 1 129 } } } 130 } 131 let cpx: i64 = px-v[o+6] 132 let cpy: i64 = py-v[o+7] 133 let cpz: i64 = pz-v[o+8] 134 let d5: i64 = abx*cpx + aby*cpy + abz*cpz 135 let d6: i64 = acx*cpx + acy*cpy + acz*cpz 136 if done == 0 { if d6 >= 0 { if d5 <= d6 { cx = v[o+6]; cy = v[o+7]; cz = v[o+8]; done = 1 } } } 137 if done == 0 { 138 let vb: i64 = d5*d2 - d1*d6 139 if vb <= 0 { if d2 >= 0 { if d6 <= 0 { 140 var den2: i64 = d2 - d6 141 if den2 == 0 { den2 = 1 } 142 let t2: i64 = d2*16384/den2 143 cx = ax + acx*t2/16384; cy = ay + acy*t2/16384; cz = az + acz*t2/16384; done = 1 144 } } } 145 } 146 if done == 0 { 147 let va: i64 = d3*d6 - d5*d4 148 if va <= 0 { if d4 - d3 >= 0 { if d5 - d6 >= 0 { 149 var den3: i64 = (d4-d3) + (d5-d6) 150 if den3 == 0 { den3 = 1 } 151 let t3: i64 = (d4-d3)*16384/den3 152 cx = v[o+3] + (v[o+6]-v[o+3])*t3/16384 153 cy = v[o+4] + (v[o+7]-v[o+4])*t3/16384 154 cz = v[o+5] + (v[o+8]-v[o+5])*t3/16384 155 done = 1 156 } } } 157 if done == 0 { 158 // range-reduce before the fixed-point: on a 100mm tri va/vb/vc ~ 3e15 and *16384 overflows i64 159 // (caught by the selftest fixture 2026-08-05); shifting all three together preserves the ratio 160 var vbn: i64 = d5*d2 - d1*d6 161 var vcn: i64 = d1*d4 - d3*d2 162 var den4: i64 = va + vbn + vcn 163 if den4 == 0 { den4 = 1 } 164 while den4 > 549755813888 { vbn = vbn/1024; vcn = vcn/1024; den4 = den4/1024; if den4 == 0 { den4 = 1 } } 165 while den4 < 0-549755813888 { vbn = vbn/1024; vcn = vcn/1024; den4 = den4/1024; if den4 == 0 { den4 = 1 } } 166 let vv: i64 = vbn*16384/den4 167 let ww: i64 = vcn*16384/den4 168 cx = ax + abx*vv/16384 + acx*ww/16384 169 cy = ay + aby*vv/16384 + acy*ww/16384 170 cz = az + abz*vv/16384 + acz*ww/16384 171 done = 1 172 } 173 } 174 let dx: i64 = px-cx 175 let dy: i64 = py-cy 176 let dz: i64 = pz-cz 177 return dx*dx + dy*dy + dz*dz 178} 179 180// ---- grid over the target mesh ---- 181static MD_GX0: i64 182static MD_GY0: i64 183static MD_GZ0: i64 184static MD_CS: i64 185func md_cell(v: i64, g0: i64) -> i64 { 186 var c: i64 = (v - g0)/MD_CS 187 if c < 0 { c = 0 } 188 if c >= MD_GRID { c = MD_GRID-1 } 189 return c 190} 191// build cell lists: cnt/start (MD_GRID^3) + refs. returns nrefs. 192func md_grid(v: *i64, nt: i64, start: *i64, refs: *i64, maxrefs: i64) -> i64 { 193 var mnx: i64 = MD_BIG 194 var mny: i64 = MD_BIG 195 var mnz: i64 = MD_BIG 196 var mxx: i64 = 0-MD_BIG 197 var mxy: i64 = 0-MD_BIG 198 var mxz: i64 = 0-MD_BIG 199 var t: i64 = 0 200 while t < nt { 201 var c: i64 = 0 202 while c < 3 { 203 let x: i64 = v[t*9+c*3] 204 let y: i64 = v[t*9+c*3+1] 205 let z: i64 = v[t*9+c*3+2] 206 if x < mnx { mnx = x } 207 if x > mxx { mxx = x } 208 if y < mny { mny = y } 209 if y > mxy { mxy = y } 210 if z < mnz { mnz = z } 211 if z > mxz { mxz = z } 212 c = c + 1 213 } 214 t = t + 1 215 } 216 var ext: i64 = mxx - mnx 217 if mxy - mny > ext { ext = mxy - mny } 218 if mxz - mnz > ext { ext = mxz - mnz } 219 MD_CS = ext/MD_GRID + 1 220 MD_GX0 = mnx 221 MD_GY0 = mny 222 MD_GZ0 = mnz 223 let nc: i64 = MD_GRID*MD_GRID*MD_GRID 224 var i: i64 = 0 225 while i <= nc { start[i] = 0; i = i + 1 } 226 // pass 1: counts (into start[c+1]) 227 t = 0 228 while t < nt { 229 // full AABB of the tri, in cells 230 var ax2: i64 = md_cell(v[t*9], MD_GX0) 231 var bx2: i64 = ax2 232 var ay2: i64 = md_cell(v[t*9+1], MD_GY0) 233 var by2: i64 = ay2 234 var az2: i64 = md_cell(v[t*9+2], MD_GZ0) 235 var bz2: i64 = az2 236 var c2: i64 = 1 237 while c2 < 3 { 238 let qx: i64 = md_cell(v[t*9+c2*3], MD_GX0) 239 let qy: i64 = md_cell(v[t*9+c2*3+1], MD_GY0) 240 let qz: i64 = md_cell(v[t*9+c2*3+2], MD_GZ0) 241 if qx < ax2 { ax2 = qx } 242 if qx > bx2 { bx2 = qx } 243 if qy < ay2 { ay2 = qy } 244 if qy > by2 { by2 = qy } 245 if qz < az2 { az2 = qz } 246 if qz > bz2 { bz2 = qz } 247 c2 = c2 + 1 248 } 249 var gx: i64 = ax2 250 while gx <= bx2 { 251 var gy: i64 = ay2 252 while gy <= by2 { 253 var gz: i64 = az2 254 while gz <= bz2 { 255 let cell: i64 = (gx*MD_GRID + gy)*MD_GRID + gz 256 start[cell+1] = start[cell+1] + 1 257 gz = gz + 1 258 } 259 gy = gy + 1 260 } 261 gx = gx + 1 262 } 263 t = t + 1 264 } 265 // prefix 266 i = 0 267 while i < nc { start[i+1] = start[i+1] + start[i]; i = i + 1 } 268 let nrefs: i64 = start[nc] 269 if nrefs > maxrefs { return 0 - 1 } 270 // pass 2: fill (cursor array) 271 let cur: *i64 = sys_mmap((nc+1)*8) as *i64 272 i = 0 273 while i <= nc { cur[i] = start[i]; i = i + 1 } 274 t = 0 275 while t < nt { 276 var ax3: i64 = md_cell(v[t*9], MD_GX0) 277 var bx3: i64 = ax3 278 var ay3: i64 = md_cell(v[t*9+1], MD_GY0) 279 var by3: i64 = ay3 280 var az3: i64 = md_cell(v[t*9+2], MD_GZ0) 281 var bz3: i64 = az3 282 var c3: i64 = 1 283 while c3 < 3 { 284 let rx: i64 = md_cell(v[t*9+c3*3], MD_GX0) 285 let ry: i64 = md_cell(v[t*9+c3*3+1], MD_GY0) 286 let rz: i64 = md_cell(v[t*9+c3*3+2], MD_GZ0) 287 if rx < ax3 { ax3 = rx } 288 if rx > bx3 { bx3 = rx } 289 if ry < ay3 { ay3 = ry } 290 if ry > by3 { by3 = ry } 291 if rz < az3 { az3 = rz } 292 if rz > bz3 { bz3 = rz } 293 c3 = c3 + 1 294 } 295 var hx: i64 = ax3 296 while hx <= bx3 { 297 var hy: i64 = ay3 298 while hy <= by3 { 299 var hz: i64 = az3 300 while hz <= bz3 { 301 let cell2: i64 = (hx*MD_GRID + hy)*MD_GRID + hz 302 refs[cur[cell2]] = t 303 cur[cell2] = cur[cell2] + 1 304 hz = hz + 1 305 } 306 hy = hy + 1 307 } 308 hx = hx + 1 309 } 310 t = t + 1 311 } 312 return nrefs 313} 314 315// nearest squared distance from p to the gridded mesh (expanding shells) 316func md_near(px: i64, py: i64, pz: i64, v: *i64, start: *i64, refs: *i64) -> i64 { 317 let cx: i64 = md_cell(px, MD_GX0) 318 let cy: i64 = md_cell(py, MD_GY0) 319 let cz: i64 = md_cell(pz, MD_GZ0) 320 var best: i64 = MD_BIG 321 var r: i64 = 0 322 var go: i64 = 1 323 while go == 1 { 324 // if the closest possible point in shell r is farther than best -> stop 325 if r > 0 { 326 let ring: i64 = (r-1)*MD_CS 327 if ring*ring > best { go = 0 } 328 } 329 if r >= MD_GRID { go = 0 } 330 if go == 1 { 331 var gx: i64 = cx-r 332 while gx <= cx+r { 333 var gy: i64 = cy-r 334 while gy <= cy+r { 335 var gz: i64 = cz-r 336 while gz <= cz+r { 337 var onshell: i64 = 0 338 if gx == cx-r { onshell = 1 } 339 if gx == cx+r { onshell = 1 } 340 if gy == cy-r { onshell = 1 } 341 if gy == cy+r { onshell = 1 } 342 if gz == cz-r { onshell = 1 } 343 if gz == cz+r { onshell = 1 } 344 if onshell == 1 { if gx >= 0 { if gx < MD_GRID { if gy >= 0 { if gy < MD_GRID { if gz >= 0 { if gz < MD_GRID { 345 let cell: i64 = (gx*MD_GRID + gy)*MD_GRID + gz 346 var k: i64 = start[cell] 347 while k < start[cell+1] { 348 let d2v: i64 = md_pt(px, py, pz, v, refs[k]*9) 349 if d2v < best { best = d2v } 350 k = k + 1 351 } 352 } } } } } } } 353 gz = gz + 1 354 } 355 gy = gy + 1 356 } 357 gx = gx + 1 358 } 359 r = r + 1 360 } 361 } 362 return best 363} 364 365// one direction: samples from src tri centroids -> nearest on dst. prints a JSON fragment. 366func md_dir(label: *u8, src: *i64, ns: i64, dst: *i64, ndst: i64, start: *i64, refs: *i64, want: i64) -> i64 { 367 var stride: i64 = ns/want 368 if stride < 1 { stride = 1 } 369 let hist: *i64 = sys_mmap(MD_HISTMAX*8) as *i64 370 var i: i64 = 0 371 while i < MD_HISTMAX { hist[i] = 0; i = i + 1 } 372 var cnt: i64 = 0 373 var sum: i64 = 0 374 var sum2: i64 = 0 375 var mx: i64 = 0 376 var in1000: i64 = 0 377 var in500: i64 = 0 378 var t: i64 = 0 379 while t < ns { 380 let o: i64 = t*9 381 let px: i64 = (src[o] + src[o+3] + src[o+6])/3 382 let py: i64 = (src[o+1] + src[o+4] + src[o+7])/3 383 let pz: i64 = (src[o+2] + src[o+5] + src[o+8])/3 384 let d2v: i64 = md_near(px, py, pz, dst, start, refs) 385 let dum: i64 = md_isqrt(d2v) * MD_UM_PER_Q 386 cnt = cnt + 1 387 sum = sum + dum 388 sum2 = sum2 + (dum/10)*(dum/10) 389 if dum > mx { mx = dum } 390 if dum <= 1000 { in1000 = in1000 + 1 } 391 if dum <= 500 { in500 = in500 + 1 } 392 var hb: i64 = dum/MD_UM_PER_Q 393 if hb >= MD_HISTMAX { hb = MD_HISTMAX-1 } 394 hist[hb] = hist[hb] + 1 395 t = t + stride 396 } 397 if cnt == 0 { return 0 - 1 } 398 // p95 from histogram 399 var acc: i64 = 0 400 var p95: i64 = 0 401 var hb2: i64 = 0 402 while hb2 < MD_HISTMAX { 403 acc = acc + hist[hb2] 404 if p95 == 0 { if acc*100 >= cnt*95 { p95 = hb2*MD_UM_PER_Q } } 405 hb2 = hb2 + 1 406 } 407 hw("\x22" as *u8); hw(label); hw("\x22:{\x22samples\x22:" as *u8); pn(cnt) 408 hw(",\x22mean_um\x22:" as *u8); pn(sum/cnt) 409 hw(",\x22rms_um\x22:" as *u8); pn(md_isqrt(sum2/cnt)*10) 410 hw(",\x22p95_um\x22:" as *u8); pn(p95) 411 hw(",\x22max_um\x22:" as *u8); pn(mx) 412 hw(",\x22within_1mm_permil\x22:" as *u8); pn(in1000*1000/cnt) 413 hw(",\x22within_0p5mm_permil\x22:" as *u8); pn(in500*1000/cnt) 414 hw("}" as *u8) 415 return 0 416} 417 418func md_run(candp: *u8, oracp: *u8, want: i64) -> i64 { 419 let cv: *i64 = sys_mmap(MD_MAXTRI*9*8 + 64) as *i64 420 let ov: *i64 = sys_mmap(MD_MAXTRI*9*8 + 64) as *i64 421 let nc: i64 = md_load(candp, cv) 422 if nc < 0 { md_refuse("candidate unreadable or not NXMSH2" as *u8); return 3 } 423 let no: i64 = md_load(oracp, ov) 424 if no < 0 { md_refuse("oracle unreadable or not NXMSH2" as *u8); return 3 } 425 // centroid align + height scale (candidate -> oracle frame), declared 426 var csx: i64 = 0 427 var csy: i64 = 0 428 var csz: i64 = 0 429 var cmny: i64 = MD_BIG 430 var cmxy: i64 = 0-MD_BIG 431 var t: i64 = 0 432 while t < nc { 433 csx = csx + cv[t*9] 434 csy = csy + cv[t*9+1] 435 csz = csz + cv[t*9+2] 436 if cv[t*9+1] < cmny { cmny = cv[t*9+1] } 437 if cv[t*9+1] > cmxy { cmxy = cv[t*9+1] } 438 t = t + 1 439 } 440 csx = csx/nc 441 csy = csy/nc 442 csz = csz/nc 443 var osx: i64 = 0 444 var osy: i64 = 0 445 var osz: i64 = 0 446 var omny: i64 = MD_BIG 447 var omxy: i64 = 0-MD_BIG 448 t = 0 449 while t < no { 450 osx = osx + ov[t*9] 451 osy = osy + ov[t*9+1] 452 osz = osz + ov[t*9+2] 453 if ov[t*9+1] < omny { omny = ov[t*9+1] } 454 if ov[t*9+1] > omxy { omxy = ov[t*9+1] } 455 t = t + 1 456 } 457 osx = osx/no 458 osy = osy/no 459 osz = osz/no 460 var ch: i64 = cmxy - cmny 461 if ch < 1 { ch = 1 } 462 var oh: i64 = omxy - omny 463 if oh < 1 { oh = 1 } 464 let scl: i64 = oh*1000/ch 465 t = 0 466 while t < nc { 467 var c: i64 = 0 468 while c < 3 { 469 cv[t*9+c*3] = (cv[t*9+c*3]-csx)*scl/1000 + osx 470 cv[t*9+c*3+1] = (cv[t*9+c*3+1]-csy)*scl/1000 + osy 471 cv[t*9+c*3+2] = (cv[t*9+c*3+2]-csz)*scl/1000 + osz 472 c = c + 1 473 } 474 t = t + 1 475 } 476 let ncell: i64 = MD_GRID*MD_GRID*MD_GRID 477 let startO: *i64 = sys_mmap((ncell+2)*8) as *i64 478 let refsO: *i64 = sys_mmap(MD_MAXTRI*24*8) as *i64 479 if md_grid(ov, no, startO, refsO, MD_MAXTRI*24) < 0 { md_refuse("oracle grid overflow" as *u8); return 4 } 480 hw("{\x22organ\x22:\x22nx_meshdist\x22,\x22quant_um\x22:10,\x22align\x22:\x22centroid+height (no rotation, v1)\x22,\x22scale_permil\x22:" as *u8); pn(scl) 481 hw(",\x22cand_tris\x22:" as *u8); pn(nc) 482 hw(",\x22oracle_tris\x22:" as *u8); pn(no) 483 hw("," as *u8) 484 if md_dir("accuracy_cand_to_oracle" as *u8, cv, nc, ov, no, startO, refsO, want) < 0 { md_refuse("no samples" as *u8); return 5 } 485 hw("," as *u8) 486 let startC: *i64 = sys_mmap((ncell+2)*8) as *i64 487 let refsC: *i64 = sys_mmap(MD_MAXTRI*24*8) as *i64 488 if md_grid(cv, nc, startC, refsC, MD_MAXTRI*24) < 0 { md_refuse("cand grid overflow" as *u8); return 4 } 489 if md_dir("completeness_oracle_to_cand" as *u8, ov, no, cv, nc, startC, refsC, want) < 0 { md_refuse("no samples" as *u8); return 5 } 490 hw(",\x22note\x22:\x22METRO/DTU convention: accuracy+completeness both required; p95=HD95; centroid sampling v1 (area-weighted = named rung)\x22}\n" as *u8) 491 return 0 492} 493 494// integer/scale -> f32 bits (the estate encoder pattern; colours are PER-MILLE by NXMSH2 convention) 495func md_enc(v: i64, scale: i64) -> i64 { 496 if v == 0 { return 0 } 497 var neg: i64 = 0 498 var m: i64 = v 499 if m < 0 { neg = 1; m = 0-m } 500 var e: i64 = 0 501 var num: i64 = m 502 var den: i64 = scale 503 while num >= den*2 { den = den*2; e = e+1 } 504 while num < den { num = num*2; e = e-1 } 505 let frac: i64 = ((num - den)*MD_M8388608)/den 506 var bits: i64 = ((e+127) << 23) | (frac & MD_M8388607) 507 if neg == 1 { bits = bits | (1<<31) } 508 return bits 509} 510func md_w32(b: *u8, o: i64, v: i64) -> i64 { 511 b[o]=(v&255) as u8; b[o+1]=((v>>8)&255) as u8; b[o+2]=((v>>16)&255) as u8; b[o+3]=((v>>24)&255) as u8 512 return 0 513} 514// ---- DEVIATION HEATMAP (the ZEISS-Inspect capability, sovereign; operator research drop 2026-08-05): 515// paint each CANDIDATE triangle by its centroid distance to the oracle, banded green->magenta, into the 516// candidate file itself -- nx_meshview already renders per-tri colour, so the picture ships through the 517// existing viewer untouched. Alignment is used for the DISTANCE only; colours land on the ORIGINAL bytes. 518func md_paint(candp: *u8, oracp: *u8, outp: *u8, band_um: i64) -> i64 { 519 if band_um < 100 { md_refuse("band under 100um is below the quantisation floor -- a sub-quantum band lies" as *u8); return 2 } 520 let fd: i64 = sys_openat_rd(candp) 521 if fd < 0 { md_refuse("candidate unreadable" as *u8); return 3 } 522 let raw: *u8 = sys_mmap(MD_CAP + 64) 523 var n: i64 = 0 524 var go: i64 = 1 525 while go == 1 { 526 let r: i64 = sys_read(fd, ((raw as i64) + n) as *u8, MD_CAP - n) 527 if r <= 0 { go = 0 } else { n = n + r } 528 if n >= MD_CAP { go = 0 } 529 } 530 sys_close(fd) 531 if n < 44 { md_refuse("candidate too small" as *u8); return 3 } 532 if raw[0] != (78 as u8) { md_refuse("candidate not NXMSH2" as *u8); return 3 } 533 let nlay: i64 = md_u32(raw, 8) 534 let nc: i64 = md_u32(raw, 12) 535 if nc <= 0 { md_refuse("candidate empty" as *u8); return 3 } 536 if nc > MD_MAXTRI { md_refuse("candidate over tri cap" as *u8); return 3 } 537 let hdr: i64 = 16 + nlay*24 538 let cv: *i64 = sys_mmap(MD_MAXTRI*9*8 + 64) as *i64 539 var t: i64 = 0 540 while t < nc { 541 var c: i64 = 0 542 while c < 9 { cv[t*9+c] = md_f32q(md_u32(raw, hdr + t*84 + c*4)); c = c + 1 } 543 t = t + 1 544 } 545 let ov: *i64 = sys_mmap(MD_MAXTRI*9*8 + 64) as *i64 546 let no: i64 = md_load(oracp, ov) 547 if no < 0 { md_refuse("oracle unreadable or not NXMSH2" as *u8); return 3 } 548 // align candidate coords to the oracle frame (distance only) 549 var csx: i64 = 0 550 var csy: i64 = 0 551 var csz: i64 = 0 552 var cmny: i64 = MD_BIG 553 var cmxy: i64 = 0-MD_BIG 554 t = 0 555 while t < nc { 556 csx = csx + cv[t*9] 557 csy = csy + cv[t*9+1] 558 csz = csz + cv[t*9+2] 559 if cv[t*9+1] < cmny { cmny = cv[t*9+1] } 560 if cv[t*9+1] > cmxy { cmxy = cv[t*9+1] } 561 t = t + 1 562 } 563 csx = csx/nc 564 csy = csy/nc 565 csz = csz/nc 566 var osx: i64 = 0 567 var osy: i64 = 0 568 var osz: i64 = 0 569 var omny: i64 = MD_BIG 570 var omxy: i64 = 0-MD_BIG 571 t = 0 572 while t < no { 573 osx = osx + ov[t*9] 574 osy = osy + ov[t*9+1] 575 osz = osz + ov[t*9+2] 576 if ov[t*9+1] < omny { omny = ov[t*9+1] } 577 if ov[t*9+1] > omxy { omxy = ov[t*9+1] } 578 t = t + 1 579 } 580 osx = osx/no 581 osy = osy/no 582 osz = osz/no 583 var ch: i64 = cmxy - cmny 584 if ch < 1 { ch = 1 } 585 var oh: i64 = omxy - omny 586 if oh < 1 { oh = 1 } 587 let scl: i64 = oh*1000/ch 588 t = 0 589 while t < nc { 590 var c2: i64 = 0 591 while c2 < 3 { 592 cv[t*9+c2*3] = (cv[t*9+c2*3]-csx)*scl/1000 + osx 593 cv[t*9+c2*3+1] = (cv[t*9+c2*3+1]-csy)*scl/1000 + osy 594 cv[t*9+c2*3+2] = (cv[t*9+c2*3+2]-csz)*scl/1000 + osz 595 c2 = c2 + 1 596 } 597 t = t + 1 598 } 599 let ncell: i64 = MD_GRID*MD_GRID*MD_GRID 600 let startO: *i64 = sys_mmap((ncell+2)*8) as *i64 601 let refsO: *i64 = sys_mmap(MD_MAXTRI*24*8) as *i64 602 if md_grid(ov, no, startO, refsO, MD_MAXTRI*24) < 0 { md_refuse("oracle grid overflow" as *u8); return 4 } 603 let hist: *i64 = sys_mmap(8*8) as *i64 604 var hb: i64 = 0 605 while hb < 5 { hist[hb] = 0; hb = hb + 1 } 606 t = 0 607 while t < nc { 608 let o9: i64 = t*9 609 let px: i64 = (cv[o9] + cv[o9+3] + cv[o9+6])/3 610 let py: i64 = (cv[o9+1] + cv[o9+4] + cv[o9+7])/3 611 let pz: i64 = (cv[o9+2] + cv[o9+5] + cv[o9+8])/3 612 let d2v: i64 = md_near(px, py, pz, ov, startO, refsO) 613 let dum: i64 = md_isqrt(d2v) * MD_UM_PER_Q 614 var band: i64 = dum/band_um 615 if band > 4 { band = 4 } 616 hist[band] = hist[band] + 1 617 var r2: i64 = 100 618 var g2: i64 = 850 619 var b2: i64 = 150 620 if band == 1 { r2 = 900; g2 = 880; b2 = 100 } 621 if band == 2 { r2 = 950; g2 = 550; b2 = 80 } 622 if band == 3 { r2 = 920; g2 = 120; b2 = 80 } 623 if band == 4 { r2 = 850; g2 = 100; b2 = 850 } 624 let co: i64 = hdr + t*84 625 md_w32(raw, co+72, md_enc(r2, 1000)) 626 md_w32(raw, co+76, md_enc(g2, 1000)) 627 md_w32(raw, co+80, md_enc(b2, 1000)) 628 t = t + 1 629 } 630 let ofd: i64 = sys_openat_wr(outp, 420) 631 if ofd < 0 { md_refuse("output unwritable" as *u8); return 6 } 632 sys_write(ofd, raw, n) 633 sys_close(ofd) 634 hw("{\x22organ\x22:\x22nx_meshdist\x22,\x22verb\x22:\x22paint\x22,\x22band_um\x22:" as *u8); pn(band_um) 635 hw(",\x22tris\x22:" as *u8); pn(nc) 636 hw(",\x22bands\x22:{\x22b0_green\x22:" as *u8); pn(hist[0]) 637 hw(",\x22b1_yellow\x22:" as *u8); pn(hist[1]) 638 hw(",\x22b2_orange\x22:" as *u8); pn(hist[2]) 639 hw(",\x22b3_red\x22:" as *u8); pn(hist[3]) 640 hw(",\x22b4_magenta\x22:" as *u8); pn(hist[4]) 641 hw("},\x22permil\x22:[" as *u8); pn(hist[0]*1000/nc) 642 hw("," as *u8); pn(hist[1]*1000/nc) 643 hw("," as *u8); pn(hist[2]*1000/nc) 644 hw("," as *u8); pn(hist[3]*1000/nc) 645 hw("," as *u8); pn(hist[4]*1000/nc) 646 hw("],\x22note\x22:\x22ZEISS-Inspect-class deviation heatmap: candidate painted by distance-to-oracle; the histogram IS the picture in numbers so they can never disagree\x22}\n" as *u8) 647 return 0 648} 649 650// ---- PER-PART TABLE (debt 1785968584): oracle layers are the donor's own named parts, so 651// completeness can be reported PER PART -- 'her left forearm vs ours: N mm'. One candidate grid, 652// then each layer's contiguous tri range is sampled independently. 653func md_parts(candp: *u8, oracp: *u8, want: i64) -> i64 { 654 let cv: *i64 = sys_mmap(MD_MAXTRI*9*8 + 64) as *i64 655 let ov: *i64 = sys_mmap(MD_MAXTRI*9*8 + 64) as *i64 656 let nc: i64 = md_load(candp, cv) 657 if nc < 0 { md_refuse("candidate unreadable or not NXMSH2" as *u8); return 3 } 658 let fd: i64 = sys_openat_rd(oracp) 659 if fd < 0 { md_refuse("oracle unreadable" as *u8); return 3 } 660 let b: *u8 = sys_mmap(MD_CAP + 64) 661 var n: i64 = 0 662 var go: i64 = 1 663 while go == 1 { 664 let r: i64 = sys_read(fd, ((b as i64) + n) as *u8, MD_CAP - n) 665 if r <= 0 { go = 0 } else { n = n + r } 666 if n >= MD_CAP { go = 0 } 667 } 668 sys_close(fd) 669 if n < 44 { md_refuse("oracle too small" as *u8); return 3 } 670 if b[0] != (78 as u8) { md_refuse("oracle not NXMSH2" as *u8); return 3 } 671 let nlay: i64 = md_u32(b, 8) 672 let no: i64 = md_u32(b, 12) 673 if no <= 0 { md_refuse("oracle empty" as *u8); return 3 } 674 if no > MD_MAXTRI { md_refuse("oracle over tri cap" as *u8); return 3 } 675 if nlay < 2 { md_refuse("oracle carries no named parts (single layer) -- run the converter with skinning" as *u8); return 6 } 676 let hdr: i64 = 16 + nlay*24 677 var t: i64 = 0 678 while t < no { 679 var c: i64 = 0 680 while c < 9 { ov[t*9+c] = md_f32q(md_u32(b, hdr + t*84 + c*4)); c = c + 1 } 681 t = t + 1 682 } 683 var csx: i64 = 0 684 var csy: i64 = 0 685 var csz: i64 = 0 686 var cmny: i64 = MD_BIG 687 var cmxy: i64 = 0-MD_BIG 688 t = 0 689 while t < nc { 690 csx = csx + cv[t*9] 691 csy = csy + cv[t*9+1] 692 csz = csz + cv[t*9+2] 693 if cv[t*9+1] < cmny { cmny = cv[t*9+1] } 694 if cv[t*9+1] > cmxy { cmxy = cv[t*9+1] } 695 t = t + 1 696 } 697 csx = csx/nc 698 csy = csy/nc 699 csz = csz/nc 700 var osx: i64 = 0 701 var osy: i64 = 0 702 var osz: i64 = 0 703 var omny: i64 = MD_BIG 704 var omxy: i64 = 0-MD_BIG 705 t = 0 706 while t < no { 707 osx = osx + ov[t*9] 708 osy = osy + ov[t*9+1] 709 osz = osz + ov[t*9+2] 710 if ov[t*9+1] < omny { omny = ov[t*9+1] } 711 if ov[t*9+1] > omxy { omxy = ov[t*9+1] } 712 t = t + 1 713 } 714 osx = osx/no 715 osy = osy/no 716 osz = osz/no 717 var ch: i64 = cmxy - cmny 718 if ch < 1 { ch = 1 } 719 var oh: i64 = omxy - omny 720 if oh < 1 { oh = 1 } 721 let scl: i64 = oh*1000/ch 722 t = 0 723 while t < nc { 724 var c2: i64 = 0 725 while c2 < 3 { 726 cv[t*9+c2*3] = (cv[t*9+c2*3]-csx)*scl/1000 + osx 727 cv[t*9+c2*3+1] = (cv[t*9+c2*3+1]-csy)*scl/1000 + osy 728 cv[t*9+c2*3+2] = (cv[t*9+c2*3+2]-csz)*scl/1000 + osz 729 c2 = c2 + 1 730 } 731 t = t + 1 732 } 733 let ncell: i64 = MD_GRID*MD_GRID*MD_GRID 734 let startC: *i64 = sys_mmap((ncell+2)*8) as *i64 735 let refsC: *i64 = sys_mmap(MD_MAXTRI*24*8) as *i64 736 if md_grid(cv, nc, startC, refsC, MD_MAXTRI*24) < 0 { md_refuse("cand grid overflow" as *u8); return 4 } 737 hw("{\x22organ\x22:\x22nx_meshdist\x22,\x22verb\x22:\x22parts\x22,\x22quant_um\x22:10,\x22scale_permil\x22:" as *u8); pn(scl) 738 hw(",\x22oracle_layers\x22:" as *u8); pn(nlay) 739 hw(",\x22parts\x22:[" as *u8) 740 var first: i64 = 1 741 var L: i64 = 0 742 while L < nlay { 743 let lb: i64 = 16 + L*24 744 let off: i64 = md_u32(b, lb+16) 745 let cnt: i64 = md_u32(b, lb+20) 746 if cnt >= 20 { if off + cnt <= no { 747 var stride: i64 = cnt/want 748 if stride < 1 { stride = 1 } 749 var sum: i64 = 0 750 var mx: i64 = 0 751 var scnt: i64 = 0 752 let hist: *i64 = sys_mmap(4096*8) as *i64 753 var hz: i64 = 0 754 while hz < 4096 { hist[hz] = 0; hz = hz + 1 } 755 var k: i64 = off 756 while k < off + cnt { 757 let o9: i64 = k*9 758 let px: i64 = (ov[o9] + ov[o9+3] + ov[o9+6])/3 759 let py: i64 = (ov[o9+1] + ov[o9+4] + ov[o9+7])/3 760 let pz: i64 = (ov[o9+2] + ov[o9+5] + ov[o9+8])/3 761 let d2v: i64 = md_near(px, py, pz, cv, startC, refsC) 762 let dum: i64 = md_isqrt(d2v) * MD_UM_PER_Q 763 sum = sum + dum 764 if dum > mx { mx = dum } 765 var hb: i64 = dum/1000 766 if hb >= 4096 { hb = 4095 } 767 hist[hb] = hist[hb] + 1 768 scnt = scnt + 1 769 k = k + stride 770 } 771 if scnt > 0 { 772 var acc: i64 = 0 773 var p95: i64 = 0 774 var hb2: i64 = 0 775 while hb2 < 4096 { 776 acc = acc + hist[hb2] 777 if p95 == 0 { if acc*100 >= scnt*95 { p95 = hb2*1000 } } 778 hb2 = hb2 + 1 779 } 780 if first == 0 { hw("," as *u8) } 781 first = 0 782 hw("{\x22part\x22:\x22" as *u8) 783 var q: i64 = 0 784 while q < 16 { let chq: i64 = b[lb+q] as i64; if chq >= 32 { if chq < 127 { sys_write(1, ((b as i64)+lb+q) as *u8, 1) } } q = q + 1 } 785 hw("\x22,\x22tris\x22:" as *u8); pn(cnt) 786 hw(",\x22samples\x22:" as *u8); pn(scnt) 787 hw(",\x22mean_um\x22:" as *u8); pn(sum/scnt) 788 hw(",\x22p95_um\x22:" as *u8); pn(p95) 789 hw(",\x22max_um\x22:" as *u8); pn(mx) 790 hw("}" as *u8) 791 } 792 } } 793 L = L + 1 794 } 795 hw("],\x22note\x22:\x22completeness per named oracle part: distance from HER part to OUR nearest surface; parts under 20 tris skipped\x22}\n" as *u8) 796 return 0 797} 798 799// ---- selftest: literal-only fixture writes (nx_cc wrong-slot bug 1785936860 workaround) ---- 800func md_fix(path: *u8, zoff_f32: i64) -> i64 { 801 let b: *u8 = sys_mmap(256) 802 b[0]=78 as u8; b[1]=88 as u8; b[2]=77 as u8; b[3]=83 as u8 803 b[4]=72 as u8; b[5]=50 as u8; b[6]=0 as u8; b[7]=0 as u8 804 b[8]=1 as u8; b[9]=0 as u8; b[10]=0 as u8; b[11]=0 as u8 805 b[12]=1 as u8; b[13]=0 as u8; b[14]=0 as u8; b[15]=0 as u8 806 var q: i64 = 0 807 while q < 16 { b[16+q] = 0 as u8; q = q + 1 } 808 b[16]=115 as u8; b[17]=107 as u8; b[18]=105 as u8; b[19]=110 as u8 809 b[32]=0 as u8; b[33]=0 as u8; b[34]=0 as u8; b[35]=0 as u8 810 b[36]=1 as u8; b[37]=0 as u8; b[38]=0 as u8; b[39]=0 as u8 811 // tri (0,0,0)(100,0,0)(0,100,v2z) mm; f32 100.0 = 0x42C80000; only v2's z comes from the arg 812 // (a TILTED PLANE survives centroid+height alignment; a pure translation or coplanar subset is nulled) 813 var o: i64 = 40 814 var k: i64 = 0 815 while k < 9 { b[o+k*4]=0 as u8; b[o+k*4+1]=0 as u8; b[o+k*4+2]=0 as u8; b[o+k*4+3]=0 as u8; k = k + 1 } 816 // v1 x=100f 817 b[o+12]=0 as u8; b[o+13]=0 as u8; b[o+14]=200 as u8; b[o+15]=66 as u8 818 // v2 y=100f, z=arg 819 b[o+28]=0 as u8; b[o+29]=0 as u8; b[o+30]=200 as u8; b[o+31]=66 as u8 820 b[o+32]=(zoff_f32&255) as u8; b[o+33]=((zoff_f32>>8)&255) as u8; b[o+34]=((zoff_f32>>16)&255) as u8; b[o+35]=((zoff_f32>>24)&255) as u8 821 // normals (0,0,1) f32 1.0=0x3F800000 x9... keep zeros (loader ignores), colors zeros 822 // layer id 823 b[o+84]=0 as u8; b[o+85]=0 as u8; b[o+86]=0 as u8; b[o+87]=0 as u8 824 let fd: i64 = sys_openat_wr(path, 420) 825 if fd < 0 { return 0 - 1 } 826 sys_write(fd, b, o+88) 827 sys_close(fd) 828 return 0 829} 830func md_selftest() -> i64 { 831 var fails: i64 = 0 832 md_fix("/tmp/md_a.nxmesh" as *u8, 0) 833 md_fix("/tmp/md_b.nxmesh" as *u8, 1112014848) 834 hw("T0 self-vs-self (expect mean ~0, floor = edge/16384 barycentric quantum):\n" as *u8) 835 if md_run("/tmp/md_a.nxmesh" as *u8, "/tmp/md_a.nxmesh" as *u8, 64) != 0 { fails = fails + 1 } 836 hw("T1 tilted plane, v2 lifted 50mm z (expect mean_um in the LOW THOUSANDS both directions):\n" as *u8) 837 if md_run("/tmp/md_a.nxmesh" as *u8, "/tmp/md_b.nxmesh" as *u8, 64) != 0 { fails = fails + 1 } 838 if md_run("/tmp/md_absent_zz.nxmesh" as *u8, "/tmp/md_a.nxmesh" as *u8, 64) == 0 { fails = fails + 1; hw("T2 FAIL absent accepted\n" as *u8) } else { hw("T2 PASS absent refused\n" as *u8) } 839 if fails == 0 { hw("MESHDIST-SELFTEST GREEN (verify T0 mean==0 and T1 mean~5000 above)\n" as *u8); return 0 } 840 hw("MESHDIST-SELFTEST RED fails=" as *u8); pn(fails); hw("\n" as *u8) 841 return 1 842} 843 844func main(argc: i64, argv: *i64) -> i64 { 845 if argc >= 2 { 846 let a1: *u8 = argv[1] as *u8 847 var m: i64 = 0 848 while a1[m] != (0 as u8) { m = m + 1 } 849 if m == 8 { 850 var ok: i64 = 1 851 let lit: *u8 = "selftest" as *u8 852 var i: i64 = 0 853 while i < 8 { if a1[i] != lit[i] { ok = 0; i = 8 } else { i = i + 1 } } 854 if ok == 1 { let rc: i64 = md_selftest(); sys_exit(rc); return rc } 855 } 856 } 857 // paint verb: deviation heatmap (chars: pa I nt vs pa R ts) 858 if a1[0] == (112 as u8) { if a1[1] == (97 as u8) { if a1[2] == (105 as u8) { 859 if argc < 5 { hw("usage: nx_meshdist paint <cand> <oracle> <out.nxmesh> [band_um=5000]\n" as *u8); sys_exit(2); return 2 } 860 var bu: i64 = 5000 861 if argc >= 6 { 862 let a5: *u8 = argv[5] as *u8 863 bu = 0 864 var q5: i64 = 0 865 while a5[q5] != (0 as u8) { bu = bu*10 + ((a5[q5] as i64) - 48); q5 = q5 + 1 } 866 } 867 let rcq: i64 = md_paint(argv[2] as *u8, argv[3] as *u8, argv[4] as *u8, bu) 868 sys_exit(rcq) 869 return rcq 870 } } } 871 // parts verb: per-named-part completeness table 872 if a1[0] == (112 as u8) { if a1[1] == (97 as u8) { if a1[2] == (114 as u8) { if a1[3] == (116 as u8) { if a1[4] == (115 as u8) { if a1[5] == (0 as u8) { 873 if argc < 4 { hw("usage: nx_meshdist parts <cand.nxmesh> <oracle-with-layers.nxmesh> [samples-per-part]\n" as *u8); sys_exit(2); return 2 } 874 var pw: i64 = 256 875 if argc >= 5 { 876 let a4: *u8 = argv[4] as *u8 877 pw = 0 878 var q4: i64 = 0 879 while a4[q4] != (0 as u8) { pw = pw*10 + ((a4[q4] as i64) - 48); q4 = q4 + 1 } 880 if pw < 16 { pw = 16 } 881 if pw > 4096 { pw = 4096 } 882 } 883 let rcp: i64 = md_parts(argv[2] as *u8, argv[3] as *u8, pw) 884 sys_exit(rcp) 885 return rcp 886 } } } } } } 887 if argc < 3 { 888 hw("usage: nx_meshdist <cand.nxmesh> <oracle.nxmesh> [samples] | parts <cand> <oracle> [n] | selftest\n" as *u8) 889 sys_exit(2) 890 return 2 891 } 892 var want: i64 = 8192 893 if argc >= 4 { 894 let a3: *u8 = argv[3] as *u8 895 want = 0 896 var q: i64 = 0 897 while a3[q] != (0 as u8) { want = want*10 + ((a3[q] as i64) - 48); q = q + 1 } 898 if want < 16 { want = 16 } 899 if want > MD_MAXSAMP { want = MD_MAXSAMP } 900 } 901 let rc2: i64 = md_run(argv[1] as *u8, argv[2] as *u8, want) 902 sys_exit(rc2) 903 return rc2 904}