Power Moments in Moment-Based OIT Using GLSL

This article provides a mathematical breakdown and GLSL implementation of power moments within Moment-Based Order-Independent Transparency (MBOIT). It covers the generation of power moment vectors during additive rendering passes, matrix conditioning to prevent numerical instability, and the mathematical reconstruction of the transmittance function via Cholesky decomposition to resolve proper alpha blending without sorting fragments.

Mathematical Formulation of Power Moments

In order-independent transparency, the goal is to evaluate the net transmittance function \(T(z)\) along a view ray up to depth \(z\):

\[T(z) = \exp\left(-\int_0^z \sigma(t) \, dt\right)\]

where \(\sigma(t)\) represents the continuous extinction or absorption coefficient. In discrete rendering, each transparent fragment \(i\) provides an optical depth (absorbance) \(a_i = -\ln(1.0 - \alpha_i)\) and a normalized depth \(z_i \in [0, 1]\).

MBOIT models the distribution of optical thickness using a moment vector of degree \(N\). The \(k\)-th power moment \(b_k\) over depth distribution is defined as:

\[b_k = \sum_{i} a_i z_i^k, \quad k \in \{0, 1, \dots, N\}\]

For a four-moment configuration (\(N = 4\)), the moment vector \(\mathbf{b} = (b_0, b_1, b_2, b_3, b_4)^T\) requires:

Zeroth moment: \(b_0 = \sum_i a_i\), which represents the total integrated optical thickness along the ray.
Higher moments: \(b_k = \sum_i a_i z_i^k\), capturing the spatial concentration of depth values.

Moment Accumulation Pass in GLSL

During the geometry pass, fragments are rendered without depth writing, and hardware additive blending (GL_SRC_ALPHA, GL_ONE with GL_FUNC_ADD) is enabled. Fragments write their absorbance-weighted power moments into floating-point framebuffer textures (typically two RGBA32F targets).

// Fragment Shader: Moment Generation
layout(location = 0) out vec4 outMoments0; // b0, b1, b2, b3
layout(location = 1) out float outMoments1; // b4

in float v_linearDepth; // Depth normalized to [0.0, 1.0]
uniform float u_alpha;

void main() {
    // Clamp alpha to avoid log(0)
    float alpha = clamp(u_alpha, 0.0, 0.999);
    float a = -log(1.0 - alpha);

    float z = clamp(v_linearDepth, 0.0, 1.0);
    float z2 = z * z;
    float z3 = z2 * z;
    float z4 = z2 * z2;

    outMoments0 = a * vec4(1.0, z, z2, z3);
    outMoments1 = a * z4;
}

Reconstruction and Transmittance Evaluation

During the resolve pass, a fragment at query depth \(z_f\) estimates the fraction of optical depth situated in front of it:

\[S(z_f) = \frac{1}{b_0} \int_0^{z_f} \sigma(t) \, dt \in [0, 1]\]

The transmittance up to depth \(z_f\) is reconstructed as:

\[T(z_f) = \exp(-b_0 \cdot S(z_f))\]

To approximate \(S(z_f)\), the normalized moments \(\mu_k = \frac{b_k}{b_0}\) define a Hankel matrix:

\[\mathbf{M} = \begin{pmatrix} 1 & \mu_1 & \mu_2 \\ \mu_1 & \mu_2 & \mu_3 \\ \mu_2 & \mu_3 & \mu_4 \end{pmatrix}\]

Matrix Regularization and Cholesky Factorization

Floating-point precision limits and steep depth gradients can cause \(\mathbf{M}\) to become non-positive definite. A regularization parameter \(\epsilon\) (typically \(10^{-7}\) to \(10^{-5}\)) is blended with an identity-like bias:

\[\mathbf{M}_{\text{reg}} = (1 - \epsilon)\mathbf{M} + \epsilon \mathbf{I}\]

Performing an \(L D L^T\) or standard Cholesky decomposition \(\mathbf{M}_{\text{reg}} = \mathbf{L} \mathbf{L}^T\) yields a lower triangular matrix \(\mathbf{L}\):

\[\mathbf{L} = \begin{pmatrix} l_{00} & 0 & 0 \\ l_{10} & l_{11} & 0 \\ l_{20} & l_{21} & l_{22} \end{pmatrix}\]

The values of \(\mathbf{L}\) are computed directly:

\(l_{00} = 1.0\)
\(l_{10} = \mu_1\)
\(l_{11} = \sqrt{\mu_2 - l_{10}^2}\)
\(l_{20} = \mu_2\)
\(l_{21} = \frac{\mu_3 - l_{20} l_{10}}{l_{11}}\)
\(l_{22} = \sqrt{\mu_4 - l_{20}^2 - l_{21}^2}\)

Evaluating the Cumulative Distribution

Using the decomposition, the vector of monomial powers at query depth \(\mathbf{q}(z_f) = (1, z_f, z_f^2)^T\) is transformed by forward substitution into vector \(\mathbf{c} = \mathbf{L}^{-1} \mathbf{q}(z_f)\):

\[c_0 = 1.0\] \[c_1 = \frac{z_f - l_{10}}{l_{11}}\] \[c_2 = \frac{z_f^2 - l_{20} - l_{21} c_1}{l_{22}}\]

The upper bound of the normalized cumulative optical depth \(S(z_f)\) is computed via the Hamburger moment bound:

\[S(z_f) \approx \frac{\sum_{j=0}^{m} c_j^2 \cdot H(z_f - z_j)}{\sum_{j=0}^{m} c_j^2}\]

In practical four-moment MBOIT implementations, this simplifies to evaluating the quadratic polynomial derived from the inverted Cholesky system to compute the step distribution:

// GLSL evaluation function for transmittance
float evaluateTransmittance(vec4 b0_to_b3, float b4, float z_f) {
    float b0 = b0_to_b3.x;
    if (b0 <= 1e-6) return 1.0; // No occlusion

    // Normalized moments
    vec4 mu = vec4(b0_to_b3.yzw, b4) / b0;

    // Regularization
    float overestimation = 0.1;
    float reg = 1e-5;
    mu = mix(mu, vec4(0.5, 1.0/3.0, 0.25, 0.2), reg);

    // Cholesky factors
    float l00 = 1.0;
    float l10 = mu.x;
    float l11 = sqrt(max(1e-6, mu.y - l10 * l10));
    float l20 = mu.y;
    float l21 = (mu.z - l20 * l10) / l11;
    float l22 = sqrt(max(1e-6, mu.w - l20 * l20 - l21 * l21));

    // Solve L * c = q
    float c0 = 1.0;
    float c1 = (z_f - l10) / l11;
    float c2 = (z_f * z_f - l20 - l21 * c1) / l22;

    // Estimate step function S(z)
    float sumSq = c0 * c0 + c1 * c1 + c2 * c2;
    float S = clamp(1.0 - (1.0 / sumSq), 0.0, 1.0);

    // Optical depth integration
    return exp(-b0 * S);
}

This mathematical evaluation enables smooth, artifact-free transmittance estimation across arbitrary overlapping geometry without relying on runtime per-fragment sorting.