I’ve had an idea floating around in my head for several months now, but evening classes, a hard drive failure, then a GPU failure prevented me from doing much about it until this weekend. First, a couple of screenshots of what I’ll be talking about.

Cornell box, I choose you!

The above screenshots are of real-time single-bounce GI in a static scene with fully dynamic lighting. There are 7182 patches, and the lighting calculations take 36ms per frame on one thread of an Intel Core 2 @2.4 GHz. The code is not optimized.

The basic idea is simple and is split into two phases.

A one-time scene compilation phase:

  • Split all surfaces in the scene into roughly equal sized patches.
  • For each patch i, build a list of all other patches visible from i along with the form factor between patches i and j:
    F_{ij} = \frac{cos \phi_i \, cos \phi_j}{\pi \, |r|^2} \, A_j
    Fij is the form factor between patches i and j
    Aj is the area of patch j
    r is the vector from i to j
    Φi is the angle between r and the normal of patch i
    Φj is the angle between r and the normal of patch j

And a per-frame lighting phase:

  • For each patch i, calculate the direct illumination.
  • For each patch i, calculate single-bounce indirect illumination from all visible patches:
  • I_i = \sum_j D_j \, F_{ij}
    Ij is the single-bounce indirect illumination for patch i
    Dj is the direct illumination for patch j

So far, so radiosity. If I understand Michal Iwanicki’s GDC presentation correctly, this is similar to the lighting tech on Milo and Kate, only they additionally project the bounce lighting into SH.

The problem with this approach is that the running time is O(N2) with the number of patches. We could work around this by making the patches quite large, running on the GPU, or both. Alternatively, we can bring the running time down to O(N.log(N)) by borrowing from Michael Bunnell’s work on dynamic AO and cluster patches into a hierarchy. I chose to perform bottom-up patch clustering similarly to the method that Miles Macklin describes in his (Almost) realtime GI blog post.

Scene compilation is now:

  • Split all surfaces in the scene into roughly equal sized patches.
  • Build a hierarchy of patches using k-means clustering.
  • For each patch i, build a list of all other patches visible from i along with the form factor between patches i and j. If a visible patch j is too far from patch i look further up the hierarchy.

And the lighting phase:

  • For each leaf patch i in the hierarchy, calculate the direct illumination.
  • Propagate the direct lighting up the hierarchy.
  • For each patch i, calculate single-bounce indirect illumination from all visible patches clusters.

Although this technique is really simple, it supports a feature set similar to that of Enlighten:

  • Global illumination reacts to changes in surface albedo.
  • Static area light sources that cast soft shadows.

That’s basically about it. There are a few of other areas I’m tempted to look into once I’ve cleaned the code up a bit:

  • Calculate directly and indirect illumination at different frequencies. This would allow scaling to much larger scenes.
  • Perform the last two lighting steps multiple times to approximate more light bounces.
  • Project the indirect illumination into SH, HL2 or the Half-Life basis.
  • Light probes for dynamic objects.

You can grab the source code from here. Expect a mess, since it’s a C++ port of a C# proof of concept with liberal use of vector and hash_map. Scene construction is particularly bad and may take a minute to run. You can fly around using WASD and left-click dragging with the mouse.

Improved light attenuation

February 10, 2011

In my previous post, I talked about the attenuation curve for a spherical light source:
f_{att} = \frac{1}{(\frac{d}{r} + 1)^2}

I had suggested applying a scale and bias to the result in order to limit a light’s influence, which is a serviceable solution, but far from ideal. Unfortunately, applying such a bias causes the gradient of the curve to become non-zero at limit of the light’s influence.

Here’s the attenuation curve for a light of radius 1.0:

And after applying a scale and bias (shown in red):

You can see that the gradient at the zero-crossing is close to, but not quite zero. This is problematic because the human eye is irritatingly sensitive to discontinuities in illumination gradients and we might easily end up with Mach bands.

I was discussing this problem with a colleague of mine, Jerome Scholler, and he came up with an excellent suggestion – to transform d in the attenuation equation by some function whose value tends to infinity as its input reaches our desired maximum distance of influence. My first thought was of using tan:
d' = 2.1tan(\frac{\pi d}{2d_{max}})
f_{att} = \frac{1}{(\frac{d'}{r} + 1)^2}

That worked well, the resulting curve has roughly the same shape as the original, while also having both a gradient and value of zero at the desired maximum distance. It does have the disadvantage of using a trig function, which isn’t so hot, so we went looking for something else. After a few minutes playing around we came up with the following rational function:
d' = \frac{d}{1-(\frac{d}{d_{max}})^2}
f_{att} = \frac{1}{(\frac{d'}{r} + 1)^2}

It’s very similar to the tan version, but may run faster, depending on your hardware.

Below are some examples of the different methods, using a light with high intensity and small influence. On the left of each is the original image, on the right is the result of a levels adjustment, which emphasizes the tail of the attenuation curve.

Disclaimer: The parameters for the analytic functions were chosen to highlight their different characteristics, not to look good.

Original ray traced reference:

Analytic scale and bias:

Analytic tan:

Analytic rational:

The graphs today were brought to you courtesy of the awesome fooplot.com.

Light attenuation

January 31, 2011

The canonical equation for point light attenuation goes something like this:
f_{att} = \frac{1}{k_c + k_ld + k_qd^2}
d = distance between the light and the surface being shaded
kc = constant attenuation factor
kl = linear attenuation factor
kq = quadratic attenuation factor

Since I first read about light attenuation in the Red Book I’ve often wondered where this equation came from and what values should actually be used for the attenuation factors, but I could never find a satisfactory explanation. Pretty much every reference to light attenuation in both books and online simply presents some variant of this equation, along with screenshots of objects being lit by lights with different attenuation factors. If you’re lucky, there’s sometimes an accompanying bit of handwaving.

Today, I did some experimentation with my path tracer and was pleasantly surprised to find a correlation between the direct illumination from a physically based spherical area light source and the point light attenuation equation.

I set up a simple scene in which to conduct the tests: a spherical area light above a diffuse plane. By setting the light’s radius and distance above the plane to different values and then sampling the direct illumination at a point on the plane directly below the light, I built up a table of attenuation values. Here’s a plot of a some of the results; the distance on the horizontal axis is that between the plane and the light’s surface, not its centre.

After looking at the results from a series of tests, it became apparent that the attenuation of a spherical light can be modeled as:
f_{att} = \frac{1}{(\frac{d}{r} + 1)^2}
d = distance between the light’s surface and the point being shaded
r = the light’s radius

Expanding this out, we get:
f_{att} = \frac{1}{1 + \frac{2}{r}d + \frac{1}{r^2}d^2}
which is the original point light attenuation equation with the following attenuation factors:
k_c = 1
k_l = \frac{2}{r}
k_q = \frac{1}{r^2}

Below are a couple of renders of four lights above a plane. The first is a ground-truth render of direct illumination calculated using Monte Carlo integration:

In this second render, direct illumination is calculated analytically using the attenuation factors derived from the light radius:

The only noticeable difference between the two is that in the second image, an area of the plane to the far left is slightly too bright due to a lack of a shadowing term.

Maybe this is old news to many people, but I was pretty happy to find out that an equation that had seemed fairly arbitrary to me for so many years actually had some physical motivation behind it. I don’t really understand why this relationship is never pointed out, not even in Foley and van Dam’s venerable tome*.

Unfortunately this attenuation model is still problematic for real-time rendering, since a light’s influence is essentially unbounded. We can, however, artificially enforce a finite influence by clipping all contributions that fall below a certain threshold. Given a spherical light of radius r and intensity Li, the illumination I at distance d is:
I = \frac{L_i}{(\frac{d}{r} + 1)^2}

Assuming we want to ignore all illumination that falls below some cutoff threshold Ic, we can solve for d to find the maximum distance of the light’s influence:
d_{max} = r(\sqrt{\frac{L_i}{I_c}}-1)

Biasing the calculated illumination by -Ic and then scaling by 1/(1-Ic) ensures that illumination drops to zero at the furthest extent, and the maximum illumination is unchanged.

Here’s the result of applying these changes with a cutoff threshold of 0.001; in the second image, areas which receive no illumination are highlighted in red:

And here’s a cutoff threshold of 0.005; if you compare to the version with no cutoff, you’ll see that the illumination is now noticeably darker:

Just to round things off, here’s a GLSL snippet for calculating the approximate direct illumination from a spherical light source. Soft shadows are left as an exercise for the reader.

vec3 DirectIllumination(vec3 P, vec3 N, vec3 lightCentre, float lightRadius, vec3 lightColour, float cutoff)
    // calculate normalized light vector and distance to sphere light surface
    float r = lightRadius;
    vec3 L = lightCentre - P;
    float distance = length(L);
    float d = max(distance - r, 0);
    L /= distance;
    // calculate basic attenuation
    float denom = d/r + 1;
    float attenuation = 1 / (denom*denom);
    // scale and bias attenuation such that:
    //   attenuation == 0 at extent of max influence
    //   attenuation == 1 when d == 0
    attenuation = (attenuation - cutoff) / (1 - cutoff);
    attenuation = max(attenuation, 0);
    float dot = max(dot(L, N), 0);
    return lightColour * dot * attenuation;

* I always felt a little sorry for Feiner and Hughes.