nx_microscan_lib.nx source
↩ module page · 366 lines · 16341 B
1// nx_microscan_lib.nx -- INVERSE MICRO-SURFACE ANALYSIS: infer a surface's micro-STRUCTURE, and its
2// DEFECTS, from the field itself. The estate already GENERATES micro-relief (nx_relief_lib owns the
3// forward law: three bands and a jittered-grid pore layer). Nothing ever read a field back. This is
4// that half, and it is deliberately the INVERSE of a generator rather than a second generator.
5//
6// WHY THIS IS A LEAF THAT KNOWS NOTHING ABOUT SKIN. Every function below takes a plain field of i64
7// samples, a width, a height and the physical size of one sample. It has no pore constant, no skin
8// constant and no relief import -- so the same analyser reads a skin height map, a machined surface,
9// or a field of WALL THICKNESSES off a mesh, and the teapot-leak question is the same call as the
10// pore question with a different quantity in the array. A ruler that imported the thing it measures
11// could not be pointed anywhere else, and its round-trip test would be circular.
12//
13// THE FOUR QUESTIONS IT ANSWERS, all as NUMBERS and never as a diagnosis:
14// 1. WHAT IS THE PERIOD of the discrete features carried here? ms_period_texels / ms_period_um
15// 2. IS THIS A LATTICE AT ALL, or a random scatter? ms_regularity_permil
16// 3. ARE TWO LAYERS BOUND OR INDEPENDENT? ms_hit_permil / ms_reg_verdict
17// 4. WHERE DOES THIS FIELD DEPART FROM ITSELF? ms_robust / ms_anom_count / ms_min_loc
18//
19// DECLARED SCOPE, because this is the half that is usually overclaimed: question 4 emits CANDIDATE
20// REGIONS WITH SCORES. It is an outlier finder over a field. It is NOT a diagnosis of anything and no
21// clinical reading may be built on it -- the resolution bar is stated in the gate's verdict note.
22//
23// license_tier: ORIGINAL No hw writes (Rule 26).
24import "nx_syscalls.nx"
25import "nx_vecmath.nx"
26
27// larger than any squared pixel distance the fields below can produce, so it is a safe minimum seed.
28const MS_D2_INF: i64 = 1000000000
29const MS_PERMIL: i64 = 1000
30const MS_MILLI: i64 = 1000
31// distance is accumulated at MILLI precision so a mean of integer distances does not quantise away
32// the very difference between a lattice and a scatter that ms_regularity_permil exists to see.
33const MS_MILLI2: i64 = 1000000
34
35// registration verdicts. UNMEASURABLE is a first-class answer: a field with no features, or a point
36// set with no points, must not be scored -- an analyser that cannot see must abstain, never acquit.
37const MS_REG_DECOUPLED: i64 = 0
38const MS_REG_PARTIAL: i64 = 1
39const MS_REG_BOUND: i64 = 2
40const MS_REG_UNMEASURABLE: i64 = 3
41// THE TWO BARS ARE DECLARED AND PRINTED, never hidden. BOUND requires the observed hit rate to beat
42// the field's OWN chance rate by this multiple; DECOUPLED is a rate within this multiple of chance.
43// Between them the honest answer is PARTIAL rather than a forced call.
44const MS_REG_K_BOUND: i64 = 4
45const MS_REG_K_CHANCE: i64 = 2
46// THE MINIMUM DENOMINATOR, AND IT IS ASYMMETRIC ON PURPOSE. Measured 2026-09-01: with 256 points
47// against a 12-permil chance rate the expected number of chance hits is THREE, and at three expected
48// events an observation of eight is ordinary noise -- so a DECOUPLED verdict there would be a claim of
49// independence made on evidence that could not have detected binding had it existed.
50// BOUND is a PRESENCE claim and one strong witness carries it. DECOUPLED is an ABSENCE claim and
51// absence needs coverage. So the guard binds the DECOUPLED path only: below this many EXPECTED chance
52// hits the honest answer is UNMEASURABLE, never DECOUPLED and never PARTIAL.
53const MS_REG_MIN_EXPECT: i64 = 10
54
55// the two ANALYTIC references regularity is read against. For a Poisson scatter the mean nearest
56// neighbour distance is half the density-derived cell edge; for a perfect lattice it is the whole
57// edge. Both are properties of the geometry, not tuned numbers, which is why this measure needs no
58// calibration constant of its own.
59const MS_REG_POISSON_PERMIL: i64 = 500
60const MS_REG_LATTICE_PERMIL: i64 = 1000
61
62func ms_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
63
64// ---------------------------------------------------------------------------------------------
65// FEATURES. A sample is a feature peak when it clears the threshold and is the maximum of its own
66// 3x3 neighbourhood under a TOTAL ORDER: strictly greater than every neighbour preceding it in scan
67// order, greater-or-equal to every one following. On a plateau -- which an integer-sampled pit floor
68// routinely is -- that counts the region EXACTLY ONCE. A naive strict-maximum test counts a plateau
69// zero times and a naive non-strict test counts it many, and both errors land directly in the
70// density that the period is derived from.
71func ms_is_peak(f: *i64, w: i64, h: i64, x: i64, y: i64, thresh: i64) -> i64 {
72 let v: i64 = f[y * w + x]
73 if v <= thresh { return 0 }
74 var ok: i64 = 1
75 var dy: i64 = 0 - 1
76 while dy <= 1 {
77 var dx: i64 = 0 - 1
78 while dx <= 1 {
79 var skip: i64 = 0
80 if dx == 0 { if dy == 0 { skip = 1 } }
81 let ax: i64 = x + dx
82 let ay: i64 = y + dy
83 if ax < 0 { skip = 1 }
84 if ax >= w { skip = 1 }
85 if ay < 0 { skip = 1 }
86 if ay >= h { skip = 1 }
87 if skip == 0 {
88 let nv: i64 = f[ay * w + ax]
89 var earlier: i64 = 0
90 if dy < 0 { earlier = 1 }
91 if dy == 0 { if dx < 0 { earlier = 1 } }
92 if earlier == 1 { if nv >= v { ok = 0 } }
93 if earlier == 0 { if nv > v { ok = 0 } }
94 }
95 dx = dx + 1
96 }
97 dy = dy + 1
98 }
99 return ok
100}
101
102// enumerate every feature peak into xs/ys. Returns the count, or -1 IF THE BUFFER FILLED -- a silent
103// cap here would turn a partial census into a smaller density and therefore a LARGER period, i.e. a
104// confident wrong answer that reads exactly like a measurement.
105func ms_features(f: *i64, w: i64, h: i64, thresh: i64, xs: *i64, ys: *i64, cap: i64) -> i64 {
106 var n: i64 = 0
107 var over: i64 = 0
108 var y: i64 = 0
109 while y < h {
110 var x: i64 = 0
111 while x < w {
112 if ms_is_peak(f, w, h, x, y, thresh) == 1 {
113 if n < cap { xs[n] = x; ys[n] = y; n = n + 1 } else { over = 1 }
114 }
115 x = x + 1
116 }
117 y = y + 1
118 }
119 if over == 1 { return 0 - 1 }
120 return n
121}
122
123// ---------------------------------------------------------------------------------------------
124// THE PERIOD, FROM DENSITY. For a field carrying one feature per lattice cell -- which is exactly
125// what a jittered-grid cellular layer is BY CONSTRUCTION -- the cell edge is sqrt(area / count), and
126// that estimate is INVARIANT TO THE JITTER. That invariance is why density is the primary estimator
127// here and a nearest-neighbour mean is not: jitter moves every neighbour distance and moves the
128// count not at all.
129func ms_period_texels(n_feat: i64, w: i64, h: i64) -> i64 {
130 if n_feat <= 0 { return 0 - 1 }
131 return vm_isqrt(w * h / n_feat)
132}
133// the same period in physical units. Computed in square microns rather than by scaling the texel
134// answer, so the integer square root is taken at the finer scale and the round trip stays exact.
135func ms_period_um(n_feat: i64, w: i64, h: i64, texel_um: i64) -> i64 {
136 if n_feat <= 0 { return 0 - 1 }
137 if texel_um <= 0 { return 0 - 1 }
138 return vm_isqrt((w * texel_um) * (h * texel_um) / n_feat)
139}
140
141// mean nearest-neighbour distance in texels at MILLI precision, or -1 when there is no pair to
142// measure. O(n^2) and deliberately so: the populations this reads are hundreds of features, and an
143// approximate neighbour structure would put an unstated error inside the one number that decides
144// whether a field is a lattice at all.
145func ms_nn_mean_milli(xs: *i64, ys: *i64, n: i64) -> i64 {
146 if n < 2 { return 0 - 1 }
147 var sum: i64 = 0
148 var i: i64 = 0
149 while i < n {
150 var best: i64 = MS_D2_INF
151 var j: i64 = 0
152 while j < n {
153 if j != i {
154 let dx: i64 = xs[j] - xs[i]
155 let dy: i64 = ys[j] - ys[i]
156 let d2: i64 = dx * dx + dy * dy
157 if d2 < best { best = d2 }
158 }
159 j = j + 1
160 }
161 sum = sum + vm_isqrt(best * MS_MILLI2)
162 i = i + 1
163 }
164 return sum / n
165}
166
167// IS THIS A LATTICE? The nearest-neighbour mean as a fraction of the density-derived cell edge.
168// Near MS_REG_POISSON_PERMIL the features are a random scatter; near MS_REG_LATTICE_PERMIL they sit
169// on a grid; a jittered grid lands between. Both ends are analytic, so this discriminates WITHOUT a
170// calibrated constant -- and it is the check that stops a period being read off a field that has no
171// period to give.
172func ms_regularity_permil(nn_mean_milli: i64, period_texels: i64) -> i64 {
173 if nn_mean_milli < 0 { return 0 - 1 }
174 if period_texels <= 0 { return 0 - 1 }
175 return nn_mean_milli / period_texels
176}
177
178// ---------------------------------------------------------------------------------------------
179// COVERAGE: the fraction of the surface the features occupy. This doubles as the CHANCE RATE for
180// registration below, which is the whole point -- the null hypothesis is measured off the same field
181// rather than assumed.
182func ms_coverage_permil(f: *i64, w: i64, h: i64, thresh: i64) -> i64 {
183 let n: i64 = w * h
184 if n <= 0 { return 0 - 1 }
185 var c: i64 = 0
186 var i: i64 = 0
187 while i < n {
188 if f[i] > thresh { c = c + 1 }
189 i = i + 1
190 }
191 return c * MS_PERMIL / n
192}
193
194// REGISTRATION: what fraction of a set of PLACED POINTS lands on a feature. Points outside the field
195// are counted in the denominator and never as hits -- dropping them would let a point set score by
196// leaving the surface.
197func ms_hit_permil(f: *i64, w: i64, h: i64, thresh: i64, px: *i64, py: *i64, n: i64) -> i64 {
198 if n <= 0 { return 0 - 1 }
199 var hit: i64 = 0
200 var i: i64 = 0
201 while i < n {
202 let x: i64 = px[i]
203 let y: i64 = py[i]
204 var inside: i64 = 1
205 if x < 0 { inside = 0 }
206 if x >= w { inside = 0 }
207 if y < 0 { inside = 0 }
208 if y >= h { inside = 0 }
209 if inside == 1 { if f[y * w + x] > thresh { hit = hit + 1 } }
210 i = i + 1
211 }
212 return hit * MS_PERMIL / n
213}
214
215// BOUND or INDEPENDENT, against the field's own chance rate. This is the shape of the question "does
216// the hair actually come out of the pore": a groom that area-samples a surface without ever
217// consulting the pore layer scores AT CHANCE, and no amount of follicle count changes that.
218//
219// n_points is REQUIRED and is not decoration -- it is the denominator the DECOUPLED claim rests on.
220// A rate comparison with three expected events cannot separate "at chance" from "twice chance", and a
221// verdict function that took only two rates would have no way to know that. This one abstains.
222func ms_reg_verdict(hit_permil: i64, chance_permil: i64, n_points: i64) -> i64 {
223 if hit_permil < 0 { return MS_REG_UNMEASURABLE }
224 if chance_permil < 0 { return MS_REG_UNMEASURABLE }
225 if n_points <= 0 { return MS_REG_UNMEASURABLE }
226 // PRESENCE FIRST, and it needs no denominator: an excess this large over the field's own chance
227 // rate is not something a small sample manufactures.
228 if hit_permil >= chance_permil * MS_REG_K_BOUND {
229 if hit_permil > chance_permil { return MS_REG_BOUND }
230 }
231 // ABSENCE NEEDS COVERAGE. Below the minimum expected count the instrument could not have SEEN
232 // binding, so it may not report the lack of it.
233 if n_points * chance_permil / MS_PERMIL < MS_REG_MIN_EXPECT { return MS_REG_UNMEASURABLE }
234 if hit_permil <= chance_permil * MS_REG_K_CHANCE { return MS_REG_DECOUPLED }
235 return MS_REG_PARTIAL
236}
237
238// ---------------------------------------------------------------------------------------------
239// DEPARTURE FROM SELF. Block means read against the MEDIAN of block means, scaled by the median
240// absolute deviation -- robust, so a handful of genuine anomalies cannot inflate the very yardstick
241// that is supposed to find them. A mean-and-standard-deviation version of this hides exactly the
242// case it is built for.
243func ms_block_means(f: *i64, w: i64, h: i64, bs: i64, out: *i64, cap: i64) -> i64 {
244 if bs <= 0 { return 0 - 1 }
245 let bw: i64 = w / bs
246 let bh: i64 = h / bs
247 if bw <= 0 { return 0 - 1 }
248 if bh <= 0 { return 0 - 1 }
249 if bw * bh > cap { return 0 - 1 }
250 var by: i64 = 0
251 while by < bh {
252 var bx: i64 = 0
253 while bx < bw {
254 var sum: i64 = 0
255 var y: i64 = 0
256 while y < bs {
257 var x: i64 = 0
258 while x < bs {
259 sum = sum + f[(by * bs + y) * w + (bx * bs + x)]
260 x = x + 1
261 }
262 y = y + 1
263 }
264 out[by * bw + bx] = sum / (bs * bs)
265 bx = bx + 1
266 }
267 by = by + 1
268 }
269 return bw * bh
270}
271
272// median and median-absolute-deviation into out2, leaving the caller's array UNTOUCHED. The sort
273// happens in a scratch buffer the caller owns: an in-place ruler that reorders its own input makes
274// every index reported afterwards refer to a different sample than the caller believes.
275func ms_robust(v: *i64, n: i64, scratch: *i64, out2: *i64) -> i64 {
276 out2[0] = 0 - 1
277 out2[1] = 0 - 1
278 if n <= 0 { return 0 - 1 }
279 var i: i64 = 0
280 while i < n { scratch[i] = v[i]; i = i + 1 }
281 var a: i64 = 1
282 while a < n {
283 let key: i64 = scratch[a]
284 var b: i64 = a - 1
285 var go: i64 = 1
286 while go == 1 {
287 if b >= 0 {
288 if scratch[b] > key { scratch[b + 1] = scratch[b]; b = b - 1 } else { go = 0 }
289 } else { go = 0 }
290 }
291 scratch[b + 1] = key
292 a = a + 1
293 }
294 let med: i64 = scratch[n / 2]
295 out2[0] = med
296 i = 0
297 while i < n { scratch[i] = ms_abs(v[i] - med); i = i + 1 }
298 a = 1
299 while a < n {
300 let key2: i64 = scratch[a]
301 var b2: i64 = a - 1
302 var go2: i64 = 1
303 while go2 == 1 {
304 if b2 >= 0 {
305 if scratch[b2] > key2 { scratch[b2 + 1] = scratch[b2]; b2 = b2 - 1 } else { go2 = 0 }
306 } else { go2 = 0 }
307 }
308 scratch[b2 + 1] = key2
309 a = a + 1
310 }
311 out2[1] = scratch[n / 2]
312 return 0
313}
314
315// candidate count at a DECLARED sensitivity k, in MILLI, so k = 3000 means three deviations.
316// DECLARED IMPRECISION, and it is a real one: when the deviation is zero the field has a single
317// dominant level and there is no scale to divide by, so any departure at all is a candidate. That is
318// stated here rather than silently returning zero, because a zero from a degenerate scale and a zero
319// from a clean field are otherwise the same number.
320func ms_anom_count(v: *i64, n: i64, med: i64, mad: i64, k_milli: i64) -> i64 {
321 if n <= 0 { return 0 - 1 }
322 var c: i64 = 0
323 var i: i64 = 0
324 while i < n {
325 let d: i64 = ms_abs(v[i] - med)
326 if mad <= 0 {
327 if d > 0 { c = c + 1 }
328 } else {
329 if d * MS_MILLI >= k_milli * mad { c = c + 1 }
330 }
331 i = i + 1
332 }
333 return c
334}
335// the strongest candidate's index, so a count is never published without the worklist that makes it
336// actionable.
337func ms_worst_idx(v: *i64, n: i64, med: i64) -> i64 {
338 if n <= 0 { return 0 - 1 }
339 var bi: i64 = 0
340 var bd: i64 = 0 - 1
341 var i: i64 = 0
342 while i < n {
343 let d: i64 = ms_abs(v[i] - med)
344 if d > bd { bd = d; bi = i }
345 i = i + 1
346 }
347 return bi
348}
349
350// THE FAIL-POINT QUERY, on a field of THICKNESSES rather than heights: the minimum and WHERE it is.
351// Same analyser, different quantity -- this is the teapot-wall question, and it composes downstream
352// of nx_wallstat, which already produces exactly such per-face thickness values on a real mesh and
353// which this estate has never once called.
354func ms_min_loc(f: *i64, n: i64, out_idx: *i64) -> i64 {
355 out_idx[0] = 0 - 1
356 if n <= 0 { return 0 - 1 }
357 var bi: i64 = 0
358 var bv: i64 = f[0]
359 var i: i64 = 1
360 while i < n {
361 if f[i] < bv { bv = f[i]; bi = i }
362 i = i + 1
363 }
364 out_idx[0] = bi
365 return bv
366}