Importance sampling for infinite area light is off by half a texel

[johguenther]

As explained in Section 14.6.5 and Figure 14.13, p.726 in PBRT v2 the importance of the environment texture should be computed by interpolating adjacent texels. Therefore InfiniteAreaLight samples the texture at the discrete coordinates uv. However, this means that the pdf in Distribution2D for uv represents the average importance for the area (uv-0.5, uv+0.5). Now, SampleContinuous selects a bin according to the pdf and samples uniformly within that bin, i.e. the returned coordinates represent [uv, uv+1.0). To account for this difference in coordinates the Sample_Li and Sample_Le functions of the infinite area light should subtract "half a texel" from uv, e.g.:

--- a/src/lights/infinite.cpp
+++ b/src/lights/infinite.cpp
@@ -98,3 +98,5 @@ Spectrum InfiniteAreaLight::Sample_Li(const Interaction &ref, const Point2f &u
     // Find $(u,v)$ sample coordinates in infinite light texture
     Float mapPdf;
     Point2f uv = distribution->SampleContinuous(u, &mapPdf);
+    uv.x -= (Float)0.5 / Lmap->Width();
+    uv.y -= (Float)0.5 / Lmap->Height();

So, Figure 14.13 is actually misleading, because for textures their values are defined at the center of each texel (at continuous coordinates uv+0.5, see Section 10.4.3, p.538): the piecewise linear function should therefore be shifted by 0.5.

The error is subtle and only gets noticed for high-contrast environment maps. For example, using a black texture with only one white texel, half of the samples are wasted (they have zero contribution). A nice way to debug these issues is to visualize the sample contributions of the infinite area light, i.e. render L = infiniteLight->Sample_Li(...) / lightPdf using the (normalized) pixel coordinates as random numbers. The more constant the resulting image, the better the importance sampling.

Some example images (before and after the correction) are attached:

Single white pixel, original:
Single white pixel, with fix:
Grace Cathedral HDRI, original:
Grace Cathedral HDRI, with fix:

Some more remarks on minor issues:

• The environment texture is resampled by MIPMap to a power-of-two resolution. To avoid aliasing artefacts the importance map should be created with that resolution as well (that's why the above patch uses Lmap->Width() and Lmap->Height()):

  --- a/src/lights/infinite.cpp
  +++ b/src/lights/infinite.cpp
  @@ -63,3 +63,3 @@ InfiniteAreaLight::InfiniteAreaLight(const Transform &LightToWorld
       // Compute scalar-valued image _img_ from environment map
  -    int width = resolution.x, height = resolution.y;
  +    int width = Lmap->Width(), height = Lmap->Height();

• Because we want to compute the average importance for the area (uv-0.5, uv+0.5) the average of sinTheta is better approximated with "sin(uv)":

  --- a/src/lights/infinite.cpp
  +++ b/src/lights/infinite.cpp
  @@ -68,5 +68,5 @@ InfiniteAreaLight::InfiniteAreaLight(const Transform &LightToWorld,
       for (int v = 0; v < height; ++v) {
           Float vp = (Float)v / (Float)height;
  -        Float sinTheta = std::sin(Pi * Float(v + .5f) / Float(height));
  +        Float sinTheta = std::sin(Pi * vp);
            for (int u = 0; u < width; ++u) {

  This leads to zero probability for the first row, which is actually a good thing: uv cannot become negative after subtracting half a texel. And sampling the north pole leads to artefacts anyway (because of texture wrapping light leaks from the south pole).

• The computation of the filter width could be skipped. The default width of zero in MIPMap::Lookup already ensures that the most detailed mipmap level is used, and the triangle filter used for bilinear interpolation computes the wanted average of the four texels.

[mmp]

Wow! Really impressive work working through those details.

I'd like to digest this carefully, do some tests (and figure out what edits
to make to the book text to reflect this); it may be a few days before this
is merged/fixed, but thank you for figuring it out!

On Mon, Mar 14, 2016 at 6:44 AM, Johannes Günther notifications@github.com
wrote:

As explained in Section 14.6.5 and Figure 14.13, p.726 in PBRT v2 the
importance of the environment texture should be computed by interpolating
adjacent texels. Therefore InfiniteAreaLight samples the texture at the
discrete coordinates uv. However, this means that the pdf in
Distribution2D for uv represents the average importance for the area (uv-0.5,
uv+0.5). Now, SampleContinuous selects a bin according to the pdf and
samples uniformly within that bin, i.e. the returned coordinates represent [uv,
uv+1.0). To account for this difference in coordinates the Sample_Li and
Sample_Le functions of the infinite area light should subtract "half a
texel" from uv, e.g.:

--- a/src/lights/infinite.cpp+++ b/src/lights/infinite.cpp@@ -98,3 +98,5 @@ Spectrum InfiniteAreaLight::Sample_Li(const Interaction &ref, const Point2f &u
// Find $(u,v)$ sample coordinates in infinite light texture
Float mapPdf;
Point2f uv = distribution->SampleContinuous(u, &mapPdf);+ uv.x -= (Float)0.5 / Lmap->Width();+ uv.y -= (Float)0.5 / Lmap->Height();

So, Figure 14.13 is actually misleading, because for textures their values
are defined at the center of each texel (at continuous coordinates uv+0.5,
see Section 10.4.3, p.538): the piecewise linear function should therefore
be shifted by 0.5.

The error is subtle and only gets noticed for high-contrast environment
maps. For example, using a black texture with only one white texel, half of
the samples are wasted (they have zero contribution). A nice way to debug
these issues is to visualize the sample contributions of the infinite area
light, i.e. render L = infiniteLight->Sample_Li(...) / lightPdf using the
(normalized) pixel coordinates as random numbers. The more constant the
resulting image, the better the importance sampling.
Some example images (before and after the correction) are attached:
Single white pixel, original: [image: single_pixel_before]
https://cloud.githubusercontent.com/assets/9035336/13745787/4309cd10-e9f0-11e5-9db2-465c36b3574b.png
Single white pixel, with fix: [image: single_pixel_after]
https://cloud.githubusercontent.com/assets/9035336/13745788/430a000a-e9f0-11e5-8e20-19b10bfa61a8.png
Grace Cathedral HDRI http://gl.ict.usc.edu/Data/HighResProbes/,
original: [image: grace-new_before]
https://cloud.githubusercontent.com/assets/9035336/13745835/7b5995ec-e9f0-11e5-9240-82fcf73834c2.png
Grace Cathedral HDRI, with fix: [image: grace-new_after]
https://cloud.githubusercontent.com/assets/9035336/13745836/7b59896c-e9f0-11e5-8786-ada4cbf959ec.png

Some more remarks on minor issues:

The environment texture is resampled by MIPMap to a power-of-two
resolution. To avoid aliasing artefacts the importance map should be
created with that resolution as well (that's why the above patch uses
Lmap->Width() and Lmap->Height()):

--- a/src/lights/infinite.cpp+++ b/src/lights/infinite.cpp@@ -63,3 +63,3 @@ InfiniteAreaLight::InfiniteAreaLight(const Transform &LightToWorld
// Compute scalar-valued image img from environment map- int width = resolution.x, height = resolution.y;+ int width = Lmap->Width(), height = Lmap->Height();

Because we want to compute the average importance for the area (uv-0.5,
uv+0.5) the average of sinTheta is better approximated with "sin(uv)":

--- a/src/lights/infinite.cpp+++ b/src/lights/infinite.cpp@@ -68,5 +68,5 @@ InfiniteAreaLight::InfiniteAreaLight(const Transform &LightToWorld,
for (int v = 0; v < height; ++v) {
Float vp = (Float)v / (Float)height;- Float sinTheta = std::sin(Pi * Float(v + .5f) / Float(height));+ Float sinTheta = std::sin(Pi * vp);
for (int u = 0; u < width; ++u) {

This leads to zero probability for the first row, which is actually a
good thing: uv cannot become negative after subtracting half a texel. And
sampling the north pole leads to artefacts anyway (because of texture
wrapping light leaks from the south pole).

• The computation of the filter width could be skipped. The default
  width of zero in MIPMap::Lookup already ensures that the most detailed
  mipmap level is used, and the triangle filter used for bilinear
  interpolation computes the wanted average of the four texels.

—
Reply to this email directly or view it on GitHub
#62.

[mmp]

Very nice find; 0.5 offset fix has been applied.

I'm not convinced about the change to the sinTheta computation, however. Consider for example a 1x1 environment map with a constant radiance value (as is created by InfiniteAreaLight if no filename is given but a spectrum is provided); in that case, the result should be a uniform radiance function in all directions, whereas with the proposed change, PDF would be zero everywhere. Similarly, I'd argue that if the environment map only had non-black pixels in the top scanline, it should still be casting illumination.

Maybe the right thing is to be computing the average value of |sin(theta)| over the continuous u range?

[johguenther]

Very valid counterexamples that sin(0)=0 is bad for the first row. I still think computing sinTheta at the discrete coordinates is a better approximation for most of the rows (sinTheta is assumed constant within a bin, so taking the value of sin at the center of the bin is slightly better than taking the value at the border). But the first row needs indeed special treatment, because that bin represents a discontinuity (both poles).
If we assume for the the first row that

• sinTheta is piecewise linear, being zero at the center of the bin and sinThetaPole=sin(pi*0.5/height) at both borders,
• the first half of the bin represents the south pole, the second half the north pole (with the texel values constant due to the wanted "clamp to edge" behavior in v),

then the importance can be approximated by

1/2 sinThetaPole texelSouth + 1/2 sinThetaPole texelNorth
= sinThetaPole 1/2 (texelSouth + texelNorth)
= sinThetaPole Lookup(u, 0),

which is what is currently implemented. Thus, sinTheta should stay as is for the first row, but could be changed for all other rows.

Still, v needs special care after subtracting half a texel in Sample_Le and Sample_Li. If v is negative, the south pole was sampled and 1.0 should be added to v.

[johguenther]

One more thing: the shift by half a texel (in the opposite direction) needs also be done in Pdf_Li and Pdf_Le before calling distribution->Pdf (and the wrap-around case when uv > 1. needs to be handled as well).

Additionally, Pdf_Le should use SphericalPhi as well to reproject negative results from atan2 – currently half of the directions return the pdf of the first bin. However, InfiniteAreaLight::Pdf_Le is not used anywhere (bdpt handles infinite lights on its own), so that bug never showed.

[mmp]

(Fixed that SphericalPhi issue in Pdf_Le--thanks!)

I've been thinking about this more and think that a different fix makes more sense. (I'd love to hear your feedback on this thinking, though!)

The issue that got me started on this was running into a bug where now, if a single-texel environment map is used, a negative PDF will often be returned. It turns out that this is because of the half-texel shift, which in turn causes a negative sin(theta) value to be computed. (Upon further digging, it seems that this sometimes happens for all environment maps.)

In turn, this got me thinking... I'd suggest that all of the InfiniteAreaLight code should only operate on continuous coordinates and shouldn't be thinking about continuous vs discrete coordinates at all. (Or, as little as possible). For the most part, we'd just like to be thinking in terms of continuous functions over the domain, sampling those functions, etc., without having to worry about how the functions are represented.

Under that assumption, almost all of the methods of InfiniteAreaLight can go back to how they were before. However, the one place where this issue does have to be addressed is in the construction of the piecewise-constant function used for importance sampling. Here is how I think that should work:

• Consider for example a 1 dimensional texture with two texels. The continuous coordinates are from [0,2], or can be scaled to [0,1]. (We'll use [0,1] for here on out.)
• The two texels are at 0.25 and 0.75, in continuous coordinates.
• The Distribution1D is specified so that the first value given is the piecewise constant function for [0, 0.5] and the second is the piecewise constant function for [0.5, 1]
• Because MIPMap::Lookup takes continuous coordinates, we want to do one lookup at 0.25 and one at 0.75. Those are exactly on the discrete texels, so we want to specify a filter width that gets over to the neighboring texel.

If I do all that and do the same visualization, the results seem as good or better. (I also doubled the resolution of the tabulated PDF, which helps as well, so it's not an equal comparison.)

Visualizations:

What do you think?

[mmp]

(The first visualization was with the initial fix, the second is with the current top of tree.)

[johguenther]

That's a nice alternative solution. Doubling the number of bins is actually crucial: it ensures that the function is linear within the bin and can thus be well approximated.

The single white pixel map now also looks even better:

Again I think the filter for the MIPMap lookup is not necessary, because bilinear interpolation suffices (essential it is now done not halfway between the texels, but with weights 0.25 and 0.75).