nx_ng_gsplat.nx source

1// nx_ng_gsplat.nx -- CAP-GAUSSIAN-SPLAT, generation 10: 3D Gaussian Splatting (Kerbl et al. 2023, arxiv 2// 2308.04079, Inria/EU). The other pillar of the radiance-field frontier: instead of an MLP queried along 3// rays (NeRF), the scene is EXPLICIT Gaussian primitives that are SPLATTED (each projects to a footprint) and 4// ALPHA-OVER composited front-to-back -- real-time, and differentiable so the primitives can be optimized. 5// 6// This rung REUSES the alpha-over compositor proven in CAP-VOLUME-RENDER (T_g = prod(1-a_h); C = sum T_g*a_g*c_g) 7// and adds the 3DGS-specific pieces: the explicit Gaussian PRIMITIVE + its SPLAT FOOTPRINT (a Gaussian falloff 8// exp(-(x-mu)^2/2sigma^2) via the autograd module's fixed-point nfa_fxexp). 1D image, Q16, no float. 9// 10// Honest gated proof (HARD evidence): N Gaussians (centers mu_g, width sigma) splat onto a P-pixel image; their 11// per-pixel opacities a_g(p)=footprint are alpha-composited (the reused compositor) -> a rendered image. Gate: 12// T1 FORWARD SPLAT -- each Gaussian's footprint PEAKS at its center and falls off (a real splat), 13// T2 DIFFERENTIABLE FIT -- GD on the Gaussian COLORS reconstructs a target render (loss->~0, colors recovered): 14// the explicit-primitive scene is differentiably optimizable (the 3DGS training loop), 15// T3 NEG-CONTROL -- ascending that gradient diverges (the gradient has teeth), 16// T4 BIT-EXACT -- train twice -> identical integer colors. 17// HONEST: 1D, fixed mu/sigma/opacity (only colors optimized); full 3DGS = 3D gaussians with optimized 18// position/covariance/opacity + tile sort + adaptive densification = the follow-on. Sovereign: 19// nx_nofloat_autograd (tape + nfa_fxexp) + syscalls, no float. license_tier: ORIGINAL expect_exit: 0 20import "nx_nofloat_autograd.nx" 21import "nx_syscalls.nx" 22const K_MAGIC_1024: i64 = 1024 23const K_MAGIC_16384: i64 = 16384 24const K_MAGIC_32768: i64 = 32768 25const K_MAGIC_49152: i64 = 49152 26const K_MAGIC_13107: i64 = 13107 27const K_MAGIC_45875: i64 = 45875 28const K_MAGIC_58982: i64 = 58982 29const K_MAGIC_40000: i64 = 40000 30const K_MAGIC_20000: i64 = 20000 31const K_MAGIC_4096: i64 = 4096 32 33const GLOG: *u8 = "knowledge/status/ng_gsplat.log" 34const Q16: i64 = 65536 35const NP: i64 = 12 // image pixels (1D) 36const NGA: i64 = 3 // gaussian primitives 37const INV2S2: i64 = 2275555 // 1/(2*sigma^2) in Q16 for sigma=0.12 (footprint width) 38const EPOCHS: i64 = 5000 39const LRQ: i64 = 4096 40 41func gp(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 42func gpn(v: i64) -> i64 { 43 let b: *u8 = sys_mmap(28); var x: i64=v 44 if x<0 { b[0]=45; sys_write(1,b,1); x=0-x } 45 if x==0 { b[0]=48; sys_write(1,b,1); return 0 } 46 var d: i64=0; var y: i64=x; while y>0 { d=d+1; y=y/10 } 47 var i: i64=d-1; y=x; while i>=0 { b[i]=(48+(y%10)) as u8; y=y/10; i=i-1 } 48 sys_write(1,b,d); return 0 49} 50func g_abs(v: i64) -> i64 { if v<0 { return 0-v } return v } 51func gqm(a: i64, b: i64) -> i64 { return (a * b) >> 16 } 52func gl_ws(fd: i64, s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(fd,s,n); return 0 } 53func gl_wn(fd: i64, v: i64) -> i64 { let b: *u8=sys_mmap(28); var m: i64=v; if m<0 {sys_write(fd,"-" as *u8,1); m=0-m} let t: *u8=sys_mmap(28); var k: i64=0; if m==0 {t[0]=48;k=1} while m>0 {t[k]=(48+(m%10)) as u8; m=m/10; k=k+1} var i: i64=0; while i<k {b[i]=t[k-1-i]; i=i+1} sys_write(fd,b,k); return 0 } 54 55// build the differentiable render: out_p = sum_g w[g][p]*c_g (w = the fixed splat+transmittance weights, c learnable) 56func gs_build(tape: *i64, vals: *i64, st: *i64, cp: *i64, w: *i64, targ: *i64, cn: *i64) -> i64 { 57 st[0]=0; st[1]=0 58 var g: i64=0 59 while g < NGA { cn[g] = nfa_leaf(tape,vals,st, 1, 1, cp, g); g = g + 1 } 60 var root: i64 = 0 - 1; var p: i64 = 0 61 while p < NP { 62 var acc: i64 = 0 - 1; var k: i64 = 0 63 while k < NGA { 64 let cm: i64 = nfa_cmul(tape,vals,st, cn[k], w[k*NP+p]) 65 if acc < 0 { acc = cm } else { acc = nfa_vadd(tape,vals,st, acc, cm) } 66 k = k + 1 67 } 68 let nt: i64 = nfa_leaf(tape,vals,st, 1, 1, targ, p) 69 let nm: i64 = nfa_mse(tape,vals,st, acc, nt) 70 if root < 0 { root = nm } else { root = nfa_vadd(tape,vals,st, root, nm) } 71 p = p + 1 72 } 73 return root 74} 75func gs_train(tape: *i64, vals: *i64, grads: *i64, st: *i64, w: *i64, targ: *i64, cout: *i64, lf: *i64, ll: *i64, sgn: i64) -> i64 { 76 let cp: *i64 = sys_mmap(NGA*8) as *i64 77 var i: i64=0; while i<NGA { cp[i]=0; i=i+1 } // learn colors from zero 78 let cn: *i64 = sys_mmap(NGA*8) as *i64; let g: *i64 = sys_mmap(NGA*8) as *i64 79 var ep: i64=0 80 while ep<EPOCHS { 81 let root: i64 = gs_build(tape,vals,st, cp,w,targ, cn) 82 nfa_backward(tape,vals,grads, st[0], root) 83 if ep==0 { *lf=nfa_val(tape,vals,root,0) } 84 *ll=nfa_val(tape,vals,root,0) 85 var k: i64=0; while k<NGA { g[k]=nfa_grad(tape,grads,cn[k],0); k=k+1 } 86 if sgn > 0 { nfa_sgd(cp, g, NGA, LRQ) } else { var j: i64=0; while j<NGA { cp[j]=cp[j]+gqm(LRQ,g[j]); j=j+1 } } // sgn<0 = ASCENT neg-control 87 ep=ep+1 88 } 89 i=0; while i<NGA { cout[i]=cp[i]; i=i+1 } 90 return 0 91} 92 93func main() -> i64 { 94 gp("nx_ng_gsplat: CAP-GAUSSIAN-SPLAT (3DGS) -- explicit gaussian primitives, splat + alpha-composite, Q16\n" as *u8) 95 let tape: *i64 = sys_mmap(K_MAGIC_1024*7*8) as *i64 96 let vals: *i64 = sys_mmap(K_MAGIC_16384*8) as *i64 97 let grads: *i64 = sys_mmap(K_MAGIC_16384*8) as *i64 98 let st: *i64 = sys_mmap(2*8) as *i64 99 100 let mu: *i64 = sys_mmap(NGA*8) as *i64 101 mu[0]=K_MAGIC_16384; mu[1]=K_MAGIC_32768; mu[2]=K_MAGIC_49152 // gaussian centers 0.25, 0.5, 0.75 102 let cstar: *i64 = sys_mmap(NGA*8) as *i64 103 cstar[0]=K_MAGIC_13107; cstar[1]=K_MAGIC_45875; cstar[2]=K_MAGIC_58982 // TRUE colors 0.2, 0.7, 0.9 104 105 // --- SPLAT: footprint[g][p] = exp(-(pos_p - mu_g)^2 / 2sigma^2) via fixed-point nfa_fxexp --- 106 let fp: *i64 = sys_mmap(NGA*NP*8) as *i64 107 var g: i64=0 108 while g < NGA { 109 var p: i64=0 110 while p < NP { 111 let pos: i64 = p * (Q16/NP) + (Q16/(2*NP)) // pixel center 112 let dlt: i64 = pos - mu[g] 113 let d2: i64 = gqm(dlt, dlt) 114 let arg: i64 = gqm(d2, INV2S2) 115 fp[g*NP+p] = nfa_fxexp(0 - arg) 116 p = p + 1 117 } 118 g = g + 1 119 } 120 // --- alpha-over composite weights (REUSE CAP-VOLUME-RENDER math) per pixel: w[g][p]=T*footprint; T*=(1-fp) --- 121 let w: *i64 = sys_mmap(NGA*NP*8) as *i64 122 var pp: i64=0 123 while pp < NP { 124 var T: i64 = Q16 125 var gg: i64=0 126 while gg < NGA { 127 let a: i64 = fp[gg*NP+pp] 128 w[gg*NP+pp] = gqm(T, a) 129 T = gqm(T, Q16 - a) 130 gg = gg + 1 131 } 132 pp = pp + 1 133 } 134 // target render from the TRUE colors 135 let targ: *i64 = sys_mmap(NP*8) as *i64 136 pp=0; while pp<NP { var s: i64=0; var k: i64=0; while k<NGA { s = s + gqm(w[k*NP+pp], cstar[k]); k=k+1 } targ[pp]=s; pp=pp+1 } 137 138 // T1 FORWARD: the middle gaussian's footprint peaks at its center pixel (p=6 ~ 0.5) and falls off at edges 139 var t1: i64=0 140 if fp[1*NP+6] > fp[1*NP+0] { if fp[1*NP+6] > K_MAGIC_40000 { if fp[1*NP+0] < K_MAGIC_20000 { t1=1 } } } 141 gp(" T1 SPLAT: mid-gaussian footprint center(p6)=" as *u8); gpn(fp[1*NP+6]); gp(" edge(p0)=" as *u8); gpn(fp[1*NP+0]); gp(" peaks-at-center=" as *u8); gpn(t1); gp("\n" as *u8) 142 143 // T2 DIFFERENTIABLE FIT: GD on colors reconstructs the target render 144 let cf: *i64 = sys_mmap(NGA*8) as *i64; let lf: *i64=sys_mmap(8) as *i64; let ll: *i64=sys_mmap(8) as *i64 145 gs_train(tape,vals,grads,st, w,targ, cf, lf,ll, 1) 146 var t2: i64=1 147 if (*ll)*20 >= (*lf) { t2=0 } 148 var gI: i64=0; while gI<NGA { if g_abs(cf[gI]-cstar[gI]) >= K_MAGIC_4096 { t2=0 } gI=gI+1 } 149 gp(" T2 FIT: SSE " as *u8); gpn(*lf); gp("->" as *u8); gpn(*ll); gp(" colors c=[" as *u8); gpn(cf[0]); gp("," as *u8); gpn(cf[1]); gp("," as *u8); gpn(cf[2]); gp("] vs c*=[13107,45875,58982] recovered=" as *u8); gpn(t2); gp("\n" as *u8) 150 151 // T3 NEG-CONTROL: ascend the gradient (bounded by EPOCHS already small lr) -> loss must grow 152 let cf2: *i64=sys_mmap(NGA*8) as *i64; let lf2: *i64=sys_mmap(8) as *i64; let ll2: *i64=sys_mmap(8) as *i64 153 gs_train(tape,vals,grads,st, w,targ, cf2, lf2,ll2, 0-1) 154 var t3: i64=0; if (*ll2) > (*lf2) { t3=1 } 155 gp(" T3 ASCENT(neg-control): SSE " as *u8); gpn(*lf2); gp("->" as *u8); gpn(*ll2); gp(" grew=" as *u8); gpn(t3); gp("\n" as *u8) 156 157 // T4 BIT-EXACT 158 let cf3: *i64=sys_mmap(NGA*8) as *i64; let lf3: *i64=sys_mmap(8) as *i64; let ll3: *i64=sys_mmap(8) as *i64 159 gs_train(tape,vals,grads,st, w,targ, cf3, lf3,ll3, 1) 160 var t4: i64=1; gI=0; while gI<NGA { if cf3[gI]!=cf[gI] { t4=0 } gI=gI+1 } 161 gp(" T4 bit-exact=" as *u8); gpn(t4); gp("\n" as *u8) 162 163 var ok: i64=1 164 if t1!=1 { ok=0 } 165 if t2!=1 { ok=0 } 166 if t3!=1 { ok=0 } 167 if t4!=1 { ok=0 } 168 let logf: i64 = sys_openat_append(GLOG, 420) 169 if logf>=0 { 170 gl_ws(logf,"NGGSPLAT authored=organ 3dgs-1d loss_first=" as *u8); gl_wn(logf,*lf); gl_ws(logf," loss_last=" as *u8); gl_wn(logf,*ll) 171 gl_ws(logf," t1=" as *u8); gl_wn(logf,t1); gl_ws(logf," t2=" as *u8); gl_wn(logf,t2); gl_ws(logf," t3=" as *u8); gl_wn(logf,t3); gl_ws(logf," t4=" as *u8); gl_wn(logf,t4) 172 if ok==1 { gl_ws(logf," verdict=GREEN\n" as *u8) } else { gl_ws(logf," verdict=RED\n" as *u8) } 173 sys_close(logf) 174 } 175 gp(" verdict=" as *u8) 176 if ok==1 { gp("GREEN (explicit gaussian primitives splat + alpha-composite + differentiable color optimization = 3DGS)\n" as *u8); sys_exit(0); return 0 } 177 gp("RED\n" as *u8) 178 sys_exit(1); return 1 179}