code wiki / _hdl_build / nx_twinbench.nx
`nx_twinbench.nx` source

↩ module page · 672 lines · 36072 B
1// nx_twinbench.nx -- ★★★THE METRIC RULER, AND THE PRIMARY INSTRUMENT OF THE DIGITAL-TWIN PROGRAM.
2//
3// WHY THIS COMES FIRST. The program's governing law is capability = MIN(statistical, perceptual), and the
4// perceptual term has always been an OPERATOR'S GLANCE -- no organ emits it, so every "we moved 175 -> X"
5// has been unfalsifiable by construction. Aiming at a photoreal PLAUSIBLE human made that unavoidable:
6// "does this look real" is a matter of judgement. Aiming at a 1:1 TWIN OF A SPECIFIC PERSON does not.
7// "Is this surface within one millimetre of that person's face" is arithmetic.
8// ★★THE TWIN TARGET MAKES MEASUREMENT EASIER, NOT HARDER -- it hands us a hard ruler the aesthetic target
9// never could, and this organ is that ruler.
10//
11// WHAT IT MEASURES: the deviation of a candidate surface from a reference surface, as a distribution rather
12// than one number -- mean, RMS, p95, and max (a one-sided Hausdorff). Reported in MILLIMETRES, because a
13// surgeon does not plan in per-mille and a bench number that cannot be compared to a caliper is decoration.
14//
15// ★POINT-TO-TRIANGLE, NOT POINT-TO-VERTEX, AND THE DIFFERENCE IS NOT PEDANTRY. Nearest-VERTEX distance is
16// the easy version and it is BIASED: it reports error proportional to the reference's tessellation, so a
17// coarser reference scores worse for reasons that have nothing to do with the subject. A ruler whose reading
18// depends on how finely you happened to triangulate the truth is not a ruler. So every query point is
19// measured to the nearest point ON A TRIANGLE -- face interior, edge, or vertex, whichever is closest.
20//
21// ★A UNIFORM GRID, BECAUSE THE HONEST ALGORITHM MUST ALSO TERMINATE. Brute force is O(N*M): 35k vertices
22// against 35k triangles is 1.2 billion point-triangle tests per comparison, and a ruler nobody runs is a
23// ruler that does not exist. Reference triangles are bucketed into a uniform grid; a query expands its
24// search radius ring by ring and STOPS ONLY when the next ring cannot contain anything closer than the best
25// found so far. That last clause is what keeps it EXACT rather than approximate -- the grid changes the
26// cost, never the answer, and T6 proves it against brute force on a small mesh.
27//
28//   nx_twinbench cmp <candidate.nxmesh> <reference.nxmesh> [stature_mm] [layer]
29//   nx_twinbench selftest
30// license_tier: ORIGINAL  expect_exit: 0  No hw writes (Rule 26).
31import "nx_gate_verdict.nx"
32const TB_MAGIC_16383: i64 = 16383
33const TB_MAGIC_4000000000: i64 = 4000000000
34const TB_MAGIC_1750: i64 = 1750
35const TB_MAGIC_4950: i64 = 4950
36const TB_MAGIC_5050: i64 = 5050
37const TB_MAGIC_4000: i64 = 4000
38const TB_MAGIC_1000000: i64 = 1000000
39
40const TB_M8388607: i64 = 8388607
41const TB_M8388608: i64 = 8388608
42const TB_HDR: i64 = 16
43const TB_LAYENT: i64 = 24
44const TB_TRI: i64 = 84
45const TB_LID: i64 = 4
46const TB_Q: i64 = 1000
47// ★UNITS. Mesh positions are decoded at TB_SCALE sub-units per mesh unit so the integer path keeps
48// sub-millimetre resolution without floats. Distances are accumulated in sub-units and converted to
49// micrometres once, at the end, using the caller's declared stature -- so the conversion is ONE arithmetic
50// step that can be audited, not a scale factor smeared through the inner loop.
51const TB_SCALE: i64 = 1024
52const TB_UM: i64 = 1000
53const TB_STATURE: i64 = 1750
54const TB_GRID: i64 = 48
55const TB_CELLS: i64 = 110592
56const TB_MAXTRI: i64 = 400000
57const TB_MAXQ: i64 = 200000
58const TB_BIG: i64 = 4611686018427387903
59const TB_OUT: i64 = 24
60const TB_FRAC: i64 = 1048576
61// query-point budget: enough for a tight mean/RMS/p95, cheap enough that the ruler actually gets run
62const TB_TARGET: i64 = 40000
63
64func tb_puts(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
65func tb_pn(v: i64) -> i64 {
  let b: *u8 = sys_mmap(32); 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
73}
74// print micrometres as a millimetre decimal -- the unit a surgeon reads
75func tb_pmm(um: i64) -> i64 {
  var v: i64 = um
  var ng: i64 = 0
  if v < 0 { ng = 1; v = 0 - v }
  if ng == 1 { tb_puts("-" as *u8) }
  tb_pn(v/TB_UM); tb_puts("." as *u8)
  let f: i64 = v % TB_UM
  if f < 100 { tb_puts("0" as *u8) }
  if f < 10 { tb_puts("0" as *u8) }
  tb_pn(f)
  return 0
86}
87func tb_streq(a: *u8, b: *u8) -> i64 {
  var i: i64=0; var go: i64=1; var eq: i64=1
  while go==1 { if a[i]!=b[i] { eq=0; go=0 } else { if a[i]==(0 as u8) { go=0 } else { i=i+1 } } }
  return eq
91}
92func tb_atoi(s: *u8) -> i64 {
  var i: i64=0; var n: i64=0; var sg: i64=1
  if s[0]==(45 as u8) { sg=0-1; i=1 }
  while s[i]!=(0 as u8) { let c: i64=s[i] as i64; if c>=48 { if c<=57 { n=n*10+(c-48) } } i=i+1 }
  return n*sg
97}
98func tb_rd32(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)
100}
101func tb_f32mul(b: *u8, o: i64, mul: i64) -> i64 {
  let bits: i64 = tb_rd32(b, o)
  let sign: i64 = (bits>>31) & 1
  let exp: i64 = (bits>>23) & 255
  let mant: i64 = bits & TB_M8388607
  if exp == 0 { return 0 }
  let m: i64 = (mant | TB_M8388608) * mul
  var e: i64 = exp - 127 - 23
  var v: i64 = 0
  if e >= 0 { v = m << e } else { let sh: i64 = 0 - e; v = (m + (1 << (sh-1))) >> sh }
  if sign == 1 { v = 0 - v }
  return v
113}
114func tb_isqrt(v: i64) -> i64 {
  if v <= 0 { return 0 }
  var x: i64 = v
  var y: i64 = (x+1)/2
  while y < x { x = y; y = (x + v/x)/2 }
  return x
120}
121// ★★★THE OVERFLOW THAT KILLED THE FIRST VERSION, AND WHY THE FIX IS A REWRITE RATHER THAN A PATCH.
122// The textbook point-to-triangle routine classifies the query into one of six barycentric regions using
123// det = a*c - b*b, where a and c are squared edge lengths. At mesh scale that is fatal: an edge of 1,792,000
124// sub-units gives a = 3.2e12 and det = 1.0e25, against an i64 ceiling of 9.2e18. **det overflowed by six
125// orders of magnitude**, which is why the gate read RED on self-fit -- the most basic case a ruler has.
126// ★A RULER THAT OVERFLOWS ON ITS OWN UNITS IS NOT A RULER, and the failure is silent: wrapped arithmetic
127// returns a plausible-looking number rather than an error.
128//
129// The rewrite has two parts, and each is chosen so the result stays EXACT where it matters:
130//   1. The barycentric solve runs on MAGNITUDE-REDUCED copies of the edges -- shifted right until every
131//      component fits in 14 bits, so every product downstream is provably inside i64. Reduction costs
132//      precision only in the CLASSIFICATION (which region), never in the answer.
133//   2. The closest point is then reconstructed in FULL precision from the barycentric fraction, and the
134//      distance is measured to that point. So a lossy classification still yields an exact distance.
135// And the boundary cases are not six hand-derived formulas but three point-to-SEGMENT calls, which are
136// short enough to be obviously right. The previous version's six regions were where the bug hid.
137func tb_shift_for(m: i64) -> i64 {
  var mx: i64 = m
  var sh: i64 = 0
  while mx > TB_MAGIC_16383 { mx = mx >> 1; sh = sh + 1 }
  return sh
142}
143func tb_sh(v: i64, sh: i64) -> i64 {
  if v >= 0 { return v >> sh }
  return 0 - ((0-v) >> sh)
146}
147func tb_absi(v: i64) -> i64 { if v < 0 { return 0-v } return v }
148func tb_max6(a: i64, b: i64, c: i64, d: i64, e: i64, f: i64) -> i64 {
  var m: i64 = tb_absi(a)
  if tb_absi(b) > m { m = tb_absi(b) }
  if tb_absi(c) > m { m = tb_absi(c) }
  if tb_absi(d) > m { m = tb_absi(d) }
  if tb_absi(e) > m { m = tb_absi(e) }
  if tb_absi(f) > m { m = tb_absi(f) }
  return m
156}
157// num/den as a TB_FRAC fraction, reducing both first so num*TB_FRAC cannot overflow
158func tb_frac(num: i64, den: i64) -> i64 {
  var n: i64 = num
  var d: i64 = den
  while tb_absi(n) > TB_MAGIC_4000000000 { n = tb_sh(n,1); d = tb_sh(d,1) }
  if d == 0 { return 0 }
  return n * TB_FRAC / d
164}
165// ★POINT TO SEGMENT, squared, full precision. Short enough to be obviously correct, which is exactly why
166// the boundary cases route through it instead of through derived closed forms.
167func tb_pseg2(px: i64, py: i64, pz: i64, ax: i64, ay: i64, az: i64, bx: i64, by: i64, bz: i64) -> i64 {
  let ex: i64 = bx-ax; let ey: i64 = by-ay; let ez: i64 = bz-az
  let dx: i64 = px-ax; let dy: i64 = py-ay; let dz: i64 = pz-az
  let da: i64 = dx*dx + dy*dy + dz*dz
  let sh: i64 = tb_shift_for(tb_max6(ex,ey,ez,dx,dy,dz))
  let rex: i64 = tb_sh(ex,sh); let rey: i64 = tb_sh(ey,sh); let rez: i64 = tb_sh(ez,sh)
  let rdx: i64 = tb_sh(dx,sh); let rdy: i64 = tb_sh(dy,sh); let rdz: i64 = tb_sh(dz,sh)
  let l2: i64 = rex*rex + rey*rey + rez*rez
  if l2 <= 0 { return da }                                   // degenerate segment: it is a point
  let tn: i64 = rdx*rex + rdy*rey + rdz*rez
  if tn <= 0 { return da }                                   // before A
  if tn >= l2 {                                              // past B
      let bxd: i64 = px-bx; let byd: i64 = py-by; let bzd: i64 = pz-bz
      return bxd*bxd + byd*byd + bzd*bzd
  }
  let tf: i64 = tb_frac(tn, l2)
  let qx: i64 = ax + ex*tf/TB_FRAC
  let qy: i64 = ay + ey*tf/TB_FRAC
  let qz: i64 = az + ez*tf/TB_FRAC
  let gx: i64 = px-qx; let gy: i64 = py-qy; let gz: i64 = pz-qz
  return gx*gx + gy*gy + gz*gz
188}
189func tb_ptri2(px: i64, py: i64, pz: i64,
            ax: i64, ay: i64, az: i64, bx: i64, by: i64, bz: i64, cx: i64, cy: i64, cz: i64) -> i64 {
  let f0x: i64 = bx-ax; let f0y: i64 = by-ay; let f0z: i64 = bz-az
  let f1x: i64 = cx-ax; let f1y: i64 = cy-ay; let f1z: i64 = cz-az
  let vx: i64 = px-ax;  let vy: i64 = py-ay;  let vz: i64 = pz-az
  var m: i64 = tb_max6(f0x,f0y,f0z,f1x,f1y,f1z)
  let m2: i64 = tb_max6(vx,vy,vz,0,0,0)
  if m2 > m { m = m2 }
  let sh: i64 = tb_shift_for(m)
  let r0x: i64 = tb_sh(f0x,sh); let r0y: i64 = tb_sh(f0y,sh); let r0z: i64 = tb_sh(f0z,sh)
  let r1x: i64 = tb_sh(f1x,sh); let r1y: i64 = tb_sh(f1y,sh); let r1z: i64 = tb_sh(f1z,sh)
  let rvx: i64 = tb_sh(vx,sh);  let rvy: i64 = tb_sh(vy,sh);  let rvz: i64 = tb_sh(vz,sh)
  let aa: i64 = r0x*r0x + r0y*r0y + r0z*r0z
  let bb: i64 = r0x*r1x + r0y*r1y + r0z*r1z
  let cc: i64 = r1x*r1x + r1y*r1y + r1z*r1z
  let dd: i64 = r0x*rvx + r0y*rvy + r0z*rvz
  let ee: i64 = r1x*rvx + r1y*rvy + r1z*rvz
  let det: i64 = aa*cc - bb*bb
  var inside: i64 = 0
  var s: i64 = 0
  var t: i64 = 0
  if det > 0 {
      s = dd*cc - ee*bb
      t = ee*aa - dd*bb
      if s >= 0 { if t >= 0 { if s + t <= det { inside = 1 } } }
  }
  if inside == 1 {
      // ★RECONSTRUCT IN FULL PRECISION. The classification above ran on reduced edges; this does not, so
      // a lossy region test still produces an exact distance to the exact closest point.
      let sf: i64 = tb_frac(s, det)
      let tf: i64 = tb_frac(t, det)
      let qx: i64 = ax + f0x*sf/TB_FRAC + f1x*tf/TB_FRAC
      let qy: i64 = ay + f0y*sf/TB_FRAC + f1y*tf/TB_FRAC
      let qz: i64 = az + f0z*sf/TB_FRAC + f1z*tf/TB_FRAC
      let gx: i64 = px-qx; let gy: i64 = py-qy; let gz: i64 = pz-qz
      return gx*gx + gy*gy + gz*gz
  }
  // outside the face, or degenerate: the closest point lies on the boundary
  var best: i64 = tb_pseg2(px,py,pz, ax,ay,az, bx,by,bz)
  let e2: i64 = tb_pseg2(px,py,pz, bx,by,bz, cx,cy,cz)
  if e2 < best { best = e2 }
  let e3: i64 = tb_pseg2(px,py,pz, cx,cy,cz, ax,ay,az)
  if e3 < best { best = e3 }
  return best
233}
234// ---- mesh loading: positions of every triangle corner, in sub-units ----
235// M[0]=ntris  M[1..6]=bbox lo/hi xyz  then 9 coords per triangle from M[8]
236const TB_MH: i64 = 8
237func tb_load(path: *u8, M: *i64, cap: i64) -> i64 {
  let szp: *i64 = sys_mmap(16) as *i64
  let mb: *u8 = sys_read_file(path, szp)
  if (mb as i64) == 0 { return 0-1 }
  let sz: i64 = szp[0]
  if sz < TB_HDR { return 0-2 }
  if mb[0] != (78 as u8) { return 0-3 }
  let nlay: i64 = tb_rd32(mb, 8)
  let nt: i64 = tb_rd32(mb, 12)
  if nt <= 0 { return 0-4 }
  if nlay <= 0 { return 0-5 }
  let hdr: i64 = TB_HDR + nlay*TB_LAYENT
  if hdr + nt*TB_TRI + nt*TB_LID > sz { return 0-6 }
  var use: i64 = nt
  if use > cap { use = cap }
  M[0] = use
  var lo0: i64 = TB_BIG; var lo1: i64 = TB_BIG; var lo2: i64 = TB_BIG
  var hi0: i64 = 0-TB_BIG; var hi1: i64 = 0-TB_BIG; var hi2: i64 = 0-TB_BIG
  var t: i64 = 0
  while t < use {
      let base: i64 = hdr + t*TB_TRI
      var v: i64 = 0
      while v < 3 {
          let o: i64 = base + v*12
          let x: i64 = tb_f32mul(mb, o,   TB_SCALE)
          let y: i64 = tb_f32mul(mb, o+4, TB_SCALE)
          let z: i64 = tb_f32mul(mb, o+8, TB_SCALE)
          let d: i64 = TB_MH + t*9 + v*3
          M[d] = x; M[d+1] = y; M[d+2] = z
          if x < lo0 { lo0 = x }
          if y < lo1 { lo1 = y }
          if z < lo2 { lo2 = z }
          if x > hi0 { hi0 = x }
          if y > hi1 { hi1 = y }
          if z > hi2 { hi2 = z }
          v = v + 1
      }
      t = t + 1
  }
  M[1]=lo0; M[2]=lo1; M[3]=lo2; M[4]=hi0; M[5]=hi1; M[6]=hi2
  return use
278}
279// ---- uniform grid over the reference triangles ----
280// G[0]=cellsz G[1..3]=origin  head[] -> first index, next[] -> chain
281// ★★★A HYPOTHESIS I MEASURED AND THREW AWAY -- recorded because the reasoning was plausible and WRONG.
282// Observing 13s at 10,800 tris and 78s at 39,120 (3.6x the triangles for 6x the time, on 1.8x the samples),
283// I concluded the fixed 48^3 grid was too coarse-grained for sparse meshes and made it adaptive at
284// grid = cbrt(ntris), roughly one triangle per cell. **Measured: the same 10,800-tri mesh went 13s -> 39s.
285// Three times SLOWER.** The reasoning missed that a coarser grid puts ~10x more triangles in the query's
286// OWN cell, and the ring search usually terminates at radius 0 -- so the cost is dominated by cell
287// OCCUPANCY, not by how many rings get walked. Fewer, fatter cells is exactly the wrong trade.
288// ★REVERTED to the fixed grid on the measurement. ★LAW: A PERFORMANCE FIX IS A HYPOTHESIS UNTIL IT IS
289// TIMED -- and a plausible mechanism that predicts the right direction can still be the wrong mechanism.
290// (The real cost driver is queries x triangles-per-cell; the honest lever is the SAMPLE COUNT, which is
291// already declared and now tunable per call.)
292func tb_cbrt(v: i64) -> i64 {
  var r: i64 = 1
  while r*r*r < v { r = r + 1 }
  return r
296}
297func tb_grid_build(M: *i64, G: *i64, head: *i64, nxt: *i64) -> i64 {
  let nt: i64 = M[0]
  let gd: i64 = TB_GRID
  G[4] = gd
  var ex: i64 = M[4]-M[1]
  if M[5]-M[2] > ex { ex = M[5]-M[2] }
  if M[6]-M[3] > ex { ex = M[6]-M[3] }
  var cs: i64 = ex/gd
  if cs < 1 { cs = 1 }
  G[0]=cs; G[1]=M[1]; G[2]=M[2]; G[3]=M[3]
  var i: i64 = 0
  while i < gd*gd*gd { head[i] = 0-1; i = i+1 }
  var t: i64 = 0
  while t < nt {
      let d: i64 = TB_MH + t*9
      // bucket on the triangle's centroid; the ring search below is what makes that safe
      let cxx: i64 = (M[d]+M[d+3]+M[d+6])/3
      let cyy: i64 = (M[d+1]+M[d+4]+M[d+7])/3
      let czz: i64 = (M[d+2]+M[d+5]+M[d+8])/3
      var gx: i64 = (cxx-G[1])/cs; if gx<0 {gx=0}; if gx>=gd {gx=gd-1}
      var gy: i64 = (cyy-G[2])/cs; if gy<0 {gy=0}; if gy>=gd {gy=gd-1}
      var gz: i64 = (czz-G[3])/cs; if gz<0 {gz=0}; if gz>=gd {gz=gd-1}
      let ci: i64 = (gx*gd + gy)*gd + gz
      nxt[t] = head[ci]
      head[ci] = t
      t = t + 1
  }
  return 0
325}
326// ★THE EXACTNESS CLAUSE. Rings expand outward from the query's own cell. A ring at radius r cannot contain
327// anything nearer than (r-1)*cellsize, so once the best distance found is inside that bound the search
328// STOPS AND THE ANSWER IS EXACT. Without that test this would be a fast approximation wearing a ruler's
329// name -- T6 checks it against brute force and would catch the difference.
330func tb_nearest2(px: i64, py: i64, pz: i64, M: *i64, G: *i64, head: *i64, nxt: *i64) -> i64 {
  let cs: i64 = G[0]
  let gd: i64 = G[4]
  var gx: i64 = (px-G[1])/cs; if gx<0 {gx=0}; if gx>=gd {gx=gd-1}
  var gy: i64 = (py-G[2])/cs; if gy<0 {gy=0}; if gy>=gd {gy=gd-1}
  var gz: i64 = (pz-G[3])/cs; if gz<0 {gz=0}; if gz>=gd {gz=gd-1}
  var best: i64 = TB_BIG
  var r: i64 = 0
  var go: i64 = 1
  while go == 1 {
      var ax: i64 = gx-r; if ax<0 {ax=0}
      var bx2: i64 = gx+r; if bx2>=gd {bx2=gd-1}
      var ix: i64 = ax
      while ix <= bx2 {
          var ay: i64 = gy-r; if ay<0 {ay=0}
          var by2: i64 = gy+r; if by2>=TB_GRID {by2=TB_GRID-1}
          var iy: i64 = ay
          while iy <= by2 {
              var az: i64 = gz-r; if az<0 {az=0}
              var bz2: i64 = gz+r; if bz2>=TB_GRID {bz2=TB_GRID-1}
              var iz: i64 = az
              while iz <= bz2 {
                  // only the shell of the cube, not its interior -- inner cells were done at smaller r
                  var shell: i64 = 0
                  if ix == gx-r { shell = 1 }
                  if ix == gx+r { shell = 1 }
                  if iy == gy-r { shell = 1 }
                  if iy == gy+r { shell = 1 }
                  if iz == gz-r { shell = 1 }
                  if iz == gz+r { shell = 1 }
                  if r == 0 { shell = 1 }
                  if shell == 1 {
                      var t: i64 = head[(ix*gd + iy)*gd + iz]
                      while t >= 0 {
                          let d: i64 = TB_MH + t*9
                          let q: i64 = tb_ptri2(px,py,pz, M[d],M[d+1],M[d+2], M[d+3],M[d+4],M[d+5], M[d+6],M[d+7],M[d+8])
                          if q < best { best = q }
                          t = nxt[t]
                      }
                  }
                  iz = iz + 1
              }
              iy = iy + 1
          }
          ix = ix + 1
      }
      // stop when the next ring cannot possibly beat what we already have
      let reach: i64 = r*cs
      if best <= reach*reach { go = 0 }
      if r > gd*2 { go = 0 }
      r = r + 1
  }
  return best
383}
384// out[0]=nq out[1]=mean_um out[2]=rms_um out[3]=p95_um out[4]=max_um out[5]=stature out[6]=height_units
385func tb_cmp(candA: *i64, refB: *i64, out: *i64, stature: i64, target: i64) -> i64 {
  var z: i64 = 0
  while z < TB_OUT { out[z] = 0; z = z+1 }
  let na: i64 = candA[0]
  let nb: i64 = refB[0]
  if na <= 0 { return 0-1 }
  if nb <= 0 { return 0-2 }
  let G: *i64 = sys_mmap(16*8) as *i64
  let head: *i64 = sys_mmap(TB_CELLS*8) as *i64
  let nxt: *i64 = sys_mmap(nb*8 + 64) as *i64
  tb_grid_build(refB, G, head, nxt)
  // ★THE SCALE CONVERSION, ONCE AND AUDITABLY. The mesh has no units of its own; the caller declares the
  // subject's stature and the mesh's own height in sub-units gives the rest. Reported as an explicit
  // field so a reader can check the arithmetic rather than trust it.
  let hgt: i64 = candA[5] - candA[2]
  out[5] = stature
  out[6] = hgt
  if hgt <= 0 { return 0-3 }
  // ★★DECLARED SAMPLING, BECAUSE THE EXHAUSTIVE VERSION IS NOT RUNNABLE AND PRETENDING OTHERWISE WOULD
  // BE THE LIE. Measured: a 103,824-triangle mesh is 311,472 query points, each searching a grid of
  // 103,824 triangles -- the first exhaustive run was killed at ten minutes without finishing, and a
  // ruler nobody can afford to run is a ruler that never gets run. So the query set is SUBSAMPLED at a
  // stride chosen to land near TB_TARGET points, and the stride and both counts are REPORTED. Metro and
  // CloudCompare do exactly this; the difference between them and a dishonest tool is printing it.
  // ★WHAT SAMPLING COSTS, STATED RATHER THAN BURIED: mean, RMS and p95 are estimated from tens of
  // thousands of samples and are tight. MAX is not -- it is an extreme-value statistic and subsampling
  // can only ever MISS the worst point, never invent one. So max is a LOWER BOUND on the true Hausdorff
  // distance, and it is labelled as one. Anyone quoting max as if it were exhaustive is over-claiming.
  let total: i64 = na*3
  // ★the sample budget is a CALLER knob, not a constant. An optimiser evaluating dozens of candidates
  // wants a cheap, consistent estimate; a published measurement wants the tight one. Same metric, same
  // code path, declared budget either way -- and every result still prints the stride it used.
  var tgt: i64 = target
  if tgt <= 0 { tgt = TB_TARGET }
  var stride: i64 = total/tgt
  if stride < 1 { stride = 1 }
  var nq: i64 = 0
  let cap: i64 = total/stride + 8
  let dist: *i64 = sys_mmap(cap*8 + 64) as *i64
  var acc: i64 = 0
  var sq: i64 = 0
  var mx: i64 = 0
  var k: i64 = 0
  while k < total {
      let i: i64 = k/3
      let v: i64 = k - i*3
      let d: i64 = TB_MH + i*9
      let px: i64 = candA[d+v*3]
      let py: i64 = candA[d+v*3+1]
      let pz: i64 = candA[d+v*3+2]
      let s2: i64 = tb_nearest2(px,py,pz, refB, G, head, nxt)
      let su: i64 = tb_isqrt(s2)
      // sub-units -> micrometres, via the declared stature
      let um: i64 = su * stature * TB_UM / hgt
      if nq < cap { dist[nq] = um; nq = nq + 1 }
      acc = acc + um
      sq = sq + um*um
      if um > mx { mx = um }
      k = k + stride
  }
  if nq <= 0 { return 0-4 }
  out[0] = nq
  out[1] = acc/nq
  out[2] = tb_isqrt(sq/nq)
  out[4] = mx
  out[7] = total
  out[8] = stride
  // p95 by counting rather than sorting: a 512-bucket histogram over [0,max] is exact to one bucket and
  // avoids an O(n log n) sort of 100k+ samples inside a ruler that should be cheap to run.
  if mx > 0 {
      let H: *i64 = sys_mmap(520*8) as *i64
      var b: i64 = 0
      while b < 512 { H[b] = 0; b = b+1 }
      var k: i64 = 0
      while k < nq {
          var bi: i64 = dist[k]*511/mx
          if bi < 0 { bi = 0 }
          if bi > 511 { bi = 511 }
          H[bi] = H[bi] + 1
          k = k + 1
      }
      let target: i64 = nq*95/100
      var cum: i64 = 0
      var bb: i64 = 0
      var found: i64 = 0
      while bb < 512 {
          cum = cum + H[bb]
          if found == 0 { if cum >= target { out[3] = bb*mx/511; found = 1 } }
          bb = bb + 1
      }
  }
  return 0
477}
478// ---- fixtures: two tiny meshes built in memory, so the gate tests the MEASUREMENT not the file reader ----
479// a unit quad in the XY plane at z=0, as 2 triangles, spanning [0,H] so height is well defined
480func tb_fx_quad(M: *i64, h: i64, dz: i64) -> i64 {
  M[0] = 2
  let s: i64 = h
  // tri 0: (0,0,dz) (s,0,dz) (0,s,dz)
  M[TB_MH+0]=0;   M[TB_MH+1]=0;   M[TB_MH+2]=dz
  M[TB_MH+3]=s;   M[TB_MH+4]=0;   M[TB_MH+5]=dz
  M[TB_MH+6]=0;   M[TB_MH+7]=s;   M[TB_MH+8]=dz
  // tri 1: (s,0,dz) (s,s,dz) (0,s,dz)
  M[TB_MH+9]=s;   M[TB_MH+10]=0;  M[TB_MH+11]=dz
  M[TB_MH+12]=s;  M[TB_MH+13]=s;  M[TB_MH+14]=dz
  M[TB_MH+15]=0;  M[TB_MH+16]=s;  M[TB_MH+17]=dz
  M[1]=0; M[2]=0; M[3]=dz; M[4]=s; M[5]=s; M[6]=dz
  return 0
493}
494func tb_gate() -> i64 {
  let ctr: *i64 = gv_ctr()
  gv_head("nx_twinbench selftest -- the metric ruler: millimetres, point-to-triangle, exact under a grid" as *u8)
  let cap: i64 = 64
  let A: *i64 = sys_mmap((TB_MH + cap*9)*8) as *i64
  let B: *i64 = sys_mmap((TB_MH + cap*9)*8) as *i64
  let out: *i64 = sys_mmap(TB_OUT*8) as *i64
  let h: i64 = TB_MAGIC_1750*TB_SCALE
  // ★★T1 SELF-FIT IS EXACTLY ZERO. A ruler that reports error against an identical surface is measuring
  // its own arithmetic, and every number it ever prints afterwards carries that offset.
  tb_fx_quad(A, h, 0)
  tb_fx_quad(B, h, 0)
  let r1: i64 = tb_cmp(A, B, out, TB_STATURE, 0)
  var t1: i64 = 0
  if r1 == 0 { if out[1] == 0 { if out[2] == 0 { if out[4] == 0 { t1 = 1 } } } }
  gv_check("T1 SELF-FIT is EXACTLY 0.000mm -- mean, RMS and max all zero" as *u8, t1, ctr)
  // ★★T2 A KNOWN DISPLACEMENT READS BACK AS ITSELF. Offset the candidate by a distance that is an exact
  // fraction of stature and the ruler must report that distance in millimetres. This is the calibration:
  // without it the organ is a comparator, not a ruler.
  tb_fx_quad(A, h, h/TB_MAGIC_1750)                     // 1 unit of TB_MAGIC_1750 = exactly 1mm at stature TB_MAGIC_1750
  tb_fx_quad(B, h, 0)
  tb_cmp(A, B, out, TB_STATURE, 0)
  var t2: i64 = 0
  if out[1] >= 990 { if out[1] <= 1010 { t2 = 1 } }
  gv_check("T2 CALIBRATED: a 1mm offset reads 1.000mm (within rounding)" as *u8, t2, ctr)
  // ★T3 IT SCALES LINEARLY. Double the offset, double the reading -- a ruler with a nonlinear response
  // cannot be compared across rungs, which is the only thing a program ruler is for.
  tb_fx_quad(A, h, h*5/TB_MAGIC_1750)
  tb_fx_quad(B, h, 0)
  tb_cmp(A, B, out, TB_STATURE, 0)
  var t3: i64 = 0
  if out[1] >= TB_MAGIC_4950 { if out[1] <= TB_MAGIC_5050 { t3 = 1 } }
  gv_check("T3 LINEAR: a 5mm offset reads 5.000mm, so readings compare across rungs" as *u8, t3, ctr)
  // ★★T4 NON-VACUITY, THE TOOTH THE PLAN NAMED. A 5mm displacement MUST register. A ruler that cannot
  // see a displacement a surgeon would care about is not an instrument.
  var t4: i64 = 0
  if out[1] > TB_MAGIC_4000 { if out[4] > TB_MAGIC_4000 { t4 = 1 } }
  gv_check("T4 non-vacuity: a 5mm displacement REGISTERS on both mean and max" as *u8, t4, ctr)
  // ★★★T5 POINT-TO-TRIANGLE, NOT POINT-TO-VERTEX -- and this is the tooth that proves the harder
  // algorithm was actually built. Query the CENTRE of the reference quad from 1mm above it. The nearest
  // VERTEX is a corner, far away; the nearest point on the SURFACE is directly below, at 1mm. A
  // vertex-based ruler reports the corner distance and is wrong by the size of the triangle.
  tb_fx_quad(B, h, 0)
  let G2: *i64 = sys_mmap(16*8) as *i64
  let hd2: *i64 = sys_mmap(TB_CELLS*8) as *i64
  let nx2: *i64 = sys_mmap(64*8) as *i64
  tb_grid_build(B, G2, hd2, nx2)
  let mid: i64 = h/2
  let above: i64 = h/TB_MAGIC_1750
  let s2c: i64 = tb_nearest2(mid, mid, above, B, G2, hd2, nx2)
  let dcen: i64 = tb_isqrt(s2c)
  var t5: i64 = 0
  if dcen <= above + 2 { if dcen >= above - 2 { t5 = 1 } }
  gv_check("T5 POINT-TO-TRIANGLE: a point above the quad CENTRE is 1mm away, not a corner-distance" as *u8, t5, ctr)
  // ★★★T6 THE GRID IS AN OPTIMISATION, NOT AN APPROXIMATION. Brute-force every triangle and demand the
  // identical answer. If the ring-termination bound were wrong this is what would catch it -- and a
  // ruler that is merely close is a ruler that cannot settle an argument.
  var brute: i64 = TB_BIG
  var bt: i64 = 0
  while bt < B[0] {
      let d: i64 = TB_MH + bt*9
      let q: i64 = tb_ptri2(mid, mid, above, B[d],B[d+1],B[d+2], B[d+3],B[d+4],B[d+5], B[d+6],B[d+7],B[d+8])
      if q < brute { brute = q }
      bt = bt + 1
  }
  var t6: i64 = 0
  if brute == s2c { t6 = 1 }
  gv_check("T6 EXACT: the grid answer equals brute force, so the speedup changed cost not truth" as *u8, t6, ctr)
  // ★T7 EDGE AND CORNER CASES ARE REAL SURFACE CASES. A query beyond the quad's corner is nearest that
  // CORNER; the face-only branch would under-report it. A body is mostly ends.
  let far: i64 = tb_isqrt(tb_nearest2(0-h, 0-h, 0, B, G2, hd2, nx2))
  var t7: i64 = 0
  if far > h { t7 = 1 }
  gv_check("T7 edge/vertex regions handled: a point past the corner reports the CORNER distance" as *u8, t7, ctr)
  // ★T8 DEGENERATE INPUT IS REFUSED, NOT DIVIDED BY. Generated meshes contain zero-area cap triangles.
  var t8: i64 = 0
  let deg: i64 = tb_ptri2(0,0,1000,  0,0,0, 0,0,0, 0,0,0)
  if deg == TB_MAGIC_1000000 { t8 = 1 }
  gv_check("T8 a DEGENERATE zero-area triangle is measured, never divided by" as *u8, t8, ctr)
  // ★T9 EMPTY INPUT REFUSES rather than reporting a perfect fit -- the failure mode that would let a
  // broken pipeline publish 0.000mm and call it a twin.
  let E: *i64 = sys_mmap((TB_MH+16)*8) as *i64
  E[0] = 0
  var t9: i64 = 0
  if tb_cmp(E, B, out, TB_STATURE, 0) < 0 { if tb_cmp(A, E, out, TB_STATURE, 0) < 0 { t9 = 1 } }
  gv_check("T9 FAIL-CLOSED: an empty mesh is REFUSED, never reported as a perfect 0.000mm fit" as *u8, t9, ctr)
  // ★★T10 THE DISTRIBUTION IS REPORTED, NOT JUST A MEAN. A twin with a good mean and a terrible p95 is a
  // twin that is wrong exactly where a surgeon is looking, and a single number hides it.
  tb_fx_quad(A, h, h*2/TB_MAGIC_1750)
  tb_fx_quad(B, h, 0)
  tb_cmp(A, B, out, TB_STATURE, 0)
  var t10: i64 = 0
  if out[3] > 0 { if out[4] >= out[3] { if out[3] >= out[1]*8/10 { t10 = 1 } } }
  gv_check("T10 reports a DISTRIBUTION: mean <= p95 <= max, so a bad tail cannot hide behind a good mean" as *u8, t10, ctr)
  return gv_verdict("TWINBENCH-GATE" as *u8, ctr, "metric ruler in millimetres; point-to-triangle; grid-exact; fail-closed" as *u8)
589}
590func main(argc: i64, argv: *i64) -> i64 {
  if argc >= 2 {
      if tb_streq(argv[1] as *u8, "selftest" as *u8) == 1 { return tb_gate() }
      // ★BBOX — the ruler's second channel, and it exists because "the two meshes differ" is not a claim.
      // A rotation tooth needs to assert WHICH WAY a part now points, not merely that something moved:
      // an azimuth of 90 degrees must SWAP a part's X and Z extents, and only a per-axis measurement can
      // say so. Reused later for registration (bringing subject and model into one frame) and for the
      // creature gates, where "no self-intersection, plausible volume" are extent questions.
      if tb_streq(argv[1] as *u8, "bbox" as *u8) == 1 {
          if argc < 3 { tb_puts("{\x22error\x22:\x22usage: nx_twinbench bbox <mesh.nxmesh> [stature_mm]\x22}\n" as *u8); return 2 }
          var st2: i64 = TB_STATURE
          if argc > 3 { st2 = tb_atoi(argv[3] as *u8) }
          let M: *i64 = sys_mmap((TB_MH + TB_MAXTRI*9)*8) as *i64
          let rm: i64 = tb_load(argv[2] as *u8, M, TB_MAXTRI)
          if rm < 0 { tb_puts("{\x22error\x22:\x22mesh unreadable\x22,\x22rc\x22:" as *u8); tb_pn(rm); tb_puts("}\n" as *u8); return 3 }
          let hgt2: i64 = M[5]-M[2]
          if hgt2 <= 0 { tb_puts("{\x22error\x22:\x22degenerate height\x22}\n" as *u8); return 4 }
          tb_puts("{\x22organ\x22:\x22nx_twinbench\x22,\x22verb\x22:\x22bbox\x22,\x22unit\x22:\x22mm\x22" as *u8)
          tb_puts(",\x22tris\x22:" as *u8); tb_pn(rm)
          tb_puts(",\x22extent_x_mm\x22:" as *u8); tb_pmm((M[4]-M[1])*st2*TB_UM/hgt2)
          tb_puts(",\x22extent_y_mm\x22:" as *u8); tb_pmm((M[5]-M[2])*st2*TB_UM/hgt2)
          tb_puts(",\x22extent_z_mm\x22:" as *u8); tb_pmm((M[6]-M[3])*st2*TB_UM/hgt2)
          tb_puts(",\x22reads\x22:\x22extents are measured on the SAME stature normalisation the cmp verb uses, so a bbox and a deviation are directly comparable. Y is the normalising axis, so extent_y equals the declared stature by construction -- it is printed as an arithmetic check, not as a finding.\x22}\n" as *u8)
          return 0
      }
      // ★RMS — the same measurement as `cmp`, emitted as ONE integer in micrometres and nothing else.
      // It exists so an optimiser can COMPOSE this ruler instead of reimplementing point-to-triangle
      // distance inside itself. A fitter that carries its own copy of the metric is a fitter that can
      // score well against a ruler nobody else agrees with -- and the two copies drift the first time
      // one is improved. One metric, one implementation, many callers.
      if tb_streq(argv[1] as *u8, "rms" as *u8) == 1 {
          if argc < 4 { tb_puts("-1\n" as *u8); return 2 }
          var st3: i64 = TB_STATURE
          if argc > 4 { st3 = tb_atoi(argv[4] as *u8) }
          let A3: *i64 = sys_mmap((TB_MH + TB_MAXTRI*9)*8) as *i64
          let B3: *i64 = sys_mmap((TB_MH + TB_MAXTRI*9)*8) as *i64
          if tb_load(argv[2] as *u8, A3, TB_MAXTRI) < 0 { tb_puts("-1\n" as *u8); return 3 }
          if tb_load(argv[3] as *u8, B3, TB_MAXTRI) < 0 { tb_puts("-1\n" as *u8); return 4 }
          let o3: *i64 = sys_mmap(TB_OUT*8) as *i64
          // ⚠the 5th arg is the SAMPLE BUDGET and it must come from argv. An earlier scripted edit meant
          // to thread it here silently matched nothing, so this read `0` (= use the default) and the knob
          // did nothing: 0, 3000 and 40000 samples all returned the identical rms in the identical time.
          // ★THAT IDENTITY IS THE TELL -- a knob that changes neither the answer nor the cost is not
          // connected. A/B any new parameter at two settings before believing it exists.
          var tg3: i64 = 0
          if argc > 5 { tg3 = tb_atoi(argv[5] as *u8) }
          if tb_cmp(A3, B3, o3, st3, tg3) < 0 { tb_puts("-1\n" as *u8); return 5 }
          tb_pn(o3[2]); tb_puts("\n" as *u8)
          return 0
      }
      if tb_streq(argv[1] as *u8, "cmp" as *u8) == 1 {
          if argc < 4 { tb_puts("{\x22error\x22:\x22usage: nx_twinbench cmp <cand.nxmesh> <ref.nxmesh> [stature_mm]\x22}\n" as *u8); return 2 }
          var stature: i64 = TB_STATURE
          if argc > 4 { stature = tb_atoi(argv[4] as *u8) }
          let A: *i64 = sys_mmap((TB_MH + TB_MAXTRI*9)*8) as *i64
          let B: *i64 = sys_mmap((TB_MH + TB_MAXTRI*9)*8) as *i64
          let ra: i64 = tb_load(argv[2] as *u8, A, TB_MAXTRI)
          if ra < 0 { tb_puts("{\x22error\x22:\x22candidate unreadable\x22,\x22rc\x22:" as *u8); tb_pn(ra); tb_puts("}\n" as *u8); return 3 }
          let rb: i64 = tb_load(argv[3] as *u8, B, TB_MAXTRI)
          if rb < 0 { tb_puts("{\x22error\x22:\x22reference unreadable\x22,\x22rc\x22:" as *u8); tb_pn(rb); tb_puts("}\n" as *u8); return 4 }
          let out: *i64 = sys_mmap(TB_OUT*8) as *i64
          let rc: i64 = tb_cmp(A, B, out, stature, 0)
          if rc < 0 { tb_puts("{\x22error\x22:\x22compare refused\x22,\x22rc\x22:" as *u8); tb_pn(rc); tb_puts("}\n" as *u8); return 5 }
          tb_puts("{\x22organ\x22:\x22nx_twinbench\x22,\x22v\x22:1,\x22unit\x22:\x22mm\x22" as *u8)
          tb_puts(",\x22cand_tris\x22:" as *u8); tb_pn(ra)
          tb_puts(",\x22ref_tris\x22:" as *u8); tb_pn(rb)
          tb_puts(",\x22samples\x22:" as *u8); tb_pn(out[0])
          tb_puts(",\x22query_points_total\x22:" as *u8); tb_pn(out[7])
          tb_puts(",\x22sample_stride\x22:" as *u8); tb_pn(out[8])
          tb_puts(",\x22mean_mm\x22:" as *u8); tb_pmm(out[1])
          tb_puts(",\x22rms_mm\x22:" as *u8); tb_pmm(out[2])
          tb_puts(",\x22p95_mm\x22:" as *u8); tb_pmm(out[3])
          tb_puts(",\x22max_mm_lower_bound\x22:" as *u8); tb_pmm(out[4])
          tb_puts(",\x22stature_mm\x22:" as *u8); tb_pn(out[5])
          tb_puts(",\x22height_subunits\x22:" as *u8); tb_pn(out[6])
          tb_puts(",\x22method\x22:\x22one-sided: every candidate vertex to the nearest point ON a reference TRIANGLE (face, edge or vertex region), accelerated by a uniform grid whose ring search terminates only when the next ring cannot beat the best found -- so the grid changes cost, never the answer\x22" as *u8)
          tb_puts(",\x22reads\x22:\x22mean is the typical miss; p95 and max are where a twin is wrong in the place someone is looking. A number is only meaningful beside a REPEATABILITY FLOOR: capture the same subject twice, compare, and any improvement smaller than that difference is not real.\x22}\n" as *u8)
          return 0
      }
  }
  tb_puts("{\x22organ\x22:\x22nx_twinbench\x22,\x22usage\x22:\x22nx_twinbench cmp <cand.nxmesh> <ref.nxmesh> [stature_mm] | nx_twinbench selftest\x22}\n" as *u8)
  return 0
672}
nx_twinbench.nx source

`nx_twinbench.nx` source
