`nx_skullsdf.nx` source

↩ module page · 2774 lines · 146229 B
1// nx_skullsdf.nx -- the skull as a FUSED SHELL WITH OPENINGS, built in the implicit path.
2//
3// ★WHY THIS EXISTS, from two render cycles of evidence rather than preference. nx_skullgen emits the
4// skull as seven LOFTED TUBES. Rendered, that is seven separate objects sharing a bounding box: the
5// nasal bones float clear of the face, the zygomatic arch is a rod passing through nothing, and the
6// braincase is a balloon with a jaw hanging near it. Moving the tubes closer helps a little and cannot
7// fix the class of problem, because A SKULL IS NOT A SET OF TUBES -- it is one continuous fused shell
8// with holes cut in it. Tubes cannot fuse and cannot have holes.
9//
10// ★THE IMPLICIT PATH DOES EXACTLY THOSE TWO THINGS AND WE ALREADY OWN IT: smooth-union joins bones into
11// one surface, subtraction cuts the orbits and the nasal aperture as real openings, and nx_isosurf's
12// surface_nets polygonizes the result into nx_trimesh. The plan reserved SDF for viscera on the grounds
13// that it cannot resolve an eyelash -- true, and irrelevant to bone, whose features are centimetres.
14// ⚠So this does NOT replace the loft path for skin. It claims the skull only.
15//
16//   nx_skullsdf <out.nxmesh> [sex 0-1000] [robust 0-1000] [grid]
17// license_tier: ORIGINAL  expect_exit: 0  No hw writes (Rule 26).
18import "nx_isosurf.nx"
19import "nx_itoa_lib.nx"   // shared MSB-first emitter (zero-alloc)
20import "nx_memfloor.nx"   // mf_admit: refuse an O(G^3) grid that would eat the host (2026-07-30 incident)
21import "nx_sdfprim_lib.nx" // the implicit primitives (ellipsoid, capsule, arch, plate, walled cylinder, ring): ONE family for every anatomy field
22import "nx_fieldfit_lib.nx" // the fit machinery (two-sided region-balanced objective, coordinate descent, prefix fit file): ONE optimiser for every field
23const SS_MAGIC_2026: i64 = 2026
24const SS_MAGIC_1024: i64 = 1024
25const SS_MAGIC_8388608: i64 = 8388608
26const SS_MAGIC_8388607: i64 = 8388607
27const SS_MAGIC_40000: i64 = 40000
28
29const SS_G: i64 = 160             // grid cells per axis
30const SS_EXT: i64 = 120           // world half-extent in MILLIMETRES. A male skull is ~185mm long and
                                // ~145mm broad, so 260mm of grid covers it with no wasted resolution
                                // (the first version sampled 600mm and spent most of the grid on air)
33// ★THE PRIMITIVES ARE AUTHORED IN AN ARBITRARY ~90-UNIT SPACE, NOT MILLIMETRES. Assuming otherwise
34// produced a skull measuring 84x75x93mm against the reference's 151x209x215 -- roughly half scale and
35// wrongly proportioned. This factor maps the authored space onto real millimetres, and the three
36// reference dimensions are named so the gate can check the result against them instead of my eye.
37// Derived from the measurement, not guessed: at 1850 the emitted skull stood 348mm against a 209mm
38// reference, so the factor is 1850 x 209/348. Each correction here is now driven by the bbox the organ
39// prints, which is the whole reason that check exists.
40// ★IDENTITY. The field below is now authored in REAL MILLIMETRES against the reference's measured box,
41// so there is no unit-space to get wrong. The previous arbitrary ~90-unit space needed a fudge factor
42// that could correct overall size but never proportion -- the skull came out 176x256x210 where the
43// reference is 151x209x215: too tall AND too shallow, which no single multiplier can fix.
44const SS_SCALE: i64 = 1000
45const SS_REF_W: i64 = 151         // reference bounding box, mm, measured off BodyParts3D: x lateral
46const SS_REF_H: i64 = 215         // HEIGHT = the file's z extent (BodyParts3D is a floor-origin z-UP body frame)
47const SS_REF_D: i64 = 209         // DEPTH = the file's y extent (antero-posterior, face at -y). Until 2026-09-18 the two
                                // were swapped, because every reader below took the file as y-up: see ss_ref_map
49const SS_TOL: i64 = 25            // mm each dimension may differ before the build is called wrong
50const SS_Q: i64 = SP_Q            // fixed point for the distance functions
51const SS_NQ: i64 = 4096           // the scale surface_nets emits gradient normals in (nx_isosurf.nx:68)
52const SS_ISO: i64 = 0
53const SS_BONE: i64 = 0xe8e0d8
54// ***REDISTANCE (argv[8]) -- REBUILD |grad f| = 1 ON THE SAMPLED GRID, the fix this file's own tissue
55// diagnosis names (see the DISABLED gradient-normalisation block): chained ss_smin leaves the field up
56// to 40 percent off a true SDF near the surface, so an iso-offset of depth d moves the skin by
57// d/|grad| instead of d (measured 3-6mm scatter), and dividing by |grad| downstream AMPLIFIES flat
58// cells (measured 62mm chins). The standard remedy is an eikonal pass over the grid we already build:
59// SEED every node adjacent to a zero crossing with its exact sub-cell distance to that crossing (so
60// the bone surface CANNOT move -- linear interpolation between two seeds reproduces the crossing at
61// exactly the raw field's own position), then FAST-SWEEP the upwind eikonal update in all 8 diagonal
62// orders until distances are true everywhere. The redistanced field is carried at SS_RDQ fixed point
63// so sub-cell seeds survive integer storage; surface_nets only ever interpolates zero crossings and
64// normalises gradients, so a uniformly scaled field polygonizes identically.
65// ADDITIVE: argv[8] absent or 0 leaves the old path byte-identical.
66const SS_RDQ: i64 = 256           // fixed-point units per millimetre in the redistanced field
67const SS_RD_INF: i64 = 76800      // 300mm x SS_RDQ -- farther than any point in the sampled box
68// ---- M2 SUB-MM RAW FIELD (debt 1786406737, 2026-08-12): the raw grid stored INTEGER-MILLIMETRE
69// distances, so every iso-crossing snapped to coarse fractions of a cell and the redistance seeds
70// inherited those mm crossings -- the onion-ring terraces that SURVIVE redistancing (proven 08-10:
71// "redistance seeds inherit the RAW field's integer-mm zero crossings; the only cure is fixed-point
72// evaluation of the primitive field itself"). The _q lanes below evaluate the SAME authored-mm
73// geometry but carry the RESULT at SS_FQ units per millimetre. SS_FQ == SS_RDQ so the raw and
74// redistanced grids share ONE fixed point and every ratio-based consumer (seed crossing t=af/(af+ag),
75// sign restore, surface_nets at SS_ISO=0) is scale-invariant by construction -- no downstream change.
76// The mm lane (ss_field/ss_field_p/ss_field_world) stays for the tuner/bench/grad verbs unchanged.
77const SS_FQ: i64 = SP_FQ
78const SS_RD_NEAR: i64 = 3         // near-band half-width, cells, for the printed |grad| A/B
79const SS_RD_LO: i64 = 800         // |grad| per-mille below this counts under (same bar as ss_grad)
80const SS_RD_HI: i64 = 1200        // |grad| per-mille above this counts over (same bar as ss_grad)
81static SS_REDIST: i64
82// ★★SMOOTH-UNION RADIUS -- AND WHY IT IS NOW ZERO. ss_smin subtracts h^2/(4k) from the result, so every
83// chained call pushes the field further negative. Chaining eleven of them means the value at a point
84// depends on how many bones are near it, and the result is NOT a distance field -- gradient magnitude
85// drifts far from 1. surface_nets places each vertex by LINEAR INTERPOLATION along a cell edge, which
86// assumes the field is locally linear in distance; feed it a distorted field and vertices land in
87// inconsistent places, giving a jagged self-intersecting shell that renders as speckle.
88// ★Plain min() is the exact union of two SDFs and cannot distort anything. It is also anatomically
89// right: bones meet at SUTURES, which are sharp, not blended.
90// Restored to a real millimetre blend. At 0 (plain union) the mandible and face floated as separate
91// lumps with visible seams; bones do meet at sutures, but a 6mm blend at 2mm sampling is what makes
92// adjacent masses read as ONE bone rather than a pile of touching balls.
93// ★2026-08-10: default 6 -> 3, the value every measured run shipped with today (redistancing repairs
94// the field distortion smin causes, so the old reason to hold 6 is gone; 3 keeps fusion w/ sharper form)
95const SS_K: i64 = 3
96// ***SS_K IS NOW RUNTIME-SETTABLE (argv[7]) SO THE SMOOTH-UNION/DISTANCE-FIDELITY TRADE CAN BE A/B'd
97// WITHOUT A REBUILD. This file already records setting K to ZERO once, for exactly the reason that
98// matters here -- plain min() is the EXACT union of two SDFs and cannot distort anything, whereas every
99// chained ss_smin subtracts h*h/(4k) and pushes the field away from a true distance function. K=0 gives
100// honest distances with visible sutures; K=6 gives one fused shell with a distorted field. Until now that
101// trade was a recompile, so nobody measured it. Default is unchanged, so every existing caller is
102// byte-identical.
103static SS_KSET: i64
104static SS_KVAL: i64
105func ss_k() -> i64 { if SS_KSET == 1 { return SS_KVAL } return SS_K }
106
107func ss_puts(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
108// MIGRATED to nxi_out (debt 1785563586). The old body mmapped 32 bytes per call and never freed it,
109// across 33 call sites. THIS ORGAN ALREADY IMPORTS nx_memfloor mf_admit after the 2026-07-30
110// incident where it took 27.7GB of 36GB and froze every seat -- the GUARD was added, the CAUSE
111// was not. nxi_out is MSB-first, allocates NOTHING, and emits identical bytes including the sign.
112func ss_pn(v: i64) -> i64 { nxi_out(v); return 0 }
113func ss_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
118}
119func ss_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
123}
124func ss_isqrt(v: i64) -> i64 { return sp_isqrt(v) }
125func ss_abs(v: i64) -> i64 { if v<0 { return 0-v } return v }
126func ss_min(a: i64, b: i64) -> i64 { if a<b { return a } return b }
127func ss_max(a: i64, b: i64) -> i64 { if a>b { return a } return b }
128
129// ---- primitives. All return a signed distance in world units: negative inside.
130// An ellipsoid's exact SDF has no closed form; the standard cheap bound (scaled radial distance times
131// the smallest semi-axis) is used, which is what every SDF modeller uses and is accurate near the
132// surface -- the only place the polygonizer looks.
133// ★★TRUE ELLIPSOID DISTANCE, not a radial bound. The previous version returned
134// (radial-1) x smallest-semi-axis, which is only correct on a sphere and drifts badly off-axis. Alone
135// that is invisible; SMOOTH-UNIONED a dozen times it produced grid-scale ripples in the field, and
136// surface nets faithfully polygonized the ripples into a lumpy self-intersecting shell that z-fought
137// against itself. ★The error was in the PRIMITIVE and only became visible under COMPOSITION.
138// This is the standard first-order approximation d = k0(k0-1)/k1, exact on a sphere and accurate near
139// the surface of an ellipsoid -- which is the only region the polygonizer samples.
140// integer Bhaskara sine: degrees in, per-mille out; exact at 0/90/180, <2 permil error between.
141// The estate's emitter already carries one for the body; the field needs its own because a
142// STRUCTURED primitive places instances around an arc and a table lookup cannot be tuned smoothly.
143// ---- rule-11 named constants: each is a DOMAIN FACT, not a dial (2026-08-14) ----
144// Bhaskara I's 7th-century sine approximation: sin(x) ~ 4x(180-x) / (40500 - x(180-x)).
145// The numerator's 4 is scaled by 1000 so the result is PER-MILLE, matching ss_isin's contract
146// above; 40500 = 5 * 90^2 belongs to the closed form itself and is not tunable.
147const SS_BHASKARA_NUM_PERMILLE: i64 = 4000
148const SS_BHASKARA_DEN: i64 = 40500
149// Trilinear interpolation carries fractions as 0..SS_LERP_ONE (Q10 fixed point).
150const SS_LERP_ONE: i64 = 1024
151// Per-mille denominator. A UNIT, not a dial: sexf/robust arrive as 0..1000 and the policy is the
152// NUMERATOR (`3*sexf/SS_PERMILLE` says three per-mille per unit of sex). Hoisted ONLY at the 23
153// sites proven by enumeration to be divisors -- 47 of the 70 occurrences of 1000 are NOT (spans,
154// coordinates, argv), and folding those in would move the skull, not rename a unit.
155const SS_PERMILLE: i64 = 1000
156// Min/max accumulators start at +/- this and are replaced by the first real sample. A
157// SENTINEL, not a limit -- no value is ever clamped to it.
158const SS_EXTENT_SENTINEL: i64 = 999999
159// |grad SDF| must be ~1000 per-mille (unit) for a true signed-distance field; the probe counts
160// samples outside a +/-20% band around unit. LO and HI are ONE fact in two halves -- nx_magic's
161// 1024 threshold flags only HI, so naming HI alone would leave the pair half-told.
162const SS_GRAD_BAND_LO: i64 = 800
163const SS_GRAD_BAND_HI: i64 = 1200
164// Marching-cubes cell budget: the argv[7] default and the floor it is clamped up to.
165const SS_MARCH_CELLS_DEFAULT: i64 = 2000
166const SS_MARCH_CELLS_MIN: i64 = 1000
167// ---- THE SKULL FIELD. Every dimension is a function of the sex/robusticity knobs, as in nx_skullgen;
168// the difference is purely that these compose by union and subtraction instead of standing apart.
169// ★AUTHORED IN MILLIMETRES, in a frame centred on the cranium: +y up, +z forward (toward the face),
170// x lateral. Every figure is a real skull dimension, so the emitted bounding box can be compared
171// directly against the reference's 151 x 209 x 215 instead of through a fudge factor.
172// ★★THE SKULL AS A PARAMETER VECTOR. Hand-placing coordinates got the bounding box to 5mm and left the
173// SHAPE at 8mm mean error, because a person adjusting numbers by eye cannot search twenty dimensions.
174// Exposing the dimensions as data lets the fit metric DRIVE them instead: `tune` runs coordinate
175// descent against the real scanned reference and keeps only changes that lower the measured error.
176const SP_N: i64 = 87
177const SS_FIT_MIN: i64 = 58     // the v1 vector length: a fit file carrying at least this many ints
                             // applies as a PREFIX (structured params added later keep their
                             // defaults until the next tune saves the full vector) -- growing the
                             // parametrization must never discard banked convergence
181const P_VW: i64 = 0            // vault half-width
182const P_VH: i64 = 1            // vault half-height
183const P_VD: i64 = 2            // vault half-depth
184const P_VY: i64 = 3            // vault centre y
185const P_VZ: i64 = 4            // vault centre z
186const P_FY: i64 = 5            // frontal centre y
187const P_FZ: i64 = 6            // frontal centre z
188const P_FH: i64 = 7            // frontal half-height
189const P_FD: i64 = 8            // frontal half-depth
190const P_MY: i64 = 9            // maxilla centre y
191const P_MZ: i64 = 10           // maxilla centre z
192const P_MW: i64 = 11           // maxilla half-width
193const P_MH: i64 = 12           // maxilla half-height
194const P_MD: i64 = 13           // maxilla half-depth
195const P_JY: i64 = 14           // mandible chin y
196const P_JZ: i64 = 15           // mandible chin z
197const P_RX: i64 = 16           // ramus x offset
198const P_RY: i64 = 17           // ramus centre y
199const P_OR: i64 = 18           // orbit radius
200const P_OX: i64 = 19           // orbit centre x
201// ***THE REST OF THE SKULL'S STATIONS, PARAMETERISED (2026-08-09, operator: "the tooling should be
202// able to see pixel by pixel why our emitter is emitting shit compared to the benchmark and get us
203// to pixel to pixel emission -- not manual"). Coordinate descent can only move what is a parameter;
204// everything below was a frozen literal, which is why the measured-wrong stations (orbit depth,
205// nasal position, arch endpoints) could never converge no matter how good the objective was.
206const P_BX: i64 = 20           // brow bar half-span
207const P_BY: i64 = 21           // brow bar y
208const P_BZ: i64 = 22           // brow bar z
209const P_ZAX: i64 = 23          // zygomatic arch face-end x
210const P_ZAY: i64 = 24          // zygomatic arch face-end y
211const P_ZAZ: i64 = 25          // zygomatic arch face-end z
212const P_ZBX: i64 = 26          // zygomatic arch ear-end x
213const P_ZBY: i64 = 27          // zygomatic arch ear-end y
214const P_ZBZ: i64 = 28          // zygomatic arch ear-end z
215const P_ZR: i64 = 29           // zygomatic arch radius
216const P_NY: i64 = 30           // nasal bones centre y
217const P_NZ: i64 = 31           // nasal bones centre z
218const P_NRX: i64 = 32          // nasal radius x
219const P_NRY: i64 = 33          // nasal radius y
220const P_NRZ: i64 = 34          // nasal radius z
221const P_JAX: i64 = 35          // jaw body chin-end x
222const P_JAZ: i64 = 36          // jaw body chin-end z offset from P_JZ
223const P_JBX: i64 = 37          // jaw body angle-end x
224const P_JBY: i64 = 38          // jaw body angle-end y offset from P_JY
225const P_JBZ: i64 = 39          // jaw body angle-end z
226const P_JR: i64 = 40           // jaw body radius
227const P_RZ: i64 = 41           // ramus centre z
228const P_RRX: i64 = 42          // ramus radius x
229const P_RRY: i64 = 43          // ramus radius y
230const P_RRZ: i64 = 44          // ramus radius z
231const P_OAY: i64 = 45          // orbit opening y
232const P_OAZ: i64 = 46          // orbit opening z
233const P_OBX: i64 = 47          // orbit apex x
234const P_OBY: i64 = 48          // orbit apex y
235const P_OBZ: i64 = 49          // orbit apex z
236const P_PAY: i64 = 50          // nasal aperture opening y
237const P_PAZ: i64 = 51          // nasal aperture opening z
238const P_PBY: i64 = 52          // nasal aperture back y
239const P_PBZ: i64 = 53          // nasal aperture back z
240const P_PR: i64 = 54           // nasal aperture radius
241const P_CRX: i64 = 55          // chin ellipsoid radius x
242const P_CRY: i64 = 56          // chin ellipsoid radius y
243const P_CRZ: i64 = 57          // chin ellipsoid radius z
244// ***THE FIRST STRUCTURED PRIMITIVE: A REPEATED ELEMENT ALONG A CURVE (operator 2026-08-09: "you
245// cant even make teeth in a jaw... how are you going to emit a tree"). A flat list of blobs cannot
246// express teeth; ~9 knobs per ARCH drive N instances of one tooth rule along an elliptical arc.
247// This construct IS the seed of tree emission: a tree is elements-along-a-curve RECURSED.
248const P_UY: i64 = 58           // upper arch: occlusal centre y
249const P_UZ: i64 = 59           // upper arch: arc centre z
250const P_UAX: i64 = 60          // upper arch: lateral semi-axis
251const P_UAZ: i64 = 61          // upper arch: anterior semi-axis
252const P_UTH: i64 = 62          // upper arch: half arc span, degrees
253const P_UN: i64 = 63           // upper arch: tooth count
254const P_UTW: i64 = 64          // upper arch: tooth width
255const P_UTHT: i64 = 65         // upper arch: tooth height
256const P_UPR: i64 = 66          // upper arch: molar widening, permil at the arc ends
257const P_LY: i64 = 67           // lower arch: occlusal centre y
258const P_LZ: i64 = 68           // lower arch: arc centre z
259const P_LAX: i64 = 69          // lower arch: lateral semi-axis
260const P_LAZ: i64 = 70          // lower arch: anterior semi-axis
261const P_LTH: i64 = 71          // lower arch: half arc span, degrees
262const P_LN: i64 = 72           // lower arch: tooth count
263const P_LTW: i64 = 73          // lower arch: tooth width
264const P_LTHT: i64 = 74         // lower arch: tooth height
265const P_LPR: i64 = 75          // lower arch: molar widening, permil at the arc ends
266// ---- AN17 primitive rungs (2026-09-19): the openings as SHAPES the ruler can read, appended so a 76-int fit file still
267// loads as a prefix. RETIRED stations (kept in the vector for file compatibility, skipped by the tuner, unread by the
268// field): P_OR P_OBX P_OBY P_OBZ (the orbit was a drilled capsule to an apex), P_PBY P_PBZ P_PR (the aperture was a
269// capsule tunnel), P_NRZ (the nasal bones are a plate with P_NTH), P_JAX (the body is one arc through the chin),
270// P_RX P_RY P_RZ P_RRY (the ramus stands between the body's angle and the arch's ear end, derived, never placed).
271const P_OSW: i64 = 76          // orbit socket semi-width x (the rim's lateral half-span)
272const P_OSH: i64 = 77          // orbit socket semi-height y
273const P_OSD: i64 = 78          // orbit socket semi-depth z (the back wall sits at P_OCZ - P_OSD, inside the bone)
274const P_OCZ: i64 = 79          // orbit socket centre z; the front pokes through the face at P_OCZ + P_OSD
275const P_PAW: i64 = 80          // piriform aperture width (rim to rim)
276const P_PAH: i64 = 81          // piriform aperture height (rhinion to nasospinale)
277const P_PAD: i64 = 82          // piriform aperture depth: the nasal floor wall sits P_PAD behind the rim at P_PAZ
278const P_NTILT: i64 = 83        // nasal plate pitch, degrees: the rhinion end forward of the nasion end
279const P_NTH: i64 = 84          // nasal plate half-thickness
280const P_AAR: i64 = 85          // alveolar process radius: the U-arch of bone the upper teeth stand in
281const P_FPR: i64 = 86          // maxillary frontal process radius (the bone either side of the aperture, up to the orbits)
282// also RETIRED: P_PAY P_PAZ -- the aperture hangs from the nasal plate's lower edge (the rhinion), derived in the field
283// the tuner's population: 1 for a station the field reads, 0 for a retired one (an evaluation spent on a dead
284// parameter is a pass that cannot move anything)
285func ss_param_live(k: i64) -> i64 {
  if k == P_PAY { return 0 }
  if k == P_PAZ { return 0 }
  if k == P_OR { return 0 }
  if k == P_OBX { return 0 }
  if k == P_OBY { return 0 }
  if k == P_OBZ { return 0 }
  if k == P_PBY { return 0 }
  if k == P_PBZ { return 0 }
  if k == P_PR { return 0 }
  if k == P_NRZ { return 0 }
  if k == P_JAX { return 0 }
  if k == P_RX { return 0 }
  if k == P_RY { return 0 }
  if k == P_RZ { return 0 }
  if k == P_RRY { return 0 }
  return 1
302}
303// ★FITTED, NOT TYPED. These are the output of `tune` -- coordinate descent against the 171k-triangle
304// scanned reference, which took the mean surface error from 8.37mm to 5.06mm in one run. Every value
305// here was accepted only because it LOWERED a measured number, which is the opposite of how the
306// previous eight rounds of hand-placed coordinates were chosen.
307func ss_defaults(P: *i64) -> i64 {
  P[P_VW]=81; P[P_VH]=71; P[P_VD]=75; P[P_VY]=24;  P[P_VZ]=0-33
  P[P_FY]=46; P[P_FZ]=15; P[P_FH]=34; P[P_FD]=33
  P[P_MY]=0-45; P[P_MZ]=58; P[P_MW]=30; P[P_MH]=13; P[P_MD]=53
  P[P_JY]=0-78; P[P_JZ]=60; P[P_RX]=57; P[P_RY]=0-46
  P[P_OR]=8; P[P_OX]=30
  // every remaining primitive station, extracted verbatim from the field's former literals so the
  // TUNER can reach them -- the orbit stations, nasal position and arch endpoints were frozen
  // literals no optimizer could move, which is why the gamed 8mm orbit could never recover
  P[P_BX]=52; P[P_BY]=22; P[P_BZ]=62
  P[P_ZAX]=52; P[P_ZAY]=6; P[P_ZAZ]=52; P[P_ZBX]=62; P[P_ZBY]=8; P[P_ZBZ]=0-30; P[P_ZR]=7
  P[P_NY]=8; P[P_NZ]=76; P[P_NRX]=9; P[P_NRY]=18; P[P_NRZ]=12
  P[P_JAX]=22; P[P_JAZ]=0-4; P[P_JBX]=46; P[P_JBY]=10; P[P_JBZ]=10; P[P_JR]=17
  P[P_RZ]=0-6; P[P_RRX]=9; P[P_RRY]=34; P[P_RRZ]=26
  P[P_OAY]=18; P[P_OAZ]=80; P[P_OBX]=20; P[P_OBY]=12; P[P_OBZ]=18
  P[P_PAY]=0-8; P[P_PAZ]=92; P[P_PBY]=0-2; P[P_PBZ]=40; P[P_PR]=13
  P[P_CRX]=26; P[P_CRY]=18; P[P_CRZ]=22
  // DENTAL ARCHES: placement DERIVES from the parent bones (never type a placement the parent
  // already knows -- the maxilla and jaw rows above are the parents); the structural dims (arch
  // radii, tooth size/count, published human dental anatomy) are STARTING values the tuner owns.
  P[P_UY]=P[P_MY]-24; P[P_UZ]=P[P_MZ]-10; P[P_UAX]=26; P[P_UAZ]=30
  P[P_UTH]=75; P[P_UN]=14; P[P_UTW]=9; P[P_UTHT]=11; P[P_UPR]=400
  P[P_LY]=P[P_JY]+16; P[P_LZ]=P[P_JZ]-12; P[P_LAX]=25; P[P_LAZ]=28
  P[P_LTH]=75; P[P_LN]=14; P[P_LTW]=9; P[P_LTHT]=11; P[P_LPR]=400
  // AN17 primitive rungs: the openings START at the oracle's own ruler readings (nx_skull_canon_gate, 2026-09-18, mm:
  // aperture 23.1 wide, 22.4 high, 10.7 deep; orbit width 28.3) and the tuner owns them from there
  P[P_OSW]=16; P[P_OSH]=14; P[P_OSD]=22; P[P_OCZ]=P[P_OAZ]-8
  P[P_PAW]=23; P[P_PAH]=22; P[P_PAD]=11
  P[P_NTILT]=25; P[P_NTH]=3; P[P_AAR]=8; P[P_FPR]=6
  return 0
337}
338
339// the millimetre lane is the SS_FQ lane divided down: ONE definition of the skull's geometry (until 2026-09-19 the two
340// lanes carried the same twenty-two primitives twice, and a primitive changed in one and not the other is half a skull)
341func ss_field_p(px: i64, py: i64, pz: i64, sexf: i64, robust: i64, P: *i64) -> i64 {
  return ss_field_p_q(px, py, pz, sexf, robust, P)/SS_FQ
343}
344// ss_field_p at SS_FQ output: same geometry, same composition order, sub-mm result. ss_smin/ss_sub
345// are scale-invariant so only the primitive outputs and the smoothing radius change units.
346func ss_field_p_q(px: i64, py: i64, pz: i64, sexf: i64, robust: i64, P: *i64) -> i64 {
  let vw: i64 = P[P_VW] - 3*sexf/SS_PERMILLE
  let vh: i64 = P[P_VH] - 2*sexf/SS_PERMILLE
  let vd: i64 = P[P_VD] - 2*sexf/SS_PERMILLE
  let brow: i64 = 10 - 3*sexf/SS_PERMILLE + 5*robust/SS_PERMILLE
  let kq: i64 = ss_k()*SS_FQ
  var d: i64 = sp_ellipsoid_q(px,py,pz, 0, P[P_VY], P[P_VZ], vw, vh, vd)
  d = sp_smin(d, sp_ellipsoid_q(px,py,pz, 0, P[P_FY], P[P_FZ], vw*78/100, P[P_FH], P[P_FD]), kq)
  d = sp_smin(d, sp_capsule_q(px,py,pz, 0-P[P_BX],P[P_BY],P[P_BZ], P[P_BX],P[P_BY],P[P_BZ], brow), kq)
  d = sp_smin(d, sp_ellipsoid_q(px,py,pz, 0, P[P_MY], P[P_MZ], P[P_MW], P[P_MH], P[P_MD]), kq)
  // (d) the alveolar process: a U-arch of bone standing on the upper dental arc; its front at the midline is the prosthion
  d = sp_smin(d, sp_arch_q(px,py,pz, P[P_UY]+P[P_UTHT]/2+P[P_AAR]*2/3, P[P_UZ], P[P_UAX], P[P_UAZ], P[P_UTH], 0, P[P_AAR]), kq)
  d = sp_smin(d, sp_ring_q(px,py,pz, P[P_UY], P[P_UZ], P[P_UAX], P[P_UAZ], P[P_UTH], P[P_UN], P[P_UTW], P[P_UTHT], P[P_UPR]), kq)
  d = sp_smin(d, sp_ring_q(px,py,pz, P[P_LY], P[P_LZ], P[P_LAX], P[P_LAZ], P[P_LTH], P[P_LN], P[P_LTW], P[P_LTHT], P[P_LPR]), kq)
  d = sp_smin(d, sp_capsule_q(px,py,pz, 0-P[P_ZAX],P[P_ZAY],P[P_ZAZ], 0-P[P_ZBX],P[P_ZBY],P[P_ZBZ], P[P_ZR]), kq)
  d = sp_smin(d, sp_capsule_q(px,py,pz,   P[P_ZAX],P[P_ZAY],P[P_ZAZ],   P[P_ZBX],P[P_ZBY],P[P_ZBZ], P[P_ZR]), kq)
  // (c) the nasal bones as a pitched PLATE from the nasion down and forward to the rhinion, never a peg
  d = sp_smin(d, sp_ellipsoid_tiltx_q(px,py,pz, 0, P[P_NY], P[P_NZ], P[P_NRX], P[P_NRY], P[P_NTH], P[P_NTILT]), kq)
  d = sp_smin(d, sp_ellipsoid_q(px,py,pz, 0, P[P_JY], P[P_JZ], P[P_CRX], P[P_CRY], P[P_CRZ]), kq)
  // (e) the mandibular body as ONE arc from angle to angle through the chin, rising P_JBY toward the angles
  d = sp_smin(d, sp_arch_q(px,py,pz, P[P_JY], P[P_JBZ], P[P_JBX], P[P_JZ]+P[P_JAZ]-P[P_JBZ], 90, P[P_JBY], P[P_JR]), kq)
  // the ramus stands between the body's angle end and the arch's ear end (the condyle under the arch root): DERIVED
  // from the two structures it joins, so the angle and the condyle cannot drift apart
  let cdy: i64 = P[P_ZBY] - 2*P[P_ZR]
  let rcy: i64 = (P[P_JY]+P[P_JBY] + cdy)/2
  let rcz: i64 = (P[P_JBZ] + P[P_ZBZ])/2
  let rcx: i64 = (P[P_JBX] + P[P_ZBX])/2
  var rhy: i64 = (cdy - (P[P_JY]+P[P_JBY]))/2 + P[P_RRX]
  if rhy < P[P_RRX] { rhy = P[P_RRX] }
  d = sp_smin(d, sp_ellipsoid_q(px,py,pz, 0-rcx, rcy, rcz, P[P_RRX], rhy, P[P_RRZ]), kq)
  d = sp_smin(d, sp_ellipsoid_q(px,py,pz,   rcx, rcy, rcz, P[P_RRX], rhy, P[P_RRZ]), kq)
  // the rhinion: the nasal plate's lower edge, where the aperture's top rim is by definition. The aperture is placed
  // FROM it (a typed aperture station drifted 20 mm in front of the bone the moment the peg became a plate: an
  // opening carved into air is not an opening)
  let nts: i64 = sp_isin(P[P_NTILT])
  let ntc: i64 = sp_icos(P[P_NTILT])
  let yrh: i64 = P[P_NY] - P[P_NRY]*ntc/SS_PERMILLE
  let zrh: i64 = P[P_NZ] + P[P_NRY]*nts/SS_PERMILLE
  // the maxilla's FRONTAL PROCESSES: the bone either side of the aperture, rising from its lateral rims to the medial
  // orbital rims, so the opening has rims to be walked and the maxilla joins the brow instead of floating under it
  let fpx: i64 = P[P_PAW]/2 + P[P_FPR]
  let fpy: i64 = yrh - P[P_PAH]/2
  let fpz: i64 = zrh - P[P_FPR]
  let fqx: i64 = P[P_OX] - P[P_OSW] - P[P_FPR]/2
  let fqy: i64 = P[P_OAY] - P[P_OSH]/2
  d = sp_smin(d, sp_capsule_q(px,py,pz, 0-fpx,fpy,fpz, 0-fqx,fqy,P[P_OCZ], P[P_FPR]), kq)
  d = sp_smin(d, sp_capsule_q(px,py,pz,   fpx,fpy,fpz,   fqx,fqy,P[P_OCZ], P[P_FPR]), kq)
  // (b) the orbits as SOCKETS: an ellipsoid carved with its back wall inside the bone and its front through the face
  d = sp_sub(d, sp_ellipsoid_q(px,py,pz, 0-P[P_OX], P[P_OAY], P[P_OCZ], P[P_OSW], P[P_OSH], P[P_OSD]))
  d = sp_sub(d, sp_ellipsoid_q(px,py,pz,   P[P_OX], P[P_OAY], P[P_OCZ], P[P_OSW], P[P_OSH], P[P_OSD]))
  // (c) the piriform aperture as an OPENING hanging from the rhinion, its nasal floor P_PAD behind the rim
  d = sp_sub(d, sp_cylz_q(px,py,pz, 0, yrh - P[P_PAH]/2, P[P_PAW]/2, P[P_PAH]/2, zrh - P[P_PAD]))
  return d
399}
400
401// THE DEFAULT PARAMETER VECTOR IS MAPPED ONCE, NOT ONCE PER SAMPLE -- the organ's single largest
402// defect, not a tuning matter. ss_field is the hot path of the whole program: the polygonizer
403// evaluates it at every one of (G+1)^3 grid corners. It used to sys_mmap(SP_N*8) on EVERY call, and
404// sys_mmap is PAGE-GRANULAR -- a 160-byte request costs a whole 4096-byte page, nothing here frees,
405// so THE FIELD LEAKED ONE PAGE PER SAMPLE.
406// MEASURED, and it reproduces the 2026-07-30 incident to the digit: G=224 is 225^3 = 11,390,625
407// samples x 4096 B = 46.6 GB of anonymous pages; the box OOM-killed it at total-vm 33.8 GB /
408// anon-rss 28.0 GB partway through. The grid array at that size is only 91 MB and cv another 90 MB,
409// so 2,966 bytes per cell was never the grid -- it was this line. It is also why G=1500 mapped 27 GB
410// and froze every other seat's builds.
411// THE BANKED DIAGNOSIS -- 'the knob allocates on the RAW argument while computing on a CAPPED one' --
412// WAS WRONG, and the memory arithmetic refutes it: the raw-argument allocation is O(G^3)x8 = 91 MB at
413// G=224, three orders of magnitude short of the observed kill. The cost was per-EVALUATION, not
414// per-cell-of-storage.
415// A static POINTER lazily mapped once is the sanctioned shape here (a static ARRAY in BSS is the
416// banked crash-on-startup hazard, and this is the idiom nx_fluid/nx_genome/nx_ghash already use). The
417// vector is immutable after ss_defaults, so ONE copy is correct for every caller and ss_field stays a
418// pure function of its arguments -- the property the tuner and the polygonizer both rely on.
419// ---- THE FIT FILE: the tuner's output is DATA the emitter loads, never a source edit. If
420// knowledge/skull_fit.dat exists and carries a COMPLETE SP_N integer vector it overrides the
421// builtin defaults; short or absent files change nothing (parsed into a temp, applied only when
422// complete, so a truncated write cannot half-poison the skull). `tune ... <savepath>` writes it,
423// every emit reports fit_src, and the optimize->emit->measure loop closes with no hand in the numbers.
424static SS_FITSRC: i64
425static SS_FITN: i64
426// the fit file through the ONE loader (nx_fieldfit_lib.ff_load_prefix): a prefix of at least SS_FIT_MIN ints applies,
427// anything shorter changes nothing; SS_FITN carries how many stations the file governed for the emit receipt
428func ss_load_fit(P: *i64) -> i64 {
  let got: i64 = ff_load_prefix("knowledge/skull_fit.dat" as *u8, P, SP_N, SS_FIT_MIN)
  if got >= SS_FIT_MIN { SS_FITN = got; return 1 }
  return 0
432}
433static SS_PDEF: i64
434func ss_pdef() -> *i64 {
  if SS_PDEF == 0 {
      SS_PDEF = sys_mmap(SP_N*8) as i64
      ss_defaults(SS_PDEF as *i64)
      SS_FITSRC = ss_load_fit(SS_PDEF as *i64)
  }
  return SS_PDEF as *i64
441}
442func ss_field(px: i64, py: i64, pz: i64, sexf: i64, robust: i64) -> i64 {
  // delegates to the parameterised field with the default vector, so every existing caller is
  // unchanged and there is exactly ONE definition of the skull's geometry
  return ss_field_p(px, py, pz, sexf, robust, ss_pdef())
446}
447
448// Evaluate the authored field at a WORLD millimetre point: map the point back into authored space,
449// then scale the returned distance forward again so the field stays a true SDF in world units.
450func ss_field_world(px: i64, py: i64, pz: i64, sexf: i64, robust: i64) -> i64 {
  let d: i64 = ss_field(px*1000/SS_SCALE, py*1000/SS_SCALE, pz*1000/SS_SCALE, sexf, robust)
  return d * SS_SCALE / 1000
453}
454// world-space field at SS_FQ units per world-unit -- the sub-mm raw lane (see the SS_FQ block)
455func ss_field_world_q(px: i64, py: i64, pz: i64, sexf: i64, robust: i64) -> i64 {
  let d: i64 = ss_field_p_q(px*1000/SS_SCALE, py*1000/SS_SCALE, pz*1000/SS_SCALE, sexf, robust, ss_pdef())
  return d * SS_SCALE / 1000
458}
459
460// ★★★SKIN AS A TRUE ENVELOPE -- THE ANSWER THE OFFSET FAILURE HANDED US. nx_skinwrap displaced every
461// bone vertex outward along its own normal by the forensic tissue depth, and it LOWERED the binding judge
462// 99 -> 88 against the cadaver skin. The reason was structural, not tuning: AN OFFSET ALONG A NORMAL
463// CANNOT BRIDGE A CAVITY, and a face is largely defined by the cavities its skin bridges -- the orbits,
464// the nasal aperture, the temporal fossa. A per-vertex offset reproduces every opening and pushes its rim
465// outward, so a socket becomes a bulging crater. That is the exact inverse of this lane's banked law that
466// relief modulates a radius and cannot CUT a hole: a normal offset cannot FILL one.
467// ★THE IMPLICIT PATH DOES IT FOR FREE, and that is why this belongs here and not in a new organ. The skull
468// is already a distance FIELD. Subtracting a depth from the distance before polygonizing moves the zero
469// level OUTWARD by that depth, and AN ISO-SURFACE AT DISTANCE d NATURALLY BRIDGES ANY OPENING NARROWER
470// THAN 2d -- which is precisely the eyelid-over-orbit and skin-over-aperture behaviour we want, arising
471// from the representation rather than from a special case.
472// ★DEPTH IS KEYED ON THE FIELD'S OWN AUTHORED STATIONS, not on a per-mille of some bounding box, because
473// the field already places every landmark in millimetres in a frame centred on the cranium. These are the
474// published forensic depths (Rhine/De Greef class) this codebase carries in knowledge/face_tissue.dat.
475// ⚠HONEST v1 LIMIT, declared here rather than discovered later: this varies with HEIGHT only. The temple
476// is 3mm while the cheekbone at nearly the same height is 7, so a height-keyed table cannot express a
477// lateral difference. The nasion dip at 3mm between brow 7 and nasal 6 IS expressible, and it is the one
478// that makes a nose read as a nose.
479static SS_TISSUE: i64
480// ★★★DENSE-ER TISSUE MAP: DEPTH IS A FUNCTION OF (HEIGHT, LATERAL DISTANCE), NOT HEIGHT ALONE.
481// ★WHY, from the 2026 literature rather than from taste: the state of the art in anatomy-guided face
482// generation (Skull-to-Face, arXiv 2403.16207; SCULPTOR skeleton-consistent parametric generator; and the
483// PLOS One automatic forensic reconstruction from DENSE STATISTICS OF SOFT TISSUE THICKNESS) does not use
484// a sparse landmark list -- it fits a DENSE MAP of tissue thickness over the skull surface. Our v1 was a
485// 9-band table keyed on height only, and I declared its limit when I shipped it: the temple is ~3mm while
486// the zygion sits at ~7mm at nearly the SAME HEIGHT, so a 1-D table structurally cannot express either
487// without corrupting the other. Adding the lateral axis is the smallest honest step toward a dense map,
488// and it is the axis that carries the largest published disagreement.
489// ★THE LATERAL STRUCTURE IS REAL ANATOMY, NOT SMOOTHING. Two rows invert with distance from the midline
490// and they are the reason this matters: at the BROW LINE the midline glabella is thick (7) while the
491// TEMPLE is the thinnest soft tissue on the whole head (3); at the JAW the midline chin is 11 but the
492// GONION is thicker still going outward (11) because the MASSETER inserts there. A height-only table gets
493// the temple and the gonion wrong in OPPOSITE directions, which is exactly the kind of error that shows up
494// as a wrong silhouette rather than a wrong number.
495// ⚠STILL NOT A DENSE MAP: this is 9 height bands x 3 lateral bands = 27 cells, against a literature
496// standard of a per-vertex statistical field. It is a strict superset of v1's information, not parity.
497const SS_LAT_MID: i64 = 25     // |x| under this is the midline column
498const SS_LAT_FAR: i64 = 55     // |x| over this is the far-lateral column (temple / zygomatic arch / gonion)
499func ss_tdepth2(y: i64, ax: i64) -> i64 {
  var col: i64 = 0
  if ax >= SS_LAT_MID { col = 1 }
  if ax >= SS_LAT_FAR { col = 2 }
  if y < 0-78 { if col==0 { return 11 } if col==1 { return 10 } return 8 }    // chin / mental
  if y < 0-60 { if col==0 { return 9 }  if col==1 { return 10 } return 11 }   // jaw -- THICKENS outward (masseter)
  if y < 0-38 { if col==0 { return 11 } if col==1 { return 9 }  return 8 }    // lower lip
  if y < 0-30 { if col==0 { return 10 } if col==1 { return 8 }  return 7 }    // upper lip
  if y < 0    { if col==0 { return 7 }  if col==1 { return 7 }  return 6 }    // maxilla / zygion
  if y < 8    { if col==0 { return 6 }  if col==1 { return 7 }  return 6 }    // nasal
  if y < 22   { if col==0 { return 3 }  if col==1 { return 6 }  return 5 }    // nasion DIP on the midline only
  if y < 30   { if col==0 { return 7 }  if col==1 { return 6 }  return 4 }    // brow: glabella thick, temple thin
  if col==0 { return 5 }
  if col==1 { return 4 }
  return 3                                                                    // temple -- thinnest on the head
514}
515// v1 kept as the midline column so the 1-D behaviour is still reachable and comparable
516func ss_tdepth(y: i64) -> i64 { return ss_tdepth2(y, 0) }
517// ***TABLE VALUES ARE LANDMARK MEASUREMENTS, NOT STEP FUNCTIONS. On the distorted field the band
518// steps were accidentally BLURRED by the field error; the redistanced field reproduces them
519// faithfully, as visible horizontal shelves at every band boundary (rendered proof:
520// /world/mh_skin_env_redist.png). The honest reading of a forensic depth table is depth AT a
521// landmark station, interpolated BETWEEN stations. Only the REDIST path consumes this, so the
522// shipped band behaviour is byte-identical.
523const SS_TD_N: i64 = 9
524static SS_TD_Y: i64
525static SS_TD_D: i64
526func ss_td_tab() -> i64 {
  if SS_TD_Y == 0 {
      SS_TD_Y = sys_mmap(SS_TD_N*8) as i64
      SS_TD_D = sys_mmap(SS_TD_N*8) as i64
      let Y: *i64 = SS_TD_Y as *i64
      let D: *i64 = SS_TD_D as *i64
      Y[0]=0-88; D[0]=11
      Y[1]=0-69; D[1]=9
      Y[2]=0-49; D[2]=11
      Y[3]=0-34; D[3]=10
      Y[4]=0-15; D[4]=7
      Y[5]=4;    D[5]=6
      Y[6]=15;   D[6]=3
      Y[7]=26;   D[7]=7
      Y[8]=45;   D[8]=5
  }
  return 0
543}
544// depth at height y in mm x SS_RDQ, piecewise-linear between the landmark stations (band centres of
545// ss_tdepth, same forensic values), clamped at both ends
546func ss_tdepth_q(y: i64) -> i64 {
  ss_td_tab()
  let Y: *i64 = SS_TD_Y as *i64
  let D: *i64 = SS_TD_D as *i64
  if y <= Y[0] { return D[0]*SS_RDQ }
  if y >= Y[SS_TD_N-1] { return D[SS_TD_N-1]*SS_RDQ }
  var i: i64 = 0
  while i < SS_TD_N-1 {
      if y <= Y[i+1] {
          let c0: i64 = Y[i]
          let c1: i64 = Y[i+1]
          return (D[i]*(c1-y) + D[i+1]*(y-c0))*SS_RDQ/(c1-c0)
      }
      i = i + 1
  }
  return D[SS_TD_N-1]*SS_RDQ
562}
563func ss_f32(v: i64, scale: i64) -> i64 {
  // integer -> IEEE-754 single bits, value = v/scale
  if v == 0 { return 0 }
  var neg: i64 = 0
  var m: i64 = v
  if m < 0 { neg = 1; m = 0 - m }
  var e: i64 = 127
  // normalise m/scale into [1,2)
  var num: i64 = m
  var den: i64 = scale
  while num >= den*2 { den = den*2; e = e + 1 }
  while num < den { num = num*2; e = e - 1 }
  // mantissa = (num/den - 1) * 2^23
  let frac: i64 = ((num - den) * SS_MAGIC_8388608) / den
  var bits: i64 = (e << 23) | (frac & SS_MAGIC_8388607)
  if neg == 1 { bits = bits | (1 << 31) }
  return bits
580}
581func ss_wr32(b: *u8, o: i64, v: i64) -> i64 {
  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
  return 0
587}
588
589// ***REALISED THICKNESS -- THE OBJECTIVE THE DENSE-MAP LITERATURE ACTUALLY NEEDS, AND IT IS FREE HERE.
590// The 2-D table failed because writing depths into a FIELD is not the same as getting them on the
591// SURFACE: surface_nets redistributes them. So the thing to build is not a better table, it is a
592// MEASUREMENT OF WHAT THE SURFACE ACTUALLY REALISED. Our skull IS a distance function, so |field| at an
593// envelope vertex IS that vertex's realised tissue thickness over the bone -- no ray cast, no nearest-
594// triangle search. Same trick ss_fit already uses against the scanned reference.
595// Simultaneously the diagnostic (how much does the field-to-surface step redistribute a requested
596// depth) and the objective a coordinate descent would minimise.
597func ss_bandix(y: i64) -> i64 {
  if y < 0-78 { return 0 }
  if y < 0-60 { return 1 }
  if y < 0-38 { return 2 }
  if y < 0-30 { return 3 }
  if y < 0    { return 4 }
  if y < 8    { return 5 }
  if y < 22   { return 6 }
  if y < 30   { return 7 }
  return 8
607}
608func ss_band_req(ix: i64) -> i64 {
  if ix == 0 { return 11 }
  if ix == 1 { return 9 }
  if ix == 2 { return 11 }
  if ix == 3 { return 10 }
  if ix == 4 { return 7 }
  if ix == 5 { return 6 }
  if ix == 6 { return 3 }
  if ix == 7 { return 7 }
  return 5
618}
619func ss_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)
621}
622func ss_r_f32(bits: i64, scale: i64) -> i64 {
  let s: i64 = (bits >> 31) & 1
  let e: i64 = (bits >> 23) & 255
  let m: i64 = bits & SS_MAGIC_8388607
  if e == 0 { return 0 }
  if e == 255 { return 0 }
  let mant: i64 = SS_MAGIC_8388608 | m
  let sh: i64 = e - 127 - 23
  var v: i64 = 0
  if sh >= 0 { if sh > 30 { return 0 } v = (mant*scale) << sh }
  else { let rs: i64 = 0 - sh; if rs > 62 { return 0 } v = (mant*scale) >> rs }
  if s == 1 { return 0 - v }
  return v
635}
636// ★SCORE AGAINST THE REAL SCANNED REFERENCE, IN MILLIMETRES -- and the metric is FREE because our skull
637// IS a distance function: |field| at a reference vertex is the error there. No registration search, no
638// closest-triangle hunt. Eight rounds of hand-placing primitives matched the bounding box to 5mm while
639// the SHAPE stayed a cluster of lumps; a bounding box cannot see shape and an eyeball cannot be
640// optimized against. This can.
641// ⚠FRAME: the reference faces -z in a full-body frame, ours faces +z centred on the cranium. The
642// transform is DERIVED from the reference's own bbox, and z is negated -- without that we would be
643// scoring our face against the back of its skull.
644// THE ORACLE'S FRAME, READ ONCE, CORRECTLY (anatomy AN17, 2026-09-18). knowledge/skull.nxmesh is BodyParts3D, and
645// BodyParts3D is a FLOOR-ORIGIN, z-UP body frame: the head sits at z 1.42..1.64 m, y runs antero-posterior with the
646// face at -y, x lateral (first triangle: x 18.7, y -156.9, z 1457.4 mm; bbox x 151, y 209, z 215 mm). Every reader in
647// this file took it as y-up and merely negated z, so the tuner fitted an upright field to a skull LYING ON ITS BACK
648// since the first tune: the trend's "vault 4 mm" was measured against that pose and knowledge/skull_fit.dat converged
649// to the egg that best matches a supine skull (vault half-height 92, a negative nasal radius). Measured by the AN17
650// canon ruler's midline profile dump: read as x, z, -y the profile from the front IS a skull (forehead receding 110 mm
651// over 126 mm to a glabella, an aperture floor 26 mm behind it, a chin 53 mm below the prosthion); read as y-up it is
652// the skull base seen from below, its 36 mm see-through the foramen magnum.
653// x_a = x - cx ; y_a = z - cz ; z_a = -(y - cy) -- a proper rotation of a quarter turn about x, applied to positions
654// (out3, integer mm at scale units) and to normals (ss_ref_map_normal). ONE mapper for every reader below.
655func ss_ref_map(buf: *u8, o: i64, cx: i64, cy: i64, cz: i64, scale: i64, out3: *i64) -> i64 {
  out3[0] = ss_r_f32(ss_rd32(buf, o), scale) - cx
  out3[1] = ss_r_f32(ss_rd32(buf, o+8), scale) - cz
  out3[2] = 0 - (ss_r_f32(ss_rd32(buf, o+4), scale) - cy)
  return 0
660}
661func ss_ref_map_normal(buf: *u8, o: i64, scale: i64, out3: *i64) -> i64 {
  out3[0] = ss_r_f32(ss_rd32(buf, o), scale)
  out3[1] = ss_r_f32(ss_rd32(buf, o+8), scale)
  out3[2] = 0 - ss_r_f32(ss_rd32(buf, o+4), scale)
  return 0
666}
667func ss_fit(path: *u8, sexf: i64, robust: i64, stride: i64) -> i64 {
  let ln: *i64 = sys_mmap(16) as *i64
  let buf: *u8 = sys_read_file(path, ln)
  if (buf as i64) == 0 { ss_puts("{\x22error\x22:\x22cannot read reference\x22}\n" as *u8); return 3 }
  let nlayer: i64 = ss_rd32(buf, 8)
  let nt: i64 = ss_rd32(buf, 12)
  let hdr: i64 = 16 + nlayer*24
  if nt <= 0 { ss_puts("{\x22error\x22:\x22no triangles\x22}\n" as *u8); return 4 }
  var lox: i64 = 0; var hix: i64 = 0
  var loy: i64 = 0; var hiy: i64 = 0
  var loz: i64 = 0; var hiz: i64 = 0
  var first: i64 = 1
  var t: i64 = 0
  while t < nt {
      let o: i64 = hdr + t*84
      var j: i64 = 0
      while j < 3 {
          let x: i64 = ss_r_f32(ss_rd32(buf,o+j*12),   1)
          let y: i64 = ss_r_f32(ss_rd32(buf,o+j*12+4), 1)
          let z: i64 = ss_r_f32(ss_rd32(buf,o+j*12+8), 1)
          if first == 1 { lox=x; hix=x; loy=y; hiy=y; loz=z; hiz=z; first=0 } else {
              if x<lox {lox=x}
              if x>hix {hix=x}
              if y<loy {loy=y}
              if y>hiy {hiy=y}
              if z<loz {loz=z}
              if z>hiz {hiz=z}
          }
          j = j + 1
      }
      t = t + 1
  }
  let cx: i64 = (lox+hix)/2
  let cy: i64 = (loy+hiy)/2
  let cz: i64 = (loz+hiz)/2
  var sum: i64 = 0
  var n: i64 = 0
  var worst: i64 = 0
  var over: i64 = 0
  t = 0
  while t < nt {
      if t % stride == 0 {
          let o: i64 = hdr + t*84
          let m3: *i64 = sys_mmap(24) as *i64
          ss_ref_map(buf, o, cx, cy, cz, 1, m3)
          let x: i64 = m3[0]
          let y: i64 = m3[1]
          let z: i64 = m3[2]
          var d: i64 = ss_field(x, y, z, sexf, robust)
          if d < 0 { d = 0 - d }
          sum = sum + d
          n = n + 1
          if d > worst { worst = d }
          if d > 10 { over = over + 1 }
      }
      t = t + 1
  }
  if n == 0 { ss_puts("{\x22error\x22:\x22no samples\x22}\n" as *u8); return 5 }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22fit\x22,\x22mean_err_mm\x22:" as *u8); ss_pn(sum/n)
  ss_puts(",\x22worst_err_mm\x22:" as *u8); ss_pn(worst)
  ss_puts(",\x22over_10mm_permil\x22:" as *u8); ss_pn(over*1000/n)
  ss_puts(",\x22samples\x22:" as *u8); ss_pn(n)
  ss_puts(",\x22ref_bbox_mm\x22:[" as *u8); ss_pn(hix-lox)
  ss_puts("," as *u8); ss_pn(hiy-loy); ss_puts("," as *u8); ss_pn(hiz-loz)
  ss_puts("]}\n" as *u8)
  return 0
733}
734
735// load reference vertices once into flat arrays, already transformed into our frame
736func ss_load_ref(path: *u8, RX: *i64, RY: *i64, RZ: *i64, cap: i64, stride: i64) -> i64 {
  let ln: *i64 = sys_mmap(16) as *i64
  let buf: *u8 = sys_read_file(path, ln)
  if (buf as i64) == 0 { return 0-1 }
  let nlayer: i64 = ss_rd32(buf, 8)
  let nt: i64 = ss_rd32(buf, 12)
  let hdr: i64 = 16 + nlayer*24
  var lox: i64 = 0; var hix: i64 = 0
  var loy: i64 = 0; var hiy: i64 = 0
  var loz: i64 = 0; var hiz: i64 = 0
  var first: i64 = 1
  var t: i64 = 0
  while t < nt {
      let o: i64 = hdr + t*84
      // all three vertices: the gapmap and voxref centre by the same bbox, and a centre read off vertex 0 alone
      // (as this loop did until 2026-09-18) registered the two readers a fraction of a cell apart
      var vj: i64 = 0
      while vj < 3 {
          let x: i64 = ss_r_f32(ss_rd32(buf,o+vj*12),   1)
          let y: i64 = ss_r_f32(ss_rd32(buf,o+vj*12+4), 1)
          let z: i64 = ss_r_f32(ss_rd32(buf,o+vj*12+8), 1)
          if first == 1 { lox=x; hix=x; loy=y; hiy=y; loz=z; hiz=z; first=0 } else {
              if x<lox {lox=x}
              if x>hix {hix=x}
              if y<loy {loy=y}
              if y>hiy {hiy=y}
              if z<loz {loz=z}
              if z>hiz {hiz=z}
          }
          vj = vj + 1
      }
      t = t + 1
  }
  let cx: i64 = (lox+hix)/2
  let cy: i64 = (loy+hiy)/2
  let cz: i64 = (loz+hiz)/2
  var n: i64 = 0
  let m3: *i64 = sys_mmap(24) as *i64
  t = 0
  while t < nt {
      if t % stride == 0 {
          if n < cap {
              let o: i64 = hdr + t*84
              ss_ref_map(buf, o, cx, cy, cz, 1, m3)
              RX[n] = m3[0]
              RY[n] = m3[1]
              RZ[n] = m3[2]
              n = n + 1
          }
      }
      t = t + 1
  }
  return n
789}
790// ★★THE OBJECTIVE IS SYMMETRIC, AND IT HAS TO BE. A one-directional score -- "how far is each REFERENCE
791// point from my surface" -- is trivially gamed by DELETING GEOMETRY: the first tuning run shrank the
792// orbits from 22mm to 8mm and flattened the maxilla from 30mm to 13mm, improving the number 40% by
793// filling in the eye sockets, because surface that does not exist cannot be wrong.
794// The reverse term fixes that by construction: sample points on OUR surface and measure how far they
795// are from the REFERENCE. A filled orbit puts our surface where the reference has a hole, and the
796// reverse term charges for it. Together they are a two-sided (Hausdorff-style) distance.
797// ⚠We have no distance function for the reference MESH, so the reverse term uses the nearest reference
798// VERTEX from a spatial grid. That is an upper bound on the true point-to-surface distance -- it can
799// overstate, never understate, which is the safe direction for a penalty.
800// the oracle box the reverse term's spatial hash covers: the scanned skull's own extent, in the field's frame
801const SG_DIM: i64 = 24            // spatial hash grid, cells per axis
802const SG_EXT: i64 = 140           // half-extent covered, mm
803const SS_TUNE_LAT: i64 = 6        // reverse-term lattice spacing, mm (the objective's sampling of OUR surface)
804const SS_GR_N: i64 = 7            // the named regions of ss_gap_region (vault, base, brow-orbits, nasal, zygomatic, mandible, maxilla)
805// ★THE LANDMARK TERM (AN17, 2026-09-19). MEASURED FIRST: with the two surface terms alone, region-balanced, the descent
806// closed the orbit sockets to 2 mm slits and deepened the aperture to a 35 mm tunnel while the objective fell 445 -> 413,
807// because a surface objective cannot see an opening. This term is the third error mean the fit library folds in: where
808// the field PLACES its features (from its own stations, closed-form) against where the ruler FOUND the oracle's
809// landmarks (knowledge/status/skull_oracle_q.dat, written by nx_skull_canon_gate from its own measurement, never typed
810// here). It is a SURROGATE of the ruler on our mesh -- the beat still runs the real ruler on the emitted mesh and gates
811// the pass on its in-band count -- and its predicted quantities are printed so the surrogate's error is readable beside
812// the ruler's reading. Quantities the ruler seeds from its own conf (orbit depth at the eye seed, the vault chord pitch)
813// are not predicted. Absent file = term OFF, announced.
814const SS_LMQ_N: i64 = 18                  // the ruler's quantity count, in skc_qname order (15 + the three of 2026-09-19)
815const SS_LMQ_PRED: i64 = 16               // the quantities this field can place from its stations
816const SS_LMQ_BYTES: i64 = 4096
817const SS_LMQ_NONE: i64 = 0 - 999999999    // "this quantity is not predicted from the stations" -- a value no landmark can reach
818static SS_LMQ: i64                        // *i64: oracle_u10 per quantity, or 0 when the file is absent
819static SS_LM_SEXF: i64
820static SS_LM_ROBUST: i64
821func ss_lm_load() -> i64 {
  if SS_LMQ != 0 { return 1 }
  let ln: *i64 = sys_mmap(16) as *i64
  let buf: *u8 = sys_read_file("knowledge/status/skull_oracle_q.dat" as *u8, ln)
  if (buf as i64) == 0 { return 0 }
  let n: i64 = ln[0]
  let T: *i64 = sys_mmap(SS_LMQ_N*8) as *i64
  // each row: <name> <oracle_u10> <band_u10>; the second field of each line is the target. Flag-exited scans: a loop
  // that clobbers its cursor to break cannot report where it stopped (the estate's cursor-sentinel law)
  var i: i64 = 0
  var row: i64 = 0
  while i < n {
      var go: i64 = 1
      while go == 1 { if i >= n { go = 0 } else { if buf[i] == (32 as u8) { go = 0 } else { i = i + 1 } } }
      if i < n { i = i + 1 }
      var neg: i64 = 0
      if i < n { if buf[i] == (45 as u8) { neg = 1; i = i + 1 } }
      var v: i64 = 0
      var any: i64 = 0
      go = 1
      while go == 1 {
          if i >= n { go = 0 } else {
              let c: i64 = buf[i] as i64
              if c >= 48 { if c <= 57 { v = v*10 + (c-48); any = 1; i = i + 1 } else { go = 0 } } else { go = 0 }
          }
      }
      if any == 1 { if row < SS_LMQ_N { if neg == 1 { v = 0 - v } T[row] = v; row = row + 1 } }
      go = 1
      while go == 1 { if i >= n { go = 0 } else { if buf[i] == (10 as u8) { go = 0 } else { i = i + 1 } } }
      if i < n { i = i + 1 }
  }
  if row < SS_LMQ_N { return 0 }
  SS_LMQ = T as i64
  return 1
855}
856// the field's own placement of the ruler's quantities, u10, from the stations (index = skc_qname order; -1 = not predicted)
857func ss_lm_predict(P: *i64, Q: *i64) -> i64 {
  let brow: i64 = 10 - 3*SS_LM_SEXF/SS_PERMILLE + 5*SS_LM_ROBUST/SS_PERMILLE
  // the ruler's glabella is the MOST ANTERIOR midline point above the aperture floor: the brow bar's front or the
  // frontal ellipsoid's front, whichever leads (measured 2026-09-19: the surrogate read the bar, the ruler the bulge)
  var glab_z: i64 = P[P_BZ] + brow
  if P[P_FZ] + P[P_FD] > glab_z { glab_z = P[P_FZ] + P[P_FD] }
  let glab_y: i64 = P[P_BY]
  let nts: i64 = sp_isin(P[P_NTILT])
  let ntc: i64 = sp_icos(P[P_NTILT])
  let nasion_y: i64 = P[P_NY] + P[P_NRY]*ntc/SS_PERMILLE
  let yrh: i64 = P[P_NY] - P[P_NRY]*ntc/SS_PERMILLE
  let nasosp_y: i64 = yrh - P[P_PAH]
  let pog_z: i64 = P[P_JZ] + P[P_CRZ]
  let gnath_y: i64 = P[P_JY] - P[P_CRY]
  let pros_z: i64 = P[P_UZ] + P[P_UAZ] + P[P_AAR]
  let pros_y: i64 = P[P_UY] + P[P_UTHT]/2 + P[P_AAR]*2/3
  let k: i64 = ss_k()
  Q[0] = (pog_z - glab_z)*10
  Q[1] = (nasion_y - nasosp_y)*10
  Q[2] = P[P_PAH]*10
  Q[3] = P[P_PAW]*10                       // the ruler's crest walk stops at the cylinder's own rim, not the process crest (measured)
  Q[4] = 2*P[P_OX]*10
  Q[5] = (glab_y - P[P_OAY])*10
  Q[6] = 2*(P[P_OSW] + k)*10
  Q[7] = 2*(P[P_ZBX] + P[P_ZR])*10
  Q[8] = 2*(P[P_JBX] + P[P_JR])*10
  Q[9] = (pros_z - glab_z)*10
  Q[10] = (nasion_y - pros_y)*10
  Q[11] = ((glab_y - nasosp_y) - (nasosp_y - gnath_y))*10
  Q[12] = P[P_PAD]*10
  Q[13] = SS_LMQ_NONE
  Q[14] = SS_LMQ_NONE
  // the three the ruler grew for the tuner's blind spots: the socket's rim height, the dip at the nasion (the plate's
  // front face at its top end against the glabella), the nasal bones' projection (the plate's slope times its length)
  let nasion_z: i64 = P[P_NZ] - P[P_NRY]*nts/SS_PERMILLE + P[P_NTH]
  let rhinion_z: i64 = P[P_NZ] + P[P_NRY]*nts/SS_PERMILLE + P[P_NTH]
  Q[15] = 2*(P[P_OSH] + k)*10
  Q[16] = (glab_z - nasion_z)*10
  Q[17] = (rhinion_z - nasion_z)*10
  return SS_LMQ_PRED
897}
898// the third error mean, centi-mm: mean |predicted - oracle| over the predicted quantities. The scratch vector is mapped
899// ONCE (this runs per descent move; a page per call is the leak this file measured in ss_field)
900static SS_LMSCR: i64
901func ss_lm_term(P: *i64) -> i64 {
  if SS_LMQ == 0 { return 0 }
  let O: *i64 = SS_LMQ as *i64
  if SS_LMSCR == 0 { SS_LMSCR = sys_mmap(SS_LMQ_N*8) as i64 }
  let Q: *i64 = SS_LMSCR as *i64
  ss_lm_predict(P, Q)
  var sum: i64 = 0
  var cnt: i64 = 0
  var q: i64 = 0
  while q < SS_LMQ_N {
      if Q[q] != SS_LMQ_NONE {
          var d: i64 = Q[q] - O[q]
          if d < 0 { d = 0 - d }
          sum = sum + d
          cnt = cnt + 1
      }
      q = q + 1
  }
  if cnt == 0 { return 0 }
  return sum*10/cnt
921}
922// ★THE TUNE COMPOSES nx_fieldfit_lib (AN17, 2026-09-19): the two-sided region-balanced objective, the coordinate descent
923// and the fit file are the library's; this organ supplies its field, its region classifier and its live-station mask
924// as function values, and keeps only what is the skull's -- the oracle reader, the parameter vector, the receipt.
925func ss_tune(path: *u8, sexf: i64, robust: i64, stride: i64, passes: i64, step0: i64, save: *u8) -> i64 {
  let cap: i64 = SS_MAGIC_40000
  let C: *FfCtx = ff_ctx_new(SP_N, SS_GR_N, cap, SG_EXT, SG_DIM, SS_TUNE_LAT, SS_FQ)
  C.field = ss_field_p_q
  C.region = ss_gap_region
  C.live = ss_param_live
  C.a = sexf
  C.b = robust
  SS_LM_SEXF = sexf
  SS_LM_ROBUST = robust
  C.extra = ss_lm_term
  C.has_extra = ss_lm_load()
  let n: i64 = ss_load_ref(path, C.RX, C.RY, C.RZ, cap, stride)
  if n <= 0 { ss_puts("{\x22error\x22:\x22cannot load reference\x22}\n" as *u8); return 3 }
  C.n = n
  ff_bin(C)
  ff_regions(C)
  let P: *i64 = sys_mmap(SP_N*8) as *i64
  ss_defaults(P)
  // iterate FROM the current fit so successive tune runs converge instead of restarting
  let seeded: i64 = ss_load_fit(P)
  let start: i64 = ff_score2(C, P)
  let best: i64 = ff_descend(C, P, passes, step0, 1)
  // persist the fit as DATA (the emitter loads it; no source edit carries these numbers)
  if (save as i64) != 0 { ff_save(save, P, SP_N) }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22tune\x22,\x22samples\x22:" as *u8); ss_pn(n)
  ss_puts(",\x22start_err_mm\x22:" as *u8); ss_pn(start/100)
  ss_puts(",\x22final_err_mm\x22:" as *u8); ss_pn(best/100)
  ss_puts(",\x22start_centi\x22:" as *u8); ss_pn(start)
  ss_puts(",\x22final_centi\x22:" as *u8); ss_pn(best)
  ss_puts(",\x22objective\x22:\x22region-balanced\x22,\x22regions\x22:" as *u8); ss_pn(SS_GR_N)
  ss_puts(",\x22live_params\x22:" as *u8); ss_pn(ff_live_count(C))
  ss_puts(",\x22seeded\x22:" as *u8); ss_pn(seeded)
  ss_puts(",\x22landmark_term\x22:" as *u8); ss_pn(C.has_extra)
  if C.has_extra == 1 {
      // the surrogate's own prediction, u10, in the ruler's quantity order (-1 = not predicted): read it beside the
      // canon gate's `q ... ours=` line to know the surrogate's error
      let LQ: *i64 = sys_mmap(SS_LMQ_N*8) as *i64
      ss_lm_predict(P, LQ)
      ss_puts(",\x22landmark_term_centi\x22:" as *u8); ss_pn(ss_lm_term(P))
      ss_puts(",\x22predicted_u10\x22:[" as *u8)
      var lq: i64 = 0
      while lq < SS_LMQ_N {
          if lq > 0 { ss_puts("," as *u8) }
          ss_pn(LQ[lq])
          lq = lq + 1
      }
      ss_puts("]" as *u8)
  }
  ss_puts(",\x22params\x22:[" as *u8)
  var q: i64 = 0
  while q < SP_N {
      if q > 0 { ss_puts("," as *u8) }
      ss_pn(P[q])
      q = q + 1
  }
  ss_puts("]}\n" as *u8)
  return 0
983}
984
985// THE `thick` VERB: report, per height band, the depth we REQUESTED versus the thickness the emitted
986// envelope ACTUALLY realised over the bone. A gap between them is the field-to-surface redistribution
987// made visible, and its mean absolute error is the number a fit would drive to zero.
988func ss_thick(path: *u8, sexf: i64, robust: i64, stride: i64) -> i64 {
  let ln: *i64 = sys_mmap(16) as *i64
  let buf: *u8 = sys_read_file(path, ln)
  if (buf as i64) == 0 { ss_puts("{\x22error\x22:\x22cannot read envelope mesh\x22}\n" as *u8); return 3 }
  let nlayer: i64 = ss_rd32(buf, 8)
  let nt: i64 = ss_rd32(buf, 12)
  if nt <= 0 { ss_puts("{\x22error\x22:\x22no triangles\x22}\n" as *u8); return 4 }
  let hdr: i64 = 16 + nlayer*24
  let cnt: *i64 = sys_mmap(9*8) as *i64
  let sm:  *i64 = sys_mmap(9*8) as *i64
  let mn:  *i64 = sys_mmap(9*8) as *i64
  let mx:  *i64 = sys_mmap(9*8) as *i64
  var b: i64 = 0
  while b < 9 { cnt[b]=0; sm[b]=0; mn[b]=SS_EXTENT_SENTINEL; mx[b]=0-SS_EXTENT_SENTINEL; b=b+1 }
  var abserr: i64 = 0
  var n: i64 = 0
  var t: i64 = 0
  while t < nt {
      if t % stride == 0 {
          var j: i64 = 0
          while j < 3 {
              let o: i64 = hdr + t*84 + j*12
              let px: i64 = ss_r_f32(ss_rd32(buf,o),   1)
              let py: i64 = ss_r_f32(ss_rd32(buf,o+4), 1)
              let pz: i64 = ss_r_f32(ss_rd32(buf,o+8), 1)
              let d: i64 = ss_field(px, py, pz, sexf, robust)
              let ix: i64 = ss_bandix(py)
              cnt[ix] = cnt[ix] + 1
              sm[ix] = sm[ix] + d
              if d < mn[ix] { mn[ix] = d }
              if d > mx[ix] { mx[ix] = d }
              var e: i64 = d - ss_band_req(ix)
              if e < 0 { e = 0 - e }
              abserr = abserr + e
              n = n + 1
              j = j + 1
          }
      }
      t = t + 1
  }
  if n == 0 { ss_puts("{\x22error\x22:\x22no samples\x22}\n" as *u8); return 5 }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22thick\x22,\x22samples\x22:" as *u8); ss_pn(n)
  ss_puts(",\x22mean_abs_err_mm\x22:" as *u8); ss_pn(abserr/n)
  ss_puts(",\x22bands\x22:[" as *u8)
  b = 0
  while b < 9 {
      if b > 0 { ss_puts("," as *u8) }
      ss_puts("{\x22band\x22:" as *u8); ss_pn(b)
      ss_puts(",\x22req\x22:" as *u8); ss_pn(ss_band_req(b))
      ss_puts(",\x22mean\x22:" as *u8)
      if cnt[b] > 0 { ss_pn(sm[b]/cnt[b]) } else { ss_pn(0-1) }
      ss_puts(",\x22min\x22:" as *u8)
      if cnt[b] > 0 { ss_pn(mn[b]) } else { ss_pn(0-1) }
      ss_puts(",\x22max\x22:" as *u8)
      if cnt[b] > 0 { ss_pn(mx[b]) } else { ss_pn(0-1) }
      ss_puts(",\x22n\x22:" as *u8); ss_pn(cnt[b])
      ss_puts("}" as *u8)
      b = b + 1
  }
  ss_puts("]}\n" as *u8)
  return 0
1049}
1050
1051// ***THE `grad` VERB -- MEASURE |grad f| DIRECTLY INSTEAD OF INFERRING IT FROM SCATTER.
1052// Every remaining hypothesis about the tissue-thickness scatter is a claim about how far this field is
1053// from a TRUE distance function, where |grad f| = 1 everywhere. Rather than argue, sample it: walk a
1054// coarse lattice, keep the points NEAR THE SURFACE (the only region that matters for an offset), and
1055// report the distribution of |grad f| in per-mille. A true SDF reads 1000.
1056// ***READ IT LIKE THIS: spread around 1000 is the offset error budget. If a requested 7mm offset lands
1057// where |grad| is 700, the surface moves 10mm; at 1400 it moves 5mm. That IS the 3-6mm scatter, and this
1058// verb says so directly instead of by inference. It also discriminates the two surviving suspects --
1059// ss_smin composition (already exonerated by K=0) versus ss_ellipsoid's first-order approximation --
1060// because you can re-run it with argv[5]=0 for a plain union and compare the SAME distribution.
1061func ss_grad(sexf: i64, robust: i64, step: i64, near: i64) -> i64 {
  var lo: i64 = SS_EXTENT_SENTINEL
  var hi: i64 = 0-SS_EXTENT_SENTINEL
  var sum: i64 = 0
  var n: i64 = 0
  var under: i64 = 0
  var over: i64 = 0
  var x: i64 = 0 - SS_EXT
  while x <= SS_EXT {
      var y: i64 = 0 - SS_EXT
      while y <= SS_EXT {
          var z: i64 = 0 - SS_EXT
          while z <= SS_EXT {
              let f: i64 = ss_field(x, y, z, sexf, robust)
              var af: i64 = f
              if af < 0 { af = 0 - af }
              if af <= near {
                  let gx: i64 = (ss_field(x+1,y,z,sexf,robust) - ss_field(x-1,y,z,sexf,robust)) * 1000 / 2
                  let gy: i64 = (ss_field(x,y+1,z,sexf,robust) - ss_field(x,y-1,z,sexf,robust)) * 1000 / 2
                  let gz: i64 = (ss_field(x,y,z+1,sexf,robust) - ss_field(x,y,z-1,sexf,robust)) * 1000 / 2
                  let g: i64 = ss_isqrt(gx*gx + gy*gy + gz*gz)
                  if g < lo { lo = g }
                  if g > hi { hi = g }
                  sum = sum + g
                  n = n + 1
                  if g < SS_GRAD_BAND_LO { under = under + 1 }
                  if g > SS_GRAD_BAND_HI { over = over + 1 }
              }
              z = z + step
          }
          y = y + step
      }
      x = x + step
  }
  if n == 0 { ss_puts("{\x22error\x22:\x22no near-surface samples\x22}\n" as *u8); return 5 }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22grad\x22,\x22samples\x22:" as *u8); ss_pn(n)
  ss_puts(",\x22true_sdf_is\x22:1000,\x22mean\x22:" as *u8); ss_pn(sum/n)
  ss_puts(",\x22min\x22:" as *u8); ss_pn(lo)
  ss_puts(",\x22max\x22:" as *u8); ss_pn(hi)
  ss_puts(",\x22under_800_permil\x22:" as *u8); ss_pn(under*1000/n)
  ss_puts(",\x22over_1200_permil\x22:" as *u8); ss_pn(over*1000/n)
  ss_puts(",\x22smooth_union_k\x22:" as *u8); ss_pn(ss_k())
  ss_puts("}\n" as *u8)
  return 0
1105}
1106
1107// ---- GAP MAP (the `gapmap` verb): WHERE and HOW MUCH we deviate from the scanned oracle ----
1108// The redistance pass makes OUR field a true SDF; the same sweep machinery run over the ORACLE'S
1109// VERTICES gives the ORACLE'S distance field on the SAME grid. Comparing the two turns "ours is
1110// wrong" into a RANKED PER-REGION WORKLIST (which primitive to fix first, by measured millimetres)
1111// and TWO HEAT MESHES the eye can read: our surface coloured by distance-to-oracle (where our
1112// geometry is wrong) and the oracle's surface coloured by distance-to-ours (what real anatomy we
1113// are MISSING -- sockets, teeth, arches glow). The oracle stays a RULER: nothing here feeds oracle
1114// geometry back into the generator.
1115// ⚠DECLARED BOUND: the oracle field is seeded from VERTICES, not triangles, so it OVERSTATES
1116// distance by at most half the oracle's vertex spacing (~1mm on the 171k-tri skull) -- safe
1117// direction for a gap report, and printed in the JSON.
1118// ***THE RESIDUAL LAYER (operator 2026-08-09: "not until you actually can emit the benchmark skull
1119// can you claim any of this"). The parametric base B is the CONTROLLABLE layer and converges near
1120// ~4.5mm -- blobs cannot carry identity. The industry decomposition is base + DISPLACEMENT: gapmap
1121// already builds B and the ingested input's own SDF O on the same grid, so the ideal correction is
1122// their pointwise difference, downsampled to a lattice and stored as DATA. Emission = B + trilinear(R)
1123// -- the input-derived residual carries identity, the base stays controllable, and a scale knob
1124// slides between the generic parametric skull (0) and the ingested benchmark (1000). R is DERIVED
1125// FROM THE INGESTED INPUT by construction: that IS the match half of the ingest->match->exceed
1126// circle, and it involves no hand-typed number anywhere.
1127const SS_RES_MAGIC: i64 = 0x4e585253      // 'NXRS'
1128const SS_RES_MAXN: i64 = 256              // lattice nodes per axis, upper bound (256^3 x 8 = 134MB;
                                        // raised 128->256 2026-08-10: a 2mm-lattice residual over a
                                        // 1mm grid blurs sub-2mm bone walls into speckle -- MEASURED
                                        // as the ghost-shell emission; identity needs lattice = grid)
1132static SS_RESP: i64                        // loaded residual values (Q units)
1133static SS_RESN: i64                        // lattice nodes per axis
1134static SS_RESORG: i64                      // lattice origin, mm
1135static SS_RESSP: i64                       // lattice spacing, mm
1136static SS_RESSC: i64                       // application scale, permil
1137// load a residual file: [magic n org sp q] then n^3 i64 values; complete-or-nothing
1138func ss_res_load(path: *u8) -> i64 {
  let ln: *i64 = sys_mmap(16) as *i64
  let buf: *u8 = sys_read_file(path, ln)
  if (buf as i64) == 0 { return 0 }
  let hdr: *i64 = buf as *i64
  if hdr[0] != SS_RES_MAGIC { return 0 }
  let n: i64 = hdr[1]
  if n < 2 { return 0 }
  if n > SS_RES_MAXN { return 0 }
  if ln[0] < 40 + n*n*n*8 { return 0 }
  SS_RESN = n
  SS_RESORG = hdr[2]
  SS_RESSP = hdr[3]
  SS_RESP = (buf as i64) + 40
  return 1
1153}
1154// trilinear sample of the residual at a world-mm point, in Q units; zero outside the lattice
1155func ss_res_at(x: i64, y: i64, z: i64) -> i64 {
  if SS_RESP == 0 { return 0 }
  let R: *i64 = SS_RESP as *i64
  let n: i64 = SS_RESN
  let sp: i64 = SS_RESSP
  var ux: i64 = x - SS_RESORG
  var uy: i64 = y - SS_RESORG
  var uz: i64 = z - SS_RESORG
  if ux < 0 { return 0 }
  if uy < 0 { return 0 }
  if uz < 0 { return 0 }
  var i0: i64 = ux/sp
  var j0: i64 = uy/sp
  var k0: i64 = uz/sp
  if i0 >= n-1 { return 0 }
  if j0 >= n-1 { return 0 }
  if k0 >= n-1 { return 0 }
  let fx: i64 = (ux - i0*sp)*SS_LERP_ONE/sp
  let fy: i64 = (uy - j0*sp)*SS_LERP_ONE/sp
  let fz: i64 = (uz - k0*sp)*SS_LERP_ONE/sp
  let c000: i64 = R[(i0*n+j0)*n+k0]
  let c100: i64 = R[((i0+1)*n+j0)*n+k0]
  let c010: i64 = R[(i0*n+(j0+1))*n+k0]
  let c110: i64 = R[((i0+1)*n+(j0+1))*n+k0]
  let c001: i64 = R[(i0*n+j0)*n+(k0+1)]
  let c101: i64 = R[((i0+1)*n+j0)*n+(k0+1)]
  let c011: i64 = R[(i0*n+(j0+1))*n+(k0+1)]
  let c111: i64 = R[((i0+1)*n+(j0+1))*n+(k0+1)]
  let x00: i64 = c000 + (c100-c000)*fx/SS_LERP_ONE
  let x10: i64 = c010 + (c110-c010)*fx/SS_LERP_ONE
  let x01: i64 = c001 + (c101-c001)*fx/SS_LERP_ONE
  let x11: i64 = c011 + (c111-c011)*fx/SS_LERP_ONE
  let y0: i64 = x00 + (x10-x00)*fy/SS_LERP_ONE
  let y1: i64 = x01 + (x11-x01)*fy/SS_LERP_ONE
  return y0 + (y1-y0)*fz/SS_LERP_ONE
1190}
1191// ***THE ORACLE FIELD, DONE RIGHT (2026-08-10, operator: "nothing still looks close to a skull do
1192// you actually know what you are doing"). Honest answer: the vertex-cloud + flood + normal-heuristic
1193// construction of O was improvised through every classic mesh-to-SDF pitfall in sequence (beads,
1194// leaks, dead signs). The robust textbook construction, which the redistance pass ALREADY uses for
1195// OUR field, is: exact SURFACE CROSSINGS + PARITY. Intersect each grid column (x-ray at fixed y,z)
1196// with the actual TRIANGLES: a 2D point-in-triangle test in (y,z) selects hits, barycentric
1197// interpolation gives the exact crossing x in Q units, the two flanking nodes seed with exact
1198// distances, and CROSSING PARITY gives inside/outside -- no normals, no flood, no beads possible.
1199// The control verb `voxref` emits surface_nets(O) ALONE so the render answers one question first:
1200// is O the skull?
1201// ***NO CAPACITY CONSTANTS (rule 11, operator 2026-08-10: "stop doing these magic number bullshit
1202// things"). Bucket granularity = the GRID CELL (derived, not chosen); insertion storage is sized by a
1203// COUNTING pass, so there is no cap to hit and nothing to silently drop; the per-column hit buffer is
1204// sized by the longest measured bucket. Every bound is the population's own, computed per run.
1205// counting pass: how many (triangle, bucket-cell) insertions will the build make?
1206func ss_vb_count(rbuf: *u8, rhdr: i64, rnt: i64, cy: i64, cz: i64, ocy: i64, ocz: i64,
               org: i64, nb: i64, bcell: i64) -> i64 {
  var ins: i64 = 0
  var t: i64 = 0
  while t < rnt {
      let o: i64 = rhdr + t*84
      var loy: i64 = 0
      var hiy: i64 = 0
      var loz: i64 = 0
      var hiz: i64 = 0
      var vj: i64 = 0
      while vj < 3 {
          let wy: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+8), 1) - cz + ocy
          let wz: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,o+vj*12+4), 1) - cy) + ocz
          if vj == 0 { loy=wy; hiy=wy; loz=wz; hiz=wz } else {
              if wy<loy {loy=wy}
              if wy>hiy {hiy=wy}
              if wz<loz {loz=wz}
              if wz>hiz {hiz=wz}
          }
          vj = vj + 1
      }
      var bj: i64 = (loy - org)/bcell
      let bj1: i64 = (hiy - org)/bcell
      while bj <= bj1 {
          if bj >= 0 { if bj < nb {
              var bk: i64 = (loz - org)/bcell
              let bk1: i64 = (hiz - org)/bcell
              while bk <= bk1 {
                  if bk >= 0 { if bk < nb { ins = ins + 1 } }
                  bk = bk + 1
              }
          } }
          bj = bj + 1
      }
      t = t + 1
  }
  return ins
1244}
1245// fill pass: identical walk; capacity is the counting pass's own number, so overflow is impossible
1246func ss_vb_build(rbuf: *u8, rhdr: i64, rnt: i64, cx: i64, cy: i64, cz: i64, ocx: i64, ocy: i64, ocz: i64,
               org: i64, nb: i64, bcell: i64, head: *i64, nxt: *i64, tref: *i64, cap: i64) -> i64 {
  var c: i64 = 0
  while c < nb*nb { head[c] = 0-1; c = c + 1 }
  var ins: i64 = 0
  var t: i64 = 0
  while t < rnt {
      let o: i64 = rhdr + t*84
      var loy: i64 = 0
      var hiy: i64 = 0
      var loz: i64 = 0
      var hiz: i64 = 0
      var vj: i64 = 0
      while vj < 3 {
          let wy: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+8), 1) - cz + ocy
          let wz: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,o+vj*12+4), 1) - cy) + ocz
          if vj == 0 { loy=wy; hiy=wy; loz=wz; hiz=wz } else {
              if wy<loy {loy=wy}
              if wy>hiy {hiy=wy}
              if wz<loz {loz=wz}
              if wz>hiz {hiz=wz}
          }
          vj = vj + 1
      }
      var bj: i64 = (loy - org)/bcell
      let bj1: i64 = (hiy - org)/bcell
      while bj <= bj1 {
          if bj >= 0 { if bj < nb {
              var bk: i64 = (loz - org)/bcell
              let bk1: i64 = (hiz - org)/bcell
              while bk <= bk1 {
                  if bk >= 0 { if bk < nb {
                      if ins < cap {
                          tref[ins] = t
                          nxt[ins] = head[bj*nb+bk]
                          head[bj*nb+bk] = ins
                          ins = ins + 1
                      }
                  } }
                  bk = bk + 1
              }
          } }
          bj = bj + 1
      }
      t = t + 1
  }
  return ins
1293}
1294func ss_gap_region(x: i64, y: i64, z: i64) -> i64 {
  // coarse NAMED boxes in the field's own authored frame (mm, +y up, +z face); priority order.
  // THE NASAL BOX IS DERIVED FROM THE FIELD'S OWN NASAL PRIMITIVES (AN17, 2026-09-18): the typed box below (|x|<20,
  // z>=55, -15<=y<30) read n_ours=0 n_oracle=0 on every trend row since 2026-08-14 -- it was authored for the default
  // nasal station (y 8, z 76) while the tuned station had drifted (y -28, z 34), and the oracle's nose, read supine,
  // was never where any box expected it. The region is now wherever the nasal-bone ellipsoid and the piriform capsule
  // stand in the CURRENT vector: y from the aperture's lowest rim to the top of the nasal bones, |x| within the wider
  // of the two radii, z no further back than the nasal ellipsoid's back face. It moves with the tune.
  let P: *i64 = ss_pdef()
  var ax: i64 = x
  if ax < 0 { ax = 0-ax }
  if y < 0-55 { return 5 }
  var nxr: i64 = P[P_NRX]
  if P[P_PAW]/2 > nxr { nxr = P[P_PAW]/2 }
  // the aperture hangs from the rhinion (the plate's lower edge), so its lowest rim is one aperture height below it
  let nylo: i64 = P[P_NY] - P[P_NRY]*sp_icos(P[P_NTILT])/SS_PERMILLE - P[P_PAH]
  let nyhi: i64 = P[P_NY] + P[P_NRY]
  let nzlo: i64 = P[P_NZ] - P[P_NRY]
  if ax <= nxr { if z >= nzlo { if y >= nylo { if y <= nyhi { return 3 } } } }
  if ax < 20 { if z >= 55 { if y >= 0-15 { if y < 30 { return 3 } } } }
  if ax < 45 { if z >= 25 { if y >= 5 { if y < 40 { return 2 } } } }
  if ax >= 42 { if y >= 0-10 { if y < 25 { if z >= 0-35 { return 4 } } } }
  if y >= 0-55 { if y < 0-20 { if z >= 25 { return 6 } } }
  if y >= 18 { return 0 }
  return 1
1319}
1320func ss_gap_regname(i: i64) -> i64 {
  if i == 0 { ss_puts("vault" as *u8); return 0 }
  if i == 1 { ss_puts("base_other" as *u8); return 0 }
  if i == 2 { ss_puts("brow_orbits" as *u8); return 0 }
  if i == 3 { ss_puts("nasal" as *u8); return 0 }
  if i == 4 { ss_puts("zygomatic" as *u8); return 0 }
  if i == 5 { ss_puts("mandible" as *u8); return 0 }
  ss_puts("maxilla" as *u8)
  return 0
1329}
1330// heat ramp: green 0mm -> yellow 3mm -> red 8mm+, written as the tri's 3 colour floats
1331func ss_gap_wrcol(b: *u8, o: i64, gq: i64) -> i64 {
  let g3: i64 = 3*SS_RDQ
  let g8: i64 = 8*SS_RDQ
  var r: i64 = 0
  var g: i64 = 750
  var bl: i64 = 250
  if gq <= g3 {
      r = 200 + 650*gq/g3
  } else {
      var t: i64 = (gq-g3)*1000/(g8-g3)
      if t > 1000 { t = 1000 }
      r = 850 + 50*t/SS_PERMILLE
      g = 750 - 500*t/SS_PERMILLE
      bl = 250 - 50*t/SS_PERMILLE
  }
  ss_wr32(b, o,   ss_f32(r,1000))
  ss_wr32(b, o+4, ss_f32(g,1000))
  ss_wr32(b, o+8, ss_f32(bl,1000))
  return 0
1350}
1351// nearest grid node index for a world-mm point; wx-org >= 0 always (org is the negative bound)
1352func ss_gap_node(w: i64, org: i64, cell: i64, n1: i64) -> i64 {
  var i: i64 = (w - org + cell/2)/cell
  if i < 0 { i = 0 }
  if i > n1-1 { i = n1-1 }
  return i
1357}
1358// write one NXMSH2 header + layer table (single "bone" layer covering nt tris)
1359func ss_gap_wrhdr(buf: *u8, nt: i64) -> i64 {
  buf[0]=78 as u8; buf[1]=88 as u8; buf[2]=77 as u8; buf[3]=83 as u8
  buf[4]=72 as u8; buf[5]=50 as u8; buf[6]=0 as u8; buf[7]=0 as u8
  ss_wr32(buf, 8, 3)
  ss_wr32(buf, 12, nt)
  var L: i64 = 0
  while L < 3 {
      let lb: i64 = 16 + L*24
      var q: i64 = 0
      while q < 16 { buf[lb+q]=0 as u8; q=q+1 }
      if L==2 { buf[lb]=98 as u8; buf[lb+1]=111 as u8; buf[lb+2]=110 as u8; buf[lb+3]=101 as u8 }
      if L==2 { ss_wr32(buf, lb+16, 0); ss_wr32(buf, lb+20, nt) } else { ss_wr32(buf, lb+16, 0); ss_wr32(buf, lb+20, 0) }
      L = L + 1
  }
  return 0
1374}
1375func ss_gapmap(ref: *u8, oursout: *u8, oraout: *u8, sexf: i64, robust: i64, cum: i64, resout: *u8, latmm: i64) -> i64 {
  var G: i64 = (SS_EXT*2*1000 + cum - 1)/cum
  if G > SS_EXT*2 { G = SS_EXT*2 }
  let cell: i64 = (SS_EXT*2 + G - 1)/G
  let org: i64 = 0 - (G*cell)/2
  let n1: i64 = G+1
  if mf_admit("nx_skullsdf gapmap grids" as *u8, n1*n1*n1*8*3) == 0 { return 8 }
  let grid: *i64 = sys_mmap(n1*n1*n1*8) as *i64
  let dist: *i64 = sys_mmap(n1*n1*n1*8) as *i64
  let odist: *i64 = sys_mmap(n1*n1*n1*8) as *i64
  let hq: i64 = cell*SS_RDQ
  // ---- OUR true SDF (fill raw, seed, sweep, sign -- the proven redistance recipe) ----
  var i: i64 = 0
  while i < n1 {
      let x: i64 = org + i*cell
      var j: i64 = 0
      while j < n1 {
          let y: i64 = org + j*cell
          var k: i64 = 0
          while k < n1 {
              let z: i64 = org + k*cell
              grid[(i*n1+j)*n1+k] = ss_field_world(x, y, z, sexf, robust)
              k = k + 1
          }
          j = j + 1
      }
      i = i + 1
  }
  var si: i64 = 0
  while si < n1 {
      var sj: i64 = 0
      while sj < n1 {
          var sk: i64 = 0
          while sk < n1 {
              let ix: i64 = (si*n1+sj)*n1+sk
              let f: i64 = grid[ix]
              var best: i64 = SS_RD_INF
              if si > 0 { best = ss_rd_try(f, grid[((si-1)*n1+sj)*n1+sk], cell, best) }
              if si < n1-1 { best = ss_rd_try(f, grid[((si+1)*n1+sj)*n1+sk], cell, best) }
              if sj > 0 { best = ss_rd_try(f, grid[(si*n1+(sj-1))*n1+sk], cell, best) }
              if sj < n1-1 { best = ss_rd_try(f, grid[(si*n1+(sj+1))*n1+sk], cell, best) }
              if sk > 0 { best = ss_rd_try(f, grid[(si*n1+sj)*n1+(sk-1)], cell, best) }
              if sk < n1-1 { best = ss_rd_try(f, grid[(si*n1+sj)*n1+(sk+1)], cell, best) }
              dist[ix] = best
              sk = sk + 1
          }
          sj = sj + 1
      }
      si = si + 1
  }
  var it: i64 = 0
  while it < 2 {
      var sw: i64 = 0
      while sw < 8 {
          var di: i64 = 1
          if (sw & 1) == 1 { di = 0-1 }
          var dj: i64 = 1
          if ((sw>>1) & 1) == 1 { dj = 0-1 }
          var dk: i64 = 1
          if ((sw>>2) & 1) == 1 { dk = 0-1 }
          let sd0: i64 = ss_rd_sweep(dist, n1, hq, di, dj, dk)
          sw = sw + 1
      }
      it = it + 1
  }
  // ---- OUR surface's own bbox centre, then sign restore ----
  // The oracle below is centred by ITS bbox; without centring OURS the same way, the authored
  // frame's off-centre z biases every gap (first run: nasal read 25mm with n_oracle=0 from
  // exactly this). And surface_nets needs the SIGNED field -- polygonizing the unsigned one
  // emits a doubled shell (first run: 177,120 tris where the skull has 97,660).
  var obx0: i64 = 0; var obx1: i64 = 0
  var oby0: i64 = 0; var oby1: i64 = 0
  var obz0: i64 = 0; var obz1: i64 = 0
  var obfirst: i64 = 1
  var bi2: i64 = 0
  while bi2 < n1 {
      var bj2: i64 = 0
      while bj2 < n1 {
          var bk2: i64 = 0
          while bk2 < n1 {
              var ad2: i64 = dist[(bi2*n1+bj2)*n1+bk2]
              if ad2 < 0 { ad2 = 0-ad2 }
              if ad2 <= hq {
                  let wx2: i64 = org + bi2*cell
                  let wy2: i64 = org + bj2*cell
                  let wz2: i64 = org + bk2*cell
                  if obfirst == 1 { obx0=wx2; obx1=wx2; oby0=wy2; oby1=wy2; obz0=wz2; obz1=wz2; obfirst=0 } else {
                      if wx2<obx0 {obx0=wx2}
                      if wx2>obx1 {obx1=wx2}
                      if wy2<oby0 {oby0=wy2}
                      if wy2>oby1 {oby1=wy2}
                      if wz2<obz0 {obz0=wz2}
                      if wz2>obz1 {obz1=wz2}
                  }
              }
              bk2 = bk2 + 1
          }
          bj2 = bj2 + 1
      }
      bi2 = bi2 + 1
  }
  let ocx: i64 = (obx0+obx1)/2
  let ocy: i64 = (oby0+oby1)/2
  let ocz: i64 = (obz0+obz1)/2
  var sgi: i64 = 0
  while sgi < n1*n1*n1 {
      if grid[sgi] <= 0 { dist[sgi] = 0 - dist[sgi] }
      sgi = sgi + 1
  }
  // ---- ORACLE distance field on the SAME grid, seeded from its vertices ----
  let ln: *i64 = sys_mmap(16) as *i64
  let rbuf: *u8 = sys_read_file(ref, ln)
  if (rbuf as i64) == 0 { ss_puts("{\x22error\x22:\x22cannot read reference\x22}\n" as *u8); return 3 }
  let rlayer: i64 = ss_rd32(rbuf, 8)
  let rnt: i64 = ss_rd32(rbuf, 12)
  let rhdr: i64 = 16 + rlayer*24
  if rnt <= 0 { ss_puts("{\x22error\x22:\x22no reference triangles\x22}\n" as *u8); return 4 }
  var lox: i64 = 0; var hix: i64 = 0
  var loy: i64 = 0; var hiy: i64 = 0
  var loz: i64 = 0; var hiz: i64 = 0
  var first: i64 = 1
  var t: i64 = 0
  while t < rnt {
      let o: i64 = rhdr + t*84
      var vj: i64 = 0
      while vj < 3 {
          let x: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12),   1)
          let y: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+4), 1)
          let z: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+8), 1)
          if first == 1 { lox=x; hix=x; loy=y; hiy=y; loz=z; hiz=z; first=0 } else {
              if x<lox {lox=x}
              if x>hix {hix=x}
              if y<loy {loy=y}
              if y>hiy {hiy=y}
              if z<loz {loz=z}
              if z>hiz {hiz=z}
          }
          vj = vj + 1
      }
      t = t + 1
  }
  let cx: i64 = (lox+hix)/2
  let cy: i64 = (loy+hiy)/2
  let cz: i64 = (loz+hiz)/2
  // ***SEEDS CARRY THE SIGN OF THE ORACLE'S OWN VERTEX NORMALS. The first version signed the whole
  // near-band negative (flood barrier), and a point-cloud distance field touches zero AT every seed
  // -- surface_nets then emitted ONE BEAD PER ORACLE VERTEX (measured: ~380K verts ~ the 171K-vertex
  // benchmark, rendered as skull-shaped filigree). A node signed by the side of the surface its
  // seeding vertex's normal says it is on crosses zero ONCE.
  if mf_admit("nx_skullsdf gapmap seed-signs" as *u8, n1*n1*n1) == 0 { return 8 }
  let osgn: *u8 = sys_mmap(n1*n1*n1)
  var oi: i64 = 0
  while oi < n1*n1*n1 { odist[oi] = SS_RD_INF; osgn[oi] = 0 as u8; oi = oi + 1 }
  t = 0
  while t < rnt {
      let o: i64 = rhdr + t*84
      var vj: i64 = 0
      while vj < 3 {
          let wx: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12),   1) - cx + ocx
          let wy: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+8), 1) - cz + ocy
          let wz: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,o+vj*12+4), 1) - cy) + ocz
          let nvx: i64 = ss_r_f32(ss_rd32(rbuf,o+36+vj*12),   1000)
          let nvy: i64 = ss_r_f32(ss_rd32(rbuf,o+36+vj*12+8), 1000)
          let nvz: i64 = 0 - ss_r_f32(ss_rd32(rbuf,o+36+vj*12+4), 1000)
          let ni: i64 = ss_gap_node(wx, org, cell, n1)
          let nj: i64 = ss_gap_node(wy, org, cell, n1)
          let nk: i64 = ss_gap_node(wz, org, cell, n1)
          let dxq: i64 = (wx - (org + ni*cell))*SS_RDQ
          let dyq: i64 = (wy - (org + nj*cell))*SS_RDQ
          let dzq: i64 = (wz - (org + nk*cell))*SS_RDQ
          let dq: i64 = ss_isqrt(dxq*dxq + dyq*dyq + dzq*dzq)
          let nix: i64 = (ni*n1+nj)*n1+nk
          if dq < odist[nix] {
              odist[nix] = dq
              // node minus vertex, dotted with the outward normal: positive = outside
              let dt: i64 = (0-dxq)*nvx + (0-dyq)*nvy + (0-dzq)*nvz
              if dt >= 0 { osgn[nix] = 1 as u8 } else { osgn[nix] = 2 as u8 }
          }
          // ***SEED THE SIX NEIGHBOURS TOO, SIGNED. MEASURED MOTIVE: the benchmark's vertices are
          // INTEGER MILLIMETRES, so at cell=1 they land EXACTLY on nodes -- offset zero, sign moot
          // (+-0 is 0: the counters read 47,529 outside / 0 inside and both residuals hashed
          // identical), and an isolated zero-node beads by construction. A signed +-cell seed on
          // each neighbour puts the crossing THROUGH the vertex node instead of islanding it.
          var nbi: i64 = 0
          while nbi < 6 {
              var oi2: i64 = 0
              var oj2: i64 = 0
              var ok2: i64 = 0
              if nbi == 0 { oi2 = 1 }
              if nbi == 1 { oi2 = 0-1 }
              if nbi == 2 { oj2 = 1 }
              if nbi == 3 { oj2 = 0-1 }
              if nbi == 4 { ok2 = 1 }
              if nbi == 5 { ok2 = 0-1 }
              let mi: i64 = ni + oi2
              let mj: i64 = nj + oj2
              let mk3: i64 = nk + ok2
              var inr: i64 = 1
              if mi < 0 { inr = 0 }
              if mi > n1-1 { inr = 0 }
              if mj < 0 { inr = 0 }
              if mj > n1-1 { inr = 0 }
              if mk3 < 0 { inr = 0 }
              if mk3 > n1-1 { inr = 0 }
              if inr == 1 {
                  let vxq: i64 = (org + mi*cell - wx)*SS_RDQ
                  let vyq: i64 = (org + mj*cell - wy)*SS_RDQ
                  let vzq: i64 = (org + mk3*cell - wz)*SS_RDQ
                  let nd2: i64 = ss_isqrt(vxq*vxq + vyq*vyq + vzq*vzq)
                  let mix: i64 = (mi*n1+mj)*n1+mk3
                  if nd2 < odist[mix] {
                      odist[mix] = nd2
                      let dt2: i64 = vxq*nvx + vyq*nvy + vzq*nvz
                      if dt2 >= 0 { osgn[mix] = 1 as u8 } else { osgn[mix] = 2 as u8 }
                  }
              }
              nbi = nbi + 1
          }
          vj = vj + 1
      }
      t = t + 1
  }
  // instrument the sign seeding itself: a fix whose counters read zero is dead code, not a fix
  var sgn_out: i64 = 0
  var sgn_in: i64 = 0
  var sgi2: i64 = 0
  while sgi2 < n1*n1*n1 {
      if osgn[sgi2] == (1 as u8) { sgn_out = sgn_out + 1 }
      if osgn[sgi2] == (2 as u8) { sgn_in = sgn_in + 1 }
      sgi2 = sgi2 + 1
  }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22seedsign\x22,\x22outside\x22:" as *u8); ss_pn(sgn_out)
  ss_puts(",\x22inside\x22:" as *u8); ss_pn(sgn_in)
  ss_puts("}\n" as *u8)
  it = 0
  while it < 2 {
      var sw2: i64 = 0
      while sw2 < 8 {
          var di2: i64 = 1
          if (sw2 & 1) == 1 { di2 = 0-1 }
          var dj2: i64 = 1
          if ((sw2>>1) & 1) == 1 { dj2 = 0-1 }
          var dk2: i64 = 1
          if ((sw2>>2) & 1) == 1 { dk2 = 0-1 }
          let sd1: i64 = ss_rd_sweep(odist, n1, hq, di2, dj2, dk2)
          sw2 = sw2 + 1
      }
      it = it + 1
  }
  // ---- per-region two-sided tallies over both near-bands ----
  let fsum: *i64 = sys_mmap(SS_GR_N*8) as *i64
  let fmax: *i64 = sys_mmap(SS_GR_N*8) as *i64
  let fn:   *i64 = sys_mmap(SS_GR_N*8) as *i64
  let rsum: *i64 = sys_mmap(SS_GR_N*8) as *i64
  let rmax: *i64 = sys_mmap(SS_GR_N*8) as *i64
  let rn:   *i64 = sys_mmap(SS_GR_N*8) as *i64
  var gi: i64 = 0
  while gi < SS_GR_N { fsum[gi]=0; fmax[gi]=0; fn[gi]=0; rsum[gi]=0; rmax[gi]=0; rn[gi]=0; gi=gi+1 }
  var pi: i64 = 0
  while pi < n1 {
      let px: i64 = org + pi*cell
      var pj: i64 = 0
      while pj < n1 {
          let py: i64 = org + pj*cell
          var pk: i64 = 0
          while pk < n1 {
              let pz: i64 = org + pk*cell
              let ix2: i64 = (pi*n1+pj)*n1+pk
              var ao: i64 = dist[ix2]
              if ao < 0 { ao = 0-ao }
              let r: i64 = ss_gap_region(px, py, pz)
              if ao <= hq {
                  let gv: i64 = odist[ix2]
                  fsum[r] = fsum[r] + gv
                  if gv > fmax[r] { fmax[r] = gv }
                  fn[r] = fn[r] + 1
              }
              if odist[ix2] <= hq {
                  rsum[r] = rsum[r] + ao
                  if ao > rmax[r] { rmax[r] = ao }
                  rn[r] = rn[r] + 1
              }
              pk = pk + 1
          }
          pj = pj + 1
      }
      pi = pi + 1
  }
  // ---- heat mesh 1: OUR surface coloured by distance-to-oracle ----
  tm_reset()
  surface_nets(dist, G, G, G, org, org, org, cell, SS_ISO, SS_BONE)
  let nt1: i64 = tm_nt()
  if nt1 > 0 {
      let h1: i64 = 16 + 3*24
      let by1: i64 = h1 + nt1*84 + nt1*4
      let b1: *u8 = sys_mmap(by1 + 64)
      ss_gap_wrhdr(b1, nt1)
      var t1: i64 = 0
      while t1 < nt1 {
          let o: i64 = h1 + t1*84
          let ia: i64 = tm_ta(t1); let ib: i64 = tm_tb(t1); let ic: i64 = tm_tc(t1)
          ss_wr32(b1, o,    ss_f32(tm_vx(ia),1)); ss_wr32(b1, o+4,  ss_f32(tm_vy(ia),1)); ss_wr32(b1, o+8,  ss_f32(tm_vz(ia),1))
          ss_wr32(b1, o+12, ss_f32(tm_vx(ib),1)); ss_wr32(b1, o+16, ss_f32(tm_vy(ib),1)); ss_wr32(b1, o+20, ss_f32(tm_vz(ib),1))
          ss_wr32(b1, o+24, ss_f32(tm_vx(ic),1)); ss_wr32(b1, o+28, ss_f32(tm_vy(ic),1)); ss_wr32(b1, o+32, ss_f32(tm_vz(ic),1))
          ss_wr32(b1, o+36, ss_f32(tm_vnx(ia),SS_NQ)); ss_wr32(b1, o+40, ss_f32(tm_vny(ia),SS_NQ)); ss_wr32(b1, o+44, ss_f32(tm_vnz(ia),SS_NQ))
          ss_wr32(b1, o+48, ss_f32(tm_vnx(ib),SS_NQ)); ss_wr32(b1, o+52, ss_f32(tm_vny(ib),SS_NQ)); ss_wr32(b1, o+56, ss_f32(tm_vnz(ib),SS_NQ))
          ss_wr32(b1, o+60, ss_f32(tm_vnx(ic),SS_NQ)); ss_wr32(b1, o+64, ss_f32(tm_vny(ic),SS_NQ)); ss_wr32(b1, o+68, ss_f32(tm_vnz(ic),SS_NQ))
          let mx: i64 = (tm_vx(ia)+tm_vx(ib)+tm_vx(ic))/3
          let my: i64 = (tm_vy(ia)+tm_vy(ib)+tm_vy(ic))/3
          let mz: i64 = (tm_vz(ia)+tm_vz(ib)+tm_vz(ic))/3
          let hi1: i64 = ss_gap_node(mx, org, cell, n1)
          let hj1: i64 = ss_gap_node(my, org, cell, n1)
          let hk1: i64 = ss_gap_node(mz, org, cell, n1)
          ss_gap_wrcol(b1, o+72, odist[(hi1*n1+hj1)*n1+hk1])
          ss_wr32(b1, h1 + nt1*84 + t1*4, 2)
          t1 = t1 + 1
      }
      let fd1: i64 = sys_openat_wr(oursout, 420)
      sys_write(fd1, b1, by1)
      sys_close(fd1)
  }
  // ---- heat mesh 2: ORACLE surface coloured by distance-to-ours (what we are MISSING) ----
  let h2: i64 = 16 + 3*24
  let by2: i64 = h2 + rnt*84 + rnt*4
  let b2: *u8 = sys_mmap(by2 + 64)
  ss_gap_wrhdr(b2, rnt)
  t = 0
  while t < rnt {
      let so: i64 = rhdr + t*84
      let dsto: i64 = h2 + t*84
      var vj: i64 = 0
      var mx2: i64 = 0
      var my2: i64 = 0
      var mz2: i64 = 0
      while vj < 3 {
          let wx: i64 = ss_r_f32(ss_rd32(rbuf,so+vj*12),   1) - cx + ocx
          let wy: i64 = ss_r_f32(ss_rd32(rbuf,so+vj*12+8), 1) - cz + ocy
          let wz: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,so+vj*12+4), 1) - cy) + ocz
          ss_wr32(b2, dsto+vj*12,   ss_f32(wx,1))
          ss_wr32(b2, dsto+vj*12+4, ss_f32(wy,1))
          ss_wr32(b2, dsto+vj*12+8, ss_f32(wz,1))
          // normals: copy with z negated (decode at mille scale, re-encode negated -- no new
          // bit-twiddling operator, and a unit normal loses nothing at 1/SS_PERMILLE resolution)
          // normals through the same quarter turn: authored ny = file nz, authored nz = -file ny
          ss_wr32(b2, dsto+36+vj*12,   ss_rd32(rbuf, so+36+vj*12))
          ss_wr32(b2, dsto+36+vj*12+4, ss_rd32(rbuf, so+36+vj*12+8))
          let nzv: i64 = ss_r_f32(ss_rd32(rbuf, so+36+vj*12+4), 1000)
          ss_wr32(b2, dsto+36+vj*12+8, ss_f32(0-nzv, 1000))
          mx2 = mx2 + wx
          my2 = my2 + wy
          mz2 = mz2 + wz
          vj = vj + 1
      }
      mx2 = mx2/3
      my2 = my2/3
      mz2 = mz2/3
      let hi2: i64 = ss_gap_node(mx2, org, cell, n1)
      let hj2: i64 = ss_gap_node(my2, org, cell, n1)
      let hk2: i64 = ss_gap_node(mz2, org, cell, n1)
      var av: i64 = dist[(hi2*n1+hj2)*n1+hk2]
      if av < 0 { av = 0-av }
      ss_gap_wrcol(b2, dsto+72, av)
      ss_wr32(b2, h2 + rnt*84 + t*4, 2)
      t = t + 1
  }
  let fd2: i64 = sys_openat_wr(oraout, 420)
  sys_write(fd2, b2, by2)
  sys_close(fd2)
  // ---- optional RESIDUAL LAYER: R = O_signed - B at lattice nodes, written as DATA ----
  if (resout as i64) != 0 {
      // sign the oracle field by flooding OUTSIDE-ness from the box boundary; the near-surface
      // band is the barrier. Enclosed cavities (sinuses) stay INSIDE, which is correct anatomy.
      if mf_admit("nx_skullsdf gapmap oracle-sign" as *u8, n1*n1*n1*8) == 0 { return 8 }
      let mk: *i64 = sys_mmap(n1*n1*n1*8) as *i64
      let thr: i64 = hq
      var fi: i64 = 0
      while fi < n1 {
          var fj: i64 = 0
          while fj < n1 {
              var fk: i64 = 0
              while fk < n1 {
                  let ixf: i64 = (fi*n1+fj)*n1+fk
                  mk[ixf] = 0
                  var onb: i64 = 0
                  if fi == 0 { onb = 1 }
                  if fi == n1-1 { onb = 1 }
                  if fj == 0 { onb = 1 }
                  if fj == n1-1 { onb = 1 }
                  if fk == 0 { onb = 1 }
                  if fk == n1-1 { onb = 1 }
                  if onb == 1 { if odist[ixf] > thr { mk[ixf] = 1 } }
                  fk = fk + 1
              }
              fj = fj + 1
          }
          fi = fi + 1
      }
      var chg: i64 = 1
      var itc: i64 = 0
      while chg == 1 {
          if itc >= 8 { chg = 0 } else {
              chg = 0
              var pi2: i64 = 0
              while pi2 < n1 {
                  var pj2: i64 = 0
                  while pj2 < n1 {
                      var pk2: i64 = 0
                      while pk2 < n1 {
                          let ixp: i64 = (pi2*n1+pj2)*n1+pk2
                          if mk[ixp] == 0 { if odist[ixp] > thr {
                              var nb: i64 = 0
                              if pi2 > 0 { if mk[((pi2-1)*n1+pj2)*n1+pk2] == 1 { nb = 1 } }
                              if pi2 < n1-1 { if mk[((pi2+1)*n1+pj2)*n1+pk2] == 1 { nb = 1 } }
                              if pj2 > 0 { if mk[(pi2*n1+(pj2-1))*n1+pk2] == 1 { nb = 1 } }
                              if pj2 < n1-1 { if mk[(pi2*n1+(pj2+1))*n1+pk2] == 1 { nb = 1 } }
                              if pk2 > 0 { if mk[(pi2*n1+pj2)*n1+(pk2-1)] == 1 { nb = 1 } }
                              if pk2 < n1-1 { if mk[(pi2*n1+pj2)*n1+(pk2+1)] == 1 { nb = 1 } }
                              if nb == 1 { mk[ixp] = 1; chg = 1 }
                          } }
                          pk2 = pk2 + 1
                      }
                      pj2 = pj2 + 1
                  }
                  pi2 = pi2 + 1
              }
              var qi2: i64 = n1-1
              while qi2 >= 0 {
                  var qj2: i64 = n1-1
                  while qj2 >= 0 {
                      var qk2: i64 = n1-1
                      while qk2 >= 0 {
                          let ixq: i64 = (qi2*n1+qj2)*n1+qk2
                          if mk[ixq] == 0 { if odist[ixq] > thr {
                              var nb2: i64 = 0
                              if qi2 > 0 { if mk[((qi2-1)*n1+qj2)*n1+qk2] == 1 { nb2 = 1 } }
                              if qi2 < n1-1 { if mk[((qi2+1)*n1+qj2)*n1+qk2] == 1 { nb2 = 1 } }
                              if qj2 > 0 { if mk[(qi2*n1+(qj2-1))*n1+qk2] == 1 { nb2 = 1 } }
                              if qj2 < n1-1 { if mk[(qi2*n1+(qj2+1))*n1+qk2] == 1 { nb2 = 1 } }
                              if qk2 > 0 { if mk[(qi2*n1+qj2)*n1+(qk2-1)] == 1 { nb2 = 1 } }
                              if qk2 < n1-1 { if mk[(qi2*n1+qj2)*n1+(qk2+1)] == 1 { nb2 = 1 } }
                              if nb2 == 1 { mk[ixq] = 1; chg = 1 }
                          } }
                          qk2 = qk2 - 1
                      }
                      qj2 = qj2 - 1
                  }
                  qi2 = qi2 - 1
              }
              itc = itc + 1
          }
      }
      var lat: i64 = latmm
      if lat < cell { lat = cell }
      var n2: i64 = (G*cell)/lat + 1
      // REFUSE, never clamp: a clamped lattice would silently cover part of the box and be read
      // as the whole -- the exact silent-cap class this organ refuses one layer down (grid finer
      // than the field). Caught in my own first version, 2026-08-09.
      if n2 > SS_RES_MAXN {
          ss_puts("{\x22error\x22:\x22residual lattice finer than the bound allows\x22,\x22requested_n\x22:" as *u8); ss_pn(n2)
          ss_puts(",\x22max_n\x22:" as *u8); ss_pn(SS_RES_MAXN)
          ss_puts(",\x22finest_allowed_lat_mm\x22:" as *u8); ss_pn((G*cell)/(SS_RES_MAXN-1) + 1)
          ss_puts("}\n" as *u8)
          return 9
      }
      let rby: i64 = 40 + n2*n2*n2*8
      let rb: *u8 = sys_mmap(rby + 64)
      let rh: *i64 = rb as *i64
      rh[0] = SS_RES_MAGIC
      rh[1] = n2
      rh[2] = org
      rh[3] = lat
      rh[4] = SS_RDQ
      let rv: *i64 = ((rb as i64) + 40) as *i64
      var li: i64 = 0
      while li < n2 {
          var gi2: i64 = li*lat/cell
          if gi2 > n1-1 { gi2 = n1-1 }
          var lj: i64 = 0
          while lj < n2 {
              var gj2: i64 = lj*lat/cell
              if gj2 > n1-1 { gj2 = n1-1 }
              var lk: i64 = 0
              while lk < n2 {
                  var gk2: i64 = lk*lat/cell
                  if gk2 > n1-1 { gk2 = n1-1 }
                  let ixr: i64 = (gi2*n1+gj2)*n1+gk2
                  var osg: i64 = odist[ixr]
                  // a normal-signed seed overrides the flood; unseeded nodes fall back to it
                  var oneg: i64 = 0
                  if osgn[ixr] == (2 as u8) { oneg = 1 }
                  if osgn[ixr] == (0 as u8) { if mk[ixr] == 0 { oneg = 1 } }
                  if oneg == 1 { osg = 0 - osg }
                  rv[(li*n2+lj)*n2+lk] = osg - dist[ixr]
                  lk = lk + 1
              }
              lj = lj + 1
          }
          li = li + 1
      }
      let rfd: i64 = sys_openat_wr(resout, 420)
      sys_write(rfd, rb, rby)
      sys_close(rfd)
      ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22fitres\x22,\x22lattice_n\x22:" as *u8); ss_pn(n2)
      ss_puts(",\x22spacing_mm\x22:" as *u8); ss_pn(lat)
      ss_puts(",\x22flood_iters\x22:" as *u8); ss_pn(itc)
      ss_puts(",\x22bytes\x22:" as *u8); ss_pn(rby)
      ss_puts(",\x22derived_from\x22:\x22the ingested reference -- this file IS the match half of the circle, benchmark evidence, not a product asset\x22}\n" as *u8)
  }
  // ---- the ranked worklist ----
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22gapmap\x22,\x22cell_mm\x22:" as *u8); ss_pn(cell)
  ss_puts(",\x22ours_centre_mm\x22:[" as *u8); ss_pn(ocx)
  ss_puts("," as *u8); ss_pn(ocy)
  ss_puts("," as *u8); ss_pn(ocz)
  ss_puts("]" as *u8)
  ss_puts(",\x22bound\x22:\x22vertex-seeded oracle field OVERSTATES gaps by at most half the oracle vertex spacing (~1mm); ours near-band = our surface, oracle near-band = its surface\x22" as *u8)
  ss_puts(",\x22regions\x22:[" as *u8)
  gi = 0
  while gi < SS_GR_N {
      if gi > 0 { ss_puts("," as *u8) }
      ss_puts("{\x22region\x22:\x22" as *u8)
      ss_gap_regname(gi)
      ss_puts("\x22,\x22ours_to_oracle_mean_mm\x22:" as *u8)
      if fn[gi] > 0 { ss_pn(fsum[gi]/fn[gi]/SS_RDQ) } else { ss_pn(0-1) }
      ss_puts(",\x22ours_to_oracle_worst_mm\x22:" as *u8); ss_pn(fmax[gi]/SS_RDQ)
      ss_puts(",\x22missing_anatomy_mean_mm\x22:" as *u8)
      if rn[gi] > 0 { ss_pn(rsum[gi]/rn[gi]/SS_RDQ) } else { ss_pn(0-1) }
      ss_puts(",\x22missing_anatomy_worst_mm\x22:" as *u8); ss_pn(rmax[gi]/SS_RDQ)
      ss_puts(",\x22n_ours\x22:" as *u8); ss_pn(fn[gi])
      ss_puts(",\x22n_oracle\x22:" as *u8); ss_pn(rn[gi])
      ss_puts("}" as *u8)
      gi = gi + 1
  }
  ss_puts("]}\n" as *u8)
  return 0
1909}
1910
1911// ---- VOXREF: the oracle voxelized by exact triangle crossings + parity (derived bounds, no caps) ----
1912func ss_voxref(ref: *u8, outp: *u8, cum: i64) -> i64 {
  var G: i64 = (SS_EXT*2*1000 + cum - 1)/cum
  if G > SS_EXT*2 { G = SS_EXT*2 }
  let cell: i64 = (SS_EXT*2 + G - 1)/G
  let org: i64 = 0 - (G*cell)/2
  let n1: i64 = G+1
  let ln: *i64 = sys_mmap(16) as *i64
  let rbuf: *u8 = sys_read_file(ref, ln)
  if (rbuf as i64) == 0 { ss_puts("{\x22error\x22:\x22cannot read reference\x22}\n" as *u8); return 3 }
  let rlayer: i64 = ss_rd32(rbuf, 8)
  let rnt: i64 = ss_rd32(rbuf, 12)
  let rhdr: i64 = 16 + rlayer*24
  if rnt <= 0 { ss_puts("{\x22error\x22:\x22no reference triangles\x22}\n" as *u8); return 4 }
  // bbox centre (mm)
  var lox: i64 = 0; var hix: i64 = 0
  var loy: i64 = 0; var hiy: i64 = 0
  var loz: i64 = 0; var hiz: i64 = 0
  var first: i64 = 1
  var t: i64 = 0
  while t < rnt {
      let o: i64 = rhdr + t*84
      var vj: i64 = 0
      while vj < 3 {
          let x: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12),   1)
          let y: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+4), 1)
          let z: i64 = ss_r_f32(ss_rd32(rbuf,o+vj*12+8), 1)
          if first == 1 { lox=x; hix=x; loy=y; hiy=y; loz=z; hiz=z; first=0 } else {
              if x<lox {lox=x}
              if x>hix {hix=x}
              if y<loy {loy=y}
              if y>hiy {hiy=y}
              if z<loz {loz=z}
              if z>hiz {hiz=z}
          }
          vj = vj + 1
      }
      t = t + 1
  }
  let cx: i64 = (lox+hix)/2
  let cy: i64 = (loy+hiy)/2
  let cz: i64 = (loz+hiz)/2
  // buckets over (y,z) AT GRID RESOLUTION -- granularity derived from the grid, storage sized by
  // the counting pass: no capacity constant exists to hit
  let nb: i64 = G + 2
  let nins: i64 = ss_vb_count(rbuf, rhdr, rnt, cy, cz, 0, 0, org, nb, cell)
  if nins <= 0 { ss_puts("{\x22error\x22:\x22no triangles land in the grid\x22}\n" as *u8); return 4 }
  if mf_admit("nx_skullsdf voxref" as *u8, n1*n1*n1*8 + nb*nb*8 + nins*16) == 0 { return 8 }
  let head: *i64 = sys_mmap(nb*nb*8) as *i64
  let nxt: *i64 = sys_mmap(nins*8) as *i64
  let tref: *i64 = sys_mmap(nins*8) as *i64
  let ins: i64 = ss_vb_build(rbuf, rhdr, rnt, cx, cy, cz, 0, 0, 0, org, nb, cell, head, nxt, tref, nins)
  if ins != nins {
      // the two passes walk identical loops; disagreement is an internal defect, never tolerable
      ss_puts("{\x22error\x22:\x22count/fill mismatch\x22,\x22counted\x22:" as *u8); ss_pn(nins)
      ss_puts(",\x22filled\x22:" as *u8); ss_pn(ins)
      ss_puts("}\n" as *u8)
      return 9
  }
  // per-column hit buffer sized by the longest measured bucket (the exact upper bound)
  var maxlist: i64 = 0
  var bc2: i64 = 0
  while bc2 < nb*nb {
      var ll: i64 = 0
      var p2: i64 = head[bc2]
      while p2 >= 0 { ll = ll + 1; p2 = nxt[p2] }
      if ll > maxlist { maxlist = ll }
      bc2 = bc2 + 1
  }
  // signed field, Q units; parity signs every node, crossings seed exact flank distances
  let fld: *i64 = sys_mmap(n1*n1*n1*8) as *i64
  let hits: *i64 = sys_mmap((maxlist+1)*8) as *i64
  let orgq: i64 = org*SS_RDQ
  let cellq: i64 = cell*SS_RDQ
  var totalhits: i64 = 0
  var repaired: i64 = 0
  var maxnd: i64 = 0
  var j: i64 = 0
  while j < n1 {
      let y0q: i64 = (org + j*cell)*SS_RDQ
      var k: i64 = 0
      while k < n1 {
          let z0q: i64 = (org + k*cell)*SS_RDQ
          let bj: i64 = j
          let bk: i64 = k
          var nh: i64 = 0
          if bj >= 0 { if bj < nb { if bk >= 0 { if bk < nb {
              var p: i64 = head[bj*nb+bk]
              while p >= 0 {
                  let tt: i64 = tref[p]
                  let o: i64 = rhdr + tt*84
                  let xa: i64 = ss_r_f32(ss_rd32(rbuf,o),    SS_RDQ) - cx*SS_RDQ
                  let ya: i64 = ss_r_f32(ss_rd32(rbuf,o+8),  SS_RDQ) - cz*SS_RDQ
                  let za: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,o+4),  SS_RDQ) - cy*SS_RDQ)
                  let xb: i64 = ss_r_f32(ss_rd32(rbuf,o+12), SS_RDQ) - cx*SS_RDQ
                  let yb: i64 = ss_r_f32(ss_rd32(rbuf,o+20), SS_RDQ) - cz*SS_RDQ
                  let zb: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,o+16), SS_RDQ) - cy*SS_RDQ)
                  let xc: i64 = ss_r_f32(ss_rd32(rbuf,o+24), SS_RDQ) - cx*SS_RDQ
                  let yc: i64 = ss_r_f32(ss_rd32(rbuf,o+32), SS_RDQ) - cz*SS_RDQ
                  let zc: i64 = 0 - (ss_r_f32(ss_rd32(rbuf,o+28), SS_RDQ) - cy*SS_RDQ)
                  // edge functions in the (y,z) plane at (y0q,z0q)
                  let e0: i64 = (yb-ya)*(z0q-za) - (zb-za)*(y0q-ya)
                  let e1: i64 = (yc-yb)*(z0q-zb) - (zc-zb)*(y0q-yb)
                  let e2: i64 = (ya-yc)*(z0q-zc) - (za-zc)*(y0q-yc)
                  var hit: i64 = 0
                  if e0 >= 0 { if e1 >= 0 { if e2 >= 0 { hit = 1 } } }
                  if e0 <= 0 { if e1 <= 0 { if e2 <= 0 { hit = 1 } } }
                  if hit == 1 {
                      let den: i64 = e0 + e1 + e2
                      if den != 0 {
                          // barycentric: e1 weights vertex a, e2 weights b, e0 weights c
                          // nh is bounded by this bucket's own length <= maxlist: no cap exists
                          let xq: i64 = (e1*xa + e2*xb + e0*xc)/den
                          hits[nh] = xq
                          nh = nh + 1
                      }
                  }
                  p = nxt[p]
              }
          } } } }
          // insertion sort + dedup within 2 Q (shared-edge doubles)
          var s1: i64 = 1
          while s1 < nh {
              let v: i64 = hits[s1]
              var s2: i64 = s1
              var go2: i64 = 1
              while go2 == 1 {
                  if s2 <= 0 { go2 = 0 } else {
                      if hits[s2-1] > v { hits[s2] = hits[s2-1]; s2 = s2 - 1 } else { go2 = 0 }
                  }
              }
              hits[s2] = v
              s1 = s1 + 1
          }
          var nd: i64 = 0
          var s3: i64 = 0
          while s3 < nh {
              var keep: i64 = 1
              // shared-edge duplicates coincide geometrically; each crossing rounds ONCE in the
              // barycentric integer division, so duplicates differ by at most 2 quanta -- the
              // tolerance is the arithmetic's own bound, not a tuned number
              if nd > 0 { if hits[s3] - hits[nd-1] <= 2 { keep = 0 } }
              if keep == 1 { hits[nd] = hits[s3]; nd = nd + 1 }
              s3 = s3 + 1
          }
          // parity repair: an ODD count means this ray passed through a hole in the scan; the
          // unpaired crossing would flip parity for the rest of the column and extrude a streak to
          // the box edge (rendered, measured). Drop the deepest crossing so parity closes; COUNT it.
          if (nd & 1) == 1 { nd = nd - 1; repaired = repaired + 1 }
          if nd > maxnd { maxnd = nd }
          totalhits = totalhits + nd
          // walk nodes: parity sign everywhere, exact seeds at crossing flanks
          var hi2: i64 = 0
          var i: i64 = 0
          while i < n1 {
              let nxq: i64 = orgq + i*cellq
              var go3: i64 = 1
              while go3 == 1 {
                  if hi2 >= nd { go3 = 0 } else {
                      if hits[hi2] <= nxq { hi2 = hi2 + 1 } else { go3 = 0 }
                  }
              }
              var fv: i64 = SS_RD_INF
              // distance to nearest crossing on this column (exact near the surface)
              if hi2 > 0 {
                  let dl: i64 = nxq - hits[hi2-1]
                  if dl < fv { fv = dl }
              }
              if hi2 < nd {
                  let dr: i64 = hits[hi2] - nxq
                  if dr < fv { fv = dr }
              }
              // parity: odd crossings passed = inside
              if (hi2 & 1) == 1 { fv = 0 - fv }
              fld[(i*n1+j)*n1+k] = fv
              i = i + 1
          }
          k = k + 1
      }
      j = j + 1
  }
  // polygonize O alone -- the control: is O the skull?
  tm_reset()
  surface_nets(fld, G, G, G, org, org, org, cell, SS_ISO, SS_BONE)
  let nv: i64 = tm_nv()
  let nt2: i64 = tm_nt()
  if nt2 <= 0 { ss_puts("{\x22error\x22:\x22no surface\x22}\n" as *u8); return 5 }
  let h1: i64 = 16 + 3*24
  let by1: i64 = h1 + nt2*84 + nt2*4
  let b1: *u8 = sys_mmap(by1 + 64)
  ss_gap_wrhdr(b1, nt2)
  var t2: i64 = 0
  while t2 < nt2 {
      let o: i64 = h1 + t2*84
      let ia: i64 = tm_ta(t2); let ib: i64 = tm_tb(t2); let ic: i64 = tm_tc(t2)
      ss_wr32(b1, o,    ss_f32(tm_vx(ia),1)); ss_wr32(b1, o+4,  ss_f32(tm_vy(ia),1)); ss_wr32(b1, o+8,  ss_f32(tm_vz(ia),1))
      ss_wr32(b1, o+12, ss_f32(tm_vx(ib),1)); ss_wr32(b1, o+16, ss_f32(tm_vy(ib),1)); ss_wr32(b1, o+20, ss_f32(tm_vz(ib),1))
      ss_wr32(b1, o+24, ss_f32(tm_vx(ic),1)); ss_wr32(b1, o+28, ss_f32(tm_vy(ic),1)); ss_wr32(b1, o+32, ss_f32(tm_vz(ic),1))
      ss_wr32(b1, o+36, ss_f32(tm_vnx(ia),SS_NQ)); ss_wr32(b1, o+40, ss_f32(tm_vny(ia),SS_NQ)); ss_wr32(b1, o+44, ss_f32(tm_vnz(ia),SS_NQ))
      ss_wr32(b1, o+48, ss_f32(tm_vnx(ib),SS_NQ)); ss_wr32(b1, o+52, ss_f32(tm_vny(ib),SS_NQ)); ss_wr32(b1, o+56, ss_f32(tm_vnz(ib),SS_NQ))
      ss_wr32(b1, o+60, ss_f32(tm_vnx(ic),SS_NQ)); ss_wr32(b1, o+64, ss_f32(tm_vny(ic),SS_NQ)); ss_wr32(b1, o+68, ss_f32(tm_vnz(ic),SS_NQ))
      ss_wr32(b1, o+72, ss_f32(910,1000)); ss_wr32(b1, o+76, ss_f32(890,1000)); ss_wr32(b1, o+80, ss_f32(845,1000))
      ss_wr32(b1, h1 + nt2*84 + t2*4, 2)
      t2 = t2 + 1
  }
  let fd1: i64 = sys_openat_wr(outp, MODE_0644)
  sys_write(fd1, b1, by1)
  sys_close(fd1)
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22voxref\x22,\x22cell_mm\x22:" as *u8); ss_pn(cell)
  ss_puts(",\x22crossings\x22:" as *u8); ss_pn(totalhits)
  ss_puts(",\x22parity_repaired_columns\x22:" as *u8); ss_pn(repaired)
  ss_puts(",\x22max_crossings_per_column\x22:" as *u8); ss_pn(maxnd)
  ss_puts(",\x22max_bucket_len\x22:" as *u8); ss_pn(maxlist)
  ss_puts(",\x22bucket_insertions\x22:" as *u8); ss_pn(ins)
  ss_puts(",\x22verts\x22:" as *u8); ss_pn(nv)
  ss_puts(",\x22tris\x22:" as *u8); ss_pn(nt2)
  ss_puts(",\x22overflow\x22:" as *u8); ss_pn(tm_ovf())
  ss_puts(",\x22method\x22:\x22exact triangle crossings per grid column + parity sign -- the control: this is O alone, no base, no residual\x22}\n" as *u8)
  return 0
2130}
2131
2132// ---- REDISTANCE HELPERS (argv[8]; see the SS_RDQ const block for why this pass exists) ----
2133// one seed candidate: if f and g straddle the iso level, the crossing sits |f|/(|f|+|g|) along the
2134// edge; the distance from this node to it, in SS_RDQ units, competes for the seed minimum
2135func ss_rd_try(f: i64, g: i64, cell: i64, best: i64) -> i64 {
  var fin: i64 = 0
  if f <= 0 { fin = 1 }
  var gin: i64 = 0
  if g <= 0 { gin = 1 }
  if fin == gin { return best }
  var af: i64 = f
  if af < 0 { af = 0-af }
  var ag: i64 = g
  if ag < 0 { ag = 0-ag }
  var den: i64 = af+ag
  if den < 1 { den = 1 }
  let t: i64 = af*cell*SS_RDQ/den
  if t < best { return t }
  return best
2150}
2151// upwind eikonal update from the per-axis neighbour minima with grid spacing h (all SS_RDQ units).
2152// Sequential form: try 1D, then 2D, then 3D. Each discriminant is positive when its case is reached
2153// (reaching 2D means a+h > b so (a-b)^2 < h^2 < 2h^2), but both are still guarded.
2154func ss_rd_solve(a0: i64, b0: i64, c0: i64, h: i64) -> i64 {
  var a: i64 = a0
  var b: i64 = b0
  var c: i64 = c0
  var tv: i64 = 0
  if a > b { tv=a; a=b; b=tv }
  if b > c { tv=b; b=c; c=tv }
  if a > b { tv=a; a=b; b=tv }
  if a >= SS_RD_INF { return SS_RD_INF }
  var d: i64 = a + h
  if d <= b { return d }
  let s2: i64 = 2*h*h - (a-b)*(a-b)
  if s2 > 0 { d = (a + b + ss_isqrt(s2))/2 }
  if d <= c { return d }
  let s: i64 = a+b+c
  let s3: i64 = s*s - 3*(a*a + b*b + c*c - h*h)
  if s3 > 0 { return (s + ss_isqrt(s3))/3 }
  return d
2172}
2173// one directional sweep; returns the largest single-node improvement so convergence is a PRINTED
2174// number rather than an assumption. Min-update only ever lowers a value, and an eikonal update from
2175// true-distance neighbours can never undercut the true distance (triangle inequality), so the exact
2176// sub-cell seeds survive every sweep without a freeze flag.
2177func ss_rd_sweep(dist: *i64, n1: i64, h: i64, di: i64, dj: i64, dk: i64) -> i64 {
  var maxd: i64 = 0
  var sa: i64 = 0
  while sa < n1 {
      var i: i64 = sa
      if di < 0 { i = n1-1-sa }
      var sb: i64 = 0
      while sb < n1 {
          var j: i64 = sb
          if dj < 0 { j = n1-1-sb }
          var sc: i64 = 0
          while sc < n1 {
              var k: i64 = sc
              if dk < 0 { k = n1-1-sc }
              let ix: i64 = (i*n1+j)*n1+k
              let cur: i64 = dist[ix]
              if cur > 0 {
                  var ax: i64 = SS_RD_INF
                  if i > 0 { ax = dist[((i-1)*n1+j)*n1+k] }
                  if i < n1-1 {
                      let wx: i64 = dist[((i+1)*n1+j)*n1+k]
                      if wx < ax { ax = wx }
                  }
                  var ay: i64 = SS_RD_INF
                  if j > 0 { ay = dist[(i*n1+(j-1))*n1+k] }
                  if j < n1-1 {
                      let wy: i64 = dist[(i*n1+(j+1))*n1+k]
                      if wy < ay { ay = wy }
                  }
                  var az: i64 = SS_RD_INF
                  if k > 0 { az = dist[(i*n1+j)*n1+(k-1)] }
                  if k < n1-1 {
                      let wz: i64 = dist[(i*n1+j)*n1+(k+1)]
                      if wz < az { az = wz }
                  }
                  let nd: i64 = ss_rd_solve(ax, ay, az, h)
                  if nd < cur {
                      dist[ix] = nd
                      let dl: i64 = cur - nd
                      if dl > maxd { maxd = dl }
                  }
              }
              sc = sc + 1
          }
          sb = sb + 1
      }
      sa = sa + 1
  }
  return maxd
2226}
2227// near-band |grad| distribution of a SAMPLED field, per-mille of a true SDF -- the same statistic
2228// ss_grad reports for the analytic field, taken here on the grid so the pass can be A/B'd in-run.
2229// q = field units per millimetre (1 for the raw grid, SS_RDQ for the redistanced one).
2230func ss_rd_gradstat(F: *i64, n1: i64, cell: i64, q: i64, tag: *u8) -> i64 {
  var lo: i64 = SS_EXTENT_SENTINEL
  var hi: i64 = 0-SS_EXTENT_SENTINEL
  var sum: i64 = 0
  var n: i64 = 0
  var under: i64 = 0
  var over: i64 = 0
  let band: i64 = SS_RD_NEAR*cell*q
  var i: i64 = 1
  while i < n1-1 {
      var j: i64 = 1
      while j < n1-1 {
          var k: i64 = 1
          while k < n1-1 {
              let f: i64 = F[(i*n1+j)*n1+k]
              var af: i64 = f
              if af < 0 { af = 0-af }
              if af <= band {
                  let gx: i64 = F[((i+1)*n1+j)*n1+k] - F[((i-1)*n1+j)*n1+k]
                  let gy: i64 = F[(i*n1+(j+1))*n1+k] - F[(i*n1+(j-1))*n1+k]
                  let gz: i64 = F[(i*n1+j)*n1+(k+1)] - F[(i*n1+j)*n1+(k-1)]
                  let g: i64 = ss_isqrt(gx*gx + gy*gy + gz*gz)*1000/(2*cell*q)
                  if g < lo { lo = g }
                  if g > hi { hi = g }
                  sum = sum + g
                  n = n + 1
                  if g < SS_RD_LO { under = under + 1 }
                  if g > SS_RD_HI { over = over + 1 }
              }
              k = k + 1
          }
          j = j + 1
      }
      i = i + 1
  }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22redist_gradstat\x22,\x22field\x22:\x22" as *u8)
  ss_puts(tag)
  ss_puts("\x22,\x22true_sdf_is\x22:1000,\x22samples\x22:" as *u8); ss_pn(n)
  if n > 0 {
      ss_puts(",\x22mean\x22:" as *u8); ss_pn(sum/n)
      ss_puts(",\x22min\x22:" as *u8); ss_pn(lo)
      ss_puts(",\x22max\x22:" as *u8); ss_pn(hi)
      ss_puts(",\x22under_800_permil\x22:" as *u8); ss_pn(under*1000/n)
      ss_puts(",\x22over_1200_permil\x22:" as *u8); ss_pn(over*1000/n)
  }
  ss_puts("}\n" as *u8)
  return 0
2277}
2278
2279func main(argc: i64, argv: *i64) -> i64 {
  if argc < 2 { ss_puts("usage: nx_skullsdf <out.nxmesh> [sex] [robust] [grid] [cell_um] [tissue_permil] [k_mm] [redist] | fit <ref> | tune <ref> | gapmap <ref> <heat_ours> <heat_oracle>\n" as *u8); return 2 }
  if ss_streq(argv[1] as *u8, "tune" as *u8) == 1 {
      if argc < 3 { ss_puts("usage: nx_skullsdf tune <ref.nxmesh> [sex] [robust] [stride] [passes] [step0] [save_dat]\n" as *u8); return 2 }
      var ts: i64 = 0
      if argc > 3 { ts = ss_atoi(argv[3] as *u8) }
      var tr: i64 = 500
      if argc > 4 { tr = ss_atoi(argv[4] as *u8) }
      var tst: i64 = 31
      if argc > 5 { tst = ss_atoi(argv[5] as *u8) }
      if tst < 1 { tst = 1 }
      var tp: i64 = 5
      if argc > 6 { tp = ss_atoi(argv[6] as *u8) }
      var tstep: i64 = 8
      if argc > 7 { tstep = ss_atoi(argv[7] as *u8) }
      var tsv: i64 = 0
      if argc > 8 { tsv = argv[8] }
      return ss_tune(argv[2] as *u8, ts, tr, tst, tp, tstep, tsv as *u8)
  }
  if ss_streq(argv[1] as *u8, "grad" as *u8) == 1 {
      var gs: i64 = 0
      if argc > 2 { gs = ss_atoi(argv[2] as *u8) }
      var gr: i64 = 500
      if argc > 3 { gr = ss_atoi(argv[3] as *u8) }
      var gp: i64 = 4
      if argc > 4 { gp = ss_atoi(argv[4] as *u8) }
      if gp < 1 { gp = 1 }
      if argc > 5 {
          SS_KVAL = ss_atoi(argv[5] as *u8)
          if SS_KVAL >= 0 { SS_KSET = 1 }
      }
      return ss_grad(gs, gr, gp, 3)
  }
  if ss_streq(argv[1] as *u8, "thick" as *u8) == 1 {
      if argc < 3 { ss_puts("usage: nx_skullsdf thick <envelope.nxmesh> [sex] [robust] [stride]\n" as *u8); return 2 }
      var ks: i64 = 0
      if argc > 3 { ks = ss_atoi(argv[3] as *u8) }
      var kr: i64 = 500
      if argc > 4 { kr = ss_atoi(argv[4] as *u8) }
      var kt: i64 = 7
      if argc > 5 { kt = ss_atoi(argv[5] as *u8) }
      if kt < 1 { kt = 1 }
      return ss_thick(argv[2] as *u8, ks, kr, kt)
  }
  if ss_streq(argv[1] as *u8, "voxref" as *u8) == 1 {
      if argc < 4 { ss_puts("usage: nx_skullsdf voxref <ref.nxmesh> <out.nxmesh> [cell_um]\n" as *u8); return 2 }
      var vc: i64 = 1000
      if argc > 4 { vc = ss_atoi(argv[4] as *u8) }
      if vc < 1000 { vc = 1000 }
      return ss_voxref(argv[2] as *u8, argv[3] as *u8, vc)
  }
  if ss_streq(argv[1] as *u8, "gapmap" as *u8) == 1 {
      if argc < 5 { ss_puts("usage: nx_skullsdf gapmap <ref.nxmesh> <heat_ours.nxmesh> <heat_oracle.nxmesh> [sex] [robust] [cell_um]\n" as *u8); return 2 }
      var ms: i64 = 0
      if argc > 5 { ms = ss_atoi(argv[5] as *u8) }
      var mr: i64 = 500
      if argc > 6 { mr = ss_atoi(argv[6] as *u8) }
      var mc: i64 = SS_MARCH_CELLS_DEFAULT
      if argc > 7 { mc = ss_atoi(argv[7] as *u8) }
      if mc < SS_MARCH_CELLS_MIN { mc = SS_MARCH_CELLS_MIN }
      var mres: i64 = 0
      if argc > 8 { mres = argv[8] }
      var mlat: i64 = 4
      if argc > 9 { mlat = ss_atoi(argv[9] as *u8) }
      return ss_gapmap(argv[2] as *u8, argv[3] as *u8, argv[4] as *u8, ms, mr, mc, mres as *u8, mlat)
  }
  if ss_streq(argv[1] as *u8, "fit" as *u8) == 1 {
      if argc < 3 { ss_puts("usage: nx_skullsdf fit <ref.nxmesh> [sex] [robust] [stride]\n" as *u8); return 2 }
      var fs: i64 = 0
      if argc > 3 { fs = ss_atoi(argv[3] as *u8) }
      var fr: i64 = 500
      if argc > 4 { fr = ss_atoi(argv[4] as *u8) }
      var st: i64 = 7
      if argc > 5 { st = ss_atoi(argv[5] as *u8) }
      if st < 1 { st = 1 }
      return ss_fit(argv[2] as *u8, fs, fr, st)
  }
  var sexf: i64 = 0
  if argc > 2 { sexf = ss_atoi(argv[2] as *u8) }
  var robust: i64 = 500
  if argc > 3 { robust = ss_atoi(argv[3] as *u8) }
  var G: i64 = SS_G
  if argc > 4 { G = ss_atoi(argv[4] as *u8) }
  // ★★CELL-SIZE-DRIVEN RESOLUTION (argv[5], MICROMETRES) -- the honest dial, additive so argv[4] still
  // means G for every existing caller.
  // WHY: `cell` below is INTEGER MILLIMETRES, so with a 320mm extent EVERY G from 161 to 319 yields
  // cell=2. Raising G in that range buys ZERO resolution while GROWING THE BOX (G=224 covers 448mm at
  // the SAME 2mm) and the allocation with it -- the dial is actively harmful in its dead zone.
  // MEASURED: G=160, 176 and 192 emit BYTE-IDENTICAL geometry (cell 2, verts 48814, tris 97660).
  // Asking for a CELL SIZE makes resolution monotonic and the box exactly the extent, with no waste.
  if argc > 5 {
      let cum: i64 = ss_atoi(argv[5] as *u8)
      if cum > 0 { G = (SS_EXT*2*1000 + cum - 1)/cum }
  }
2373
  // ★TISSUE ENVELOPE (argv[6], per-mille of the published depths). ADDITIVE by construction: absent or 0
  // means SS_TISSUE stays 0 and the bone surface is emitted exactly as before.
  if argc > 6 {
      let tp: i64 = ss_atoi(argv[6] as *u8)
      if tp > 0 { SS_TISSUE = tp }
  }
  // SMOOTH-UNION RADIUS (argv[7], mm). Set 0 for a TRUE distance field with visible sutures.
  if argc > 7 {
      SS_KVAL = ss_atoi(argv[7] as *u8)
      if SS_KVAL >= 0 { SS_KSET = 1 }
  }
  // REDISTANCE (argv[8]). ★DEFAULT NOW 1 (2026-08-10, shipped winners: a plain call gets the true
  // SDF; pass 0 explicitly for the legacy raw field).
  SS_REDIST = 1
  if argc > 8 {
      SS_REDIST = ss_atoi(argv[8] as *u8)
  }
  // RESIDUAL (argv[9] scale permil, argv[10] path): the input-derived identity layer on top of the
  // parametric base. ★DEFAULT NOW 250 PERMIL (2026-08-10, operator: a bug to not ship): measured ladder
  // 0->73, 250->125 with silhouettes held and the onion-ring class dead by eye, 1000->bead cloud.
  // Explicit argv[9] overrides (0 disables). A missing/invalid residual file FAILS OPEN to the
  // parametric base on the DEFAULT path -- only an EXPLICIT request errors loudly.
  if argc <= 9 { SS_RESSC = 250 }
  if argc > 9 { SS_RESSC = ss_atoi(argv[9] as *u8) }
  if SS_RESSC > 0 {
      var rp: *u8 = "knowledge/skull_res.dat" as *u8
      if argc > 10 { rp = argv[10] as *u8 }
      if ss_res_load(rp) == 0 {
          if argc > 9 {
              ss_puts("{\x22error\x22:\x22residual requested but the file did not load (magic/size/bounds)\x22}\n" as *u8)
              return 6
          }
          SS_RESSC = 0
      }
      if SS_RESSC > 0 { if SS_REDIST != 1 {
          if argc > 9 {
              ss_puts("{\x22error\x22:\x22residual requires redist=1 (argv[8]): the base must be a true SDF\x22}\n" as *u8)
              return 6
          }
          SS_RESSC = 0
      }}
  }
2416
  // ⚠CEILING, not truncation. SS_EXT*2/G with 260mm over 144 cells gives cell=1, so the grid covered
  // only 144mm and CLIPPED the skull at its own boundary -- which is why width and depth both came
  // back as 142: they were not measurements of the skull, they were measurements of the box.
  let cell: i64 = (SS_EXT*2 + G - 1)/G
  if cell*G < SS_EXT*2 { ss_puts("{\x22error\x22:\x22grid does not cover the extent\x22}\n" as *u8); return 7 }
  // THE DIAL'S OTHER DEAD ZONE, AND IT IS THE SAME DEFECT ONE LEVEL DOWN. `cell` is an INTEGER
  // NUMBER OF MILLIMETRES, so the finest resolution this field can be sampled at is 1mm, reached at
  // G = SS_EXT*2. Every G ABOVE that still yields cell=1 -- but G*cell then EXCEEDS the extent, so the
  // grid silently covers a LARGER BOX at the SAME resolution: G=480 sampled 480mm instead of 240mm,
  // costing 8x the memory and 8x the time for ZERO extra detail, while still printing a `grid` field
  // that implies the caller got what they asked for.
  // That is exactly the failure the micrometre dial above was added to cure, so silently tolerating
  // it here would leave the cure with the disease inside it. REFUSE, and NAME the finest cell that is
  // actually expressible -- an honest non-answer beats a wrong mesh (the nx_memfloor law, applied to
  // over-resolution rather than over-allocation).
  if G > SS_EXT*2 {
      ss_puts("{\x22error\x22:\x22grid finer than the field can be sampled\x22,\x22requested_grid\x22:" as *u8); ss_pn(G)
      ss_puts(",\x22max_grid\x22:" as *u8); ss_pn(SS_EXT*2)
      ss_puts(",\x22finest_cell_mm\x22:1,\x22why\x22:\x22cell size is integer millimetres, so G above 2*SS_EXT grows the sampled box without adding resolution\x22}\n" as *u8)
      return 9
  }
  // ⚠CENTRE THE GRID ON THE SPAN IT ACTUALLY HAS. Ceiling division makes G*cell LARGER than 2*EXT, so
  // an origin of -EXT put the extra span entirely on the positive side and clipped the mandible off
  // the bottom -- which made the measured height depend on the grid rather than on the skull.
  let org: i64 = 0 - (G*cell)/2
  let n1: i64 = G+1
  // UPPER BOUND. The check above guards G being too SMALL; nothing guarded it being too LARGE, and
  // this allocation is O(G^3). On 2026-07-30 a G of ~1500 mapped 1501^3*8 = 27.06 GB on a 36 GB host
  // in under three minutes: swap hit 98%, load hit 41, and the mgmt API began refusing builds for
  // every other seat. seq1547 was the same shape in nx_ssdf and took mgmt, tools and sshd down.
  // Ask BEFORE mapping, and REFUSE rather than silently shrinking -- a caller who asked for 1500 and
  // quietly got 160 would get a wrong mesh, which is worse than an honest non-answer.
  if mf_admit("nx_skullsdf sdf grid" as *u8, n1*n1*n1*8) == 0 { return 8 }
  let grid: *i64 = sys_mmap(n1*n1*n1*8) as *i64
  var i: i64 = 0
  while i < n1 {
      let x: i64 = org + i*cell
      var j: i64 = 0
      while j < n1 {
          let y: i64 = org + j*cell
          var k: i64 = 0
          while k < n1 {
              let z: i64 = org + k*cell
              // M2: the raw grid now carries SS_FQ units/mm -- seeds, signs and surface_nets are
              // ratio/sign-based so nothing downstream changes; the tissue subtract below becomes
              // unit-CONSISTENT (ss_tdepth is already mm x SS_RDQ and SS_FQ == SS_RDQ).
              var fv: i64 = ss_field_world_q(x, y, z, sexf, robust)
              // SS_TISSUE = 0 leaves the bone iso-surface EXACTLY as it was, so every existing caller
              // is byte-identical; > 0 moves the zero level outward by the forensic depth at this height.
              // depth now varies with BOTH height and lateral distance (|x|), per the dense-map
              // literature; ss_tdepth2 at ax=0 reproduces the old midline column exactly.
              // ★★★MEASURED AND REVERTED TO THE MIDLINE COLUMN. The 2-D (height x lateral) table above is
              // RETAINED UNUSED, the GX-34 emitLimbTube precedent, because it MADE THE FACE WORSE on the
              // judge built to catch exactly this. A/B at grid 80 vs the cadaver skin oracle:
              // variance headline 128 -> 129 (noise), but PLACEMENT 135 -> 123 with recall 306 -> 278 and
              // precision unchanged at 444. For positional work this lane trusts PLACEMENT over variance.
              // ★★AND THE REASON IS STRUCTURAL, worth more than the revert: the midline depths in ss_tdepth2
              // are IDENTICAL to v1 by construction, yet the MIDLINE PROFILE STILL CHANGED. In an
              // ISO-SURFACE formulation the depths are NOT independent per-vertex knobs -- subtracting a
              // different depth off-midline moves the field in a neighbourhood, and surface_nets
              // interpolates that across cells, so a LATERAL edit propagates onto the MIDLINE SURFACE.
              // ★That is the difference from the offset formulation, where every vertex moved alone, and it
              // means the literature's DENSE TISSUE MAP cannot be transplanted cell-by-cell onto an SDF
              // envelope: the map must be fitted AGAINST THE RESULTING SURFACE, not applied to the field.
              // ***THE PER-POINT SUBTRACTION IS THE ONE THAT MEASURES BEST. I moved this into a
              // gradient-normalised second pass and it was WORSE BY A LOT -- see the disabled pass below.
              // under REDIST the offset moves to a separate pass AFTER redistancing, so the raw
              // bone field survives for seeding; without it this line is the exact old behaviour
              if SS_TISSUE > 0 { if SS_REDIST == 0 { fv = fv - ss_tdepth(y)*SS_TISSUE/SS_PERMILLE } }
              grid[(i*n1+j)*n1+k] = fv
              k = k + 1
          }
          j = j + 1
      }
      i = i + 1
  }
  // ***GRADIENT-NORMALISED TISSUE OFFSET -- THE FIX THE MEASUREMENT DEMANDED, NOT THE ONE I EXPECTED.
  // MEASURED: requesting a depth d realised 3-6mm of SCATTER around it, and that scatter barely moved
  // when the cell went 3mm -> 2mm (total spread 37 -> 34, 8 percent for 2.25x the cells). So it is NOT a
  // sampling artifact. THE CAUSE IS IN THIS FILE'S OWN HEADER: ss_smin subtracts h*h/(4k) on every
  // chained call, and 'the result is NOT a distance field -- gradient magnitude drifts far from 1'.
  // SUBTRACTING d FROM A FIELD WHOSE GRADIENT MAGNITUDE IS NOT 1 MOVES THE SURFACE BY d/|grad f|, NOT d.
  // With ~16 chained unions the drift is worst exactly where primitives overlap, which is why the nasal
  // and upper-lip bands scattered most and why everything came out systematically THIN (chained smin
  // pushes the field negative, so |grad| > 1 and d/|grad| lands short).
  // ***THE REMEDY IS THE STANDARD SDF ONE: divide the requested depth by the LOCAL GRADIENT MAGNITUDE,
  // estimated by central differences on the grid we have already built. It costs one extra pass and no
  // extra field evaluations. This also subsumes the -1mm bias, which was the same cause, not a second.
  // Gradient is carried in PER-MILLE so a true distance field reads 1000 and integer division stays exact.
  // ***DISABLED AND RETAINED (GX-34 precedent). MEASURED AT 3mm, normalised vs per-point:
  // chin req 11 realised 14 with range 8..62 (was 10, range 9..13); cheek 7 -> 8 range 4..24 (was 6,
  // range 5..8); mean_abs_err 0 -> 2. DIVIDING BY |grad f| MADE IT MUCH WORSE.
  // ***WHY, and it sharpens the diagnosis rather than excusing it: the gradient does not merely drift
  // ABOVE 1, it drifts BOTH WAYS, and where it drifts LOW the division AMPLIFIES -- a 62mm chin came
  // from a near-flat cell. My guard only caught g < 0.1; a g of 0.2 still multiplies the requested
  // depth by five. So the field is not slightly distorted, it is distorted enough that neither naive
  // remedy works: IGNORING |grad| scatters by 3-6mm, DIVIDING BY IT scatters by 50.
  // ***THE REAL FIX IS TO MAKE THE FIELD A TRUE SDF, not to compensate downstream: either drop SS_K to 0
  // (this file records having done exactly that once, for exactly this reason, before restoring 6mm for
  // appearance) and accept sutured seams, or REDISTANCE the grid -- an eikonal/fast-sweeping pass that
  // rebuilds |grad| = 1 everywhere before any offset is taken. Redistancing is the standard answer and
  // it operates on the grid we already build.
  if SS_TISSUE < 0 {
      let gm: *i64 = sys_mmap(n1*n1*n1*8) as *i64
      var gi: i64 = 0
      while gi < n1 {
          var gj: i64 = 0
          while gj < n1 {
              var gk: i64 = 0
              while gk < n1 {
                  var ai: i64 = gi+1
                  if ai > n1-1 { ai = n1-1 }
                  var bi: i64 = gi-1
                  if bi < 0 { bi = 0 }
                  var aj: i64 = gj+1
                  if aj > n1-1 { aj = n1-1 }
                  var bj: i64 = gj-1
                  if bj < 0 { bj = 0 }
                  var ak: i64 = gk+1
                  if ak > n1-1 { ak = n1-1 }
                  var bk: i64 = gk-1
                  if bk < 0 { bk = 0 }
                  let dx: i64 = (grid[(ai*n1+gj)*n1+gk] - grid[(bi*n1+gj)*n1+gk]) * 1000 / ((ai-bi)*cell*SS_FQ)
                  let dy: i64 = (grid[(gi*n1+aj)*n1+gk] - grid[(gi*n1+bj)*n1+gk]) * 1000 / ((aj-bj)*cell*SS_FQ)
                  let dz: i64 = (grid[(gi*n1+gj)*n1+ak] - grid[(gi*n1+gj)*n1+bk]) * 1000 / ((ak-bk)*cell*SS_FQ)
                  var g: i64 = ss_isqrt(dx*dx + dy*dy + dz*dz)
                  // a degenerate gradient (flat field, e.g. deep inside a solid) carries no direction
                  // information; fall back to the unnormalised depth rather than dividing by ~0
                  if g < 100 { g = 1000 }
                  gm[(gi*n1+gj)*n1+gk] = g
                  gk = gk + 1
              }
              gj = gj + 1
          }
          gi = gi + 1
      }
      gi = 0
      while gi < n1 {
          let gx: i64 = org + gi*cell
          var gj: i64 = 0
          while gj < n1 {
              let gy: i64 = org + gj*cell
              var gk: i64 = 0
              while gk < n1 {
                  let ix: i64 = (gi*n1+gj)*n1+gk
                  let d: i64 = ss_tdepth(gy)*SS_TISSUE/SS_PERMILLE
                  grid[ix] = grid[ix] - d*1000/gm[ix]
                  gk = gk + 1
              }
              gj = gj + 1
          }
          gi = gi + 1
      }
  }
  // ---- REDISTANCE PASS (argv[8]; mechanism at the SS_RDQ const block) ----
  var FLD: *i64 = grid
  if SS_REDIST == 1 {
      // its own array: grid keeps the raw signed field (for seeds, signs and the printed A/B),
      // dist receives true unsigned distances at SS_RDQ fixed point. Ask before mapping (memfloor law).
      if mf_admit("nx_skullsdf redistance grid" as *u8, n1*n1*n1*8) == 0 { return 8 }
      let dist: *i64 = sys_mmap(n1*n1*n1*8) as *i64
      let hq: i64 = cell*SS_RDQ
      var seeded: i64 = 0
      var si: i64 = 0
      while si < n1 {
          var sj: i64 = 0
          while sj < n1 {
              var sk: i64 = 0
              while sk < n1 {
                  let ix: i64 = (si*n1+sj)*n1+sk
                  let f: i64 = grid[ix]
                  var best: i64 = SS_RD_INF
                  if si > 0 { best = ss_rd_try(f, grid[((si-1)*n1+sj)*n1+sk], cell, best) }
                  if si < n1-1 { best = ss_rd_try(f, grid[((si+1)*n1+sj)*n1+sk], cell, best) }
                  if sj > 0 { best = ss_rd_try(f, grid[(si*n1+(sj-1))*n1+sk], cell, best) }
                  if sj < n1-1 { best = ss_rd_try(f, grid[(si*n1+(sj+1))*n1+sk], cell, best) }
                  if sk > 0 { best = ss_rd_try(f, grid[(si*n1+sj)*n1+(sk-1)], cell, best) }
                  if sk < n1-1 { best = ss_rd_try(f, grid[(si*n1+sj)*n1+(sk+1)], cell, best) }
                  dist[ix] = best
                  if best < SS_RD_INF { seeded = seeded + 1 }
                  sk = sk + 1
              }
              sj = sj + 1
          }
          si = si + 1
      }
      // 8 diagonal sweep orders, run twice; the second set's largest improvement is printed as the
      // convergence evidence (0 = fully converged after the first set)
      var lastd: i64 = 0
      var it: i64 = 0
      while it < 2 {
          var setd: i64 = 0
          var sw: i64 = 0
          while sw < 8 {
              var di: i64 = 1
              if (sw & 1) == 1 { di = 0-1 }
              var dj: i64 = 1
              if ((sw>>1) & 1) == 1 { dj = 0-1 }
              var dk: i64 = 1
              if ((sw>>2) & 1) == 1 { dk = 0-1 }
              let sd: i64 = ss_rd_sweep(dist, n1, hq, di, dj, dk)
              if sd > setd { setd = sd }
              sw = sw + 1
          }
          lastd = setd
          it = it + 1
      }
      // restore the sign from the raw field BEFORE measuring: on an unsigned field a central
      // difference straddling the surface reads ~0 and the statistic would convict its own input
      var pi: i64 = 0
      while pi < n1 {
          var pj: i64 = 0
          while pj < n1 {
              var pk: i64 = 0
              while pk < n1 {
                  let ix2: i64 = (pi*n1+pj)*n1+pk
                  if grid[ix2] <= 0 { dist[ix2] = 0 - dist[ix2] }
                  pk = pk + 1
              }
              pj = pj + 1
          }
          pi = pi + 1
      }
      // the in-run A/B: same statistic, same band, raw grid versus redistanced grid
      ss_rd_gradstat(grid, n1, cell, SS_FQ, "raw" as *u8)
      ss_rd_gradstat(dist, n1, cell, SS_RDQ, "redistanced" as *u8)
      ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22verb\x22:\x22redist\x22,\x22seeded\x22:" as *u8); ss_pn(seeded)
      ss_puts(",\x22last_sweep_delta_q\x22:" as *u8); ss_pn(lastd)
      ss_puts(",\x22q_per_mm\x22:" as *u8); ss_pn(SS_RDQ)
      ss_puts("}\n" as *u8)
      // the identity layer rides ON TOP of the honest base: F = B + trilinear(R) x scale
      if SS_RESSC > 0 {
          var ri: i64 = 0
          while ri < n1 {
              let rx: i64 = org + ri*cell
              var rj: i64 = 0
              while rj < n1 {
                  let ry: i64 = org + rj*cell
                  var rk: i64 = 0
                  while rk < n1 {
                      let rix: i64 = (ri*n1+rj)*n1+rk
                      dist[rix] = dist[rix] + ss_res_at(rx, ry, org + rk*cell)*SS_RESSC/SS_PERMILLE
                      rk = rk + 1
                  }
                  rj = rj + 1
              }
              ri = ri + 1
          }
      }
      // tissue offset on the TRUE distance field: an iso-offset of depth d now moves the surface by
      // exactly d, which is the whole point of the pass
      if SS_TISSUE > 0 {
          var ti: i64 = 0
          while ti < n1 {
              var tj: i64 = 0
              while tj < n1 {
                  let ty: i64 = org + tj*cell
                  // interpolated landmark depths (ss_tdepth_q) rather than the band steps: on an
                  // honest field the steps render as literal shelves; see the function's comment
                  let toff: i64 = ss_tdepth_q(ty)*SS_TISSUE/SS_PERMILLE
                  var tk: i64 = 0
                  while tk < n1 {
                      let ix3: i64 = (ti*n1+tj)*n1+tk
                      dist[ix3] = dist[ix3] - toff
                      tk = tk + 1
                  }
                  tj = tj + 1
              }
              ti = ti + 1
          }
      }
      FLD = dist
  }
  tm_reset()
  surface_nets(FLD, G, G, G, org, org, org, cell, SS_ISO, SS_BONE)
  let nv: i64 = tm_nv()
  let nt: i64 = tm_nt()
  if nt <= 0 { ss_puts("{\x22error\x22:\x22no surface -- field never crosses the iso level\x22}\n" as *u8); return 5 }
2695
  // ---- write NXMSH2, one layer ("bone"), so the existing viewer and benches read it unchanged
  let hdr: i64 = 16 + 3*24
  let bytes: i64 = hdr + nt*84 + nt*4
  let buf: *u8 = sys_mmap(bytes + 64)
  buf[0]=78 as u8; buf[1]=88 as u8; buf[2]=77 as u8; buf[3]=83 as u8
  buf[4]=72 as u8; buf[5]=50 as u8; buf[6]=0 as u8; buf[7]=0 as u8
  ss_wr32(buf, 8, 3); ss_wr32(buf, 12, nt)
  var L: i64 = 0
  while L < 3 {
      let lb: i64 = 16 + L*24
      var q: i64 = 0
      while q < 16 { buf[lb+q]=0 as u8; q=q+1 }
      if L==2 { buf[lb]=98 as u8; buf[lb+1]=111 as u8; buf[lb+2]=110 as u8; buf[lb+3]=101 as u8 }
      if L==2 { ss_wr32(buf, lb+16, 0); ss_wr32(buf, lb+20, nt) } else { ss_wr32(buf, lb+16, 0); ss_wr32(buf, lb+20, 0) }
      L = L + 1
  }
  var t: i64 = 0
  while t < nt {
      let o: i64 = hdr + t*84
      let ia: i64 = tm_ta(t); let ib: i64 = tm_tb(t); let ic: i64 = tm_tc(t)
      ss_wr32(buf, o,    ss_f32(tm_vx(ia),1)); ss_wr32(buf, o+4,  ss_f32(tm_vy(ia),1)); ss_wr32(buf, o+8,  ss_f32(tm_vz(ia),1))
      ss_wr32(buf, o+12, ss_f32(tm_vx(ib),1)); ss_wr32(buf, o+16, ss_f32(tm_vy(ib),1)); ss_wr32(buf, o+20, ss_f32(tm_vz(ib),1))
      ss_wr32(buf, o+24, ss_f32(tm_vx(ic),1)); ss_wr32(buf, o+28, ss_f32(tm_vy(ic),1)); ss_wr32(buf, o+32, ss_f32(tm_vz(ic),1))
      // ⚠surface_nets writes gradient normals at 4096 scale (nx_isosurf.nx:68), NOT the Q14 16384 the
      // mesh format uses elsewhere. Reading them as Q14 made every normal a QUARTER length, so the
      // lambert term collapsed and the whole surface rendered dark and mottled -- which I spent three
      // passes blaming on the ellipsoid, the rasterizer and the z-buffer in turn.
      ss_wr32(buf, o+36, ss_f32(tm_vnx(ia),SS_NQ)); ss_wr32(buf, o+40, ss_f32(tm_vny(ia),SS_NQ)); ss_wr32(buf, o+44, ss_f32(tm_vnz(ia),SS_NQ))
      ss_wr32(buf, o+48, ss_f32(tm_vnx(ib),SS_NQ)); ss_wr32(buf, o+52, ss_f32(tm_vny(ib),SS_NQ)); ss_wr32(buf, o+56, ss_f32(tm_vnz(ib),SS_NQ))
      ss_wr32(buf, o+60, ss_f32(tm_vnx(ic),SS_NQ)); ss_wr32(buf, o+64, ss_f32(tm_vny(ic),SS_NQ)); ss_wr32(buf, o+68, ss_f32(tm_vnz(ic),SS_NQ))
      ss_wr32(buf, o+72, ss_f32(910,1000)); ss_wr32(buf, o+76, ss_f32(890,1000)); ss_wr32(buf, o+80, ss_f32(845,1000))
      ss_wr32(buf, hdr + nt*84 + t*4, 2)
      t = t + 1
  }
  let fd: i64 = sys_openat_wr(argv[1] as *u8, 420)
  sys_write(fd, buf, bytes)
  sys_close(fd)
2733
  // ★MEASURE THE RESULT AGAINST THE REFERENCE, in the organ, every run. The scale defect shipped to a
  // published page because nothing compared the emitted skull to the 151x209x215mm it is meant to be.
  var lox: i64 = 0; var hix: i64 = 0
  var loy: i64 = 0; var hiy: i64 = 0
  var loz: i64 = 0; var hiz: i64 = 0
  var vi: i64 = 0
  while vi < nv {
      let X: i64 = tm_vx(vi); let Y: i64 = tm_vy(vi); let Z: i64 = tm_vz(vi)
      if vi == 0 { lox=X; hix=X; loy=Y; hiy=Y; loz=Z; hiz=Z } else {
          if X < lox { lox = X }
          if X > hix { hix = X }
          if Y < loy { loy = Y }
          if Y > hiy { hiy = Y }
          if Z < loz { loz = Z }
          if Z > hiz { hiz = Z }
      }
      vi = vi + 1
  }
  let bw: i64 = hix-lox
  let bh: i64 = hiy-loy
  let bd: i64 = hiz-loz
  var fit: i64 = 1
  if ss_abs(bw - SS_REF_W) > SS_TOL { fit = 0 }
  if ss_abs(bh - SS_REF_H) > SS_TOL { fit = 0 }
  if ss_abs(bd - SS_REF_D) > SS_TOL { fit = 0 }
  ss_puts("{\x22organ\x22:\x22nx_skullsdf\x22,\x22bbox_mm\x22:[" as *u8); ss_pn(bw)
  ss_puts("," as *u8); ss_pn(bh); ss_puts("," as *u8); ss_pn(bd)
  ss_puts("],\x22reference_mm\x22:[" as *u8); ss_pn(SS_REF_W)
  ss_puts("," as *u8); ss_pn(SS_REF_H); ss_puts("," as *u8); ss_pn(SS_REF_D)
  ss_puts("],\x22dims_within_tolerance\x22:" as *u8); ss_pn(fit)
  ss_puts(",\x22grid\x22:" as *u8); ss_pn(G)
  ss_puts(",\x22cell\x22:" as *u8); ss_pn(cell)
  ss_puts(",\x22verts\x22:" as *u8); ss_pn(nv)
  ss_puts(",\x22tris\x22:" as *u8); ss_pn(nt)
  ss_puts(",\x22overflow\x22:" as *u8); ss_pn(tm_ovf())
  ss_puts(",\x22fit_src\x22:" as *u8); ss_pn(SS_FITSRC)
  ss_puts(",\x22fit_n\x22:" as *u8); ss_pn(SS_FITN)
  ss_puts(",\x22res_scale\x22:" as *u8); ss_pn(SS_RESSC)
  ss_puts(",\x22method\x22:\x22smooth-union of bone primitives, orbits and nasal aperture SUBTRACTED as real openings, polygonized by surface_nets\x22}\n" as *u8)
  return 0
2774}
nx_skullsdf.nx source

`nx_skullsdf.nx` source
