This text brings a number of concepts related to PRT together, trying to turn them into a straight A-to-Z derivation of the basic technique (diffuse, no bounce, no scattering). ”Straight” here does not mean simple; this is not really a tutorial. Depending on how much you already know about the subject, this is a math cheat sheet, a quick refresher, or a starting point to learn what's involved.

Note: for worry-free reading and printing, you can also use this pdf file which is exactly the same thing. The html below is autogenerated so it may or may not look very good in your browser.



Let's start with the rendering equation, which tells us that radiance leaving a point in some direction is what that point is emitting, plus any incoming radiance reflected in that direction:



We can simplify this if we assume that all the light is coming from some external environment (ie. a sky or environment map), meaning L_e is zero everywhere. Further, a diffuse-only scene simplifies the BRDF to just a (normalized) reflectivity:





The star over theta indicates that the outgoing direction doesn’t matter, since the reflection is diffuse; but at the same time it's important to keep in mind that the radiance does have a direction, since it's radiance and not radiosity.

Simplifying to a lightsource that's infinitely far away, and looking at direct light only for now (ie. no indirect bounces or other transport), means that the incoming radiance is independent of the location x. The indirect radiance at x simply becomes the radiance arriving from the direction Ψ, provided that x can see the infinitely far away light in that direction:




At this point, the incoming radiance is approximated by spherical harmonics; the function L is rewritten as a weighted sum of these basis functions yi:



where i = l(l + 1) + m and n = order². Substituting and shuffling a bit gives:



It follows that if we can compute that integral offline for every point x, then we can recover an approximate radiance leaving x by calculating the dot product with those weights, and scaling.

So how do we compute the integral? Via Monte Carlo sampling, which says that the Expected Value E of a function can be approximated by averaging a bunch of samples, taken according to the probability distribution function (pdf) used in evaluating E:



It follows that if we estimate the expected value of f∕p, instead of f, we end up estimating the integral of f, rather than its expected value:



This is known as importance sampling. The integral we want in equation 8 is the product of two parts: the spherical harmonic y, and the cosine weighted visibility term. We could importance sample the harmonic, but then we would have to evaluate the visibility term for every sample, for every harmonic. Since the visibility term is the expensive part (it's where the raytracing happens), the other way around is more efficient: importance sample the cosine part, and use the sample to update the estimate for every harmonic in parallel.

Sampling according to the cosine lobe requires a pdf that ideally follows these values (unnormalized for now):



The sine term appears to compensate for the fact that an equal jump along φ, will cover less distance as θ gets closer to azimuth, resulting in a smaller solid angle.

To normalize this pdf properly, and also to efficiently generate samples accordingly, we need the cumulative distribution function (cdf) P :



We substitute the full range to figure out the normalization factor for p:



Dividing p by this, and in turn dividing f by p, means that we simply multiply f by π, eliminating the divide-by-pi coming out of normalizing the diffuse BRDF.

To generate samples according to p∕π, invert the separable P :



where the ζ’s are uniform numbers chosen via stratified sampling, low discrepancy sampling, etc.

Now remember that the point of choosing samples via p was so that we can divide those terms away in the monte carlo estimator:



So at runtime, the radiance becomes:



Now that we have the transfer coefficients t, how do we figure out the weights w for a given lighting setup? Spherical harmonic functions are orthonormal, and real valued. The former means we can process them one by one, and not worry about normalization. The latter means they are their own complex conjugate. The result is we simply integrate the light with the harmonic:



Some lights have been defined as a sum of spherical harmonics straight out
of the box, for example sky models.

References:

• Advanced Global Illumination; Philip Dutré, Kavita Bala, Philippe Bekaert


• Principles of Digital Image Synthesis; Andrew Glassner