nx_gsplat4d_lib.nx source
↩ module page · 356 lines · 16221 B
1// nx_gsplat4d_lib.nx -- ★4D SPATIOTEMPORAL GAUSSIAN SPLATTING: the temporal axis the estate did not have.
2//
3// WHAT WAS ALREADY HERE (measured 2026-08-14 three independent ways: nx_capsearch over 1033 registered
4// tools, a coverage-complete grep of 23,023 sources, and reading the source): nx_gsplat is a real, gated,
5// sovereign 3D splatter -- isotropic blobs AND anisotropic surfels whose 2D covariance Sigma' is computed
6// EXACTLY by projecting the principal axes through the camera (the Jacobian, applied as a finite
7// projection rather than a symbolic J R S S^T R^T J^T), a conic power test against an integer exp-LUT, a
8// counting sort by depth, and front-to-back compositing C = sum c_i a_i prod(1-a_j). That is the
9// rasterizer. It is NOT re-implemented here.
10//
11// WHAT WAS ABSENT: time. Every Gaussian in the estate was static. This organ adds the 4th axis and
12// nothing else -- it deforms state and hands the result to gs_render_aniso. One rasterizer, still.
13//
14// THE KERNEL is the standard 4DGS temporal radial basis, integer-exact:
15// x_t = x_0 + v * (t - t_mu) * exp(-(t - t_mu)^2 / (2 * t_sigma^2))
16// ★AND THE ENVELOPE COSTS NOTHING NEW: the incumbent's exp-LUT is explut[k] = GFXA*exp(-k/(2*GLUTU)), so
17// indexing it at k = GLUTU * dt^2/sigma^2 yields exp(-dt^2/(2 sigma^2)) EXACTLY -- the Gaussian envelope,
18// from the table the renderer already builds. A second exp table would have been a duplicate ruler.
19//
20// ★★NO RE-DECLARED CONSTANTS. The fixed-point scale, the LUT resolution and the record stride are FACTS
21// OF THE RASTERIZER and are read from it (gs_fxa / gs_lutu / gs_stride_aniso). An earlier revision of this
22// file hardcoded 256, 16 and 12; nx_magic at threshold 2 found them. They are not cosmetic: each was a
23// SECOND RULER, and the failure mode of a second ruler is that the two disagree silently, compile
24// cleanly, and render something plausible. There is no symptom to notice.
25//
26// ★★THE ARENA COUNTS ITSELF. Every allocation goes through g4_alloc, which accumulates the byte total as
27// a side effect of allocating -- so g4_bytes() cannot drift from reality. The previous revision multiplied
28// a HAND-WRITTEN COUNT of the field arrays; adding a field and forgetting to bump it would have silently
29// under-reported the resource envelope, which is a fabricated number wearing the shape of a measurement.
30//
31// LAYOUT: Structure-of-Arrays, one array per field -- a field sweep touches contiguous memory instead of
32// striding over a record. Arrays are lazily allocated once against a declared capacity.
33//
34// ⚠DELIBERATE, DECLARED DIFFERENCES FROM A LITERAL FIELD-LIST TRANSCRIPTION -- each is a defect avoided:
35// - POSITION uses the signed ramp x envelope above (the kernel as specified).
36// - RADIUS uses the envelope ALONE, not the ramp. A signed ramp on a radius crosses zero and inverts
37// the splat, which the rasterizer cannot represent. A breathing surfel is the honest analogue.
38// - NORMAL uses ramp x envelope and is then RENORMALISED, so the surfel tips without its length
39// drifting. It is also the G-buffer normal the lighting rung will read.
40// - QUATERNION and 3-axis SCALE are NOT stored. The incumbent surfel is a DISK (normal + in-plane
41// radius): an in-plane spin is invisible on a disk and a 3-axis scale has nowhere to land. Storing
42// them would be fields no consumer reads. They arrive WITH the elliptical-surfel rung.
43// - SPHERICAL HARMONICS beyond DC are NOT stored: gs_render_aniso reads one flat colour, so bands 1..3
44// would be dead weight until the rasterizer evaluates view direction. ⚠AND THE SPEC AS COMMONLY
45// WRITTEN IS INCONSISTENT: 3rd-order SH is 16 coefficients PER COLOUR CHANNEL = 48, not 16.
46// license_tier: ORIGINAL No hw writes (Rule 26).
47import "nx_syscalls.nx"
48import "nx_gsplat.nx"
49
50// t and t_mu are counted in TICKS. This is a UNIT DEFINITION (like choosing milliseconds), not a tuned
51// threshold -- it fixes what one unit of time means for the velocity term, and it is exported so callers
52// and gates read it instead of assuming a number.
53const G4_TSCALE: i64 = 1000
54// A prefetch cache needs the frame in use PLUS at least one ahead, or it is not a prefetch cache at all.
55// This is a structural floor that follows from the definition, not a tuning knob: depth is CHOSEN BY THE
56// CALLER, because how far ahead to run is a property of the workload, not of this code.
57const G4_RING_MIN: i64 = 2
58const G4_ERRFD: i64 = 2
59const G4_W64: i64 = 8 // the machine word this SoA is built from
60const G4_E_CAP: i64 = 4
61const G4_E_SIGMA: i64 = 3
62const G4_E_RANGE: i64 = 5
63const G4_E_RING: i64 = 6
64
65static G4_CAP: i64
66static G4_N: i64
67static G4_RINGN: i64
68static G4_ALLOC: i64
69static G4_X: i64
70static G4_Y: i64
71static G4_Z: i64
72static G4_NX: i64
73static G4_NY: i64
74static G4_NZ: i64
75static G4_RT: i64
76static G4_CR: i64
77static G4_CG: i64
78static G4_CB: i64
79static G4_OP: i64
80static G4_TMU: i64
81static G4_TSIG: i64
82static G4_VX: i64
83static G4_VY: i64
84static G4_VZ: i64
85static G4_DRT: i64
86static G4_DNX: i64
87static G4_DNY: i64
88static G4_DNZ: i64
89static G4_EXPLUT: i64
90static G4_EXPN: i64
91static G4_RBUF: i64
92static G4_RT_T: i64
93static G4_RVALID: i64
94static G4_RHEAD: i64
95
96func g4_slen(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } return n }
97func g4_err(s: *u8) -> i64 { sys_write(G4_ERRFD, s, g4_slen(s)); return 0 }
98
99// ★EVERY allocation goes through here, so the byte total is accumulated AS A SIDE EFFECT OF ALLOCATING.
100// You cannot add an array to this organ without its cost appearing in g4_bytes(). That is the difference
101// between a measurement and a number somebody remembered to update.
102func g4_alloc(n: i64) -> i64 {
103 G4_ALLOC = G4_ALLOC + n
104 return sys_mmap(n) as i64
105}
106
107// allocate the SoA arena for cap gaussians with a caller-chosen prefetch depth.
108func g4_init(cap: i64, ring: i64) -> i64 {
109 if cap < 1 { g4_err("GSPLAT4D-REFUSE capacity must be positive\n" as *u8); return 0 - G4_E_CAP }
110 if ring < G4_RING_MIN { g4_err("GSPLAT4D-REFUSE prefetch depth must hold the frame in use plus at least one ahead\n" as *u8); return 0 - G4_E_RING }
111 if G4_CAP >= cap { if G4_RINGN == ring { G4_N = 0; return 0 } }
112 G4_ALLOC = 0
113 let w: i64 = cap*G4_W64
114 G4_X = g4_alloc(w)
115 G4_Y = g4_alloc(w)
116 G4_Z = g4_alloc(w)
117 G4_NX = g4_alloc(w)
118 G4_NY = g4_alloc(w)
119 G4_NZ = g4_alloc(w)
120 G4_RT = g4_alloc(w)
121 G4_CR = g4_alloc(w)
122 G4_CG = g4_alloc(w)
123 G4_CB = g4_alloc(w)
124 G4_OP = g4_alloc(w)
125 G4_TMU = g4_alloc(w)
126 G4_TSIG = g4_alloc(w)
127 G4_VX = g4_alloc(w)
128 G4_VY = g4_alloc(w)
129 G4_VZ = g4_alloc(w)
130 G4_DRT = g4_alloc(w)
131 G4_DNX = g4_alloc(w)
132 G4_DNY = g4_alloc(w)
133 G4_DNZ = g4_alloc(w)
134 G4_EXPN = gs_expn()
135 G4_EXPLUT = g4_alloc(G4_EXPN*G4_W64)
136 gs_build_explut(G4_EXPLUT as *i64)
137 G4_RBUF = g4_alloc(ring*cap*gs_stride_aniso()*G4_W64)
138 G4_RT_T = g4_alloc(ring*G4_W64)
139 G4_RVALID = g4_alloc(ring*G4_W64)
140 G4_CAP = cap
141 G4_RINGN = ring
142 G4_N = 0
143 G4_RHEAD = 0
144 return 0
145}
146func g4_cap() -> i64 { return G4_CAP }
147func g4_n() -> i64 { return G4_N }
148func g4_bytes() -> i64 { return G4_ALLOC }
149func g4_ring_depth() -> i64 { return G4_RINGN }
150func g4_tscale() -> i64 { return G4_TSCALE }
151func g4_ring_min() -> i64 { return G4_RING_MIN }
152
153// declare one Gaussian's REST state -- the same tuple gs_set_aniso takes, so a static scene ported here
154// is a field-for-field move with no reinterpretation.
155func g4_set_rest(i: i64, x: i64, y: i64, z: i64, nx: i64, ny: i64, nz: i64, rt: i64, r: i64, g: i64, b: i64, op: i64) -> i64 {
156 if i < 0 { return 0 - G4_E_RANGE }
157 if i >= G4_CAP { g4_err("GSPLAT4D-REFUSE index past declared capacity\n" as *u8); return 0 - G4_E_RANGE }
158 let ax: *i64 = G4_X as *i64
159 let ay: *i64 = G4_Y as *i64
160 let az: *i64 = G4_Z as *i64
161 let bx: *i64 = G4_NX as *i64
162 let by: *i64 = G4_NY as *i64
163 let bz: *i64 = G4_NZ as *i64
164 let ar: *i64 = G4_RT as *i64
165 let cr: *i64 = G4_CR as *i64
166 let cg: *i64 = G4_CG as *i64
167 let cb: *i64 = G4_CB as *i64
168 let ao: *i64 = G4_OP as *i64
169 ax[i]=x; ay[i]=y; az[i]=z; bx[i]=nx; by[i]=ny; bz[i]=nz; ar[i]=rt; cr[i]=r; cg[i]=g; cb[i]=b; ao[i]=op
170 // a Gaussian with no temporal row is STATIC BY DEFAULT: a unit-width window and zero velocity, so the
171 // envelope exists but moves nothing. A static scene therefore bakes identically at every t.
172 let tm: *i64 = G4_TMU as *i64
173 let ts: *i64 = G4_TSIG as *i64
174 let vx: *i64 = G4_VX as *i64
175 let vy: *i64 = G4_VY as *i64
176 let vz: *i64 = G4_VZ as *i64
177 let dr: *i64 = G4_DRT as *i64
178 let dx: *i64 = G4_DNX as *i64
179 let dy: *i64 = G4_DNY as *i64
180 let dz: *i64 = G4_DNZ as *i64
181 tm[i]=0; ts[i]=1; vx[i]=0; vy[i]=0; vz[i]=0; dr[i]=0; dx[i]=0; dy[i]=0; dz[i]=0
182 if i >= G4_N { G4_N = i + 1 }
183 return 0
184}
185
186// declare the temporal window. sigma <= 0 is REFUSED: it is a division by zero in the envelope, and a
187// silently-clamped sigma would render a plausible frame that means nothing.
188func g4_set_time(i: i64, tmu: i64, tsig: i64) -> i64 {
189 if i < 0 { return 0 - G4_E_RANGE }
190 if i >= G4_N { g4_err("GSPLAT4D-REFUSE temporal row names an undeclared gaussian\n" as *u8); return 0 - G4_E_RANGE }
191 if tsig <= 0 { g4_err("GSPLAT4D-REFUSE t_sigma must be positive -- a zero-width temporal window is a division by zero, not a still frame\n" as *u8); return 0 - G4_E_SIGMA }
192 let tm: *i64 = G4_TMU as *i64
193 let ts: *i64 = G4_TSIG as *i64
194 tm[i] = tmu
195 ts[i] = tsig
196 return 0
197}
198
199// declare the deformation: linear velocity, radius pulse amplitude, and normal tip rate.
200func g4_set_motion(i: i64, vx: i64, vy: i64, vz: i64, drt: i64, dnx: i64, dny: i64, dnz: i64) -> i64 {
201 if i < 0 { return 0 - G4_E_RANGE }
202 if i >= G4_N { g4_err("GSPLAT4D-REFUSE motion row names an undeclared gaussian\n" as *u8); return 0 - G4_E_RANGE }
203 let ax: *i64 = G4_VX as *i64
204 let ay: *i64 = G4_VY as *i64
205 let az: *i64 = G4_VZ as *i64
206 let ar: *i64 = G4_DRT as *i64
207 let bx: *i64 = G4_DNX as *i64
208 let by: *i64 = G4_DNY as *i64
209 let bz: *i64 = G4_DNZ as *i64
210 ax[i]=vx; ay[i]=vy; az[i]=vz; ar[i]=drt; bx[i]=dnx; by[i]=dny; bz[i]=dnz
211 return 0
212}
213
214// ★THE TEMPORAL ENVELOPE. exp(-(t-t_mu)^2 / (2 sigma^2)) read straight out of the RASTERIZER'S OWN LUT:
215// explut[k] = gs_fxa()*exp(-k/(2*gs_lutu())), so k = gs_lutu()*dt^2/sigma^2 gives exactly the envelope.
216// Past the end of the table the weight is 0 -- a hard temporal horizon that is a PROPERTY OF THE TABLE,
217// derivable as g4_horizon(), not a number chosen here.
218func g4_envelope(i: i64, t: i64) -> i64 {
219 let tm: *i64 = G4_TMU as *i64
220 let ts: *i64 = G4_TSIG as *i64
221 let sg: i64 = ts[i]
222 if sg <= 0 { return 0 }
223 var dt: i64 = t - tm[i]
224 if dt < 0 { dt = 0 - dt }
225 let k: i64 = dt*dt*gs_lutu()/(sg*sg)
226 if k >= G4_EXPN { return 0 }
227 let el: *i64 = G4_EXPLUT as *i64
228 return el[k]
229}
230// the full-weight value IS the table's first entry, by definition. Never re-stated as a literal.
231func g4_wfull() -> i64 { let el: *i64 = G4_EXPLUT as *i64; return el[0] }
232// the table itself, so a referee can check the envelope against the ENTRY ITS INDEX FORMULA NAMES rather
233// than against a hand-picked tolerance band. The thing that can actually be wrong here is the index
234// arithmetic, and only an exact comparison against the table tests it.
235func g4_lut(k: i64) -> i64 {
236 if k < 0 { return 0 }
237 if k >= G4_EXPN { return 0 }
238 let el: *i64 = G4_EXPLUT as *i64
239 return el[k]
240}
241// how many sigma out the table still carries weight: k < EXPN and k = lutu*(dt/sigma)^2, so the horizon is
242// sqrt(EXPN/lutu) sigma. DERIVED from the table, so it tracks any change to it.
243func g4_horizon_sigma() -> i64 { return gs_isqrt(G4_EXPN/gs_lutu()) }
244
245// ★BAKE the whole cloud at time t into the record gs_render_aniso consumes. This is the ONLY bridge
246// between 4D state and pixels, and it deliberately produces the incumbent's format so that the
247// rasterizer, its sort, its conic test and its blend are reused unchanged.
248func g4_bake(t: i64, out: *i64) -> i64 {
249 if G4_N < 1 { return 0 }
250 let ax: *i64 = G4_X as *i64
251 let ay: *i64 = G4_Y as *i64
252 let az: *i64 = G4_Z as *i64
253 let bx: *i64 = G4_NX as *i64
254 let by: *i64 = G4_NY as *i64
255 let bz: *i64 = G4_NZ as *i64
256 let ar: *i64 = G4_RT as *i64
257 let cr: *i64 = G4_CR as *i64
258 let cg: *i64 = G4_CG as *i64
259 let cb: *i64 = G4_CB as *i64
260 let ao: *i64 = G4_OP as *i64
261 let tm: *i64 = G4_TMU as *i64
262 let vx: *i64 = G4_VX as *i64
263 let vy: *i64 = G4_VY as *i64
264 let vz: *i64 = G4_VZ as *i64
265 let dr: *i64 = G4_DRT as *i64
266 let dnx: *i64 = G4_DNX as *i64
267 let dny: *i64 = G4_DNY as *i64
268 let dnz: *i64 = G4_DNZ as *i64
269 let fxa: i64 = gs_fxa()
270 let den: i64 = G4_TSCALE * fxa
271 var i: i64 = 0
272 while i < G4_N {
273 let w: i64 = g4_envelope(i, t)
274 let dts: i64 = t - tm[i]
275 let px: i64 = ax[i] + vx[i]*dts*w/den
276 let py: i64 = ay[i] + vy[i]*dts*w/den
277 let pz: i64 = az[i] + vz[i]*dts*w/den
278 var rt: i64 = ar[i] + dr[i]*w/fxa
279 if rt < 1 { rt = 1 }
280 var nx2: i64 = bx[i] + dnx[i]*dts*w/den
281 var ny2: i64 = by[i] + dny[i]*dts*w/den
282 var nz2: i64 = bz[i] + dnz[i]*dts*w/den
283 let nl: i64 = gs_isqrt(nx2*nx2 + ny2*ny2 + nz2*nz2)
284 if nl > 0 {
285 nx2 = nx2*fxa/nl
286 ny2 = ny2*fxa/nl
287 nz2 = nz2*fxa/nl
288 }
289 gs_set_aniso(out, i, px, py, pz, nx2, ny2, nz2, rt, cr[i], cg[i], cb[i], ao[i])
290 i = i + 1
291 }
292 return G4_N
293}
294
295// ===== ★THE BOUNDED FRAME CACHE =====
296// A caller-chosen number of slots; each holds one fully baked frame and the t it was baked at. prime()
297// fills the ring ahead of playback, find() answers whether a frame is already resident, advance()
298// recycles the OLDEST slot. Cost is bounded and reported by g4_bytes(); it never grows during playback.
299//
300// ⚠WHAT THIS IS NOT, STATED PLAINLY RATHER THAN IMPLIED: it is not lock-free, and it does not claim to be.
301// This runtime is single-threaded and deterministic, so there is no competing writer for a lock-free
302// structure to protect against, and asserting that property about a one-threaded program would be a
303// fabricated guarantee. Likewise there is NO GPU, NO host-visible VRAM mapping and NO PCIe pacing here:
304// the estate rasterizes on the CPU and no GPU binding exists to stream into. Named gap, named blocker.
305func g4_ring_init() -> i64 {
306 let rv: *i64 = G4_RVALID as *i64
307 var s: i64 = 0
308 while s < G4_RINGN { rv[s] = 0; s = s + 1 }
309 G4_RHEAD = 0
310 return 0
311}
312func g4_ring_frame(s: i64) -> *i64 {
313 return ((G4_RBUF as i64) + s*G4_CAP*gs_stride_aniso()*G4_W64) as *i64
314}
315func g4_ring_time(s: i64) -> i64 { let rt: *i64 = G4_RT_T as *i64; return rt[s] }
316func g4_ring_valid(s: i64) -> i64 { let rv: *i64 = G4_RVALID as *i64; return rv[s] }
317func g4_ring_head() -> i64 { return G4_RHEAD }
318
319func g4_ring_fill(s: i64, t: i64) -> i64 {
320 let rt: *i64 = G4_RT_T as *i64
321 let rv: *i64 = G4_RVALID as *i64
322 g4_bake(t, g4_ring_frame(s))
323 rt[s] = t
324 rv[s] = 1
325 return 0
326}
327func g4_ring_prime(t0: i64, step: i64) -> i64 {
328 g4_ring_init()
329 var s: i64 = 0
330 while s < G4_RINGN { g4_ring_fill(s, t0 + s*step); s = s + 1 }
331 G4_RHEAD = 0
332 return 0
333}
334// is time t already resident? returns the slot, or -1. A cache that cannot say MISS is a cache that lies.
335func g4_ring_find(t: i64) -> i64 {
336 let rt: *i64 = G4_RT_T as *i64
337 let rv: *i64 = G4_RVALID as *i64
338 var s: i64 = 0
339 while s < G4_RINGN {
340 if rv[s] == 1 { if rt[s] == t { return s } }
341 s = s + 1
342 }
343 return 0 - 1
344}
345// recycle the OLDEST slot to the frame one step beyond the newest. Bounded: head walks modulo the depth
346// and no allocation happens here at all.
347func g4_ring_advance(step: i64) -> i64 {
348 let rt: *i64 = G4_RT_T as *i64
349 var newest: i64 = rt[0]
350 var s: i64 = 1
351 while s < G4_RINGN { if rt[s] > newest { newest = rt[s] } s = s + 1 }
352 let victim: i64 = G4_RHEAD
353 g4_ring_fill(victim, newest + step)
354 G4_RHEAD = (G4_RHEAD + 1) % G4_RINGN
355 return victim
356}