Variable GSPLAT_VERTEX_SHADERConst
GSPLAT_VERTEX_SHADER: "\n precision highp float;\n\n \nbool isInvalidFloat(float v) {\n return isnan(v) || isinf(v);\n}\n\nfloat sanitizePositive(float v, float fallback) {\n return (isInvalidFloat(v) || v <= 0.0) ? fallback : v;\n}\n\nfloat sanitizeNonNegative(float v, float fallback) {\n return (isInvalidFloat(v) || v < 0.0) ? fallback : v;\n}\n\n// Per-element opacity sanitizer: NaN/Inf route to the 1.0 opaque\n// identity (corruption stays LOUD), finite values clamp to [0, 1]\n// (alpha is opacity, never HDR — Python pins the range at write; this\n// guards hand-crafted zarr). The clamp keeps the zero boundary\n// CONTINUOUS (a -1e-4 epsilon vanishes like +0.0 renders, instead of\n// jumping to full opacity) and keeps the value mediump-varying-safe.\nfloat sanitizeAlpha(float v) {\n return isInvalidFloat(v) ? 1.0 : clamp(v, 0.0, 1.0);\n}\n\n \nfloat perspectiveNearFade(int isOrtho, float viewZ, float nearCull) {\n if (isOrtho == 1) return 1.0;\n if (viewZ >= 0.0) return 0.0;\n return smoothstep(nearCull, nearCull * 2.0, -viewZ);\n}\n\n \nint luxarIsOrthoProjection() {\n return projectionMatrix[3][3] > 0.5 ? 1 : 0;\n}\nfloat luxarProjectionSizeScale() {\n return abs(projectionMatrix[1][1]);\n}\n\n\n // Quad corner attribute (static geometry)\n in vec2 aQuadCorner; // (-1,-1), (1,-1), (-1,1), (1,1)\n\n // Draw-slot → storage-slot mapping, double-buffered so a new ordering\n // swaps atomically (declaration + luxarSortedIndex() in glsl-lib).\n // Uint32Array attributes → bound via vertexAttribIPointer, matching\n // the uint declarations.\n \n// Element index split into two 16-bit halves, low in .x and high in .y.\n// The pick pass carries the index through an RGBA32F buffer, and float32\n// has a 24-bit mantissa — so a single float channel cannot represent\n// consecutive indices past 16,777,216, while a node's capacity reaches\n// 2^25 on a 32768-texel device. Both halves are <= 65535, hence exact,\n// and the pick decoder recombines them (see picking-system/pick-render.ts).\n// Kept in INT space: doing the split on a float would already have lost\n// the bit it is meant to preserve.\nvec2 luxarElementIdSplit(uint i) {\n return vec2(float(i & 0xFFFFu), float(i >> 16u));\n}\n\nin uint aSortedIndex;\nin uint aSortedIndexB;\nuniform int uSortedIndexSlot;\n\nuint luxarSortedIndex() {\n return uSortedIndexSlot == 1 ? aSortedIndexB : aSortedIndex;\n}\n\n// The STORAGE slot's id parts — the id the rest of the pipeline (loaders,\n// selection) addresses elements by, not the transient draw slot. Mesh reads\n// the shared split directly instead, off gl_VertexID (spec §6.5).\nvec2 luxarElementIdParts() {\n return luxarElementIdSplit(luxarSortedIndex());\n}\n// Projected-density thinning (scene/density-guard.ts): the fraction of this\n// node's elements to DROP, chosen per STORAGE index with a deterministic\n// integer hash so the kept subset is stable under depth re-sorting, identical\n// across the visual and picking passes, and spatially uniform (storage order\n// is Hilbert/BSP-coherent, so a prefix would be a hole). A material that does\n// not set the uniform reads 0 and drops nothing.\nuniform float uDensityDrop;\nbool luxarDensityDropped() {\n if (uDensityDrop <= 0.0) return false;\n uint h = luxarSortedIndex();\n h ^= h >> 16u;\n h *= 0x7feb352du;\n h ^= h >> 15u;\n h *= 0x846ca68bu;\n h ^= h >> 16u;\n return float(h) * (1.0 / 4294967296.0) < uDensityDrop;\n}\n\n\n // Splat data texture: RGBA32F, 4 texels/splat (see\n // rendering/element-texture-layout.ts for the texel layout).\n uniform highp sampler2D uSplatTex;\n\n // Uniforms (modelViewMatrix and projectionMatrix are built-in THREE.js uniforms)\n uniform vec2 uResolution;\n uniform float uTruncate; // Truncation radius (in sigmas)\n uniform float uRayIntegralFactor; // Shifted Gaussian ray integral factor\n uniform int uProjectionMode; // 0 = sum (ray-integral: additive/luminous/volumetric), 1 = peak (2D-projected surfaces: max/normal/opaque)\n uniform float uNearCull; // Near cull distance (scene-scale-aware)\n uniform float uMaxExtentFactor; // Max projected extent as fraction of viewport before fade\n uniform float uCov2DDilation; // 2D-covariance low-pass dilation in CSS px² (3DGS anti-aliasing)\n uniform float uPixelRatio;\n uniform int uLabelColorMode;\n uniform int uLabelFilterIndex;\n\n // Colormap uniforms (only active when USE_COLORMAP is defined)\n #ifdef USE_COLORMAP\n uniform sampler2D uColormapTex; // 256x1 LUT texture\n uniform float uScalarMin; // Scalar range minimum (display-range window)\n uniform float uScalarScale; // 1.0 / (max - min)\n uniform mediump float uInvGamma; // Gamma applied to the VALUE, pre-LUT\n #endif\n\n // Varyings to fragment - all per-instance varyings use \"flat\" (no interpolation needed)\n // OPTIMIZATION: flat qualifier skips GPU interpolation hardware for constant values\n flat out mediump vec3 vColor;\n flat out mediump float vAmplitude2D;\n flat out mediump float vAlpha; // per-splat opacity (texel3.y; 1.0 for RGB data)\n // These need highp for screen-space calculations\n // OPTIMIZATION: vL2D stores [1/L00, L10, 1/L11] to replace fragment divisions with multiplications\n flat out highp vec3 vL2D; // 2D Cholesky packed as [invL00, L10, invL11]\n flat out highp vec2 vCenterScreen; // Splat center in screen pixels\n\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 // invalidFloat is an alias for the shared isInvalidFloat helper in glsl-lib.\n bool invalidFloat(float v) {\n return isInvalidFloat(v);\n }\n\n bool invalidCov2D(mat2 S) {\n return invalidFloat(S[0][0]) || invalidFloat(S[0][1]) || invalidFloat(S[1][0]) || invalidFloat(S[1][1]);\n }\n\n // Compute 2D Cholesky from 2D covariance (symmetric positive definite)\n // OPTIMIZATION: Returns [1/L00, L10, 1/L11] for faster fragment shader (MUL instead of DIV)\n vec3 cholesky2x2(mat2 S) {\n float s00 = max(S[0][0], 1e-6);\n float s10 = invalidFloat(S[1][0]) ? 0.0 : S[1][0];\n float s11 = invalidFloat(S[1][1]) ? 1e-6 : S[1][1];\n float L00 = sqrt(s00);\n float invL00 = 1.0 / L00;\n float L10 = s10 * invL00; // Use reciprocal here too\n float L11 = sqrt(max(s11 - L10 * L10, 1e-6));\n float invL11 = 1.0 / L11;\n return vec3(invL00, L10, invL11); // Pack reciprocals for fragment shader\n }\n\n vec3 categoricalColor(float index) {\n return 0.25 + 0.75 * fract(index * vec3(0.61803398875, 0.38196601125, 0.75487766625));\n }\n\n void main() {\n // === Splat-texture fetch prologue ===\n // Four texelFetch reads reconstruct the per-splat values into\n // the exact local names the math below has always used — zero\n // changes downstream of this block. The width is a multiple of\n // 4 (element-texture-layout.ts), so a splat's 4 texels share one\n // row and only x advances.\n // Projected-density thinning (density-guard): drop this instance\n // before any texel fetch; the rasterizer discards a z=-2 vertex.\n if (luxarDensityDropped()) {\n gl_Position = vec4(0.0, 0.0, -2.0, 1.0);\n return;\n }\n int splatBase = int(luxarSortedIndex()) * 4;\n int splatTexW = LUXAR_SPLAT_TEX_W;\n ivec2 texel0 = ivec2(splatBase % splatTexW, splatBase / splatTexW);\n vec4 splatT0 = texelFetch(uSplatTex, texel0, 0);\n vec4 splatT1 = texelFetch(uSplatTex, ivec2(texel0.x + 1, texel0.y), 0);\n vec4 splatT2 = texelFetch(uSplatTex, ivec2(texel0.x + 2, texel0.y), 0);\n vec4 splatT3 = texelFetch(uSplatTex, ivec2(texel0.x + 3, texel0.y), 0);\n vec3 aCenter = splatT0.xyz; // 3D center (after nD slicing)\n float aAmplitude = splatT0.w; // Already attenuated by hidden dims\n vec2 aCholesky01 = splatT1.xy; // [L00, L10]\n vec2 aCholesky23 = splatT1.zw; // [L11, L20]\n vec2 aCholesky45 = splatT2.xy; // [L21, L22]\n vec3 aColor = vec3(splatT2.zw, splatT3.x);\n float aAlpha = splatT3.y; // per-splat opacity (1.0 when the dataset is RGB)\n float aLabelIndex = splatT3.z;\n if (uLabelFilterIndex > 0 && int(aLabelIndex + 0.5) != uLabelFilterIndex) {\n gl_Position = vec4(0.0, 0.0, -2.0, 1.0);\n return;\n }\n\n // Transform center to camera space\n vec4 centerCam4 = modelViewMatrix * vec4(aCenter, 1.0);\n vec3 centerCam = centerCam4.xyz;\n\n // Clip-space centre through the projection THIS draw uses. The screen\n // centre, the covariance Jacobian, the coverage extent and the ortho\n // branch all come from it and from P, so a splat is placed and sized\n // for whatever camera three draws with: a cube-capture face (fov -90\n // flips P), a zoomed or asymmetric frustum, an embedder's camera.\n // CPU mirror + tests: projection-math.ts.\n vec4 centerClip = projectionMatrix * centerCam4;\n int isOrtho = luxarIsOrthoProjection();\n float invW = 1.0 / centerClip.w;\n\n // === Unified near handling (shared perspectiveNearFade helper;\n // point + line shaders use the same) ===\n // Perspective: behind-camera splats fade to 0 (subsumes the old\n // standalone centerCam.z >= 0 reject) and the near-plane\n // approach fades across [uNearCull, 2*uNearCull] — prevents the\n // 1/z Jacobian singularity. Ortho: fade = 1; a behind-camera\n // splat falls through to NDC clipping, which drops it (ortho\n // near > 0 in this viewer), and no 1/z is consumed on the\n // ortho path.\n // uNearCull is scene-bounds-scaled (diagonal * 0.001); the\n // 1e-20 floor only guards uNearCull == 0 (degenerate\n // smoothstep). An absolute 1e-4 floor overrode the\n // scene-relative value on tiny-unit scenes and faded out the\n // whole scene. (Consequence: surviving zDepth is only bounded\n // by ~uNearCull, so the unguarded 1/zDepth below can get large\n // on a sub-camera-plane splat — J/Sigma2D then go non-finite\n // and the invalidCov2D reject drops the splat safely.)\n float depthFade = perspectiveNearFade(isOrtho, centerCam.z, max(uNearCull, 1e-20));\n if (depthFade < 0.01) {\n gl_Position = vec4(0.0, 0.0, -2.0, 1.0);\n return;\n }\n\n // Transform Cholesky to camera space (rotation only)\n mat3 R = mat3(modelViewMatrix);\n mat3 L3D = unpackCholesky3D(aCholesky01, aCholesky23, aCholesky45);\n mat3 L_cam = R * L3D;\n mat3 Sigma_cam = L_cam * transpose(L_cam);\n\n // Positive depth (camera looks down -Z); reused by fades, Jacobian, and projection\n float zDepth = -centerCam.z;\n\n // === Screen-coverage safety guard (independent of 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 (below) kicks in.\n // This avoids visible hard edges from clamped quads. Applies in BOTH\n // projections (ortho projected size is depth-independent, divisor 1)\n // and is computed UNCONDITIONALLY: below maxExtent*0.5 the smoothstep\n // is 0 and the fade is a no-op, so no size gate is needed. (The former\n // absolute maxLateralVar > 0.01 gate — a perf leftover from the\n // pre-#51 two-stage near cull — skipped the fade for splats with\n // spatial sigma < 0.1 world units while the extent clamp still\n // applied, leaving hard-edged clamped rectangles on deep-zoomed\n // tiny-sigma / nm-unit-scale scenes.) The 1e-20 floors match\n // the TSL twin's expressions exactly and are pure\n // div-by-zero/sqrt guards, NOT scale floors: maxLateralVar is a\n // WORLD-unit² variance, so the old absolute 1e-8 floor inflated\n // valid tiny-unit variances (sigma ~ 1e-7 => var ~ 1e-14) up to\n // sqrt(1e-8) = 1e-4 world units — projectedExtent exploded and\n // coverageFade culled EVERY splat in the scene. zDepth is\n // likewise bounded below by the scene-relative near fade.\n float halfResY = 0.5 * uResolution.y;\n float coverageFade;\n {\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)) * uTruncate / extentDivisor;\n float maxExtent = max(uResolution.x, uResolution.y) * uMaxExtentFactor;\n coverageFade = 1.0 - smoothstep(maxExtent * 0.5, maxExtent, projectedExtent);\n if (coverageFade < 0.01) {\n gl_Position = vec4(0.0, 0.0, -2.0, 1.0);\n return;\n }\n }\n\n // Combined fade: most restrictive wins (both use smoothstep → no popping)\n float nearFade = min(depthFade, coverageFade);\n\n // Projection Jacobian at the splat centre, general form (valid for any\n // P): J[k] = res/2 * (P[k].xy / w - clip.xy * P[k].w / w^2). For\n // three's symmetric perspective P it is the classic\n // [[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\n // (not a shared fx = fy), so a camera aspect that differs from the\n // buffer aspect is honoured instead of assumed away.\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 // Project covariance to 2D: Σ_2D = J · Σ_cam · Jᵀ\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\n // 2D low-pass dilation (standard 3DGS anti-aliasing): widen the diagonal\n // so every splat covers at least ~1px. Guarantees near-degenerate\n // (edge-on/flat) splats render as a soft ellipse instead of a razor-thin\n // sub-pixel spike. Applied before the Cholesky + eigen extent below so\n // the fragment footprint and the quad stay consistent. Diagonal only —\n // adding to the off-diagonal would rotate/shear the ellipse.\n //\n // ENERGY COMPENSATION (Mip-Splatting): widening the footprint without\n // touching the peak CREATES light — a 2D Gaussian's screen-integrated\n // brightness is 2*pi*peak*sqrt(det Sigma2D), so dilation inflates it by\n // sqrt(detDilated/detRaw), i.e. (sigma_px^2 + d)/sigma_px^2 for an\n // isotropic splat; the compensation multiplier below is the reciprocal,\n // sqrt(detRaw/detDilated). The inflation diverges as the splat shrinks on screen\n // (measured 1.43x at 0.84 px, 3.75x at 0.33 px), so a scene silently\n // brightened as the camera pulled back, and a lifted points->gsplat LOD\n // ladder could not match its Points level in ANY mode. Points already\n // compensate their own sub-pixel widening (the sizeScale^2 term in\n // materials/point/shader-glsl.ts); gsplats now do too.\n float detRaw2D = Sigma2D[0][0] * Sigma2D[1][1] - Sigma2D[0][1] * Sigma2D[1][0];\n // Keep the historical framebuffer-pixel low-pass below 1× render scale.\n float dilationPixelRatio = max(uPixelRatio, 1.0);\n float cov2DDilation = uCov2DDilation * dilationPixelRatio * dilationPixelRatio;\n Sigma2D[0][0] += cov2DDilation;\n Sigma2D[1][1] += cov2DDilation;\n float detDilated2D = Sigma2D[0][0] * Sigma2D[1][1] - Sigma2D[0][1] * Sigma2D[1][0];\n float dilationCompensation = sqrt(max(detRaw2D, 0.0) / max(detDilated2D, 1e-12));\n\n if (invalidCov2D(Sigma2D) || invalidFloat(aAmplitude)) {\n gl_Position = vec4(0.0, 0.0, -2.0, 1.0);\n return;\n }\n\n // Projection mode determines amplitude calculation:\n // - Sum projection (uProjectionMode = 0): Integrate Gaussian along ray → ray boost\n // - Max projection (uProjectionMode = 1): Use peak Gaussian value → no boost\n //\n // OPTIMIZATION: Use branch instead of branchless mix() to skip expensive operations\n // (normalize, sqrt, exp) when in max mode. Warps are typically coherent on this uniform.\n float sigmaRay = 1.0; // Default for max mode (no ray integration)\n if (uProjectionMode == 0) {\n // Sum projection: compute ray-integral standard deviation.\n //\n // The line-integral of an anisotropic Gaussian along ray\n // direction r has 1D std-dev sigma_line = 1 / sqrt(rᵀ Σ⁻¹ r),\n // not sqrt(rᵀ Σ r). The two only agree when r is aligned with\n // a covariance eigenvector or Σ is isotropic.\n //\n // Implementation: compute Σ_cam⁻¹ via the closed-form 3×3\n // inverse and clamp to a minimum determinant. The shader is\n // already paying for a covariance matrix-vector product, so\n // a one-off explicit inverse is a small constant factor.\n vec3 rayDir = (isOrtho == 1) ? vec3(0.0, 0.0, -1.0) : normalize(centerCam);\n // Cofactor expansion for 3x3 inverse. Σ_cam is symmetric SPD,\n // so the inverse is symmetric SPD too.\n // SCALE-FREE inversion: normalize Σ_cam by its mean diagonal\n // variance s = trace/3 first. det(Σ) is world-units⁶ — on a\n // tiny-unit scene (sigma ~ 1e-7 => det ~ 1e-42) it\n // underflows float32 (GPUs flush denormals to zero) and the\n // absolute 1e-12 clamp turned Σ⁻¹ into garbage (sum-mode\n // brightness off by many orders of magnitude); huge-unit\n // scenes overflow the same way. With Σn = Σ/s the\n // determinant and ray quadratic are O(1) at ANY scene\n // scale, so the 1e-12 / 1e-8 floors below act as pure\n // scale-free CONDITION-NUMBER guards. Σ⁻¹ = Σn⁻¹ / s, so\n // sigmaRay = sqrt(s / quadN) restores the world-unit\n // result exactly. The 1e-30 floor on s only guards an\n // all-zero (degenerate) covariance.\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\n // well-conditioned splat at any scale); un-normalize via\n // sqrt(sTrace): sigmaRay = 1/sqrt(rᵀ Σ⁻¹ r) = sqrt(s/quadN).\n float quad = max(rayDir.x * prx + rayDir.y * pry + rayDir.z * prz, 1e-8);\n sigmaRay = inversesqrt(quad) * sqrt(sTrace);\n // Shifted Gaussian ray integral: sqrt(2π)·erf(T/√2) - 2·T·exp(-0.5·T²)\n // Precomputed in TypeScript as uRayIntegralFactor (≈2.433 for T=3)\n float rayIntegrationBoost = sigmaRay * uRayIntegralFactor; // voxelSpacing = 1.0\n // dilationCompensation keeps the screen-integrated light invariant\n // under the 2D low-pass above (see its derivation there). Sum\n // projection only: this branch's quantity IS that integral, so\n // conserving it is exactly right. The peak branch below reports a\n // peak, not an integral, and is left alone.\n vAmplitude2D = aAmplitude * rayIntegrationBoost * nearFade * dilationCompensation;\n } else {\n // Max projection: no boost needed\n vAmplitude2D = aAmplitude * nearFade;\n }\n\n // Compute 2D Cholesky for fragment shader\n vL2D = cholesky2x2(Sigma2D);\n\n // Eigenvalues of Σ_2D for quad extents (oriented quad)\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 float lambda1 = max(0.5 * (trace + sqrtDisc), 1e-6);\n float lambda2 = max(0.5 * (trace - sqrtDisc), 1e-6);\n\n // Eigenvector for major axis (for oriented quad)\n vec2 majorAxis;\n if (abs(Sigma2D[0][1]) > 1e-6) {\n majorAxis = normalize(vec2(lambda1 - Sigma2D[1][1], Sigma2D[0][1]));\n } else {\n // Near-diagonal covariance: pick axis with larger variance\n majorAxis = (Sigma2D[0][0] >= Sigma2D[1][1]) ? vec2(1.0, 0.0) : vec2(0.0, 1.0);\n }\n vec2 minorAxis = vec2(-majorAxis.y, majorAxis.x);\n\n // Quad extents: truncation radius × sqrt(eigenvalue)\n // Shifted Gaussian: effectiveTruncate = uTruncate\n float extent1 = uTruncate * sqrt(lambda1);\n float extent2 = uTruncate * sqrt(lambda2);\n\n // Clamp quad extents so no splat exceeds uMaxExtentFactor × viewport.\n // The amplitude fade (nearFade) handles the visual transition smoothly;\n // this clamp prevents the rasterizer from shading oversized quads.\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\n // Screen centre in pixels from the clip-space centre (gl_FragCoord\n // convention: origin at the viewport's bottom-left corner).\n vCenterScreen = (centerClip.xy * invW * 0.5 + 0.5) * uResolution;\n\n // Expand quad vertex in screen space (oriented)\n vec2 quadOffset = aQuadCorner.x * majorAxis * extent1\n + aQuadCorner.y * minorAxis * extent2;\n vec2 screenPos = vCenterScreen + quadOffset;\n\n // Convert screen pixels to NDC (xy only)\n vec2 ndcXY = (screenPos / uResolution) * 2.0 - 1.0;\n\n // Pass through color — either from vertex attribute or colormap LUT.\n // Colormap mode: display range (uScalarMin/uScalarScale) and gamma\n // operate on the scalar VALUE (here the amplitude) before the LUT\n // lookup, not on the resulting color. See fragment-shader note.\n if (uLabelColorMode == 1 && aLabelIndex > 0.0) {\n vColor = categoricalColor(aLabelIndex);\n } else {\n #ifdef USE_COLORMAP\n float t = clamp((aAmplitude - uScalarMin) * uScalarScale, 0.0, 1.0);\n #ifndef LUXAR_GAMMA_ONE\n t = pow(t, uInvGamma); // gamma on the value, pre-LUT\n #endif\n vColor = texture(uColormapTex, vec2(t, 0.5)).rgb;\n #else\n vColor = aColor;\n #endif\n }\n // Per-splat opacity rides regardless of color source (in colormap\n // mode an RGBA dataset keeps its alpha; RGB data carries 1.0).\n // Sanitized: alpha is load-bearing in EVERY mode (linear\n // contribution scale) and maps into optical depth under\n // volumetric, where a NaN/Inf poisons τ past the discard into\n // NaN pixels — and a huge finite value would blow out the\n // linear folds (or overflow the mediump varying). Python\n // validation pins alpha to [0, 1] at write; this guards\n // hand-crafted zarr. NaN/Inf → the 1.0 opaque identity (loud);\n // finite values clamp to [0, 1] (a negative epsilon vanishes\n // continuously instead of flipping opaque). The point twin\n // does the same.\n vAlpha = sanitizeAlpha(aAlpha);\n\n // Depth through the same projection (the clip-space centre above).\n float ndcZ = centerClip.z / centerClip.w;\n\n // Output final clip position\n gl_Position = vec4(ndcXY, ndcZ, 1.0);\n }\n " = ...