NICO | GALOPPO || COMP 238 3: Irradiance Gradients

Note: click any of the images to see them at full resolution.

Method and implementation
Results
Discussion
Code/References

Method and implementation

Irradiance gradients [1] introduce a new method for improving the accuracy of a diffuse interreflection calculation. With normal irradiance caching, irradiance values at positions between cached positions are interpolated by a normal weighted sum. The irradiance gradients paper shows that cubic interpolation of irradiance based on gradient information is muchmore accurate than a standard linear interpolation. The gradients predict the change in irradiance as a function of position (translational gradients) and surface orientation (rotational gradients). The basic idea is to take advantage of all the information contained in a hemisphere sampling (such as direction, distance to surface, ...).

I have started from pbrt's IrradianceCache interpolator implementation, to produce my own new IrradianceGradients interpolator plugin. The first major change was the hemisphere sampling pattern. Pbrt uses a cosine-weighted hemisphere sampling; I had to implement a stratified hemisphere sampling pattern to be able to correctly compute gradients between the strata (it uses the dimensions of the strata). I sample the hemisphere at spherical coordinates (&theta_j, &phi_k), and calculate the indirect radiance L_j,k in that direction.

The indirect irradiance E is then simply the Monte Carlo sum of the samples, and is then stored in an octree datastructure:

I created a basic test scene with a sphere that receives only indirect lighting from the floor and nearby walls to emphasize the difference among the irradiance caching algorithms. The following image shows a distant view. A plane separates the sphere and the spotlight on the ceiling. Most of my results will be shown in grayscale (non-colored lights, walls and objects) although my implementation does not require this. As a matter of fact, I show a colored result at the end of this page.

The test scene. The horizontal plane blocks the sphere from direct lighting coming from the spotlight on the ceiling.

With plain irradiance caching using basic weighted linear interpolation as described above, the scene looks as in the following images. Here, I compare pbrt cosine-weighted sampling to the stratified hemisphere sampling. Note the discontinuities at the borders of the regions of influence of the cached samples. This causes the typical blotchiness.

Plain irradiance caching, cosine-weighted versus stratified hemisphere sampling. (512 samples on the hemisphere, 0.1 error acceptance margin)

Gradient calculation and interpolation

In a similar approach to the paper's, I opted to represent the irradiance gradient as the combination of two separate three-dimensional vectors, or rather two multi-component three-dimensional data types. I created a SpectrumGradient class as a natural extension to pbrt's Spectrum class, in order to account for multiple-component radiances (ie. colors). This was necessary because the computation of the gradients requires basic linear algebra between these data types.

Hence, each octree entry now also contains gradient information in a translational gradient SpectrumGradient tGrad and a rotational gradient SpectrumGradient rGrad. The first encodes the change in irradiance as we move over the surface, the second encodes the change as we rotate the hemisphere. The next section will go more into detail on how I computed these gradients. Basic linear weighted interpolation is now extended to cubic interpolation, using the stored irradiance and irradiance gradient information with the following equation:

The weights used are the same as in the original pbrt irradiance cache implementation, which differs slightly from the one in the irradiance gradients paper, but they seem to work well.

Note that the irradiance gradients are not used to determine the spacing of the samples, because that could bias the calculation, since areas that just happened to have small irradiance gradients would be sampled at low density, even though there could still be sudden changes in the irradiance value due to nearby surfaces.

With the creation of the SpectrumGradient class, I effectively performed the gradient calculation separately for each sample wavelength instead of combining into a spectral average. This is beneficial for example in a 'cornell box' type of setup, where there can be sudden changes in color value, for example where two nearby samples lie on a red wall and on a green wall respectively. In this case, even though the spectral average (luminosity) might not change too much, the individual components of the radiance does.

Results

We'll look at the effect of plugging in the rotational and translational gradients into the previously defined cubic interpolation equation separately. I provided to possibility to do this in the scene definition file.

Rotational gradients

First, let's look at the change in irradiance when we rotate the hemisphere. Because the contribution of each hemisphere sample affects the total irradiance magnitude proportional to the cosine projections of the contribution, the rotational gradient can be calculated by summing up the differentials of the cosines. Because of the implicit cosine in the sampling process, tangents show up in the final formula (sine over cosine):

v_k is the base plane unit vector in the φ_k - π/2 direction. Note that I subtract π/2 instead of adding as in the paper, because pbrt uses a left-handed coordinate system, whereas I *suspect* the authors of the paper use a right-handed system.

The following two images show the effect of plugging in the rotational gradient into the cubic interpolation equation only. Here, I used a higher error acceptance threshold, in order to emphasize the effect. By doing this, fewer radiance values are generated and stored into the cache, allowing for interpolation over a wider spatial range instead. Hence, the 'blotches' become larger, but it nicely shows the improvement of using the gradients.

No gradients used vs. rotational gradients. Note the discontinuities between the influence blotches disappear almost completely on the sphere (both images with 512 hemisphere samples and 0.2 error acceptance margin).

Note that the areas where the normal of the interpolated point differs from the normal of the cached point (ie. on the sphere) benefit greatly from this improved interpolation. As can be seen on the floor surface, the areas with constant normals obviously don't. This follows trivially from the interpolation equation. Compare the following image to the first two. It was generated with the same parameters (lower error acceptance threshold, smaller patches), but with rotational gradients added. Again, the major benefit is on the sphere's surface.

Rotational gradients added to the setup of the first section (with lower acceptance error threshold 0.1, hence the smaller blotches)

Here are some extra images, a top view of the sphere. The camera is located between the horizontal plane and the sphere.

No gradients versus rotational gradients in the interpolation. Top view. Error threshold 0.2

Translational gradients added

To compute the translational irradiance gradient, we consider how the projected solid angle of each of the hemisphere cells is affected by a translation of the hemisphere center. It's the change in projected solid angle, when multiplied by the radiance of the cell that will give the change in that cell's contribution. The rate of change is determined in part by the distance to the contributing surface. In fact, it is always the distance to the closer surface that determines the rate of change in occulusion, since the relative motion of a foreground surface is greater than that of a background surface. This observation explains the Min operators in the final equation. The equation consists of a running sum of two terms, the first term dictates the change for translation in the direction parallel to φ_k, the second term dictates the change for the direction perpendicular to φ. Please refer to the paper for further details.

Note that there seems to be a problem with the formula. I believe that the unit vector u_k in the φ_k direction should be inside the inner sum. Even though it seems that u_k only depends on k it must be noted that this is not the case, because it is defined to be the base plane unit vector in the φ_k direction, which is different for every sample (j,k).

Even more so, I presume the authors were using a right-handed system. Again, I subtract π/2 instead of adding it to compute v_{k_-}, and flip the sign of the first term in the running sum. The following pair of images shows that the addition of the translational gradient component improves the interpolation. Note the smoother shadow beneath the sphere. Nevertheless, there are artefacts at the right side of the shadow, which are not there in the image without the translational gradient component.

Rotational gradients only (left) vs. rotational plus translational gradients added to the interpolation (added). Note the smoother on the left, but also the artefacts on the right side of the shadow. (512 samples of the hemisphere, 0.2 error acceptance threshold)

There are several possibilities for the artefacts that I am observing, and probably it's a combination of factors. First and foremost, I suspect there is still a bug in the code. I tried to flip signs (as described above) consistently in such a fashion that it would flip only one axis of the coordinate system (ie. going from right-handed to left-handed system), but it could very well be that I have not taken all transformations into account. I went over the code several times, but could not find such an instance.

Secondly, there might be some incompatibilities with the algoritm and the specific renderer that it is being plugged into (in this case pbrt) as opposed to the author's Radiance. I will discuss these in the following section.

Discussion

Scanline renderer
After talking to several people at the department that have implemented irradiance gradients in the past (at Rhythm&Hues Studios), it became clear that the algorithm isn't really suited for scanline renderers (as is pbrt). Even more, it has been reported by others on the web [2]. The effect can be observed in my images by looking at the sphere. Indeed, I noted that the rotational gradients did a much better job at interpolation between the cache sites, but you can still see the vertical pattern coming through. Indeed, that is because, in a scanline tracer, no cache sites could have been generated 'behind' the current pixel being traced.
A quick but dirty solution could be to do a cache site generation (no shading) run of the image pixels first as a preprocess, and then the second shading pass as before. Better approaches are a uniformly distributed random tracer, or a PSR (pixel-selected-ray) tracer. Unfortunately, due to time constraints I couldn't try out either approach.
World space vs. local space
The equations in the paper freely mix surface-space quantities and vectors (u_k, v_k, φ_k, θ_j with world-space vectors (translational and rotational gradients). Obviously, the interpolation is done in world space, whereas the paper defines the gradient calculation in (local) surface-space. One must beware of that face and make sure to transform from local to world space somewhere down the line, even if the paper omits such details. I solved this by computing the gradients in local space, and then use the local surface geometry to convert to world space.
Handedness
Again, the paper doesn't mention the handedness being used, even though I am pretty convinced that it does matter in the equations. I believe the paper uses a right-handed coordinate system.
Errors in the summations
Beware of the inner sum in the equation for the translational gradient, as I noted in the previous section.
Efficiency
The paper mentions this, and my experiments confirm that the computational overhead to plain irradiance caching is marginal. For a very small cost, the visual quality improves substantially. Not only the computational cost, but also the additional memory requirements are very small, because the number of cache sites remains constant, and we only need to store 2 extra 3-dimensional vectors (or c times that quantity for c components or wavelengths).

To end with, I'll throw in some colored images, even though I believe the color tends to smoothen out the benefits as well as the artefacts. Note that the sphere was still gray, but the left wall is red and the right wall is green.

The red left wall of the box and the green right wall cast off colored light on the sphere by indirect illumination coming from the white spotlight on the ceiling. As before, the sphere is obstructed from direct lighting by the plane in between the sphere and the spotlight.

References and code

The code can be found here. It includes the implementation of a IrradianceGradients plugin in irradiancegradients.cpp, and my basic setup scene file irrtest.pbrt.

Irradiance Gradients by Greg Ward and Paul Heckbert.
Discussion on irradiance gradients limitations in scanline raytracers.