nx_sdfprim_lib.nx source
↩ module page · 216 lines · 11530 B
1// nx_sdfprim_lib.nx -- IMPLICIT PRIMITIVES FOR PARAMETRIC ANATOMY FIELDS: integer-millimetre geometry in, a signed
2// distance at SP_FQ units per millimetre out (negative = inside), composed by smooth union and subtraction.
3// Extracted VERBATIM from nx_skullsdf (anatomy AN17, 2026-09-19; the emitted skull mesh is byte-identical before and
4// after the move, proven by sha256) so an arm, a foot, a muscle belly or a soft-tissue layer composes the SAME
5// ellipsoid, capsule, swept arch, pitched plate, walled cylinder and element ring the skull is built from -- and a
6// primitive fixed here is fixed for every field at once instead of in one organ.
7//
8// SIBLING FAMILIES, DECLARED (nx_spendgate listed them; this is the adjudication, not a silent fourth copy):
9// nx_sdf.nx the Q14-millimetre CSG algebra (sphere, box, cylinder, torus, rounded box, smooth booleans),
10// nx_sdfrender.nx the renderer's op-table evaluator (sdf_ellip, sdf_ellip_rot, sdf_lipsweep ...) shaded per pixel.
11// They differ in UNITS (Q14 mm vs integer mm in / SP_FQ out) and in EVALUATION MODEL (op tables vs a composed function
12// the polygonizer and the tuner call millions of times). Consolidating the three under one evaluator is a declared
13// ecosystem rung; until it lands, this file is the family a POLYGONIZED, TUNED field uses.
14//
15// UNITS. Points and geometry are integer millimetres. Every *_q primitive returns distance * SP_FQ; sp_smin's k is in
16// the same output units (a caller passes k_mm * SP_FQ). The dimensionless ratios inside the ellipsoid estimator are
17// carried at SP_Q, and its k1 term at SP_Q squared: at plain SP_Q that term truncates to ZERO for a small offset on a
18// large axis (10*1024/100^2 = 1), which quantised the distance into concentric steps and left the polygonized surface
19// ringed and mottled (measured in the skull lane before the extraction).
20//
21// THE ELLIPSOID IS A FIRST-ORDER ESTIMATOR, ON PURPOSE -- MEASURED, NOT PREFERRED. d = k0 (k0 - 1) / k1 with
22// k0 = |q/r| and k1 = |q/r^2| is exact on a sphere and accurate near an ellipsoid's surface. A Newton refinement
23// (walk d0 along the unit normal g/|g|, re-estimate, add) was tried in the skull lane and LOST at every damping:
24// field fidelity by the `grad` verb, total deviation outside +/-20 percent of a true SDF, FIRST-ORDER 402 permil,
25// undamped Newton 550, half-damped 425; undamped also pushed the mean |grad| from 7 percent LOW to 13 percent HIGH
26// (927 -> 1126) and the max from 1870 to 2121. The reason is worth more than the attempt: the refinement re-evaluates
27// the SAME first-order formula at the stepped point, so it inherits that formula's bias -- ITERATING AN APPROXIMATION
28// WITH ITSELF CONVERGES TO THE APPROXIMATION'S FIXED POINT, NOT TO THE TRUTH. A real fix needs an exact ellipsoid
29// residual (bisection on the Lagrange multiplier), never another pass of the same estimator.
30// LIB, no main. license_tier: ORIGINAL. No hw writes (Rule 26).
31import "nx_syscalls.nx"
32
33const SP_Q: i64 = 1024 // fixed point for the dimensionless ratios inside the estimators
34const SP_FQ: i64 = 256 // output units per millimetre: the sub-millimetre lane
35const SP_PERMILLE: i64 = 1000 // the sine tables' unit, and the unit of a per-mille policy knob
36const SP_LERP_ONE: i64 = 1024 // a capsule's parameter t, Q10
37const SP_INF_MM: i64 = 300 // farther than any point in a 240 mm box: where an element-min starts
38const SP_ARCH_SEG: i64 = 12 // capsule segments a swept arch is polygonised into
39// Bhaskara I's 7th-century sine approximation: sin(x) ~ 4x(180-x) / (40500 - x(180-x)). The numerator's 4 is scaled by
40// 1000 so the result is PER-MILLE; 40500 = 5 * 90^2 belongs to the closed form itself and is not tunable.
41const SP_BHASKARA_NUM_PERMILLE: i64 = 4000
42const SP_BHASKARA_DEN: i64 = 40500
43// element rings: the count and the arc span are clamped so an optimizer cannot request a degenerate arc
44const SP_RING_N_MIN: i64 = 4
45const SP_RING_N_MAX: i64 = 20
46const SP_RING_TH_MIN: i64 = 20
47const SP_RING_TH_MAX: i64 = 88
48const SP_ARCH_TH_MIN: i64 = 10
49const SP_ARCH_TH_MAX: i64 = 90
50
51func sp_isqrt(v: i64) -> i64 { if v<=0 { return 0 } var x: i64=v; var y: i64=(x+1)/2; while y<x { x=y; y=(x+v/x)/2 } return x }
52func sp_abs(v: i64) -> i64 { if v<0 { return 0-v } return v }
53func sp_min(a: i64, b: i64) -> i64 { if a<b { return a } return b }
54func sp_max(a: i64, b: i64) -> i64 { if a>b { return a } return b }
55
56// integer sine: degrees in, per-mille out; exact at 0/90/180, under 2 permil error between
57func sp_isin(a: i64) -> i64 {
58 var t: i64 = a
59 var sg: i64 = 1
60 if t < 0 { sg = 0-1; t = 0-t }
61 if t > 180 { t = 180 }
62 let q: i64 = t*(180-t)
63 let s: i64 = SP_BHASKARA_NUM_PERMILLE*q/(SP_BHASKARA_DEN - q)
64 if sg < 0 { return 0-s }
65 return s
66}
67func sp_icos(a: i64) -> i64 {
68 var t: i64 = a
69 if t < 0 { t = 0-t }
70 return sp_isin(90 - t)
71}
72
73// ellipsoid centred (cx,cy,cz) with semi-axes (rx,ry,rz): the first-order estimator above, output at SP_FQ
74func sp_ellipsoid_q(px: i64, py: i64, pz: i64, cx: i64, cy: i64, cz: i64, rx: i64, ry: i64, rz: i64) -> i64 {
75 var ax: i64 = rx; var ay: i64 = ry; var az: i64 = rz
76 if ax < 1 { ax = 1 }
77 if ay < 1 { ay = 1 }
78 if az < 1 { az = 1 }
79 let qx: i64 = px-cx
80 let qy: i64 = py-cy
81 let qz: i64 = pz-cz
82 let n0x: i64 = qx*SP_Q/ax
83 let n0y: i64 = qy*SP_Q/ay
84 let n0z: i64 = qz*SP_Q/az
85 let k0: i64 = sp_isqrt(n0x*n0x + n0y*n0y + n0z*n0z)
86 if k0 == 0 { return (0 - sp_min(sp_min(ax,ay),az))*SP_FQ }
87 let n1x: i64 = qx*SP_Q*SP_Q/(ax*ax)
88 let n1y: i64 = qy*SP_Q*SP_Q/(ay*ay)
89 let n1z: i64 = qz*SP_Q*SP_Q/(az*az)
90 var k1: i64 = sp_isqrt(n1x*n1x + n1y*n1y + n1z*n1z)
91 if k1 < 1 { k1 = 1 }
92 return k0*(k0-SP_Q)*SP_FQ/k1
93}
94// capsule between (ax,ay,az) and (bx,by,bz) with radius r: the squared radial distance is scaled by SP_FQ^2 BEFORE the
95// integer root, so the root itself carries sub-millimetre precision (an isqrt of a mm^2 quantity can never return finer
96// than 1 mm). The right primitive for an ARCH -- which is what the zygomatic actually is
97func sp_capsule_q(px: i64, py: i64, pz: i64, ax: i64, ay: i64, az: i64, bx: i64, by: i64, bz: i64, r: i64) -> i64 {
98 let ex: i64 = bx-ax; let ey: i64 = by-ay; let ez: i64 = bz-az
99 let den: i64 = ex*ex + ey*ey + ez*ez
100 var t: i64 = 0
101 if den > 0 {
102 t = ((px-ax)*ex + (py-ay)*ey + (pz-az)*ez) * SP_LERP_ONE / den
103 if t < 0 { t = 0 }
104 if t > SP_LERP_ONE { t = SP_LERP_ONE }
105 }
106 let qx: i64 = ax + ex*t/SP_LERP_ONE
107 let qy: i64 = ay + ey*t/SP_LERP_ONE
108 let qz: i64 = az + ez*t/SP_LERP_ONE
109 let dx: i64 = px-qx; let dy: i64 = py-qy; let dz: i64 = pz-qz
110 return sp_isqrt((dx*dx + dy*dy + dz*dz)*SP_FQ*SP_FQ) - r*SP_FQ
111}
112// smooth union: the operation that makes separate bones ONE SURFACE instead of touching objects; k in output units
113func sp_smin(a: i64, b: i64, k: i64) -> i64 {
114 if k <= 0 { return sp_min(a,b) }
115 var h: i64 = k - sp_abs(a-b)
116 if h < 0 { h = 0 }
117 return sp_min(a,b) - h*h/(4*k)
118}
119// subtraction: carve B out of A. This is how an orbit becomes a real OPENING rather than a dent
120func sp_sub(a: i64, b: i64) -> i64 { return sp_max(a, 0-b) }
121
122// THE STRUCTURED PRIMITIVE: N ellipsoid elements along an elliptical arc in the xz plane at height ty (a dental arch:
123// tooth width tw, height tht, tprof widening the arc-end elements over the centre ones, per-mille). The same construct
124// recursed -- elements along a curve, each element itself carrying elements -- is branch-on-trunk, leaf-on-branch.
125// best starts at SP_INF_MM*SP_FQ*SP_FQ because an infinity a measurement can reach is not an infinity
126func sp_ring_q(px: i64, py: i64, pz: i64, ty: i64, tz0: i64, ax: i64, az: i64, thmax: i64, n: i64, tw: i64, tht: i64, tprof: i64) -> i64 {
127 var nn: i64 = n
128 if nn < SP_RING_N_MIN { nn = SP_RING_N_MIN }
129 if nn > SP_RING_N_MAX { nn = SP_RING_N_MAX }
130 var tm2: i64 = thmax
131 if tm2 < SP_RING_TH_MIN { tm2 = SP_RING_TH_MIN }
132 if tm2 > SP_RING_TH_MAX { tm2 = SP_RING_TH_MAX }
133 var best: i64 = SP_INF_MM*SP_FQ*SP_FQ
134 var k: i64 = 0
135 while k < nn {
136 let thk: i64 = 0-tm2 + (2*tm2*k)/(nn-1)
137 let s: i64 = sp_isin(thk)
138 let c: i64 = sp_icos(thk)
139 let cx: i64 = ax*s/SP_PERMILLE
140 let cz: i64 = tz0 + az*c/SP_PERMILLE
141 var athk: i64 = thk
142 if athk < 0 { athk = 0-athk }
143 let wk: i64 = tw + tw*tprof*athk/(tm2*1000)
144 let d: i64 = sp_ellipsoid_q(px,py,pz, cx, ty, cz, wk/2+1, tht/2+1, wk/2+1)
145 if d < best { best = d }
146 k = k + 1
147 }
148 return best
149}
150// a capsule swept along an ELLIPTICAL ARC in the xz plane: a polyline of SP_ARCH_SEG capsule segments over theta in
151// [-thmax, +thmax] degrees, x = ax sin, z = cz + az cos, y = cy + ylift*|sin| (the ends rise toward the angles). An
152// alveolar process, a mandibular body, a rib, a pelvic brim
153func sp_arch_q(px: i64, py: i64, pz: i64, cy: i64, cz: i64, ax: i64, az: i64, thmax: i64, ylift: i64, r: i64) -> i64 {
154 var tm: i64 = thmax
155 if tm < SP_ARCH_TH_MIN { tm = SP_ARCH_TH_MIN }
156 if tm > SP_ARCH_TH_MAX { tm = SP_ARCH_TH_MAX }
157 var best: i64 = SP_INF_MM*SP_FQ*SP_FQ
158 var lx: i64 = 0
159 var ly: i64 = 0
160 var lz: i64 = 0
161 var k: i64 = 0
162 while k <= SP_ARCH_SEG {
163 let th: i64 = 0-tm + (2*tm*k)/SP_ARCH_SEG
164 let s: i64 = sp_isin(th)
165 let c: i64 = sp_icos(th)
166 var sa: i64 = s
167 if sa < 0 { sa = 0-sa }
168 let x: i64 = ax*s/SP_PERMILLE
169 let z: i64 = cz + az*c/SP_PERMILLE
170 let y: i64 = cy + ylift*sa/SP_PERMILLE
171 if k > 0 {
172 let d: i64 = sp_capsule_q(px,py,pz, lx,ly,lz, x,y,z, r)
173 if d < best { best = d }
174 }
175 lx = x
176 ly = y
177 lz = z
178 k = k + 1
179 }
180 return best
181}
182// an ellipsoid pitched about the x axis by tilt degrees (positive = the lower end comes forward), evaluated by rotating
183// the point into the ellipsoid's own frame: a nasal plate, a scapular blade, a sloping tendon sheet
184func sp_ellipsoid_tiltx_q(px: i64, py: i64, pz: i64, cx: i64, cy: i64, cz: i64, rx: i64, ry: i64, rz: i64, tilt: i64) -> i64 {
185 let s: i64 = sp_isin(tilt)
186 let c: i64 = sp_icos(tilt)
187 let qy: i64 = py - cy
188 let qz: i64 = pz - cz
189 let ry2: i64 = (qy*c - qz*s)/SP_PERMILLE
190 let rz2: i64 = (qy*s + qz*c)/SP_PERMILLE
191 return sp_ellipsoid_q(px - cx, ry2, rz2, 0, 0, 0, rx, ry, rz)
192}
193// an elliptic cylinder along z with a flat BACK wall at zback and no front wall -- an OPENING with a floor: seen from the
194// front its rim is the ellipse, its floor is the wall, and the solid behind the wall stays solid. The 2D ellipse distance
195// is the same first-order estimator (exact on a circle); the wall is a half-space
196func sp_cylz_q(px: i64, py: i64, pz: i64, cx: i64, cy: i64, rx: i64, ry: i64, zback: i64) -> i64 {
197 var ax: i64 = rx
198 var ay: i64 = ry
199 if ax < 1 { ax = 1 }
200 if ay < 1 { ay = 1 }
201 let qx: i64 = px-cx
202 let qy: i64 = py-cy
203 let n0x: i64 = qx*SP_Q/ax
204 let n0y: i64 = qy*SP_Q/ay
205 let k0: i64 = sp_isqrt(n0x*n0x + n0y*n0y)
206 var d2: i64 = (0 - sp_min(ax,ay))*SP_FQ
207 if k0 > 0 {
208 let n1x: i64 = qx*SP_Q*SP_Q/(ax*ax)
209 let n1y: i64 = qy*SP_Q*SP_Q/(ay*ay)
210 var k1: i64 = sp_isqrt(n1x*n1x + n1y*n1y)
211 if k1 < 1 { k1 = 1 }
212 d2 = k0*(k0-SP_Q)*SP_FQ/k1
213 }
214 let dz: i64 = (zback - pz)*SP_FQ
215 return sp_max(d2, dz)
216}