The Grand Picture
Most introductions and implementations of Precomputed Radiance Transfer will deal fairly well with the easiest use case: The double product integral for diffuse reflections. Beyond this scope however, few show how proper rotation of frequency-space encoded lighting works, and even fewer dive into the problem of view-dependency in PRT. Both of these are related, as they rely on the ability to transform one set of coefficents into another: Rotating a vector creates another vector, and view dependent reflection transforms incident light represented as coefficients into coefficients of reflected light. From linear algebra, we know how to deal with such a scenario already: To transform one vector into another we need a matrix, so in terms of a lighting scenario we need to talk about Matrix Radiance Transfer and the Triple Product Integral.
Function transforms
In the last post on function transforms I took a quick look at the convolution theorem, which one can roughly describe as the ability to shortcut an integration over the product of two functions, i.e., a convolution, with a dot product in the frequency domain. Typically, the convolution theorem is introduced with the Fourier basis, but other function bases can be used as well.
This holds true for any function basis that is orthonormal. To test whether a set of functions forms an orthonormal basis, one needs to do two things: Make sure that any two basis functions and integrate to 0 if , and to 1 if .
This works completely analogous to a vector basis.
The Triple Product Integral
So lets assume we have two functions and and we want to compute the function product of both as a new function .
Assume further that we have projected and into the frequency domain of the function basis . This means that we have the coefficients for both at hand.
It may be possible to do the same operation in Equation - multiplying two functions together - in the frequency space of the basis that we chose. How would we be able to determine the coefficients of ? We can start the usual way by integrating it with the basis function.
But since we already have the coefficients for both and , we can first replace by Equation and then replace further with both Equation and .
We can rearrange some things: All the terms that are not dependent on can be moved out of the integral.
Now that is an interesting looking integral! Three basis functions in one, and at first sight the whole process looks quite familiar too. It is almost identical to what we did to derive the convolution theorem, now we only need to get rid of this integral and we’re done.
The term is called the triple product integral for more or less obvious reasons, and is a tripling coefficient. Its smaller sibling, the double product integral, is reducable to a Kroenecker Delta if the basis functions form an orthonormal basis. But for , things are not so easy. In fact, there is no general analytical solution for any arbitrary function basis.
Computing tripling coefficients
Luckily, for the two most common bases used for PRT, Haar-Wavelets and Spherical Harmonics, an analytical formula to compute their tripling coefficents exists.
Tripling coefficients of Spherical Harmonics can be expressed through a relationship with Clebsch-Gordan coefficients, which in turn can be expressed by Wigner 3j symbols. Below is an older implementation of mine, but as far as I’m aware the GNU Scientific Library contains one as well.
template<typename T>
inline T wigner_3j(int j1, int j2, int j3,
int m1, int m2, int m3)
{
assert(std::abs(m1) <= j1 &&
"wigner_3j: m1 is out of bounds");
assert(std::abs(m2) <= j2 &&
"wigner_3j: m2 is out of bounds");
assert(std::abs(m3) <= j3 &&
"wigner_3j: m3 is out of bounds");
if (!triangle_inequality(j1, j2, j3)) return 0;
if (m1+m2+m3 != 0) return 0;
if (std::abs(m1) > j1) return 0;
if (std::abs(m2) > j2) return 0;
if (std::abs(m3) > j3) return 0;
int t1 = j2 - m1 - j3;
int t2 = j1 + m2 - j3;
int t3 = j1 + j2 - j3;
int t4 = j1 - m1;
int t5 = j2 + m2;
int tmin = std::max(0, std::max(t1, t2));
int tmax = std::min(t3, std::min(t4, t5));
if (tmin > tmax)
return 0;
T wigner = 0;
auto f = factorial<T>;
for (int t = tmin; t <= tmax; ++t)
wigner = wigner + static_cast<T>(std::pow(-1.0, t))/
(f<T>(t) * f<T>(t-t1) * f<T>(t-t2) *
f<T>(t3-t) * f<T>(t4-t) * f<T>(t5-t));
wigner = wigner * static_cast<T>(std::pow(-1.0, j1-j2-m3)) *
static_cast<T>(std::sqrt(
f<T>(j1+j2-j3) * f<T>(j1-j2+j3) *
f<T>(-j1+j2+j3) / f<T>(j1+j2+j3+1) *
f<T>(j1+m1) * f<T>(j1-m1) *
f<T>(j2+m2) * f<T>(j2-m2) *
f<T>(j3+m3) * f<T>(j3-m3)
));
return wigner;
}With the tripling coefficients at hand, we can create a function product of aribtrary functions encoded in the same frequency domain. A major downside of regular PRT is that the BRDF, visibility and Lambert factor are all represented simply as one transfer function, which therefore gets encoded as one transfer vector. The Lambert factor however is a very low-frequency signal, whilst visibility and BRDF may contain high-frequency peaks. Packing them all together means choosing either to produce a lot of waste coefficients or compromising on the quality of the representation. Ng et al. get around this by decoupling visibility and the BRDF in the transfer function again with a triple product composition.
Another use case presented in Precomputed Shadow Fields for Dynamic Scenes wraps an object in a shell which, at several sample points, has coefficients encoding a visibility function from the point on that shell. If the shell collides with another PRT-object, coefficients of both objects visibility can be combined to compute dynamic shadowing between the two, getting rid of some of the rigidness requirements that PRT imposes on the scene. This paper is naturally enhanced in Precomputed Radiance Transfer Field for Rendering Interreflections in Dynamic Scenes, where the authors not just encode visibility, but indirect diffuse transfer between objects as well.
Transforming coefficients
Returning to Equation we can easily construct in another way.
Instead of using the tripling coefficient tensor and , we can also use a matrix (one way to construct it would be to multiply and ). This matrix will transform coefficients into coefficients , which represent the function in the basis that we chose.
The concept of a matrix to transform a coefficient vector is sometimes mentioned in PRT tutorials for signal rotation. Rather than having an environment map , rotate it with and then resample and recompute the Spherical Harmonic coefficients for it, we can likewise rotate the coefficents with a special frequency-space rotation matrix . In essence, we would produce a new set of coefficients that simply match those of the rotated environment map in image space.
Here is an old implementation of mine of one such matrix to do the job. The input is a regular 3x3 rotation matrix, but the result is a matrix that needs to be as big as the coefficient vector you want to rotate. In case of 4-band Spherical Harmonic vector for instance, the result will be a 16x16 Spherical Harmonic rotation matrix.
Matrix Radiance Transfer
In the diffuse case of PRT, we only need to compute a single number/color, because diffuse surfaces reflect the same intensity in all directions. The double product integral fits the job, because the convolution acts like a blur-filter that adds up all the light hitting the surface and reflects a single value. But for non-diffuse surfaces like metals this outgoing radiance varies in different directions, so we need more than just one number; we need a function that represents all light reflected into all different directions. We need a coefficient vector of the reflected light that we generate from the incident light coefficient vector.
For diffuse PRT, we turn a transfer function of just the incidenct light direction into the transfer coefficient vector. Now however we must turn a transfer function of an incident and outgoing direction into a transfer matrix. But how to construct the matrix for a function that has two variables? We can do so by imagining that we have a special set of transfer coefficients for each outgoing direction , rather than just a single set for all of them. Where before we would turn a transfer function into a set of coefficients, we now basically convert it into a set of functions which return the -th coefficient of the reflected light for an outgoing direction . To do this, we first integrate along all incident directions .
We then integrate each resulting function , this time over all outgoing directions .
In the previous article, the double integral product was used to compute the outgoing reflected light from the incident light coefficients and the surface transfer coefficients . We transferred incident light into single value (or, if you want to think about it this way, into function which equals a constant value), which is perfect for diffuse reflection. Now with Matrix Radiance Transfer a matrix instead represents the transfer of incident to exit light.
Where we used to have a dot product of two vectors, we now produce a reflection coefficent vector from incident light, and then reconstruct the light reflected into a specific direction by multiplying it with the basis .
Conclusion
Like regular vectors in 3D, coefficient vectors can likewise be transformed with matrices. This is useful when a function - encoded as such a vector - should be transformed into a different function, such as unoccluded to occluded light, normal to rotated environment, incident to exit radiance or when simply decoupling one transfer function into several parts and then multiplying them back together again.
References
- Wikipedia, Clebsch-Gordon Coefficients as Wigner-3j expression
- Free Software Foundation, GNU Scientific Library
- Ng et al., Triple Product Wavelet Integrals for All-Frequency Relighting
- Zhou et al., Precomputed Shadow Fields for Dynamic Scenes
- Pan et al., Precomputed Radiance Transfer Field for Rendering Interreflections in Dynamic Scenes
- Tobias Alexander Franke, SH rotation matrix implementation
- Lehtinen and Kautz, Matrix Radiance Transfer