The 2D Gaussian: Evaluation & Bounds

In the Slang shader, each Gaussian is a struct holding its center, covariance matrix, color, and opacity. The two key methods are eval() — which evaluates the Gaussian's contribution at a UV coordinate — and bounds() — which computes an oriented bounding box for culling.

struct Gaussian2D : IDifferentiable
{
    float2 center;
    float2x2 sigma;   // covariance matrix
    float3 color;
    float opacity;

    [Differentiable]
    float4 eval(float2 uv)
    {
        float2x2 invCov = inverse(sigma);
        float2 diff = uv - center;
        float power = -0.5f * (
            diff.x * diff.x * invCov[0][0] +
            diff.y * diff.y * invCov[1][1] +
            diff.x * diff.y * invCov[0][1] +
            diff.y * diff.x * invCov[1][0]);
        float weight = min(.99f, opacity * exp(power));
        return float4(color, weight);
    }

    OBB bounds()
    {
        float a = sigma[0][0], b = sigma[0][1];
        float c = sigma[1][0], d = sigma[1][1];
        float n_stddev = 4.f;

        if (abs(b) < 1e-6 || abs(c) < 1e-6) {
            float2x2 eigenvectors = float2x2(float2(1,0), float2(0,1));
            float2 scale = float2(sqrt(a), sqrt(d));
            return OBB(center, eigenvectors, scale * n_stddev);
        } else {
            float trace = a + d;
            float det = a * d - b * c;
            float lambda1 = 0.5 * (trace + sqrt(trace*trace - 4*det));
            float lambda2 = 0.5 * (trace - sqrt(trace*trace - 4*det));
            float2x2 eigenvectors;
            eigenvectors[0] = normalize(float2(lambda1 - d, c));
            eigenvectors[1] = normalize(float2(b, lambda2 - a));
            float2 scale = float2(sqrt(lambda1), sqrt(lambda2));
            return OBB(center, eigenvectors, scale * n_stddev);
        }
    }
}

A 2D Gaussian is fully defined by a center (mean) μ and a covariance matrix Σ. The unnormalized value at point p is:

f(p) = exp(-0.5 * (p - μ)ᵀ Σ⁻¹ (p - μ))

The covariance matrix Σ is a 2×2 symmetric matrix:

Σ = [[σ_xx,  σ_xy],
     [σ_yx,  σ_yy]]   where σ_xy = σ_yx

The diagonal entries σ_xx and σ_yy control how spread out the Gaussian is along the x and y axes respectively. The off-diagonal entry σ_xy controls the tilt — a positive value tilts the blob toward 45°, a negative value toward −45°. When σ_xy = 0, the Gaussian is axis-aligned.

The power computation in eval() expands the matrix-vector quadratic form manually into scalar operations:

(p-μ)ᵀ Σ⁻¹ (p-μ) = diff.x² * invCov[0][0]
                   + diff.y² * invCov[1][1]
                   + diff.x * diff.y * invCov[0][1]
                   + diff.y * diff.x * invCov[1][0]

Because Σ is symmetric, its inverse is also symmetric: invCov[0][1] = invCov[1][0]. The last two terms therefore double up the cross contribution — this is equivalent to 2 * diff.x * diff.y * invCov[0][1].

The weight returned by eval() is min(0.99, opacity * exp(power)) — clamped below 1 to avoid fully occluding pixels in a single splat.

Different covariance matrices produce very different Gaussian shapes. A valid covariance matrix must be positive semi-definite — all eigenvalues must be ≥ 0. This ensures the Gaussian's "spread" is always real-valued and non-negative in every direction.

The five archetypes below cover the main shapes you'll encounter. Notice how the off-diagonal entry rotates the principal axes away from the x/y grid:

The eval() method calls inverse(sigma) to get the precision matrix before computing the quadratic form. For a 2×2 matrix the formula is:

inverse([[a, b], [c, d]]) = (1 / (ad - bc)) * [[d, -b], [-c, a]]

The inverse covariance (also called the precision matrix) determines how "tight" or "loose" the Gaussian is:

A large variance (large diagonal in Σ) gives small precision values — the exponent decays slowly, giving a broad, diffuse blob. A small variance gives large precision values — the exponent decays quickly, giving a sharp, narrow spike.

The matrix display below the Demo 1 canvas updates in real time as you move the sliders, letting you see exactly how Σ⁻¹ changes with the parameters.

[Differentiable]
float2x2 inverse(float2x2 mat) {
    float det = determinant(mat);
    return float2x2(
        mat[1][1] / det, -mat[0][1] / det,
       -mat[1][0] / det,  mat[0][0] / det);
}

Because inverse() is marked [Differentiable], the autograd system can backpropagate gradients through it — including through the determinant division. This is what makes it possible to train the covariance parameters directly via gradient descent on a rendering loss.

During rendering, Gaussians are processed tile by tile. Before computing any pixels, we need to know which tiles a Gaussian can possibly affect. A Gaussian with very small weight at some pixel contributes nothing — so we can skip it entirely. The bounds() method computes an Oriented Bounding Box (OBB) that encloses the region where the Gaussian's weight exceeds a meaningful threshold.

The key insight: the Gaussian's shape is exactly described by the eigenvectors and eigenvalues of Σ. The eigenvectors give the principal axes (directions of maximum/minimum spread), and the square roots of the eigenvalues give the standard deviations along those axes.

Eigenvalues of Σ:

trace = σ_xx + σ_yy
det   = σ_xx * σ_yy - σ_xy²

λ₁ = (trace + √(trace² - 4·det)) / 2
λ₂ = (trace - √(trace² - 4·det)) / 2

Eigenvectors give the principal axes (directions of maximum/minimum spread):

v₁ = normalize((λ₁ - σ_yy, σ_yx))
v₂ = normalize((σ_xy, λ₂ - σ_xx))

OBB scale = n_stddev * (√λ₁, √λ₂) where n_stddev = 4 — the box extends 4 standard deviations in each principal direction. At 4σ, exp(-0.5 * 16) ≈ 0.00034, meaning the Gaussian has decayed to well under 0.1% of its peak value. This captures >99.99% of the Gaussian's mass while keeping tiles tiny.

The diagonal shortcut: when |σ_xy| < 1e-6, the covariance is already diagonal — the eigenvectors are exactly the x and y unit vectors, so the OBB aligns with the canvas axes and no eigendecomposition is needed.

During optimisation, Gaussian parameters are stored as raw unconstrained floats in a differentiable buffer. The load() function maps these raw values into the valid parameter ranges using smoothStep — a smooth Hermite interpolation that approaches 0 at x → −∞ and 1 at x → +∞, with zero derivative at both ends. This prevents the optimiser from getting gradient signals that push parameters out of valid ranges.

The mapping pipeline for each parameter:

• center.x / center.y — raw float → smoothStep → [0, 1] UV space
• sigma[0][0] / sigma[1][1] — scaled raw float → smoothStep → small padding (0.005) added to avoid singularities in the matrix inverse
• anisotropy — raw float → smoothStep → mapped to (−0.825, +0.825), then sigma[0][1] = sigma[1][0] = sqrt(sxx * syy) * aniso ensures the matrix stays positive semi-definite
• color.r / g / b — raw floats → smoothStep → [0, 1] per channel
• opacity — raw float → smoothStep → [0, 1]

Gaussian2D load(DiffTensorView buf, int idx) {
    Gaussian2D g;
    // Center: smoothstep maps any real -> [0,1]
    g.center.x = smoothstep(buf[idx * 9 + 0]);
    g.center.y = smoothstep(buf[idx * 9 + 1]);

    // Diagonals: scale raw value, smoothstep, add padding
    g.sigma[0][0] = smoothstep(buf[idx * 9 + 2] * 0.1) * 0.05 + 0.005;
    g.sigma[1][1] = smoothstep(buf[idx * 9 + 3] * 0.1) * 0.05 + 0.005;

    // Anisotropy: smoothstep -> [-0.825, +0.825]
    float aniso = (smoothstep(buf[idx * 9 + 4]) - 0.5) * 1.65;
    g.sigma[0][1] = g.sigma[1][0] =
        sqrt(g.sigma[0][0] * g.sigma[1][1]) * aniso;

    // Color channels
    g.color.r = smoothstep(buf[idx * 9 + 5]);
    g.color.g = smoothstep(buf[idx * 9 + 6]);
    g.color.b = smoothstep(buf[idx * 9 + 7]);

    // Opacity
    g.opacity  = smoothstep(buf[idx * 9 + 8]);
    return g;
}

The interactive visualiser below lets you explore how raw buffer values flow through the smoothStep mapping into final Gaussian parameters, and see the resulting Gaussian in real time.