Variable GLSL_GSPLAT_PROJECTIONConst
GLSL_GSPLAT_PROJECTION: "\n // Unpack 3D Cholesky to matrix (column-major order for GLSL mat3)\n // Packed order: [L00, L10, L11, L20, L21, L22]\n // c01 = [L00, L10], c23 = [L11, L20], c45 = [L21, L22]\n mat3 unpackCholesky3D(vec2 c01, vec2 c23, vec2 c45) {\n return mat3(\n c01.x, c01.y, c23.y, // Column 0: [L00, L10, L20]\n 0.0, c23.x, c45.x, // Column 1: [0, L11, L21]\n 0.0, 0.0, c45.y // Column 2: [0, 0, L22]\n );\n }\n\n // Camera-space covariance: rotate the Cholesky factor by the model-view\n // rotation, Σ_cam = (R·L)(R·L)ᵀ.\n mat3 gsplatSigmaCam(vec2 c01, vec2 c23, vec2 c45) {\n mat3 R = mat3(modelViewMatrix);\n mat3 L3D = unpackCholesky3D(c01, c23, c45);\n mat3 L_cam = R * L3D;\n return L_cam * transpose(L_cam);\n }\n\n // A Σ_2D entry that is NaN/Inf would poison the eigendecomposition; both\n // stages reject such a splat.\n bool invalidCov2D(mat2 S) {\n return isInvalidFloat(S[0][0]) || isInvalidFloat(S[0][1]) || isInvalidFloat(S[1][0]) || isInvalidFloat(S[1][1]);\n }\n\n // === Screen-coverage safety fade (independent of the depth fade) ===\n // Prevents GPU overload from splats whose projected quad is too large.\n // Fade starts at 50% of the limit and reaches ~0% AT the limit, so the\n // amplitude is negligible before the extent clamp kicks in. Applies in\n // BOTH projections (ortho projected size is depth-independent, divisor 1)\n // and is computed UNCONDITIONALLY: below maxExtent*0.5 the smoothstep is 0\n // and the fade is a no-op, so no size gate is needed. (A former absolute\n // maxLateralVar > 0.01 gate skipped the fade for spatial sigma < 0.1 world\n // units while the extent clamp still applied, leaving hard-edged clamped\n // rectangles on deep-zoomed tiny-unit scenes.) The 1e-20 floors are pure\n // div-by-zero/sqrt guards, NOT scale floors: maxLateralVar is a WORLD-unit²\n // variance, and an absolute 1e-8 floor once coverage-culled EVERY splat of\n // a tiny-unit scene. zDepth is bounded below by the scene-relative near\n // fade. truncate is the DRAW pass's T (the pick stage passes its\n // uCoverageTruncate, synced from the visual material).\n float gsplatCoverageFade(mat3 Sigma_cam, float zDepth, int isOrtho, float truncate) {\n float halfResY = 0.5 * uResolution.y;\n float maxLateralVar = max(Sigma_cam[0][0], max(Sigma_cam[1][1], Sigma_cam[2][2]));\n float extentDivisor = (isOrtho == 1) ? 1.0 : max(zDepth, 1e-20);\n float projectedExtent = (halfResY * luxarProjectionSizeScale()) * sqrt(max(maxLateralVar, 1e-20)) * truncate / extentDivisor;\n float maxExtent = max(uResolution.x, uResolution.y) * uMaxExtentFactor;\n return 1.0 - smoothstep(maxExtent * 0.5, maxExtent, projectedExtent);\n }\n\n // Σ_2D = J · Σ_cam · Jᵀ with the projection Jacobian at the splat centre,\n // general form (valid for any P): J[k] = res/2 * (P[k].xy / w -\n // clip.xy * P[k].w / w^2). For three's symmetric perspective P it is the\n // classic [[fx/z, 0], [0, fy/z], [fx*x/z^2, fy*y/z^2]]; for ortho (w = 1,\n // P[k].w = 0) it is [[fx, 0], [0, fy], [0, 0]]. The x terms use P00 (not a\n // shared fx = fy), so a camera aspect that differs from the buffer aspect\n // is honoured instead of assumed away. Before the 2D low-pass dilation.\n mat2 gsplatProjectedCovariance(mat3 Sigma_cam, vec4 centerClip, float invW) {\n vec2 halfRes = 0.5 * uResolution;\n vec2 clipTerm = centerClip.xy * (invW * invW);\n mat3x2 J;\n J[0] = halfRes * (projectionMatrix[0].xy * invW - clipTerm * projectionMatrix[0].w);\n J[1] = halfRes * (projectionMatrix[1].xy * invW - clipTerm * projectionMatrix[1].w);\n J[2] = halfRes * (projectionMatrix[2].xy * invW - clipTerm * projectionMatrix[2].w);\n\n // Compute J * Sigma_cam first\n vec2 JS0 = J[0] * Sigma_cam[0][0] + J[1] * Sigma_cam[0][1] + J[2] * Sigma_cam[0][2];\n vec2 JS1 = J[0] * Sigma_cam[1][0] + J[1] * Sigma_cam[1][1] + J[2] * Sigma_cam[1][2];\n vec2 JS2 = J[0] * Sigma_cam[2][0] + J[1] * Sigma_cam[2][1] + J[2] * Sigma_cam[2][2];\n\n // Sigma2D = (J * Sigma_cam) * J^T\n // For M = J*Sigma (cols JS0, JS1, JS2), compute M * J^T:\n // (M*J^T)[i,j] = sum_k M[i,k] * J[j,k] = sum_k JS_k[i] * J[k][j]\n mat2 Sigma2D;\n Sigma2D[0][0] = JS0.x * J[0].x + JS1.x * J[1].x + JS2.x * J[2].x;\n Sigma2D[1][0] = JS0.x * J[0].y + JS1.x * J[1].y + JS2.x * J[2].y;\n Sigma2D[0][1] = Sigma2D[1][0]; // Symmetric (J*S*J^T preserves symmetry)\n Sigma2D[1][1] = JS0.y * J[0].y + JS1.y * J[1].y + JS2.y * J[2].y;\n return Sigma2D;\n }\n\n // Sum projection: the 1D std-dev of the splat along the view ray,\n // sigma_line = 1 / sqrt(rᵀ Σ⁻¹ r) — not sqrt(rᵀ Σ r); the two only agree\n // when r is aligned with a covariance eigenvector or Σ is isotropic.\n // SCALE-FREE inversion: Σ_cam is normalized by its mean diagonal variance\n // s = trace/3 first. det(Σ) is world-units⁶ — on a tiny-unit scene (sigma\n // ~ 1e-7 => det ~ 1e-42) it underflows float32 (GPUs flush denormals to\n // zero) and an absolute 1e-12 clamp turned Σ⁻¹ into garbage; huge-unit\n // scenes overflow the same way. With Σn = Σ/s the determinant and ray\n // quadratic are O(1) at ANY scene scale, so the 1e-12 / 1e-8 floors below\n // act as pure scale-free CONDITION-NUMBER guards. Σ⁻¹ = Σn⁻¹ / s, so\n // sigmaRay = sqrt(s / quadN) restores the world-unit result exactly. The\n // 1e-30 floor on s only guards an all-zero (degenerate) covariance.\n float gsplatRaySigma(mat3 Sigma_cam, vec3 rayDir) {\n float sTrace = max((Sigma_cam[0][0] + Sigma_cam[1][1] + Sigma_cam[2][2]) * (1.0 / 3.0), 1e-30);\n float invS = 1.0 / sTrace;\n float a = Sigma_cam[0][0] * invS;\n float b = Sigma_cam[0][1] * invS;\n float c = Sigma_cam[0][2] * invS;\n float d = Sigma_cam[1][1] * invS;\n float e = Sigma_cam[1][2] * invS;\n float f = Sigma_cam[2][2] * invS;\n // det(Σn) for 3x3 symmetric — clamped against numerical singularity.\n float detSigma = a * (d * f - e * e) - b * (b * f - c * e) + c * (b * e - c * d);\n float invDet = 1.0 / max(detSigma, 1e-12);\n // Cofactors of the inverse (symmetric).\n float i00 = (d * f - e * e) * invDet;\n float i11 = (a * f - c * c) * invDet;\n float i22 = (a * d - b * b) * invDet;\n float i01 = -(b * f - c * e) * invDet;\n float i02 = (b * e - c * d) * invDet;\n float i12 = -(a * e - b * c) * invDet;\n // r' = Σ⁻¹ r ; precision quadratic = rᵀ Σ⁻¹ r\n float prx = i00 * rayDir.x + i01 * rayDir.y + i02 * rayDir.z;\n float pry = i01 * rayDir.x + i11 * rayDir.y + i12 * rayDir.z;\n float prz = i02 * rayDir.x + i12 * rayDir.y + i22 * rayDir.z;\n // quad is rᵀ Σn⁻¹ r (normalized space, O(1) for a well-conditioned\n // splat at any scale); un-normalize via sqrt(sTrace).\n float quad = max(rayDir.x * prx + rayDir.y * pry + rayDir.z * prz, 1e-8);\n return inversesqrt(quad) * sqrt(sTrace);\n }\n\n // Eigenvalues of the (dilated) Σ_2D, (lambda1, lambda2), for the oriented quad.\n vec2 gsplatEigenvalues(mat2 Sigma2D) {\n float trace = Sigma2D[0][0] + Sigma2D[1][1];\n float det = Sigma2D[0][0] * Sigma2D[1][1] - Sigma2D[0][1] * Sigma2D[1][0];\n float disc = max(trace * trace - 4.0 * det, 0.0); // Clamp for numerical stability\n float sqrtDisc = sqrt(disc);\n return vec2(max(0.5 * (trace + sqrtDisc), 1e-6), max(0.5 * (trace - sqrtDisc), 1e-6));\n }\n\n // Eigenvector of the major axis (for the oriented quad).\n vec2 gsplatMajorAxis(mat2 Sigma2D, float lambda1) {\n if (abs(Sigma2D[0][1]) > 1e-6) {\n return normalize(vec2(lambda1 - Sigma2D[1][1], Sigma2D[0][1]));\n }\n // Near-diagonal covariance: pick axis with larger variance\n return (Sigma2D[0][0] >= Sigma2D[1][1]) ? vec2(1.0, 0.0) : vec2(0.0, 1.0);\n }\n\n // Quad half extents truncate × sqrt(eigenvalue), clamped so no splat\n // exceeds uMaxExtentFactor × viewport. The amplitude fade (coverage)\n // handles the visual transition smoothly; this clamp prevents the\n // rasterizer from shading oversized quads.\n vec2 gsplatClampedExtents(float truncate, vec2 lambdas) {\n float extent1 = truncate * sqrt(lambdas.x);\n float extent2 = truncate * sqrt(lambdas.y);\n float maxExtentPx = max(uResolution.x, uResolution.y) * uMaxExtentFactor;\n float largestExtent = max(extent1, extent2);\n if (largestExtent > maxExtentPx) {\n float clampScale = maxExtentPx / largestExtent;\n extent1 *= clampScale;\n extent2 *= clampScale;\n }\n return vec2(extent1, extent2);\n }\n" = ...