Converting SH Radiance to Irradiance

In the last post we explained an alternative, simpler way of defining the Spherical Harmonic basis functions. We are now going to use this definition of the L0 and L1 basis functions to investigate a classic lighting use case: converting radiance to irradiance.

Radiance is the incoming light at a point. Since this is a spherical function, it can be approximated by a truncated SH series, and this approximation can be computed by sampling (in, say, your favourite path tracer or real-time global illumination system). However, for shading geometry we need to know outgoing light, also known as irradiance.

For the purposes of this post we assume the perfectly diffuse Lambertian BRDF, which is a constant function. Therefore irradiance is the following hemispherical integral:

    \[ I(\vec{\boldsymbol{n}})=\frac{1}{\pi}\int_H\vec{\boldsymbol{n}}\!\cdot\!\vec{\boldsymbol{\omega}} R(\vec{\boldsymbol{\omega}}) d\Omega \]

where \vec{\boldsymbol{n}}\!\cdot\!\vec{\boldsymbol{\omega}} is the geometry term.

Computing the SH approximation we can replace R(\vec{\boldsymbol{\omega}}) by R_{sh}(\vec{\boldsymbol{\omega}}) and compute the SH coefficients of independently:

    \[ I_0=\frac{1}{\pi}\int_H\vec{\boldsymbol{n}}\!\cdot\!\vec{\boldsymbol{\omega}}\, R_0 \, d\Omega = R_0 \]

    \[ \vec{\boldsymbol{I}}_1=\frac{1}{\pi}\int_H \vec{\boldsymbol{n}}\!\cdot\!\vec{\boldsymbol{\omega}} \, 3 \vec{\boldsymbol{R}}_1\!\cdot\!\vec{\boldsymbol{\omega}} \, d\Omega =2\vec{\boldsymbol{R}}_1 \]

    \[ I_2^{ij}=\frac{1}{\pi}\int_H \vec{\boldsymbol{n}}\!\cdot\!\vec{\boldsymbol{\omega}} \, \tfrac{15}{2}\omega_i \omega_j R_2^{ij} =\tfrac{15}{8}R_2^{ij} \]

Putting it all together, we find that the irradiance reconstruction is:

    \[ I_{sh}(\vec{\boldsymbol{\omega}})=R_0+2\vec{\boldsymbol{R}}_1\!\cdot\!\vec{\boldsymbol{\omega}}+\tfrac{15}{8} \omega_i \omega_j R_2^{ij}+\ldots \]

Here the radiance reconstruction constants have been absorbed into the irradiance computation. This means we have only one set of constants to apply to convert the measured radiance coefficients into irradiance.

…and we’re done?

For quite a while the Enlighten SH solver generated L1 irradiance using this linear reconstruction. We considered it an optimisation to fold the conversion into the solver itself. However it’s pretty easy to see this has serious problems. One nice feature of our definition for L1 is that we know 0 \le |\vec{\boldsymbol{R}}_1| \le R_0. The lower bound is attained when the lighting is completely symmetrical or ambient (with no preferred direction), while the upper bound occurs if the incoming light comes from a single direction. However, the factor of two in the linear radiance-to-irradiance conversion means that for very strongly directional lighting environments, our approximation for irradiance will be negative in some direction. This is clearly undesirable since irradiance can never be negative in reality, and causes obvious visual artefacts in practice.

Improving L1 irradiance with a non-linear model

To improve this situation, we eventually unbaked the conversion constants in the solver (so it once again generated SH radiance), and applied a more complex conversion to irradiance in shader code. The rest of this post contains the motivation and derivation of this improved L1 irradiance model.

We already know that the formula given above is the best possible linear conversion to irradiance. Therefore we need to look for a non-linear alternative. To narrow this down we’re first going to find a better polynomial approximation for the pure directional case |\vec{\boldsymbol{R}}_1| = R_0. Once we’re happy with that we can generalise as |\vec{\boldsymbol{R}}_1| / R_0 varies. Finally, we need to ensure the energy and dynamic range of the outgoing irradiance are correct for our measured radiance values.

Pure directional irradiance

For an ideal point light source with energy 1 in direction \vec{\boldsymbol{d}}, Lambertian irradiance is given by:

    \[ I(\vec{\boldsymbol{n}}) = \begin{cases} 4 (\vec{\boldsymbol{d}} . \vec{\boldsymbol{n}}) & \text{if}\ \vec{\boldsymbol{d}} . \vec{\boldsymbol{n}} \ge 0 \\ 0 & \text{otherwise} \end{cases} \]

This has a minimum value of 0 (attained on a whole hemisphere), and a maximum value of 4. Compare this to the linear reconstruction for L1 irradiance:

    \[ I_{lin}(\vec{\boldsymbol{n}}) = 1 + 2 (\vec{\boldsymbol{d}}.\vec{\boldsymbol{n}}) \]

This has a minimum of -1 and a maximum of 3. Although this integrates to the correct energy as written, when we use it in practice in a shader we’ll clamp negative values to 0, with the net result that the total energy of the function will be too high. Therefore the linear reconstruction gives us the wrong energy, and the wrong dynamic range (incorrect maximum value).

Best smooth approximation

To find our non-linear approximation in the pure directional case, let’s first specify that the approximation is smooth. This makes sense since we want it to be one endpoint of a smoothly varying overall model. To satisfy this, we will try to fit a polynomial function f (which will automatically be smooth) in the variable x = \vec{\boldsymbol{d}}.\vec{\boldsymbol{n}}.

There are constraints on this function that we already know: the total energy (integral under the curve) should be 2, the minimum value 0 and the maximum value 4.

Additionally, we are going to specify derivative conditions at -1 (the point opposite the light source). The true irradiance is 0 on the whole back-facing hemisphere. We can’t attain that with a polynomial, but we can at least specify that f'(-1) = f''(-1) = 0.

This gives us five constraints to satisfy, so we can fit a unique quartic ax^4 + bx^3 + cx^2 + dx + e. Solving for the coefficients gives us:

    \[ f(x) = \frac{1}{2} (1 + x)^3$ \]

Generalise as length of L1 vector varies

We now know what to do at the two extreme cases. When |\vec{\boldsymbol{R}}_1 | = R_0 we can plug in our non-linear approximation. Alternatively, if |\vec{\boldsymbol{R}}_1 | = 0 we have no directional variation at all, so irradiance is simply a constant function. But what should be do in between the two extremes?

For each value in between there is a range of possibilities, from a linear model to something “as non-linear” as our extreme model. It’s reasonable to suppose that our model should become more linear as it becomes more ambient (less directional). So let’s define the model power to be:

    \[ p = 1 + 2 \frac{| \vec{\boldsymbol{R}}_1 |}{R_0} \]

This varies between 1 and 3 based on the length of the L1 vector. Our final model will then be based on (1+x)^p where x = \hat{\boldsymbol{R}_1}.\vec{\boldsymbol{n}}, that is the dot product with the normalised direction of the L1 vector.

At this point it’s useful to introduce a new variable q = \frac{1}{2}(1 + x), so the model is of the form a + bq^p where a and b are constants we need to find.

We have two constraints to apply to find a and b: total energy, and dynamic range.

We are assuming that the total energy is 2 (only the length of L1 is varying, not L0). Therefore we want to find b so that \int_{-1}^{1} bq^p dx = 2. It follows that b = p + 1. Our model is then a linear interpolation between the ambient and directionally-varying terms:

    \[ a + (1-a)(p+1)q^p \]

The correct dynamic range (maximum – minimum) varies linearly with the length of L1 and equals 4 \frac{|\vec{\boldsymbol{R}}_1 |}{R_0}. The dynamic range of the model is (1-a)(p+1). Substituting in the definition of p and solving for a gives:

    \[ a = (1 - \frac{|\vec{\boldsymbol{R}}_1 |}{R_0}) / (1 + \frac{|\vec{\boldsymbol{R}}_1 |}{R_0}) \]

Final improved L1 irradiance model

Putting it all together, the final model is:

    \[ I_{sh}(\vec{\boldsymbol{n}}) = R_0 (a + (1 - a)(p + 1)q^p ) \]

where

    \[ q = \frac{1}{2}(1 + \hat{\boldsymbol{R}_1}.\vec{\boldsymbol{n}}) \]

    \[ p = 1 + 2 \frac{| \vec{\boldsymbol{R}}_1 |}{R_0} \]

    \[ a = (1 - \frac{|\vec{\boldsymbol{R}}_1 |}{R_0}) / (1 + \frac{|\vec{\boldsymbol{R}}_1 |}{R_0}) \]

Conclusion

We’ve discussed how to convert SH radiance (incoming light) into SH irradiance (outgoing bounced light), assuming a diffuse Lambertian BRDF. The natural, linear best-fit method generates functions which can attain negative values, which in the case of light is an obviously undesirable property.

In the case of L1 we fix this with a non-linear reconstruction, fitting the model to the measured direction, energy and dynamic range. This model gives very good “bang for the buck”: it produces decent shading quality, despite only using a small number of coefficients (four per colour channel).

Unfortunately, it’s frustratingly unclear how to extend this non-linear reconstruction naturally to higher SH bands. Instead, the most practical current approach is to remove negative values by “de-ringing”; Peter-Pike Sloan presented a robust numerical method for doing this at SIGGRAPH Asia 2017.

Leave a Reply

Your email address will not be published. Required fields are marked *