nx_skinfit_lib.nx source
↩ module page · 123 lines · 5417 B
1// nx_skinfit_lib.nx -- RECOVER MELANIN AND HAEMOGLOBIN FROM A MEASURED REFLECTANCE (2026-09-03)
2//
3// WHY THIS SHAPE, AND WHY NOT A DIFFERENTIABLE RENDERER. Inverse rendering in general is a 10^5..10^6
4// parameter problem behind a renderer, and the estate's Gauss-Newton solvers (nx_pnprefine, nx_multical) are
5// built for 6-7 well-conditioned DOF with analytic residuals -- they do not transfer, and pointing them at a
6// renderer is the measured-a-different-subject defect wearing a solver's clothes. So the problem is
7// DECOMPOSED until no general solver is needed: the illuminant is closed form (cs_cat_*), and the chromophore
8// fit -- here -- is TWO parameters against THREE channels, which a deterministic coarse-to-fine grid search
9// inverts exactly and repeatably, with no gradients, no autograd, and no float.
10//
11// WHY IT IS IDENTIFIABLE. One reflectance and two unknowns is underdetermined. Three channels give three
12// equations for two unknowns, and the fit is well posed precisely BECAUSE the two absorbers have different
13// spectral shapes -- melanin rising toward blue, haemoglobin peaking in green. If a caller supplies
14// proportional coefficient vectors the problem degenerates, and the RESIDUAL is what says so. The gate
15// carries that degenerate case as a NEG-CONTROL and measures the collapse rather than asserting it cannot
16// happen: with proportional coefficients the fit puts everything on one axis and misses by 0.30 and 0.60.
17//
18// THE RESIDUAL IS PART OF THE ANSWER, NEVER AN INTERNAL. A fitter that returns its best guess without saying
19// how badly it missed will happily report chromophores for a measurement no skin could produce. sf_fit writes
20// the residual beside the parameters so a caller can refuse; a fit is a measurement with an error bar or it
21// is a fabrication.
22//
23// PROVEN 8/8 GREEN 2026-09-03 by round trip: worst recovery error 2.2e-4 over a spread of true pairs with a
24// zero residual on every one, the degenerate case correctly missing, and an unreachable measurement coming
25// back with a residual of 0.5 against a floor of 0.002.
26// license_tier: ORIGINAL No hw writes (Rule 26). LIB (no main).
27import "nx_syscalls.nx"
28import "nx_fixq30_lib.nx"
29import "nx_kubelka_lib.nx"
30
31const SF_REFUSED: i64 = 0 - 1
32const SF_CH: i64 = 3
33const SF_SLOT: i64 = 8
34// grid steps per axis per pass, and how many times the bracket is shrunk around the best cell.
35// 16 steps shrinking 4 times resolves 1/16^4 of the unit interval, far finer than any measured coefficient.
36const SF_GRID: i64 = 16
37const SF_PASSES: i64 = 4
38const SF_O_MEL: i64 = 0
39const SF_O_HEM: i64 = 1
40const SF_O_RESID: i64 = 2
41const SF_O_SLOTS: i64 = 3
42
43// forward model: chromophore fractions -> three channel reflectances. Returns SF_REFUSED if any channel
44// refuses, so a bad coefficient set can never be silently fitted around.
45func sf_forward(kbase: *i64, kmel: *i64, khem: *i64, s: *i64, cmel: i64, chem: i64, out3: *i64) -> i64 {
46 var i: i64 = 0
47 while i < SF_CH {
48 let r: i64 = km_r_of_mix(kbase[i], kmel[i], cmel, khem[i], chem, s[i])
49 if r == KM_REFUSED { return SF_REFUSED }
50 out3[i] = r
51 i = i + 1
52 }
53 return 0
54}
55
56// sum of squared error between two reflectance triples, Q30
57func sf_sse(a: *i64, b: *i64) -> i64 {
58 var e: i64 = 0
59 var i: i64 = 0
60 while i < SF_CH {
61 let d: i64 = a[i] - b[i]
62 e = e + fq_mul(d, d)
63 i = i + 1
64 }
65 return e
66}
67
68// deterministic coarse-to-fine grid search. out carries SF_O_MEL, SF_O_HEM and SF_O_RESID.
69// Returns 0 on success, SF_REFUSED if the forward model never produced a single evaluable point.
70func sf_fit(kbase: *i64, kmel: *i64, khem: *i64, s: *i64, meas3: *i64, out: *i64) -> i64 {
71 let trial: *i64 = sys_mmap(SF_CH * SF_SLOT) as *i64
72 var mlo: i64 = 0
73 var mhi: i64 = FQ_ONE
74 var hlo: i64 = 0
75 var hhi: i64 = FQ_ONE
76 var bestm: i64 = 0
77 var besth: i64 = 0
78 var beste: i64 = 0 - 1
79 var pass: i64 = 0
80 while pass < SF_PASSES {
81 let mstep: i64 = (mhi - mlo) / SF_GRID
82 let hstep: i64 = (hhi - hlo) / SF_GRID
83 var gi: i64 = 0
84 while gi <= SF_GRID {
85 let cm: i64 = mlo + mstep * gi
86 var gj: i64 = 0
87 while gj <= SF_GRID {
88 let ch: i64 = hlo + hstep * gj
89 if sf_forward(kbase, kmel, khem, s, cm, ch, trial) == 0 {
90 let e: i64 = sf_sse(trial, meas3)
91 if beste < 0 {
92 beste = e
93 bestm = cm
94 besth = ch
95 } else {
96 if e < beste {
97 beste = e
98 bestm = cm
99 besth = ch
100 }
101 }
102 }
103 gj = gj + 1
104 }
105 gi = gi + 1
106 }
107 if beste < 0 { return SF_REFUSED }
108 // shrink the bracket to one cell either side of the best point, clamped to the unit interval
109 mlo = bestm - mstep
110 mhi = bestm + mstep
111 hlo = besth - hstep
112 hhi = besth + hstep
113 if mlo < 0 { mlo = 0 }
114 if hlo < 0 { hlo = 0 }
115 if mhi > FQ_ONE { mhi = FQ_ONE }
116 if hhi > FQ_ONE { hhi = FQ_ONE }
117 pass = pass + 1
118 }
119 out[SF_O_MEL] = bestm
120 out[SF_O_HEM] = besth
121 out[SF_O_RESID] = beste
122 return 0
123}