nx_gsplat.nx source

1// nx_gsplat.nx -- ★SOVEREIGN 3D GAUSSIAN SPLATTING (the frontier rendering rep for photoreal avatars, confirmed 2// #1 by the 2024-2025 SOTA scan: GaussianAvatars/VRGaussianAvatar/HumanSplat all use it). A scene = a cloud of 3// 3D Gaussians {position, scale(σ), colour, opacity}; render = PROJECT each to a 2D splat, DEPTH-SORT, and 4// FRONT-TO-BACK alpha-composite. ALL INTEGER (fx1024 positions, fx256 opacity/transmittance, an integer exp-LUT 5// for the Gaussian falloff), deterministic, VM-vettable. ★The RENDERER needs NO trained weights (only a 6// generator would) -- so it is fully sovereign-buildable TODAY. v0 = ISOTROPIC Gaussians (spherical σ -> a 7// circular screen splat); anisotropic 3D covariance (the Jacobian projection) + SH view-dependent colour are the 8// next rungs. Camera convention MATCHES nx_sdfrender (yaw orbit, FOCAL 586) so splats align with our SDF/mesh. 9// license_tier: ORIGINAL 10import "nx_syscalls.nx" 11import "nx_itrig.nx" 12import "nx_vecmath.nx" 13const K_MAGIC_4096: i64 = 4096 14const K_MAGIC_3970: i64 = 3970 15const K_MAGIC_2000000000: i64 = 2000000000 16const K_MAGIC_65536: i64 = 65536 17const K_MAGIC_999999: i64 = 999999 18 19const GW: i64 = 512 20const GH: i64 = 384 21const GHW: i64 = 256 22const GHH: i64 = 192 23const GFOCAL: i64 = 586 24const GFX: i64 = 1024 // model-unit fixed point (matches sdf) 25const GEXPN: i64 = 176 // exp-LUT entries: k=0..175 -> u=k/16 in [0,11) (11σ² = past 3σ) 26const GNB: i64 = 2048 // depth-sort buckets 27 28// ★THE FIXED-POINT AND LAYOUT FACTS THIS RENDERER DEFINES -- EXPORTED so no consumer re-declares them. 29// A consumer that hardcodes 256, 16 or 12 becomes a SECOND RULER: change the value here and the copy 30// silently disagrees while still compiling and still producing a plausible-looking frame. That failure has 31// no symptom, which is exactly why the numbers have to leave the consumers and live at their one source. 32const GFXA: i64 = 256 // opacity / normal / transmittance fixed point 33const GLUTU: i64 = 16 // exp-LUT sub-steps per unit of u: explut[k] = GFXA*exp(-k/(2*GLUTU)) 34const GSTA: i64 = 12 // stride of the ANISOTROPIC record gs_render_aniso consumes 35const GSTI: i64 = 8 // stride of the ISOTROPIC record gs_render consumes 36// the model-unit fixed point. Exported for the same reason GFXA is: a consumer that re-declares 1024 37// becomes a SECOND RULER for world scale, and the two disagree silently while still rendering a frame. 38func gs_fx() -> i64 { return GFX } 39func gs_fxa() -> i64 { return GFXA } 40func gs_lutu() -> i64 { return GLUTU } 41func gs_stride_aniso() -> i64 { return GSTA } 42func gs_stride_iso() -> i64 { return GSTI } 43 44// RETIRED ONTO THE SHARED OWNER 2026-08-24: this was a private copy of the same Newton floor-sqrt; vm_isqrt is 45// gate-proven exact over 20,000 inputs. The alias keeps every call site untouched. 46func gs_isqrt(v: i64) -> i64 { return vm_isqrt(v) } 47 48// exp-LUT: explut[k] = round(256 * exp(-0.5 * k/16)), built by an integer recurrence (factor exp(-0.5/16) 49// ~= 3970/4096 in fx4096, then scaled to fx256). No float, deterministic. 50func gs_build_explut(explut: *i64) -> i64 { 51 var v: i64 = K_MAGIC_4096 // fx4096, exp(0)=1 52 var k: i64 = 0 53 while k < GEXPN { 54 explut[k] = v * 256 / K_MAGIC_4096 // -> fx256 55 v = v * K_MAGIC_3970 / K_MAGIC_4096 // *= exp(-0.5/16) 56 k = k + 1 57 } 58 return 0 59} 60 61// clear the per-pixel accumulators: acc = 0 (fx256 colour sum), trans = 256 (fx256 transmittance = 1.0) 62func gs_clear_at(acc: *i64, trans: *i64, npx: i64) -> i64 { 63 var i: i64 = 0 64 while i < npx { acc[i*3] = 0; acc[i*3+1] = 0; acc[i*3+2] = 0; trans[i] = GFXA; i = i + 1 } 65 return 0 66} 67func gs_clear(acc: *i64, trans: *i64) -> i64 { return gs_clear_at(acc, trans, GW * GH) } 68 69// render the Gaussian cloud. gauss = ng*8 i64 {x,y,z,scale,r,g,b,opacity(fx256)}; camera (yaw it4096, camz units). 70// scratch: acc(GW*GH*3), trans(GW*GH), depth(ng), sx(ng), sy(ng), sig(ng), order(ng), count(GNB+1), explut(GEXPN). 71// writes fb (GW*GH packed rgb). Returns the number of splatted (visible) gaussians. 72func gs_render(gauss: *i64, ng: i64, yaw: i64, camz: i64, fb: *i64, acc: *i64, trans: *i64, depth: *i64, sxb: *i64, syb: *i64, sigb: *i64, order: *i64, count: *i64, explut: *i64, bgr: i64, bgg: i64, bgb: i64) -> i64 { 73 let sy4: i64 = it_sin4096(yaw) 74 let cy4: i64 = it_cos4096(yaw) 75 let R: i64 = camz * GFX 76 // --- project all gaussians; mark culled with depth = -1 --- 77 var dmin: i64 = K_MAGIC_2000000000 78 var dmax: i64 = 0 - K_MAGIC_2000000000 79 var vis: i64 = 0 80 var i: i64 = 0 81 while i < ng { 82 let x: i64 = gauss[i*8] 83 let y: i64 = gauss[i*8+1] 84 let z: i64 = gauss[i*8+2] 85 let sc: i64 = gauss[i*8+3] 86 let cz: i64 = (x * sy4 + z * cy4) / K_MAGIC_4096 + R // camera-forward depth 87 if cz > 64 { 88 let cx: i64 = (x * cy4 - z * sy4) / K_MAGIC_4096 89 let sx: i64 = GHW + GFOCAL * cx / cz 90 let sy: i64 = GHH - GFOCAL * y / cz 91 var sig: i64 = GFOCAL * sc / cz // screen-space sigma (px) 92 if sig < 1 { sig = 1 } 93 depth[i] = cz; sxb[i] = sx; syb[i] = sy; sigb[i] = sig 94 if cz < dmin { dmin = cz } 95 if cz > dmax { dmax = cz } 96 vis = vis + 1 97 } else { depth[i] = 0 - 1 } 98 i = i + 1 99 } 100 if vis == 0 { return 0 } 101 // --- counting sort visible gaussians by depth (near->far = ascending cz) --- 102 let span: i64 = dmax - dmin + 1 103 var b: i64 = 0 104 while b <= GNB { count[b] = 0; b = b + 1 } 105 i = 0 106 while i < ng { if depth[i] >= 0 { var bk: i64 = (depth[i] - dmin) * GNB / span; if bk < 0 { bk = 0 } if bk >= GNB { bk = GNB - 1 } count[bk] = count[bk] + 1 } i = i + 1 } 107 var acc2: i64 = 0 // prefix sum -> bucket start offsets 108 b = 0 109 while b < GNB { let c: i64 = count[b]; count[b] = acc2; acc2 = acc2 + c; b = b + 1 } 110 i = 0 111 while i < ng { if depth[i] >= 0 { var bk: i64 = (depth[i] - dmin) * GNB / span; if bk < 0 { bk = 0 } if bk >= GNB { bk = GNB - 1 } order[count[bk]] = i; count[bk] = count[bk] + 1 } i = i + 1 } 112 // --- clear + composite front-to-back --- 113 gs_clear(acc, trans) 114 var oi: i64 = 0 115 while oi < vis { 116 let g: i64 = order[oi] 117 let sx: i64 = sxb[g] 118 let sy: i64 = syb[g] 119 let sig: i64 = sigb[g] 120 let rad: i64 = 3 * sig 121 let s2: i64 = sig * sig 122 let gr: i64 = gauss[g*8+4] 123 let gg: i64 = gauss[g*8+5] 124 let gb: i64 = gauss[g*8+6] 125 let op: i64 = gauss[g*8+7] 126 var py: i64 = sy - rad 127 if py < 0 { py = 0 } 128 var pye: i64 = sy + rad 129 if pye >= GH { pye = GH - 1 } 130 while py <= pye { 131 let dy: i64 = py - sy 132 var px: i64 = sx - rad 133 if px < 0 { px = 0 } 134 var pxe: i64 = sx + rad 135 if pxe >= GW { pxe = GW - 1 } 136 while px <= pxe { 137 let dx: i64 = px - sx 138 let d2: i64 = dx * dx + dy * dy 139 var k: i64 = d2 * 16 / s2 // u=d2/sig2 in 1/16ths 140 if k < GEXPN { 141 let pix: i64 = py * GW + px 142 let tr: i64 = trans[pix] 143 if tr > 1 { 144 let alpha: i64 = op * explut[k] / 256 // fx256 splat opacity at this pixel 145 let contrib: i64 = tr * alpha / 256 // transmittance-weighted (fx256) 146 acc[pix*3] = acc[pix*3] + gr * contrib 147 acc[pix*3+1] = acc[pix*3+1] + gg * contrib 148 acc[pix*3+2] = acc[pix*3+2] + gb * contrib 149 trans[pix] = tr - tr * alpha / 256 // *= (1-alpha) 150 } 151 } 152 px = px + 1 153 } 154 py = py + 1 155 } 156 oi = oi + 1 157 } 158 // --- resolve: colour = accum/256 + background * leftover transmittance --- 159 i = 0 160 while i < GW * GH { 161 let tr: i64 = trans[i] 162 var r: i64 = acc[i*3] / 256 + bgr * tr / 256 163 var gg2: i64 = acc[i*3+1] / 256 + bgg * tr / 256 164 var bb: i64 = acc[i*3+2] / 256 + bgb * tr / 256 165 if r > 255 { r = 255 } 166 if gg2 > 255 { gg2 = 255 } 167 if bb > 255 { bb = 255 } 168 fb[i] = r + gg2 * 256 + bb * K_MAGIC_65536 169 i = i + 1 170 } 171 return vis 172} 173 174// set one isotropic gaussian 175func gs_set(gauss: *i64, i: i64, x: i64, y: i64, z: i64, sc: i64, r: i64, g: i64, b: i64, op: i64) -> i64 { 176 gauss[i*8]=x; gauss[i*8+1]=y; gauss[i*8+2]=z; gauss[i*8+3]=sc; gauss[i*8+4]=r; gauss[i*8+5]=g; gauss[i*8+6]=b; gauss[i*8+7]=op 177 return 0 178} 179 180// ===== ★ANISOTROPIC SURFACE SPLATTING (the defining 3DGS feature: oriented ELLIPSE splats, not blobs) ===== 181// each gaussian = a DISK in the tangent plane (a "surfel") with normal n + in-plane radius rtan -> projects to a 182// 2D ELLIPSE (EWA-style). This aligns splats to the surface -> smooth coverage between mesh verts (kills the 183// isotropic-blob grid/dotty look). 12 i64/gaussian: {x,y,z, nx,ny,nz(fx256 unit), rtan, r,g,b, op, _spare}. 184func gs_set_aniso(gauss: *i64, i: i64, x: i64, y: i64, z: i64, nx: i64, ny: i64, nz: i64, rtan: i64, r: i64, g: i64, b: i64, op: i64) -> i64 { 185 gauss[i*12]=x; gauss[i*12+1]=y; gauss[i*12+2]=z; gauss[i*12+3]=nx; gauss[i*12+4]=ny; gauss[i*12+5]=nz 186 gauss[i*12+6]=rtan; gauss[i*12+7]=r; gauss[i*12+8]=g; gauss[i*12+9]=b; gauss[i*12+10]=op; gauss[i*12+11]=0 187 return 0 188} 189// ★★RESOLUTION IS A PARAMETER, NOT A COMPILE-TIME FACT (widened 2026-08-23). 190// Until today every stage read the module consts GW/GH/GHW/GHH/GFOCAL, so this rasterizer -- AND 191// EVERY "second renderer" the three-stage decomposition above exists to enable (nx_gsplat_tile_lib 192// is one, measured: it calls gs_project_sort/gs_blend_rect/gs_resolve) -- could only ever produce 193// 512x384 = 196,608 pixels. A beauty-tier frame is 1920x1080 (10.6x those pixels) or 3840x2160 194// (42x). THE CEILING WAS NEVER A GRAPHICS LIMIT; IT WAS FIVE CONSTANTS READ FROM MODULE SCOPE. 195// The inner stages were ALREADY viewport-agnostic -- gs_blend_rect takes stride_w, gs_resolve takes 196// npx -- so only the projection and the entry points were nailed down. The _at variants below take 197// the viewport explicitly; the historic names are now thin wrappers passing the historic values, so 198// every existing caller renders BYTE-IDENTICALLY and the widening is neutral BY CONSTRUCTION. 199// 200// ★FOCAL IS DERIVED, NEVER PICKED. gs_focal_for(h) scales the incumbent focal by h/GH, which holds 201// tan(vfov/2) = (h/2)/focal EXACTLY equal to the historic (GH/2)/GFOCAL -- the SAME field of view at 202// any height. Picking a focal instead would silently change the LENS, which is a different decision 203// and has to be made on purpose. Half-extents are DERIVED from the viewport (hw=w/2, hh=h/2) at the 204// entry point rather than passed alongside it, so a caller cannot hand in a centre that disagrees 205// with its own width -- a second ruler this file already refuses everywhere else. 206func gs_focal_for(h: i64) -> i64 { return GFOCAL * h / GH } 207func gs_proj_x_at(x: i64, y: i64, z: i64, sy4: i64, cy4: i64, R: i64, hw: i64, focal: i64) -> i64 { 208 let cz: i64 = (x*sy4 + z*cy4)/K_MAGIC_4096 + R 209 if cz <= 64 { return 0 - K_MAGIC_999999 } 210 let cx: i64 = (x*cy4 - z*sy4)/K_MAGIC_4096 211 return hw + focal*cx/cz 212} 213func gs_proj_y_at(x: i64, y: i64, z: i64, sy4: i64, cy4: i64, R: i64, hh: i64, focal: i64) -> i64 { 214 let cz: i64 = (x*sy4 + z*cy4)/K_MAGIC_4096 + R 215 if cz <= 64 { return 0 - K_MAGIC_999999 } 216 return hh - focal*y/cz 217} 218func gs_proj_x(x: i64, y: i64, z: i64, sy4: i64, cy4: i64, R: i64) -> i64 { 219 return gs_proj_x_at(x, y, z, sy4, cy4, R, GHW, GFOCAL) 220} 221func gs_proj_y(x: i64, y: i64, z: i64, sy4: i64, cy4: i64, R: i64) -> i64 { 222 return gs_proj_y_at(x, y, z, sy4, cy4, R, GHH, GFOCAL) 223} 224// ★★THE RASTERIZER, DECOMPOSED INTO ITS THREE STAGES so a second renderer (tiled, foveated, stereo) can 225// be built WITHOUT retyping any of them. A retyped projection or blend would drift silently: two renderers 226// would disagree by rounding nobody could see in a frame. These are the one implementation of each stage. 227const GNEARZ: i64 = 64 // near-clip in camera depth units: nearer than this is culled 228const GUPTHRESH: i64 = 210 // |n.y| (fx256) above which the up-vector must be swapped to avoid 229 // a degenerate cross product when the normal IS the up axis 230const GSIGMA_CUT: i64 = 3 // AABB half-extent in sigma; the LUT carries weight to sqrt(GEXPN/GLUTU) 231 232// STAGE 1: project every gaussian to a screen position + 2D conic, and depth-sort the visible ones. 233// depth[i] = -1 marks culled. Returns the visible count; order[0..vis) is near-to-far. 234func gs_project_sort_at(gauss: *i64, ng: i64, yaw: i64, camz: i64, depth: *i64, sxb: *i64, syb: *i64, pa: *i64, pb: *i64, pc: *i64, pdet: *i64, order: *i64, count: *i64, hw: i64, hh: i64, focal: i64) -> i64 { 235 let sy4: i64 = it_sin4096(yaw) 236 let cy4: i64 = it_cos4096(yaw) 237 let R: i64 = camz * GFX 238 var dmin: i64 = K_MAGIC_2000000000 239 var dmax: i64 = 0 - K_MAGIC_2000000000 240 var vis: i64 = 0 241 var i: i64 = 0 242 while i < ng { 243 let x: i64 = gauss[i*GSTA] 244 let y: i64 = gauss[i*GSTA+1] 245 let z: i64 = gauss[i*GSTA+2] 246 let cz: i64 = (x*sy4 + z*cy4)/K_MAGIC_4096 + R 247 if cz > GNEARZ { 248 let nx: i64 = gauss[i*GSTA+3] 249 let ny: i64 = gauss[i*GSTA+4] 250 let nz: i64 = gauss[i*GSTA+5] 251 let rt: i64 = gauss[i*GSTA+6] 252 var ux: i64 = 0 253 var uy: i64 = GFXA 254 var uz: i64 = 0 255 if ny > GUPTHRESH { ux = GFXA; uy = 0; uz = 0 } 256 if ny < 0 - GUPTHRESH { ux = GFXA; uy = 0; uz = 0 } 257 var t1x: i64 = (ny*uz - nz*uy)/GFXA 258 var t1y: i64 = (nz*ux - nx*uz)/GFXA 259 var t1z: i64 = (nx*uy - ny*ux)/GFXA 260 let l1: i64 = gs_isqrt(t1x*t1x + t1y*t1y + t1z*t1z) 261 if l1 > 0 { t1x = t1x*GFXA/l1; t1y = t1y*GFXA/l1; t1z = t1z*GFXA/l1 } 262 var t2x: i64 = (ny*t1z - nz*t1y)/GFXA 263 var t2y: i64 = (nz*t1x - nx*t1z)/GFXA 264 var t2z: i64 = (nx*t1y - ny*t1x)/GFXA 265 let l2: i64 = gs_isqrt(t2x*t2x + t2y*t2y + t2z*t2z) 266 if l2 > 0 { t2x = t2x*GFXA/l2; t2y = t2y*GFXA/l2; t2z = t2z*GFXA/l2 } 267 let a1x: i64 = t1x*rt/GFXA 268 let a1y: i64 = t1y*rt/GFXA 269 let a1z: i64 = t1z*rt/GFXA 270 let a2x: i64 = t2x*rt/GFXA 271 let a2y: i64 = t2y*rt/GFXA 272 let a2z: i64 = t2z*rt/GFXA 273 let sx: i64 = gs_proj_x_at(x, y, z, sy4, cy4, R, hw, focal) 274 let sy: i64 = gs_proj_y_at(x, y, z, sy4, cy4, R, hh, focal) 275 let e1x: i64 = gs_proj_x_at(x+a1x, y+a1y, z+a1z, sy4, cy4, R, hw, focal) - sx 276 let e1y: i64 = gs_proj_y_at(x+a1x, y+a1y, z+a1z, sy4, cy4, R, hh, focal) - sy 277 let e2x: i64 = gs_proj_x_at(x+a2x, y+a2y, z+a2z, sy4, cy4, R, hw, focal) - sx 278 let e2y: i64 = gs_proj_y_at(x+a2x, y+a2y, z+a2z, sy4, cy4, R, hh, focal) - sy 279 var ca: i64 = e1x*e1x + e2x*e2x + 1 280 let cb: i64 = e1x*e1y + e2x*e2y 281 var cc: i64 = e1y*e1y + e2y*e2y + 1 282 var det: i64 = ca*cc - cb*cb 283 if det < 1 { det = 1 } 284 depth[i] = cz; sxb[i] = sx; syb[i] = sy; pa[i] = ca; pb[i] = cb; pc[i] = cc; pdet[i] = det 285 if cz < dmin { dmin = cz } 286 if cz > dmax { dmax = cz } 287 vis = vis + 1 288 } else { depth[i] = 0 - 1 } 289 i = i + 1 290 } 291 if vis == 0 { return 0 } 292 let span: i64 = dmax - dmin + 1 293 var b: i64 = 0 294 while b <= GNB { count[b] = 0; b = b + 1 } 295 i = 0 296 while i < ng { if depth[i] >= 0 { var bk: i64 = (depth[i]-dmin)*GNB/span; if bk<0{bk=0} if bk>=GNB{bk=GNB-1} count[bk]=count[bk]+1 } i=i+1 } 297 var pre: i64 = 0 298 b = 0 299 while b < GNB { let c: i64 = count[b]; count[b] = pre; pre = pre + c; b = b + 1 } 300 i = 0 301 while i < ng { if depth[i] >= 0 { var bk: i64 = (depth[i]-dmin)*GNB/span; if bk<0{bk=0} if bk>=GNB{bk=GNB-1} order[count[bk]]=i; count[bk]=count[bk]+1 } i=i+1 } 302 return vis 303} 304// The historic entry point: the viewport IS the module default. Both call sites of this stage keep 305// their signature and their exact output (this file's gs_render_aniso and nx_gsplat_tile_lib's tiled 306// walk -- the complete population, corpus_complete=1 over 59,540 files); only a caller that WANTS a 307// bigger frame passes one. One implementation, two entry widths -- never a second sort. 308func gs_project_sort(gauss: *i64, ng: i64, yaw: i64, camz: i64, depth: *i64, sxb: *i64, syb: *i64, pa: *i64, pb: *i64, pc: *i64, pdet: *i64, order: *i64, count: *i64) -> i64 { 309 return gs_project_sort_at(gauss, ng, yaw, camz, depth, sxb, syb, pa, pb, pc, pdet, order, count, GHW, GHH, GFOCAL) 310} 311 312// the screen-space AABB half-extents of a splat's conic. Derived from the conic, not chosen. 313func gs_splat_rx(ca: i64) -> i64 { return GSIGMA_CUT * gs_isqrt(ca) } 314func gs_splat_ry(cc: i64) -> i64 { return GSIGMA_CUT * gs_isqrt(cc) } 315 316// STAGE 2: blend one splat into an ALREADY-CLIPPED pixel rectangle, front-to-back. 317// ★Returns the number of CONIC EVALUATIONS performed -- the cost unit. Any renderer built on this reports 318// its work in the same currency, so two renderers can be compared without a stopwatch or a guess. 319// satcnt[0] is INCREMENTED once for each pixel that transitions from live to saturated inside this call. 320// A tiled renderer uses it to stop feeding a region that can no longer change; the full-screen policy 321// simply ignores the number. Counting it HERE, at the only place a pixel can saturate, is what makes the 322// two policies share one blend instead of one growing a private copy with an early-out bolted on. 323func gs_blend_rect(gauss: *i64, g: i64, sx: i64, sy: i64, ca: i64, cb: i64, cc: i64, det: i64, x0: i64, x1: i64, y0: i64, y1: i64, acc: *i64, trans: *i64, explut: *i64, stride_w: i64, satcnt: *i64) -> i64 { 324 let gr: i64 = gauss[g*GSTA+7] 325 let gg: i64 = gauss[g*GSTA+8] 326 let gb: i64 = gauss[g*GSTA+9] 327 let op: i64 = gauss[g*GSTA+10] 328 var evals: i64 = 0 329 var py: i64 = y0 330 while py <= y1 { 331 let dy: i64 = py - sy 332 var px: i64 = x0 333 while px <= x1 { 334 let dx: i64 = px - sx 335 // power = d^T Sigma2d^-1 d = (cc dx^2 - 2 cb dx dy + ca dy^2)/det ; LUT idx k = power*GLUTU 336 let pnum: i64 = cc*dx*dx - 2*cb*dx*dy + ca*dy*dy 337 let k: i64 = pnum * GLUTU / det 338 evals = evals + 1 339 if k < GEXPN { if k >= 0 { 340 let pix: i64 = py*stride_w + px 341 let tr: i64 = trans[pix] 342 if tr > 1 { 343 let alpha: i64 = op * explut[k] / GFXA 344 let contrib: i64 = tr * alpha / GFXA 345 acc[pix*3] = acc[pix*3] + gr*contrib 346 acc[pix*3+1] = acc[pix*3+1] + gg*contrib 347 acc[pix*3+2] = acc[pix*3+2] + gb*contrib 348 let ntr: i64 = tr - tr*alpha/GFXA 349 trans[pix] = ntr 350 if ntr <= 1 { satcnt[0] = satcnt[0] + 1 } 351 } 352 } } 353 px = px + 1 354 } 355 py = py + 1 356 } 357 return evals 358} 359 360// STAGE 3: resolve accumulated colour + leftover transmittance against the background. 361func gs_resolve(fb: *i64, acc: *i64, trans: *i64, npx: i64, bgr: i64, bgg: i64, bgb: i64) -> i64 { 362 var i: i64 = 0 363 while i < npx { 364 let tr: i64 = trans[i] 365 var r: i64 = acc[i*3]/GFXA + bgr*tr/GFXA 366 var g2: i64 = acc[i*3+1]/GFXA + bgg*tr/GFXA 367 var bb: i64 = acc[i*3+2]/GFXA + bgb*tr/GFXA 368 if r>255{r=255} if g2>255{g2=255} if bb>255{bb=255} 369 fb[i] = r + g2*256 + bb*K_MAGIC_65536 370 i = i + 1 371 } 372 return 0 373} 374 375// anisotropic render: pa/pb/pc/pdet = per-gaussian 2D-covariance scratch (ng each). Same sort + front-to-back 376// composite as the isotropic path; only the per-splat kernel is an oriented ellipse. 377// ★NOW COMPOSED FROM THE THREE STAGES ABOVE -- this function is the FULL-SCREEN policy over them, and the 378// tiled renderer is a different policy over the SAME stages. Neither can drift from the other. 379func gs_render_aniso_at(gauss: *i64, ng: i64, yaw: i64, camz: i64, fb: *i64, acc: *i64, trans: *i64, depth: *i64, sxb: *i64, syb: *i64, pa: *i64, pb: *i64, pc: *i64, pdet: *i64, order: *i64, count: *i64, explut: *i64, bgr: i64, bgg: i64, bgb: i64, w: i64, h: i64, focal: i64) -> i64 { 380 // hw/hh are DERIVED from the viewport, never passed beside it: a centre that disagrees with its 381 // own width is a second ruler, and this file refuses those everywhere else. 382 let hw: i64 = w / 2 383 let hh: i64 = h / 2 384 let vis: i64 = gs_project_sort_at(gauss, ng, yaw, camz, depth, sxb, syb, pa, pb, pc, pdet, order, count, hw, hh, focal) 385 if vis == 0 { return 0 } 386 gs_clear_at(acc, trans, w * h) 387 var oi: i64 = 0 388 while oi < vis { 389 let g: i64 = order[oi] 390 let sx: i64 = sxb[g] 391 let sy: i64 = syb[g] 392 let ca: i64 = pa[g] 393 let cb: i64 = pb[g] 394 let cc: i64 = pc[g] 395 let det: i64 = pdet[g] 396 let rx: i64 = gs_splat_rx(ca) 397 let ry: i64 = gs_splat_ry(cc) 398 var y0: i64 = sy - ry 399 if y0 < 0 { y0 = 0 } 400 var y1: i64 = sy + ry 401 if y1 >= h { y1 = h - 1 } 402 var x0: i64 = sx - rx 403 if x0 < 0 { x0 = 0 } 404 var x1: i64 = sx + rx 405 if x1 >= w { x1 = w - 1 } 406 // the full-screen policy has no early-out, so it discards the saturation count into the slot the 407 // sort scratch already reserves past its own range (callers allocate GNB+2). No extra allocation 408 // per frame, which a long-running renderer would otherwise pay a page for. 409 gs_blend_rect(gauss, g, sx, sy, ca, cb, cc, det, x0, x1, y0, y1, acc, trans, explut, w, ((count as i64) + (GNB+1)*8) as *i64) 410 oi = oi + 1 411 } 412 gs_resolve(fb, acc, trans, w*h, bgr, bgg, bgb) 413 return vis 414} 415// The historic full-screen entry point, unchanged for every existing caller: the module viewport, the 416// module focal, byte-identical output. A caller wanting a beauty-tier frame calls gs_render_aniso_at 417// with its own w/h and gs_focal_for(h) -- same lens, more pixels. SCRATCH IS THE CALLER'S CONTRACT: 418// acc needs w*h*3 and trans w*h i64, so a wider viewport must allocate wider buffers; the renderer 419// cannot check that for you, and the sizes are stated here rather than discovered by a fault. 420func gs_render_aniso(gauss: *i64, ng: i64, yaw: i64, camz: i64, fb: *i64, acc: *i64, trans: *i64, depth: *i64, sxb: *i64, syb: *i64, pa: *i64, pb: *i64, pc: *i64, pdet: *i64, order: *i64, count: *i64, explut: *i64, bgr: i64, bgg: i64, bgb: i64) -> i64 { 421 return gs_render_aniso_at(gauss, ng, yaw, camz, fb, acc, trans, depth, sxb, syb, pa, pb, pc, pdet, order, count, explut, bgr, bgg, bgb, GW, GH, GFOCAL) 422} 423// ===== ★GAUSSIAN-PARAMETER OPTIMIZER (the generator/"training" side = the real photoreal lift) ===== 424// v0 = DIFFERENTIABLE COLOUR fit: a splatted pixel is LINEAR in the Gaussian colours (pixel = Σ w_i·c_i), so the 425// image-loss gradient wrt each colour is ANALYTIC (dL/dc_i = Σ_px (render-target)·w_i). One gradient-descent 426// step: forward-render, accumulate the per-Gaussian colour gradient over its splat footprint, take a 427// least-squares-preconditioned step (÷ Σw²). Integer, deterministic. This is the 3DGS training mechanism (colour 428// channel); position/covariance optimization + adaptive density = the next rungs. Returns the L1 image loss. 429func gs_color_grad_step(gauss: *i64, ng: i64, target: *i64, yaw: i64, camz: i64, fb: *i64, acc: *i64, trans: *i64, depth: *i64, sxb: *i64, syb: *i64, sigb: *i64, order: *i64, count: *i64, explut: *i64, gradbuf: *i64, wnorm: *i64) -> i64 { 430 gs_render(gauss, ng, yaw, camz, fb, acc, trans, depth, sxb, syb, sigb, order, count, explut, 26, 28, 44) 431 // L1 loss + zero the grad accumulators 432 var loss: i64 = 0 433 var i: i64 = 0 434 while i < GW*GH { 435 var dr: i64 = (fb[i]&255) - (target[i]&255); if dr<0 {dr=0-dr} 436 var dg: i64 = ((fb[i]>>8)&255) - ((target[i]>>8)&255); if dg<0 {dg=0-dg} 437 var db: i64 = ((fb[i]>>16)&255) - ((target[i]>>16)&255); if db<0 {db=0-db} 438 loss = loss + dr + dg + db 439 i = i + 1 440 } 441 i = 0 442 while i < ng { gradbuf[i*3]=0; gradbuf[i*3+1]=0; gradbuf[i*3+2]=0; wnorm[i]=0; i=i+1 } 443 // accumulate per-gaussian colour gradient over its splat footprint (uses the projections gs_render just filled) 444 i = 0 445 while i < ng { 446 if depth[i] >= 0 { 447 let sx: i64 = sxb[i] 448 let sy: i64 = syb[i] 449 let sig: i64 = sigb[i] 450 let rad: i64 = 3*sig 451 let s2: i64 = sig*sig 452 let op: i64 = gauss[i*8+7] 453 var py: i64 = sy-rad 454 if py<0 {py=0} 455 var pye: i64 = sy+rad 456 if pye>=GH {pye=GH-1} 457 while py <= pye { 458 let dy: i64 = py-sy 459 var px: i64 = sx-rad 460 if px<0 {px=0} 461 var pxe: i64 = sx+rad 462 if pxe>=GW {pxe=GW-1} 463 while px <= pxe { 464 let dx: i64 = px-sx 465 let d2: i64 = dx*dx+dy*dy 466 let k: i64 = d2*16/s2 467 if k < GEXPN { 468 let pix: i64 = py*GW+px 469 let w: i64 = op*explut[k]/256 // this gaussian's splat weight (fx256) 470 gradbuf[i*3] = gradbuf[i*3] + ((fb[pix]&255) - (target[pix]&255)) * w 471 gradbuf[i*3+1] = gradbuf[i*3+1] + (((fb[pix]>>8)&255) - ((target[pix]>>8)&255)) * w 472 gradbuf[i*3+2] = gradbuf[i*3+2] + (((fb[pix]>>16)&255) - ((target[pix]>>16)&255)) * w 473 wnorm[i] = wnorm[i] + w*w/256 474 } 475 px = px + 1 476 } 477 py = py + 1 478 } 479 // least-squares-preconditioned step (lr 0.75): c -= 3*grad/(4*wnorm) 480 let nd: i64 = wnorm[i]*4 + 1 481 var nr: i64 = gauss[i*8+4] - gradbuf[i*3]*3/nd 482 var nge: i64 = gauss[i*8+5] - gradbuf[i*3+1]*3/nd 483 var nb: i64 = gauss[i*8+6] - gradbuf[i*3+2]*3/nd 484 if nr<0 {nr=0} if nr>255 {nr=255} 485 if nge<0 {nge=0} if nge>255 {nge=255} 486 if nb<0 {nb=0} if nb>255 {nb=255} 487 gauss[i*8+4]=nr; gauss[i*8+5]=nge; gauss[i*8+6]=nb 488 } 489 i = i + 1 490 } 491 return loss 492} 493 494// ===== ★POSITION GRADIENT: the leg that turns a RENDERER into a RECONSTRUCTOR ===== 495// The colour step above is analytic because a pixel is LINEAR in colour. Position is NOT: moving a 496// gaussian changes WHICH pixels it covers and by HOW MUCH. But the Gaussian carries its own derivative 497// in closed form, so this needs neither finite differences nor a second render: 498// w = op * exp(-d^2 / (2 sigma^2)) => dw/dsx = w * dx / sigma^2, dx = px - sx 499// A pixel's contribution from gaussian i is c * w/GFXA, so per channel 500// J_ch = c_ch * (w/GFXA) * dx / sigma^2 501// and the Gauss-Newton step is the ratio of two sums the same pixel walk already visits: 502// dL/dsx = SUM resid_ch * J_ch curv = SUM J_ch^2 dsx = -dL/dsx / curv 503// ★THE STEP SIZE IS DERIVED, NOT TUNED. Gauss-Newton fixes the magnitude; the only free choice is 504// damping, and it is set to 1/2 because the model is NONLINEAR in position (unlike colour, where the 505// incumbent can afford 3/4 precisely because the model is linear there). The screen correction is 506// returned to world units through the projection Jacobian INVERTED, dx_world = dsx * cz/GFOCAL -- 507// world-units-per-pixel AT THAT GAUSSIAN'S DEPTH, read from the camera rather than picked. 508// ⚠DECLARED APPROXIMATION, the same one the colour step makes and names: the gradient ignores 509// transmittance, i.e. it treats each gaussian's footprint independently of what occludes it. That is 510// exact for a non-overlapping scene and an approximation elsewhere; it is stated here rather than 511// discovered later, and it is why the fitting gate's fixture is separated in depth. 512func gs_pos_grad_step(gauss: *i64, ng: i64, target: *i64, yaw: i64, camz: i64, fb: *i64, acc: *i64, trans: *i64, depth: *i64, sxb: *i64, syb: *i64, sigb: *i64, order: *i64, count: *i64, explut: *i64, gradbuf: *i64, wnorm: *i64) -> i64 { 513 gs_render(gauss, ng, yaw, camz, fb, acc, trans, depth, sxb, syb, sigb, order, count, explut, 26, 28, 44) 514 var loss: i64 = 0 515 var i: i64 = 0 516 while i < GW*GH { 517 var dr: i64 = (fb[i]&255) - (target[i]&255); if dr<0 {dr=0-dr} 518 var dg: i64 = ((fb[i]>>8)&255) - ((target[i]>>8)&255); if dg<0 {dg=0-dg} 519 var db: i64 = ((fb[i]>>16)&255) - ((target[i]>>16)&255); if db<0 {db=0-db} 520 loss = loss + dr + dg + db 521 i = i + 1 522 } 523 i = 0 524 while i < ng { gradbuf[i*3]=0; gradbuf[i*3+1]=0; gradbuf[i*3+2]=0; wnorm[i]=0; i=i+1 } 525 let sy4: i64 = it_sin4096(yaw) 526 let cy4: i64 = it_cos4096(yaw) 527 let R: i64 = camz * GFX 528 i = 0 529 while i < ng { 530 if depth[i] >= 0 { 531 let sx: i64 = sxb[i] 532 let sy: i64 = syb[i] 533 let sig: i64 = sigb[i] 534 let rad: i64 = 3*sig 535 let s2: i64 = sig*sig 536 let op: i64 = gauss[i*8+7] 537 let cr: i64 = gauss[i*8+4] 538 let cg: i64 = gauss[i*8+5] 539 let cb: i64 = gauss[i*8+6] 540 let cm2: i64 = cr*cr + cg*cg + cb*cb 541 var py: i64 = sy-rad 542 if py<0 {py=0} 543 var pye: i64 = sy+rad 544 if pye>=GH {pye=GH-1} 545 while py <= pye { 546 let dy: i64 = py-sy 547 var px: i64 = sx-rad 548 if px<0 {px=0} 549 var pxe: i64 = sx+rad 550 if pxe>=GW {pxe=GW-1} 551 while px <= pxe { 552 let dx: i64 = px-sx 553 let d2: i64 = dx*dx+dy*dy 554 let k: i64 = d2*GLUTU/s2 555 if k < GEXPN { 556 let pix: i64 = py*GW+px 557 let w: i64 = op*explut[k]/GFXA 558 let rc: i64 = ((fb[pix]&255) - (target[pix]&255))*cr + (((fb[pix]>>8)&255) - ((target[pix]>>8)&255))*cg + (((fb[pix]>>16)&255) - ((target[pix]>>16)&255))*cb 559 let ufx: i64 = w*dx*GFXA/s2 560 let ufy: i64 = w*dy*GFXA/s2 561 gradbuf[i*3] = gradbuf[i*3] + rc*ufx 562 gradbuf[i*3+1] = gradbuf[i*3+1] + rc*ufy 563 gradbuf[i*3+2] = gradbuf[i*3+2] + cm2*ufx*ufx/K_MAGIC_65536 564 wnorm[i] = wnorm[i] + cm2*ufy*ufy/K_MAGIC_65536 565 } 566 px = px + 1 567 } 568 py = py + 1 569 } 570 // world step through the INVERTED projection Jacobian, damped by 2 (see header). 571 let cz: i64 = (gauss[i*8]*sy4 + gauss[i*8+2]*cy4)/K_MAGIC_4096 + R 572 if cz > GNEARZ { 573 let cxd: i64 = gradbuf[i*3+2]*2*GFOCAL + 1 574 let cyd: i64 = wnorm[i]*2*GFOCAL + 1 575 gauss[i*8] = gauss[i*8] - gradbuf[i*3]*cz/cxd 576 gauss[i*8+1] = gauss[i*8+1] + gradbuf[i*3+1]*cz/cyd 577 } 578 } 579 i = i + 1 580 } 581 return loss 582} 583 584// ===== ★ANISOTROPIC POSITION GRADIENT: fitting on the path that carries QUALITY ===== 585// gs_pos_grad_step fits the ISOTROPIC path (stride GSTI). That is the blob path. Real surfels live on 586// the ANISOTROPIC path (stride GSTA), and until now it could be RENDERED but never FITTED -- a 587// reconstructor that stops exactly where the representation starts to matter. 588// ★THE DERIVATIVE IS THE RENDERER'S OWN ARITHMETIC, DIFFERENTIATED -- not a re-derivation that could 589// drift from it. gs_blend_rect evaluates the conic as 590// pnum = cc*dx^2 - 2*cb*dx*dy + ca*dy^2, u = pnum/det, w = op*exp(-u/2) 591// so, with dx = px - sx (hence d(dx)/dsx = -1): 592// du/dsx = -2*(cc*dx - cb*dy)/det => dw/dsx = w*(cc*dx - cb*dy)/det 593// du/dsy = -2*(ca*dy - cb*dx)/det => dw/dsy = w*(ca*dy - cb*dx)/det 594// ★AND IT REDUCES TO THE ISOTROPIC CASE EXACTLY, which is the check that it is the same mathematics and 595// not a second one: with ca = cc = sigma^2 and cb = 0, det = sigma^4, dw/dsx = w*dx/sigma^2 -- the 596// isotropic formula above, recovered term for term. 597// ⚠DECLARED PRECISION: the per-pixel derivative is accumulated as an integer of order 10..100, so it 598// carries ~1 percent truncation. That is a magnitude error in a DIRECTION, and Gauss-Newton re-measures 599// the direction every step, so it costs iterations and not correctness. Stated rather than discovered. 600// ⚠DECLARED APPROXIMATION, INHERITED AND UNCHANGED: like both incumbent steps, this ignores 601// transmittance -- exact for a non-overlapping scene, an approximation under occlusion. It is why the 602// fitting fixtures separate gaussians in depth, and it is the thing to revisit FIRST when adaptive 603// density starts producing overlapping splats. 604// ★★TRANSMITTANCE IS NOW A PARAMETER, SO ITS COST IS MEASURABLE INSTEAD OF ASSERTED (2026-08-23). 605// Both incumbent gradient steps DROP transmittance: they weight every splat's footprint as if nothing 606// occluded it. That is exact for a non-overlapping scene and an approximation elsewhere -- and 607// ADAPTIVE DENSITY PRODUCES OVERLAP BY CONSTRUCTION, so densifying against the approximation would 608// tune the fit through a ruler its own output invalidates (the same shape as judging beauty at 512x384). 609// So the approximation and the exact direct form live in ONE implementation, selected by usetrans: 610// usetrans=0 -> T is treated as full (GFXA). Byte-identical to the incumbent behaviour. 611// usetrans=1 -> T is the REAL front-to-back transmittance at that pixel before this splat. 612// usetrans=2 -> direct + INDIRECT. ⚠⚠QUARANTINED, DO NOT SHIP. Measured 2026-08-23: it is not 613// merely non-beneficial, it is NUMERICALLY UNSTABLE. Single-view understated it (100 614// against 79, which reads as a mild loss); under MULTI-VIEW the same code diverged to 615// 46,128,546,235 against a 52-pixel floor -- ten orders of magnitude, ending far 616// worse than it began. The mathematics in the branch is believed right (it reduces to 617// the direct form and its units were checked twice); the instability is unfound, and 618// the likely suspects are the S_i subtraction losing significance when acc_final and 619// the prefix are close, and the 1/(GFXA-alpha) divisor near full opacity where the 620// clamp at 1 is doing real work. LEFT IN, NOT DELETED, because the arm is the only 621// way to re-measure it -- but nothing may depend on it until it is bounded. 622// A splat's contribution to a pixel is c*T*alpha/GFXA^2, so the exact direct derivative just carries the 623// T/GFXA factor: J_ch = c_ch * T * vs / GFXA^2. With T=GFXA that is c_ch*vs/GFXA -- the incumbent form, 624// recovered exactly, which is why usetrans=0 is a true control and not an approximation of the control. 625// ⚠THE INDIRECT TERM IS STILL ABSENT AND IS NAMED HERE: moving splat i also changes how much light 626// reaches every splat BEHIND it (dT_j/dparam_i for j after i). That needs a back-to-front suffix pass 627// and is the next rung; what lands here is the direct term, which is the dominant one and the one that 628// makes an occluded splat stop pulling as hard as an unoccluded one. 629// ⚠AND THE WALK ORDER IS NOW LOAD-BEARING: transmittance only means anything front-to-back, so the 630// accumulation walks order[] (the renderer's own depth sort) rather than gaussian index. 631// ⚠⚠ARITY CEILING, HIT AND HONOURED. The compiler REFUSED this function at 25 parameters: "a call with 632// more arguments than operand slots would silently truncate -- refusing is the only honest answer" 633// (cap 23, already raised 4 -> 8 -> 16 -> 24 slots; capability=call-arity-ceiling, tracked on 634// /compare/lang). That is a REAL ceiling with a stated rationale, not a phony cap: the alternative is a 635// silently truncated call. Its named fix is a parameter struct, which is the right long-term shape and 636// is owed. Two parameters were removed instead, and BOTH removals are improvements in their own right: 637// focal -- now DERIVED as gs_focal_for(vh) rather than passed, so caller and callee can no longer 638// disagree about the lens. That disagreement is exactly the second-ruler defect this file 639// refuses everywhere else, so the ceiling forced the better design. 640// tcur -- the transmittance scratch reuses acc, which is DEAD inside this function from the moment 641// the render returns (the gradient reads fb and target, never acc) and is npx*3 >= npx. 642func gs_pos_grad_step_aniso(gauss: *i64, ng: i64, target: *i64, yaw: i64, camz: i64, fb: *i64, acc: *i64, trans: *i64, depth: *i64, sxb: *i64, syb: *i64, pa: *i64, pb: *i64, pc: *i64, pdet: *i64, order: *i64, count: *i64, explut: *i64, gradbuf: *i64, wnorm: *i64, vw: i64, vh: i64, usetrans: i64) -> i64 { 643 let focal: i64 = gs_focal_for(vh) 644 // ★SCRATCH LAYOUT (usetrans=2 needs the suffix colour, so the buffers are re-cast): 645 // trans[0..npx) -- RUNNING transmittance. Its post-render value is dead here, so the walk 646 // re-uses it rather than demanding a parameter the arity ceiling forbids. 647 // acc[0..3npx) -- acc_final, the render's accumulated colour. PRESERVED (v1 aliased tcur 648 // onto acc and destroyed exactly this; the indirect term needs it). 649 // acc[3npx..6npx) -- P, the running PREFIX colour. Caller must size acc npx*6 for usetrans=2. 650 // ★★AND THIS IS WHY NO BACKWARD PASS IS NEEDED. The textbook backward pass reconstructs T by 651 // dividing it out, which in integers divides by (GFXA-alpha) and loses precision exactly where 652 // alpha is large. The suffix is available in the SAME forward walk by arithmetic instead: 653 // S_i = acc_final - P_(i-1) - (this splat's own contribution) 654 // one subtraction, no reconstruction, no division by a near-zero, and it cannot drift from the 655 // renderer because acc_final IS the renderer's own output. 656 let tcur: *i64 = trans 657 let pbase: i64 = vw*vh*3 658 let vis: i64 = gs_render_aniso_at(gauss, ng, yaw, camz, fb, acc, trans, depth, sxb, syb, pa, pb, pc, pdet, order, count, explut, 26, 28, 44, vw, vh, focal) 659 let npx: i64 = vw * vh 660 var loss: i64 = 0 661 var i: i64 = 0 662 while i < npx { 663 var dr: i64 = (fb[i]&255) - (target[i]&255); if dr<0 {dr=0-dr} 664 var dg: i64 = ((fb[i]>>8)&255) - ((target[i]>>8)&255); if dg<0 {dg=0-dg} 665 var db: i64 = ((fb[i]>>16)&255) - ((target[i]>>16)&255); if db<0 {db=0-db} 666 loss = loss + dr + dg + db 667 i = i + 1 668 } 669 i = 0 670 while i < ng { gradbuf[i*3]=0; gradbuf[i*3+1]=0; gradbuf[i*3+2]=0; wnorm[i]=0; i=i+1 } 671 // transmittance starts FULL and is attenuated exactly as the renderer attenuates it, so the value 672 // this walk sees at a pixel is the value that pixel actually had when the splat was composited. 673 i = 0 674 while i < npx { tcur[i] = GFXA; i = i + 1 } 675 // the prefix starts empty: before any splat composites, nothing is in front of anything. 676 if usetrans == 2 { i = 0; while i < npx*3 { acc[pbase+i] = 0; i = i + 1 } } 677 // the prefix starts empty: before any splat composites, nothing is in front of anything. 678 if usetrans == 2 { i = 0; while i < npx*3 { acc[pbase+i] = 0; i = i + 1 } } 679 let sy4: i64 = it_sin4096(yaw) 680 let cy4: i64 = it_cos4096(yaw) 681 let R: i64 = camz * GFX 682 var oi: i64 = 0 683 while oi < vis { 684 let g: i64 = order[oi] 685 if depth[g] >= 0 { 686 let sx: i64 = sxb[g] 687 let sy: i64 = syb[g] 688 let ca: i64 = pa[g] 689 let cb: i64 = pb[g] 690 let cc: i64 = pc[g] 691 let det: i64 = pdet[g] 692 let op: i64 = gauss[g*GSTA+10] 693 let cr: i64 = gauss[g*GSTA+7] 694 let cg: i64 = gauss[g*GSTA+8] 695 let cbl: i64 = gauss[g*GSTA+9] 696 let cm2: i64 = cr*cr + cg*cg + cbl*cbl 697 let rx: i64 = gs_splat_rx(ca) 698 let ry: i64 = gs_splat_ry(cc) 699 var py: i64 = sy-ry 700 if py<0 {py=0} 701 var pye: i64 = sy+ry 702 if pye>=vh {pye=vh-1} 703 while py <= pye { 704 let dy: i64 = py-sy 705 var px: i64 = sx-rx 706 if px<0 {px=0} 707 var pxe: i64 = sx+rx 708 if pxe>=vw {pxe=vw-1} 709 while px <= pxe { 710 let dx: i64 = px-sx 711 let pnum: i64 = cc*dx*dx - 2*cb*dx*dy + ca*dy*dy 712 let k: i64 = pnum * GLUTU / det 713 if k < GEXPN { if k >= 0 { 714 let pix: i64 = py*vw+px 715 let w: i64 = op*explut[k]/GFXA 716 let rr: i64 = (fb[pix]&255) - (target[pix]&255) 717 let rg: i64 = ((fb[pix]>>8)&255) - ((target[pix]>>8)&255) 718 let rb: i64 = ((fb[pix]>>16)&255) - ((target[pix]>>16)&255) 719 let rc: i64 = rr*cr + rg*cg + rb*cbl 720 // T is the light still available at this pixel when THIS splat composites. 721 // usetrans=0 pins it full, which reproduces the incumbent arithmetic exactly. 722 var tw: i64 = GFXA 723 if usetrans > 0 { tw = tcur[pix] } 724 let dadx: i64 = w*(cc*dx - cb*dy)/det 725 let dady: i64 = w*(ca*dy - cb*dx)/det 726 if usetrans == 2 { 727 // ★THE INDIRECT TERM. Moving this splat changes the light reaching EVERY splat 728 // behind it: dT_j/dalpha_i = -T_j/(GFXA-alpha_i), so summing over j>i gives 729 // dC/dalpha_i = c_i*T_i/GFXA - S_i/(GFXA-alpha_i) 730 // with S_i the colour accumulated BEHIND this splat. S_i needs no backward 731 // pass: it is acc_final minus the running prefix minus this splat's own 732 // contribution -- three known quantities in a forward walk. 733 let cn: i64 = tw*w/GFXA 734 let ownr: i64 = cr*cn 735 let owng: i64 = cg*cn 736 let ownb: i64 = cbl*cn 737 let sr: i64 = acc[pix*3] - acc[pbase+pix*3] - ownr 738 let sg: i64 = acc[pix*3+1] - acc[pbase+pix*3+1] - owng 739 let sb: i64 = acc[pix*3+2] - acc[pbase+pix*3+2] - ownb 740 // a fully opaque splat leaves nothing behind it to disturb; clamping the 741 // divisor at 1 is the honest floor, not a tuned epsilon. 742 var dn: i64 = GFXA - w 743 if dn < 1 { dn = 1 } 744 // ⚠UNITS. v1 wrote sr/dn and the arm came out WORSE than direct-only (81 vs 79), 745 // which is what sent me to the algebra instead of to the tooth. In pixel units 746 // d(pixel)/d(alpha_i) = [ c_i*T_i - S_i/(GFXA-alpha_i) ] / GFXA^2 747 // and this accumulator carries the DIRECT part as c*T/GFXA, i.e. already one 748 // factor of GFXA down. The suffix therefore needs the SAME division or it 749 // enters 256x oversized -- a term that dominates the very gradient it corrects. 750 let er: i64 = cr*tw/GFXA - sr/dn/GFXA 751 let eg: i64 = cg*tw/GFXA - sg/dn/GFXA 752 let eb: i64 = cbl*tw/GFXA - sb/dn/GFXA 753 let coef: i64 = rr*er + rg*eg + rb*eb 754 let em2: i64 = er*er + eg*eg + eb*eb 755 gradbuf[g*3] = gradbuf[g*3] + coef*dadx 756 gradbuf[g*3+1] = gradbuf[g*3+1] + coef*dady 757 gradbuf[g*3+2] = gradbuf[g*3+2] + (em2*dadx/GFXA/GFXA)*dadx 758 wnorm[g] = wnorm[g] + (em2*dady/GFXA/GFXA)*dady 759 acc[pbase+pix*3] = acc[pbase+pix*3] + ownr 760 acc[pbase+pix*3+1] = acc[pbase+pix*3+1] + owng 761 acc[pbase+pix*3+2] = acc[pbase+pix*3+2] + ownb 762 } else { 763 let vsx: i64 = dadx*tw/GFXA 764 let vsy: i64 = dady*tw/GFXA 765 gradbuf[g*3] = gradbuf[g*3] + rc*vsx 766 gradbuf[g*3+1] = gradbuf[g*3+1] + rc*vsy 767 // squared term reduced BEFORE squaring so a near-degenerate conic cannot overflow 768 gradbuf[g*3+2] = gradbuf[g*3+2] + (cm2*vsx/GFXA)*vsx/GFXA 769 wnorm[g] = wnorm[g] + (cm2*vsy/GFXA)*vsy/GFXA 770 } 771 // attenuate for the splats behind this one, mirroring gs_blend_rect exactly 772 let tr: i64 = tcur[pix] 773 tcur[pix] = tr - tr*w/GFXA 774 } } 775 px = px + 1 776 } 777 py = py + 1 778 } 779 let cz: i64 = (gauss[g*GSTA]*sy4 + gauss[g*GSTA+2]*cy4)/K_MAGIC_4096 + R 780 if cz > GNEARZ { 781 let dnx: i64 = gradbuf[g*3+2]*2*GFXA*focal + 1 782 let dny: i64 = wnorm[g]*2*GFXA*focal + 1 783 gauss[g*GSTA] = gauss[g*GSTA] - gradbuf[g*3]*cz/dnx 784 gauss[g*GSTA+1] = gauss[g*GSTA+1] + gradbuf[g*3+1]*cz/dny 785 } 786 } 787 oi = oi + 1 788 } 789 return loss 790} 791 792func gs_w() -> i64 { return GW } 793func gs_h() -> i64 { return GH } 794func gs_expn() -> i64 { return GEXPN } 795func gs_nb() -> i64 { return GNB } 796 797// ===== ★MESH -> GAUSSIANS: THE INGEST HALF OF THE SPLAT LANE (2026-09-04) ===== 798// The estate could RENDER splats, FIT them by integer gradient descent, ADAPT their density and BIND 799// them to a driving cage -- but nothing could turn an authored or bought mesh INTO them, so a sculpt 800// could not enter this representation at all. This is that missing edge, and it deliberately emits 801// ANISOTROPIC SURFELS rather than isotropic blobs: a triangle KNOWS its normal exactly, so throwing 802// that away and re-deriving orientation later would be inventing work and losing precision. The file's 803// own aniso comment says why -- oriented disks give smooth coverage between verts instead of the dotty 804// blob look. 805// 806// ★THE ORACLE, AND WHY THIS IS NOT A PORT OF IT: EA SEED's Mesh2Splat converts glb to 3DGS in 807// milliseconds by exploiting the rasterizer's interpolator. It is BSD-3-Clause plus an EA marks clause 808// and it stays OUTSIDE our build path as a bar, never a dependency. It also emits ONE splat set per 809// model; the exceed axis declared on procgen.matrix before any run is PER-PART emission, which this 810// signature supports by construction because the caller passes one triangle range at a time. 811// 812// ★EVERY QUANTITY IS DERIVED, NONE IS PICKED (Rule 11). Sample count comes from DENSITY x AREA, not a 813// constant. The surfel radius is the radius whose disk area equals the surface each sample owns, so 814// splats just cover the face at any density -- a picked radius would gap at low density and overdraw at 815// high. Sampling is a stratified lattice folded into the triangle by reflection, so it is deterministic 816// and has no RNG to seed. 817// 818// ⚠NO SILENT CAP: stats[] reports what was REQUESTED beside what was WRITTEN, and truncation is decided 819// by comparing them rather than by a flag set mid-loop, so a buffer that happens to fill exactly is not 820// misreported as truncated. A cap reached in silence becomes a measurement nobody knows is partial. 821const GTRI_VS: i64 = 6 // per-vertex stride inside a triangle record: x y z r g b 822const GTRI_STRIDE: i64 = 18 // 3 * GTRI_VS -- the whole triangle record 823const GTRI_CO: i64 = 3 // colour offset inside one vertex record 824const GFM_HALF: i64 = 2 // the cross-product magnitude is TWICE the triangle area 825const GFM_MIN_SIDE: i64 = 1 // a face never vanishes: every triangle yields at least one surfel 826const GFM_MIN_RTAN: i64 = 1 // a zero-radius surfel would render as nothing at all 827const GFM_STATS: i64 = 4 // slots the caller must provide 828const GFM_ST_VISITED: i64 = 0 829const GFM_ST_TRUNC: i64 = 1 830const GFM_ST_WANT: i64 = 2 831const GFM_ST_CLAMP: i64 = 3 832 833func gs_from_mesh(tris: *i64, ntri: i64, gauss: *i64, cap: i64, density: i64, opacity: i64, stats: *i64) -> i64 { 834 stats[GFM_ST_VISITED] = 0 835 stats[GFM_ST_TRUNC] = 0 836 stats[GFM_ST_WANT] = 0 837 stats[GFM_ST_CLAMP] = 0 838 var ng: i64 = 0 839 var t: i64 = 0 840 while t < ntri { 841 let b: i64 = t * GTRI_STRIDE 842 let b1: i64 = b + GTRI_VS 843 let b2: i64 = b + GTRI_VS * 2 844 let e1x: i64 = tris[b1] - tris[b] 845 let e1y: i64 = tris[b1+1] - tris[b+1] 846 let e1z: i64 = tris[b1+2] - tris[b+2] 847 let e2x: i64 = tris[b2] - tris[b] 848 let e2y: i64 = tris[b2+1] - tris[b+1] 849 let e2z: i64 = tris[b2+2] - tris[b+2] 850 let cnx: i64 = e1y*e2z - e1z*e2y 851 let cny: i64 = e1z*e2x - e1x*e2z 852 let cnz: i64 = e1x*e2y - e1y*e2x 853 let clen: i64 = gs_isqrt(cnx*cnx + cny*cny + cnz*cnz) 854 if clen > 0 { 855 var side: i64 = gs_isqrt(density * clen / (GFM_HALF * GFX * GFX)) 856 if side < GFM_MIN_SIDE { side = GFM_MIN_SIDE; stats[GFM_ST_CLAMP] = stats[GFM_ST_CLAMP] + 1 } 857 let nsamp: i64 = side * side 858 var rtan: i64 = gs_isqrt(clen / (GFM_HALF * nsamp)) 859 if rtan < GFM_MIN_RTAN { rtan = GFM_MIN_RTAN } 860 let unx: i64 = cnx * GFXA / clen 861 let uny: i64 = cny * GFXA / clen 862 let unz: i64 = cnz * GFXA / clen 863 stats[GFM_ST_WANT] = stats[GFM_ST_WANT] + nsamp 864 var i: i64 = 0 865 while i < side { 866 var j: i64 = 0 867 while j < side { 868 if ng < cap { 869 var u: i64 = (2*i+1) * GFX / (2*side) 870 var v: i64 = (2*j+1) * GFX / (2*side) 871 if u + v > GFX { u = GFX - u; v = GFX - v } 872 let w: i64 = GFX - u - v 873 let px: i64 = tris[b] + (e1x*u + e2x*v)/GFX 874 let py: i64 = tris[b+1] + (e1y*u + e2y*v)/GFX 875 let pz: i64 = tris[b+2] + (e1z*u + e2z*v)/GFX 876 let cr: i64 = (tris[b+GTRI_CO]*w + tris[b1+GTRI_CO]*u + tris[b2+GTRI_CO]*v)/GFX 877 let cg: i64 = (tris[b+GTRI_CO+1]*w + tris[b1+GTRI_CO+1]*u + tris[b2+GTRI_CO+1]*v)/GFX 878 let cb: i64 = (tris[b+GTRI_CO+2]*w + tris[b1+GTRI_CO+2]*u + tris[b2+GTRI_CO+2]*v)/GFX 879 gs_set_aniso(gauss, ng, px, py, pz, unx, uny, unz, rtan, cr, cg, cb, opacity) 880 ng = ng + 1 881 } 882 j = j + 1 883 } 884 i = i + 1 885 } 886 } 887 stats[GFM_ST_VISITED] = stats[GFM_ST_VISITED] + 1 888 t = t + 1 889 } 890 if stats[GFM_ST_WANT] > ng { stats[GFM_ST_TRUNC] = 1 } 891 return ng 892} 893 894func gs_tri_stride() -> i64 { return GTRI_STRIDE } 895func gs_fm_stats() -> i64 { return GFM_STATS }