Stochastic Light Culling in Vulkan

Contents

• Introduction
• Problem Statement
• Stochastic Fall-off Function
• Probability Function
• Light Influence Range
• Variance Control
• Constructing Equation for \(\alpha_i\)
• Implementation
• References

Introduction

Consider a scene with \(N\) light sources: \(L_0\), \(L_1\), …, \(L_{N-1}\). Let \(x\) be a shading point, and \(\hat{w}\) be the view direction. The radiance \(L(x, \hat{w})\) is defined as follows

where:

• \(x_i\) is the position of the light \(L_i\)
• \(\hat{n}\) is the surface normal at \(x\)
• \(\lVert x_i - x \rVert\) is the distance between \(x_i\) and \(x\)
• \(I_i(-\hat{w}_{i}^{\prime})\) is the radiant intensity of \(L_i\)
• \(ρ(x, \hat{w}, \hat{w}_{i}^{\prime})\) is the bidirectional reflectance distribution function (BRDF)
• \(V(x, x_i)\) is the visibility between \(x\) and \(x_i\) to represent shadow. It is a measure of how occluded \(x\) is from the light sources.
• \(f(\lVert x_i - x \rVert)\) is the fall-off function. It measures the attenuation of the light as it travels away from the light source. In other words, it models how the contribution of the light source decreases, as the distance between the shading point and the light increases.
• \(\hat{w}_{i}^{\prime}\) is the direction of the light \(L_i\), where

In a traditional rendering pipeline, also known as the Forward Rendering Pipeline, equation \(1\) is executed by the fragment shader to compute the radiance at each shading point. The problem with this setup, however, is that it does not scale well with the number of lights. If \(N\) grows, the execution time of the fragment shader increases, and the pipeline will fail to maintain adequate frame-rate, and hence why, the idea of light culling was introduced.

Light culling is an optimization technique often used in rendering pipelines. The main idea behind this method is to simply reduce the number of light sources being evaluated for each shading point. Instead of evaluating \(N\) light sources, we evaluate only a subset of them. This is done by assuming a limited range of influence \(r_i\) for each light source \(L_i\). If a shading point falls outside that influence range, the light \(L_i\) is ignored when computing the radiance at \(x\).

A high-level overview of the light culling algorithm is as follows:

The screen is split into a 2D grid of tiles.

A frustum is constructed for each tile.

A bounding sphere is built for each light source. The center of the sphere is the light position, and the radius of the sphere is the influence range of the light.

Intersection is tested between all the light bounding spheres and all tile frustums.

If there is an intersection, the index of the light source is appended to the tile's private list of lights.

Pixels within this tile will only evaluate the light sources whose indices are in this tile's list

However, the assumption of limited influence range of each light made by the light culling is physically inaccurate. Physical-based lights have infinite range of influence. Regardless of the distance \(d\) between the shading point \(x\) and the light source \(L_i\), \(x\) should still receive some radiance contribution from \(L_i\). This is typically modelled using the fall-off function \(f(d)\)

The physical inaccuracy of light culling is apparent in the scenes it produces. Rendering pipelines that use light culling produce scenes that are darker than the original. This is especially true if the influence range \(r_i\) of the light is reduced aggressively. However, despite the physical inaccuracy it introduces, light culling is still commonly used in rendering pipelines. This is due to the performance benefits it provides for scenes with large number of light sources. Although we always aim to produce realistic scenes in computer graphics, we tend to favor performance over realism. A scene with high frame-rate that is not physically accurate is preferred over a scene that is physically accurate, but non-interactive or has slow frame-rate. That’s why light culling is still used nowadays in rendering pipelines despite its shortcomings.

Now, the question is:

Is there a way we can "fix" these darkened bits in the scenes caused by light culling, while at the same time maintain the performance benefits it provides ?

The answer is:

Yes
The method is explained in the 2016 paper "Stochastic Light Culling" by Y.Tokuyoshi, T.Harada.

In this blog post, we will explain the method, derive its equations and formulas, and provide an implementation in C++.

Problem Statement

The truncated fall-off function \(f^{\prime}(d)\) used by light culling is defined as follows:

where \(r_i\) is the influence range of the light \(L_i\).

As shown, if the shading point \(x\) is within the influence range \(r_i\), the function \(f^{\prime}(d)\) applies the physically-based fall-off function \(f(d)\). However, if the shading point \(x\) is beyond the range \(r_i\), the light \(L_i\) is completely ignored. This explains why scenes produced by light culling has darkened regions.

In fact, if we computed the bias \(E\) of \(f^{\prime}(d)\), we will find it biased:

In other words, inside the cutoff radius \(r_i\), the fall-off function \(f^{\prime}(d)\) is unbiased. Otherwise, the estimator underestimates the light contribution significantly by setting it to zero. This lines up with the results we are seeing. If \(f^{\prime}(d)\) was truly unbiased, the scene would have looked identical to the one produced by the traditional rendering pipeline, which is not the case. This is the problem we are trying to solve.

The solution the paper is introducing is to replace the truncated fall-off function \(f^{\prime}(d)\) with an unbiased stochastic fall-off function \(\hat{f}(d)\). This way, we get to produce scenes that are almost identical to the physically-based one, yet keep the performance benefits of light culling.

Stochastic Fall-off Function

Let \(p_i(d)\) be a probability function, and \(ξ_i \in [0, 1)\) be a uniform random number for each light source \(L_i\). The stochastic fall-off function \(\hat{f}(d)\) approximated using Russian roulette is given by

Let’s take a step back and explain this function more clearly, and how exactly we derived it.

Consider, for example, the physically-based fall-off function \(f(d)\) defined in equation \(\). Every light \(L_i\) is contributing exactly the amount of \(\frac{1}{d_i^2}\) to the shading point \(x\). With the introduction of the truncated fall-off function \(f^{\prime}(d)\) in light culling, only a subset of the lights is contributing the amount of \(\frac{1}{d_i^2}\), decreasing the total contribution that the shading point \(x\) should physically receive. The key idea of the stochastic fall-off function is to artifically increase the radiance contribution provided by the surviving lights to account for the contribution that would have been provided by the culled lights, keeping the total contribution that the shading point \(x\) is receiving practically the same.

So, we define the probability function \(p_i(d) \in [0, 1]\) as the probability that light \(L_i\) is evaluated at the shading point \(x\). For every light \(L_i\), we flip a weighted coin \(ξ_i\). If \(p_i(d) > ξ_i\), we evaluate the light \(L_i\), and artifically inflate its contribution by a factor of \(\frac{1}{p_i(d)}\). Otherwise, we ignore it.

Suppose for example that \(f(d) = 40\) for the light \(L_i\), and \(p_i(d) = 0.25\). This means that \(25\%\) of the time, \(L_i\) survives, and the other \(75\%\) of the time, it is ignored. The average contribution of light \(L_i\) is

What we do then here is that we inflate the radiance contribution of light \(L_i\) in the times it is evaluated with a factor of \(\frac{1}{p_i(d)}\), such that

The average contribution becomes equal to the physically-accurate value of \(f(d) = 40\).

The stochastic fall-off function \(\hat{f}(d)\) will be used in the main shader of the rendering pipeline to compute the attenutation of the light intensity, in place of the truncated fall-off function \(f^{\prime}(d)\) used by light culling.

Probability Function

In the previous section, we defined the stochastic fall-off function \(\hat{f}(d)\), but we did not define the probability function \(p_i(d)\). This is because it did not matter. Regardless of how \(p_i(d)\) is defined, \(\hat{f}(d)\) is still unbiased. The definition of \(\hat{f}(d)\) is not dependent on \(p_i(d)\) at all.

Some probability functions, however, are better than other in this context. For example, we could set

With this choice of \(p_i(d)\), the stochastic fall-off function \(\hat{f}(d)\) is still unbiased. However, we would be evaluating 50% of all light sources for each shading point \(x\), regardless of the distance between them. This does not align with our light culling pipeline, since there is no relationship between \(p_i(d)\) and the distance \(d\) from the shading point \(x\).

We need to define a probability function \(p_i(d)\) that is proportional to \(f(d)\);

This way, lights that are far from the shading point \(x\) gets low probability of survival, and lights close to \(x\) gets high probability of survival, and are evaluated.

The paper defines the probability function \(p_i(d)\) as follows:

where \(\alpha_i\) is a constant value to control variance.

Since \(\frac{f(d)}{p_i{d}} = max(\alpha_i, f(d))\), we can update \(\hat{f}(d)\) as follows:

Light's Influence Range

According to the stochastic fall-off function \(\hat{f}(d)\), if the probability \(p_i(d) > ξ_i\), the light \(L_i\) is ignored. We need to compute the influence range \(r_i\) of the light \(L_i\), such that it is set at the point in which the light contribution becomes zero. This point is at the distance \(d\) at which

That is,

Since \(ξ \in [0, 1)\), then \(ξ < 1\). This means that

Variance Control

One of the variables we have not discussed till now is \(\alpha_i\). We briefly mentioned that it controls the variance \(σ^2(\hat{f}(d))\) of the stochastic fall-off function \(\hat{f}(d)\), but in this section, we will explain how. The equation for computing the variance \(σ^2(X)\) of an estimator \(X\) is defined as follows:

We will start by computing the expected value \(E\left[\hat{f}(d)\right]\) of the estimator \(\hat{f}(d)\).

Next, we will define the estimator \(\hat{f}^2(d)\) as follows:

The expected value \(E\) of \(\hat{f}^2(d)\); \(E\left[\hat{f}^2(d)\right]\) is defined as follows:

Hence, the equation for computing the varince \(σ^2(\hat{f}(d))\) of the estimator \(\hat{f}(d)\), is defined as follows:

By substituting \(p_i(d) = \frac{f(d)}{\alpha_i}\), we get

As shown in equation \(7\), \(\alpha_i\) is proportional to the variance \(σ^2(\hat{f}(d))\). This means that if \(\alpha_i\) increases, the variance \(σ^2(\hat{f}(d))\) also increases, and vice versa. This realization lines up accurately with the framework we discussed so far. Suppose, for instance, we increased \(\alpha_i\). This will cause \(p_i(d)\) to decrease, and, in turn, more lights will fail the condition of \(p_i(d) > ξ_i\) and get ignored. The influence range \(r_i\) of the lights will also decrease causing more aggressive culling by the pipeline. Moreover, the radiance contribution of the surviving lights; namely,

will also increase. This means that the values that the estimator \(\hat{f}(d)\) is producing are either zero, or very high values, leading to the high variance. Conversely, if we decreased \(\alpha_i\), the number of culled lights will decrease, leading to low variance, and in turn, better rendering quality, but at the expense of worse performance.

Regardless of the level of variance, the estimator is still unbiased. The large fluctuations of data in estimators with high variance, however, are not always favourable to have. Suppose, for example, we evaluate the stochastic fall-off \(\hat{f}(d)\) every frame with a newly generated random number. There will be large data variations between every frame, causing noise in the rendered scene.

Luckily, in our setup, we generate a random number only once at the begining, evaluate \(\hat{f}(d)\) then, and re-use the same data every frame. As a result, our rendered scene will not have any flickering. However, there will be big discrepancy between each program run. One program run could have decent rendering quality, and the next could be awful. The key idea here is to determine how to pick the value of \(\alpha_i\) that balances the competing factors of rendering quality and performance. The paper defines an equation that calculates \(\alpha_i\) based on a maximum error value \(\epsilon_{max}\) to be tolerated by the user.

Constructing Equation for \(\alpha_i\)

Recall equation \(1\). We will replace the fall-off function \(f(\lVert x_i - x \rVert)\) with our stochastic fall-off function \(\hat{f}(\lVert x_i - x \rVert)\). Namely,

Next, let’s reduce equation \(8\) by only considering the radiance contribution of one light source:

We want to compute the radiance error \(\epsilon_i(x, \hat{w})\) inflicted by the term \(\hat{f}(\lVert x_i - x \rVert)\). To do this, we will replace the term \(\hat{f}(\lVert x_i - x \rVert)\) with the standard deviation \(σ(\lVert x_i - x \rVert)\), which is a measure of how far the value produced by the estimator \(\hat{f}(\lVert x_i - x \rVert)\) is from the expected value \(E\left[\hat{f}(d)\right]\). That is,

We will compute the error bound \(\epsilon_{max}\) by assuming the maximum value for each term of equation \(9\). That is,

• The maximum value for the visibility variable \(V(x, x_i)\) is \(1\).
• The maximum value for the dot product \((\hat{w}_{i}^{\prime} \cdot \hat{n})\) is \(1\).
• We will assume the usage of the Lambertian model for the BRDF function \(ρ(x, \hat{w}, \hat{w}_{i}^{\prime})\). Hence, \(ρ(x, \hat{w}, \hat{w}_{i}^{\prime}) = \frac{1}{\pi}\)
• We will compute the maximum possible value for the standard deviation \(σ(\hat{f}(d))\); \(σ_{max}(\hat{f}(d))\) by, first, computing the maximum value for the variance \(σ_{max}^2(\hat{f}(d))\), and then applying the equation \(σ(X) = \sqrt(σ^2)\). The maximum value \(σ_{max}(\hat{f}(d))\) can be computed by calculating the derivative of equation \(7\) and setting it to \(0\). Namely,

By substituting the maximum values in equation \(9\), the resultant equation for the error bound \(\epsilon_{max}\) is defined as follows:

Therefore, the equation for computing \(\alpha_i\) is given as:

Implementation

The following is a function that computes the uniform random number \(\epsilon_i\) that is generated for each light source \(L_i\):

float GenerateEpsilon()
{
    std::mt19937 rng(std::random_device{}());
    std::uniform_real_distribution<float> dist01(0.0f, 1.0f);
    
    return dist01(rng);
}

The following is a function that computes \(\alpha_i\) using equation \(11\), where \(errorMax\) is a user-define error bound.

float ComputeAlpha(float errorMax, glm::vec3 lightColor, float lightIntensity)
{
    glm::vec3 I = lightIntensity * lightColor;
    float maxLightIntensityComponent = std::max({I.x, I.y, I.z});
    
    return ((2.0f * PI * errorMax) / (maxLightIntensityComponent));
}

The following is a function that computes the bounding sphere of the light source. The center of the sphere is the position of the light, and the radius is computed using equation \(6\):

glm::vec4 ComputeBoundingSphere(float alpha, float epsilon)
{
    glm::vec3 center = light.config.position;
    float radius = (1.0f / (sqrt(alpha * epsilon)));
    
    return glm::vec4(center, radius);
}

In the main shader, the attenutation of the light source is computed using the stochastic fall-off function \(\hat{f(d)}\) defined in equation \(5\):

float ComputeAttenutation(float dist, float alpha)
{
    float f = 1.0 / (dist * dist);
        
    return max(alpha, f);
}

References

[1] \(\) \(\) Y.Tokuyoshi, T.Harada. Stochastic Light Culling. 2016.