nx_mmdev_lib.nx source
↩ module page · 380 lines · 16161 B
1// nx_mmdev_lib.nx -- THE MILLIMETRE SURFACE-DEVIATION RULER (PG4, procgen.plan, 2026-08-24).
2//
3// WHY: every realism claim on the human lane is graded by nx_bodybench in PERMIL of silhouette and normal
4// agreement. A cosmetic surgeon buys none of that. The clinical unit is millimetres of predicted-vs-actual
5// surface deviation, and the published soft-tissue prediction bar sits near 1-2 mm mean error. Without a
6// ruler in that unit every fidelity claim after this is aesthetic opinion wearing a number.
7//
8// WHAT IT MEASURES: for EVERY vertex of mesh A, the exact Euclidean distance to the nearest point on the
9// SURFACE of mesh B (point-to-triangle, Ericson 5.1.5, all integer), then mean, p95 and max. Full
10// population, never a sample: `verts` is printed and the gate binds every assertion to it.
11//
12// UNIT DERIVED, NOT PICKED: coordinates are held in tenths of a millimetre. The barycentric test forms
13// products of two dot products; each dot is bounded by 3 * extent^2. For an extent of 2.5 m at 0.1 mm
14// that is 3 * 25000^2 = 1.9e9 per dot and 3.5e18 per product, under the 9.2e18 i64 ceiling. A finer unit
15// overflows a whole-body subject; a coarser one is too blunt for a 1 mm bar. The extent bound is CHECKED
16// at load and the ruler REFUSES by name above it -- wrong in the direction of doing nothing.
17//
18// ACCELERATION WITHOUT SAMPLING: B's triangles go into a uniform grid whose resolution is derived from the
19// triangle count (about four triangles per cell). A query walks expanding rings of cells and stops when
20// the ring's inner face is provably farther than the best hit so far. Exact, never approximate.
21//
22// DECLARED IMPRECISIONS: (1) A's vertices are visited per triangle-vertex, so a vertex shared by k faces
23// weighs k times -- the mean is valence-weighted, which favours dense regions; printed as `tri_verts`.
24// (2) p95 is read from a histogram with 0.1 mm bins, so it is exact to the bin. (3) One-sided A->B;
25// the gate runs both directions and reports both. license_tier: ORIGINAL
26
27import "nx_syscalls.nx"
28import "nx_vecmath.nx"
29import "nx_nxmesh_lib.nx"
30
31// 2.5 m in tenths of a millimetre -- the overflow-derived extent bound above.
32const MD_EXTENT_MAX: i64 = 25000
33// Histogram: one bin per unit (0.1 mm) up to 500 mm; anything beyond lands in the last bin and is
34// counted so p95 can never be reported below a deviation that fell off the end.
35const MD_HIST_BINS: i64 = 5000
36const MD_TARGET_TRIS_PER_CELL: i64 = 4
37const MD_P95_NUM: i64 = 95
38const MD_P95_DEN: i64 = 100
39
40// Result slots.
41const MD_R_VERTS: i64 = 0
42const MD_R_TRIS_B: i64 = 1
43const MD_R_SUM: i64 = 2
44const MD_R_MAX: i64 = 3
45const MD_R_P95: i64 = 4
46const MD_R_CELLS_AXIS: i64 = 5
47const MD_R_OVERFLOW_BIN: i64 = 6
48const MD_R_REFUSED: i64 = 7 // 0 ok, 1 extent bound, 2 unreadable, 3 empty
49const MD_R_N: i64 = 8
50
51func md_res() -> *i64 {
52 let r: *i64 = sys_mmap(MD_R_N * 8) as *i64
53 var i: i64 = 0
54 while i < MD_R_N { r[i] = 0; i = i + 1 }
55 return r
56}
57
58// ---- point-to-triangle squared distance, Ericson 5.1.5, integer -------------------------------
59func md_dist2_pt_tri(px: i64, py: i64, pz: i64,
60 ax: i64, ay: i64, az: i64,
61 bx: i64, by: i64, bz: i64,
62 cx: i64, cy: i64, cz: i64) -> i64 {
63 let abx: i64 = bx - ax; let aby: i64 = by - ay; let abz: i64 = bz - az
64 let acx: i64 = cx - ax; let acy: i64 = cy - ay; let acz: i64 = cz - az
65 let apx: i64 = px - ax; let apy: i64 = py - ay; let apz: i64 = pz - az
66 let d1: i64 = abx*apx + aby*apy + abz*apz
67 let d2: i64 = acx*apx + acy*apy + acz*apz
68 var qx: i64 = ax; var qy: i64 = ay; var qz: i64 = az
69 var done: i64 = 0
70 if d1 <= 0 { if d2 <= 0 { done = 1 } }
71 if done == 0 {
72 let bpx: i64 = px - bx; let bpy: i64 = py - by; let bpz: i64 = pz - bz
73 let d3: i64 = abx*bpx + aby*bpy + abz*bpz
74 let d4: i64 = acx*bpx + acy*bpy + acz*bpz
75 if d3 >= 0 { if d4 <= d3 { qx = bx; qy = by; qz = bz; done = 1 } }
76 if done == 0 {
77 let vc: i64 = d1*d4 - d3*d2
78 if vc <= 0 { if d1 >= 0 { if d3 <= 0 {
79 let den: i64 = d1 - d3
80 if den != 0 {
81 qx = ax + abx*d1/den; qy = ay + aby*d1/den; qz = az + abz*d1/den
82 }
83 done = 1
84 } } }
85 }
86 if done == 0 {
87 let cpx: i64 = px - cx; let cpy: i64 = py - cy; let cpz: i64 = pz - cz
88 let d5: i64 = abx*cpx + aby*cpy + abz*cpz
89 let d6: i64 = acx*cpx + acy*cpy + acz*cpz
90 if d6 >= 0 { if d5 <= d6 { qx = cx; qy = cy; qz = cz; done = 1 } }
91 if done == 0 {
92 let vb: i64 = d5*d2 - d1*d6
93 if vb <= 0 { if d2 >= 0 { if d6 <= 0 {
94 let den: i64 = d2 - d6
95 if den != 0 {
96 qx = ax + acx*d2/den; qy = ay + acy*d2/den; qz = az + acz*d2/den
97 }
98 done = 1
99 } } }
100 }
101 if done == 0 {
102 let va: i64 = d3*d6 - d5*d4
103 if va <= 0 { if (d4 - d3) >= 0 { if (d5 - d6) >= 0 {
104 let den: i64 = (d4 - d3) + (d5 - d6)
105 if den != 0 {
106 let w: i64 = d4 - d3
107 qx = bx + (cx - bx)*w/den; qy = by + (cy - by)*w/den; qz = bz + (cz - bz)*w/den
108 }
109 done = 1
110 } } }
111 }
112 if done == 0 {
113 // interior: barycentric (v, w) over the face
114 let den: i64 = va + vb + vc
115 if den != 0 {
116 qx = ax + (abx*vb + acx*vc)/den
117 qy = ay + (aby*vb + acy*vc)/den
118 qz = az + (abz*vb + acz*vc)/den
119 }
120 done = 1
121 }
122 }
123 }
124 let dx: i64 = px - qx; let dy: i64 = py - qy; let dz: i64 = pz - qz
125 return dx*dx + dy*dy + dz*dz
126}
127
128// ---- the grid over B ---------------------------------------------------------------------------
129// gp layout: [0..2] min xyz, [3] cell size, [4] n per axis, [5] head base, [6] next base, [7] tri base,
130// [8] ntris, [9] entries
131const MD_G_MIN: i64 = 0
132const MD_G_CS: i64 = 3
133const MD_G_N: i64 = 4
134const MD_G_HEAD: i64 = 5
135const MD_G_NEXT: i64 = 6
136const MD_G_TRI: i64 = 7
137const MD_G_NTRI: i64 = 8
138const MD_G_ENTRIES: i64 = 9
139// DECODED COORDINATE CACHE, added 2026-08-26 for SPEED ONLY -- the numbers this ruler returns are
140// unchanged and nx_mmdev_gate's exact known answers are what prove that.
141// MEASURED: md_query called nm_coord NINE TIMES PER POINT-TRIANGLE TEST, and each nm_coord is a
142// float32 decode. On a real subject (43,386 reference triangles, 362,112 query points) that is roughly
143// 3e8 float decodes in the inner loop, and a full-population paint did not return in twenty minutes.
144// The estate already contained the answer: nx_meshdist's md_load pre-decodes its whole mesh into an
145// i64 array ONCE and completes the comparable work in about a minute. Two designs for one job, and the
146// slower one was the extracted lib everything new is supposed to compose.
147// The cache is filled by the bounds pass md_grid_build ALREADY runs, so it costs one allocation and no
148// extra reads: nt*9 words, 3.1 MB for a whole human body.
149const MD_G_COORD: i64 = 10
150const MD_G_WORDS: i64 = 11
151
152func md_cell_of(v: i64, mn: i64, cs: i64, n: i64) -> i64 {
153 var c: i64 = (v - mn) / cs
154 if c < 0 { c = 0 }
155 if c >= n { c = n - 1 }
156 return c
157}
158
159// Returns 0 ok, 1 extent refused, 3 empty. Fills gp.
160func md_grid_build(mb: *u8, gp: *i64) -> i64 {
161 let nt: i64 = nm_ntris(mb)
162 if nt <= 0 { return 3 }
163 let tb: i64 = nm_tri_base(mb)
164 gp[MD_G_TRI] = tb
165 gp[MD_G_NTRI] = nt
166 // ONE decode pass over the reference mesh fills the bounds AND the coordinate cache together, so
167 // the cache costs an allocation and not a single extra float32 decode.
168 let co: *i64 = sys_mmap(nt * 9 * 8 + 64) as *i64
169 gp[MD_G_COORD] = co as i64
170 var mnx: i64 = nm_coord(mb, tb, 0, 0, 0); var mxx: i64 = mnx
171 var mny: i64 = nm_coord(mb, tb, 0, 0, 1); var mxy: i64 = mny
172 var mnz: i64 = nm_coord(mb, tb, 0, 0, 2); var mxz: i64 = mnz
173 var t: i64 = 0
174 while t < nt {
175 var v: i64 = 0
176 while v < 3 {
177 let x: i64 = nm_coord(mb, tb, t, v, 0)
178 let y: i64 = nm_coord(mb, tb, t, v, 1)
179 let z: i64 = nm_coord(mb, tb, t, v, 2)
180 co[t*9 + v*3] = x
181 co[t*9 + v*3 + 1] = y
182 co[t*9 + v*3 + 2] = z
183 if x < mnx { mnx = x } if x > mxx { mxx = x }
184 if y < mny { mny = y } if y > mxy { mxy = y }
185 if z < mnz { mnz = z } if z > mxz { mxz = z }
186 v = v + 1
187 }
188 t = t + 1
189 }
190 var ext: i64 = mxx - mnx
191 if mxy - mny > ext { ext = mxy - mny }
192 if mxz - mnz > ext { ext = mxz - mnz }
193 if ext > MD_EXTENT_MAX { return 1 }
194 if ext < 1 { ext = 1 }
195 // n per axis DERIVED: smallest n with n^3 * target >= ntris
196 var n: i64 = 1
197 while n * n * n * MD_TARGET_TRIS_PER_CELL < nt { n = n + 1 }
198 var cs: i64 = ext / n
199 if cs < 1 { cs = 1 }
200 // cells must cover the extent after integer division
201 while cs * n < ext + 1 { cs = cs + 1 }
202 gp[MD_G_MIN] = mnx; gp[MD_G_MIN+1] = mny; gp[MD_G_MIN+2] = mnz
203 gp[MD_G_CS] = cs
204 gp[MD_G_N] = n
205 let ncell: i64 = n * n * n
206 let head: *i64 = sys_mmap(ncell * 8) as *i64
207 var i: i64 = 0
208 while i < ncell { head[i] = 0 - 1; i = i + 1 }
209 // pass 1: count entries (a triangle enters every cell its AABB touches)
210 var entries: i64 = 0
211 t = 0
212 while t < nt {
213 var lo0: i64 = n; var hi0: i64 = 0 - 1; var lo1: i64 = n; var hi1: i64 = 0 - 1; var lo2: i64 = n; var hi2: i64 = 0 - 1
214 var v: i64 = 0
215 while v < 3 {
216 let c0: i64 = md_cell_of(co[t*9 + v*3], mnx, cs, n)
217 let c1: i64 = md_cell_of(co[t*9 + v*3 + 1], mny, cs, n)
218 let c2: i64 = md_cell_of(co[t*9 + v*3 + 2], mnz, cs, n)
219 if c0 < lo0 { lo0 = c0 } if c0 > hi0 { hi0 = c0 }
220 if c1 < lo1 { lo1 = c1 } if c1 > hi1 { hi1 = c1 }
221 if c2 < lo2 { lo2 = c2 } if c2 > hi2 { hi2 = c2 }
222 v = v + 1
223 }
224 entries = entries + (hi0 - lo0 + 1) * (hi1 - lo1 + 1) * (hi2 - lo2 + 1)
225 t = t + 1
226 }
227 // entry records: tri index and next pointer, two words each
228 let nxt: *i64 = sys_mmap(entries * 16 + 16) as *i64
229 var e: i64 = 0
230 t = 0
231 while t < nt {
232 var lo0: i64 = n; var hi0: i64 = 0 - 1; var lo1: i64 = n; var hi1: i64 = 0 - 1; var lo2: i64 = n; var hi2: i64 = 0 - 1
233 var v: i64 = 0
234 while v < 3 {
235 let c0: i64 = md_cell_of(co[t*9 + v*3], mnx, cs, n)
236 let c1: i64 = md_cell_of(co[t*9 + v*3 + 1], mny, cs, n)
237 let c2: i64 = md_cell_of(co[t*9 + v*3 + 2], mnz, cs, n)
238 if c0 < lo0 { lo0 = c0 } if c0 > hi0 { hi0 = c0 }
239 if c1 < lo1 { lo1 = c1 } if c1 > hi1 { hi1 = c1 }
240 if c2 < lo2 { lo2 = c2 } if c2 > hi2 { hi2 = c2 }
241 v = v + 1
242 }
243 var a: i64 = lo0
244 while a <= hi0 {
245 var b2: i64 = lo1
246 while b2 <= hi1 {
247 var c: i64 = lo2
248 while c <= hi2 {
249 let cell: i64 = (a * n + b2) * n + c
250 nxt[e*2] = t
251 nxt[e*2+1] = head[cell]
252 head[cell] = e
253 e = e + 1
254 c = c + 1
255 }
256 b2 = b2 + 1
257 }
258 a = a + 1
259 }
260 t = t + 1
261 }
262 gp[MD_G_HEAD] = head as i64
263 gp[MD_G_NEXT] = nxt as i64
264 gp[MD_G_ENTRIES] = e
265 return 0
266}
267
268// Exact nearest squared distance from point p to B's surface via the grid.
269func md_query(mb: *u8, gp: *i64, px: i64, py: i64, pz: i64) -> i64 {
270 let n: i64 = gp[MD_G_N]
271 let cs: i64 = gp[MD_G_CS]
272 let head: *i64 = gp[MD_G_HEAD] as *i64
273 let nxt: *i64 = gp[MD_G_NEXT] as *i64
274 // the decoded reference coordinates, not the raw buffer: see MD_G_COORD for why
275 let co: *i64 = gp[MD_G_COORD] as *i64
276 let c0: i64 = md_cell_of(px, gp[MD_G_MIN], cs, n)
277 let c1: i64 = md_cell_of(py, gp[MD_G_MIN+1], cs, n)
278 let c2: i64 = md_cell_of(pz, gp[MD_G_MIN+2], cs, n)
279 var best: i64 = 0 - 1
280 var r: i64 = 0
281 var go: i64 = 1
282 while go == 1 {
283 // provable lower bound for ring r: (r-1) whole cells of clearance
284 if r > 0 { if best >= 0 {
285 let lb: i64 = (r - 1) * cs
286 if lb * lb > best { go = 0 }
287 } }
288 if r > n { go = 0 }
289 if go == 1 {
290 var a: i64 = c0 - r
291 while a <= c0 + r {
292 var b2: i64 = c1 - r
293 while b2 <= c1 + r {
294 var c: i64 = c2 - r
295 while c <= c2 + r {
296 var shell: i64 = 0
297 if a == c0 - r { shell = 1 } if a == c0 + r { shell = 1 }
298 if b2 == c1 - r { shell = 1 } if b2 == c1 + r { shell = 1 }
299 if c == c2 - r { shell = 1 } if c == c2 + r { shell = 1 }
300 if r == 0 { shell = 1 }
301 if shell == 1 { if a >= 0 { if a < n { if b2 >= 0 { if b2 < n { if c >= 0 { if c < n {
302 var e: i64 = head[(a * n + b2) * n + c]
303 while e >= 0 {
304 let t: i64 = nxt[e*2]
305 let b9: i64 = t*9
306 let d2: i64 = md_dist2_pt_tri(px, py, pz,
307 co[b9], co[b9+1], co[b9+2],
308 co[b9+3], co[b9+4], co[b9+5],
309 co[b9+6], co[b9+7], co[b9+8])
310 if best < 0 { best = d2 }
311 if d2 < best { best = d2 }
312 e = nxt[e*2+1]
313 }
314 } } } } } } }
315 c = c + 1
316 }
317 b2 = b2 + 1
318 }
319 a = a + 1
320 }
321 r = r + 1
322 }
323 }
324 if best < 0 { return 0 }
325 return best
326}
327
328// A -> B deviation over EVERY triangle-vertex of A. Fills res. Returns the refusal code (0 ok).
329func md_measure(ma: *u8, mb: *u8, res: *i64) -> i64 {
330 let gp: *i64 = sys_mmap(MD_G_WORDS * 8) as *i64
331 let grc: i64 = md_grid_build(mb, gp)
332 if grc != 0 { res[MD_R_REFUSED] = grc; return grc }
333 res[MD_R_TRIS_B] = gp[MD_G_NTRI]
334 res[MD_R_CELLS_AXIS] = gp[MD_G_N]
335 let nta: i64 = nm_ntris(ma)
336 if nta <= 0 { res[MD_R_REFUSED] = 3; return 3 }
337 let tba: i64 = nm_tri_base(ma)
338 let hist: *i64 = sys_mmap(MD_HIST_BINS * 8) as *i64
339 var i: i64 = 0
340 while i < MD_HIST_BINS { hist[i] = 0; i = i + 1 }
341 var sum: i64 = 0
342 var mx: i64 = 0
343 var cnt: i64 = 0
344 var t: i64 = 0
345 while t < nta {
346 var v: i64 = 0
347 while v < 3 {
348 let d2: i64 = md_query(mb, gp, nm_coord(ma, tba, t, v, 0), nm_coord(ma, tba, t, v, 1), nm_coord(ma, tba, t, v, 2))
349 let d: i64 = vm_isqrt(d2)
350 sum = sum + d
351 if d > mx { mx = d }
352 var bin: i64 = d
353 if bin >= MD_HIST_BINS { bin = MD_HIST_BINS - 1; res[MD_R_OVERFLOW_BIN] = res[MD_R_OVERFLOW_BIN] + 1 }
354 hist[bin] = hist[bin] + 1
355 cnt = cnt + 1
356 v = v + 1
357 }
358 t = t + 1
359 }
360 res[MD_R_VERTS] = cnt
361 res[MD_R_SUM] = sum
362 res[MD_R_MAX] = mx
363 // p95 from the histogram: the first bin at which the cumulative count reaches 95 percent
364 let target: i64 = cnt * MD_P95_NUM / MD_P95_DEN
365 var acc: i64 = 0
366 var p95: i64 = MD_HIST_BINS - 1
367 i = 0
368 var found: i64 = 0
369 while i < MD_HIST_BINS {
370 if found == 0 {
371 acc = acc + hist[i]
372 if acc >= target { p95 = i; found = 1 }
373 }
374 i = i + 1
375 }
376 res[MD_R_P95] = p95
377 return 0
378}
379
380func md_mean(res: *i64) -> i64 { if res[MD_R_VERTS] == 0 { return 0 } return res[MD_R_SUM] / res[MD_R_VERTS] }