Where that #!@&?$^*%! PI comes from – The Burning Basis Vector

Another “educational” post:

In computer graphics kindergarten, we were all taught that diffuse is:

$L_d = k_d(N \cdot L)$

Where $k_d$ is a constant that determines how bright the diffuse response is.

In computer graphics grade-school, we then learn about the rendering equation, and we’re given:

$L_{\omega_o} = \int_\Omega L_i f(\omega_o,\omega_i) \cos(\theta_i)\,d\omega$

Where $f(\omega_o,\omega_i)$ is the BRDF (bidirectional reflectance distribution function). For the diffuse case, the BRDF is :

$f(\omega_o,\omega_i) = \frac{r}{\pi}$

where r is albedo

This $\pi$ is the source of much confusion. The albedo is, by definition, the fraction of incident light that gets reflected from the surface. So how come the BRDF isn’t just r?

In this post, I’m going to attempt to give an intuitive answer to this question. The target audience is anybody who’s been puzzled by this question for as long as I was. I’ll be using diffuse here for simplicity, but pretty much every BRDF you see in the literature has got a factor of $\pi$ in it’s normalization term someplace, for the same reasons.

The short answer is this: Integrating intensity over cosine-weighted solid angle gives an answer which is off by a factor of $\pi$ . For the long version, I need to pretend I understand calculus. If you are a CS freshman: pay attention in math class. You’d be surprised how relevant that stuff becomes later in your career.

About Integrals

Lets imagine we have a greyscale image. If we want to compute its average value, we know how do do that, just sum over pixels and divide:

$A=\frac{\displaystyle\sum pixel(i,j)}{N_{pixels}}$

I find that a good, intuitive mental model for an integral is an average in the limit. If we could somehow average over all infinity points on the continuous domain, the integral is the answer that we would get. If we happen to have a closed-form expression for the function, we can often derive the integral analytically, but in graphics, we never do. Luckily, we get can arbitrarily close to the right answer just by taking a crapload of samples and averaging them appropriately. Another important aspect of an integral is that it is measuring a density. In our example above, if our image represented some commodity like light or happiness, then the integral over the domain is the amount of white stuff per unit area.

If we liked, we could average over scan-lines. We take the average in each line, and then average the averages, like this:

$A=\frac{ \displaystyle\sum_j\frac{\displaystyle\sum_i pixel(i,j)}{Width}}{Height}$

It’s easy to see that we’ll get exactly the same answer. This is what a multi-dimensional integral is (kinda,sorta)

The Rendering Equation is Not What You Think It Is

Here’s that rendering equation again:

$L_{\omega_o} = \int_\Omega L_i f(\omega_o,\omega_i) \cos(\theta_i)\,d\omega$

Since we’re limiting ourselves to the diffuse case, the BRDF is a constant, so we can yank it out:

$L_{\omega_o} = f\int_\Omega L_i\cos(\theta_i)\,d\omega$

One important fact which graphics literature rarely mentions is that the integration over solid angle can and should be expressed as a multi-dimensional integral in spherical coordinates:

$L_{\omega_o} =f \int\limits_{\phi=-\pi}^{\pi} \int\limits_{\theta=0}^{\frac{\pi}{2}} L_i \cos(\theta)sin(\theta)\,d\theta\,d\phi$

You may not have seen this $sin(\theta)$ term before, but its the source of all the trouble.

The rendering equation is trying to average something over solid angle, which the same as averaging over the surface area on a unit sphere above our shading point. It turns out that the spherical coodinate system is skewed. If you hold $\theta$ fixed, and sweep the $\phi$ angle around the pole, you sweep out a great circle, as shown below.

The circumference of the circle increases with $\theta$ . These little circles are analogous to the scanlines in our earlier example. If we imagine cutting up the sphere into circles and assigning the same number of ‘pixels’ to each, we end up with sample points that are closer together at the top of the sphere than at the bottom. In order to compute the correct density, we now need to do a weighted average, so that each sample’s contribution is proportional to the amount of surface area it covers. This is what the $sin(\theta)$ term is doing.

Now we have everything we need to explain that $\pi$ away.

Finding the PI

Let’s make things simple and build a test case that we can actually solve in closed form. We’ll assume that the incoming light intensity is constant over the entire hemisphere of directions. Now we have:

$L_{\omega_o} =fL_i \int\limits_{\phi=-\pi}^{\pi} \int\limits_{\theta=0}^{\frac{\pi}{2}}\cos(\theta)sin(\theta)\,d\theta\,d\phi$

Time for some mathematical slight-of-hand. Had I paid more attention in calculus class, this part would have been much easier to write. Observe that:

$\int\limits_{\theta=0}^{\frac{\pi}{2}} cos(\theta)sin(\theta) = \frac{1}{2} \int\limits_{\theta=0}^{\frac{\pi}{2}}2cos(\theta)sin(\theta)$

We need the antiderivative of $2cos(\theta)sin(\theta)$ which we can get by abusing the product rule, which states:

$D( f(x)g(x) ) = f$

Let $f(x)=g(x)=sin(x)$ . Now:

$D( sin(x)sin(x) ) = cos(x)sin(x) + sin(x)cos(x) = 2sin(x)cos(x)$

So:
$\int\limits_{\theta=0}^{\frac{\pi}{2}} cos(\theta)sin(\theta) = \frac{1}{2} \int\limits_{\theta=0}^{\frac{\pi}{2}}2cos(\theta)sin(\theta) = \frac{1}{2} [ sin^2(\frac{\pi}{2}) - sin^2(0) ]= \frac{1}{2}$

Now we just plug that in the outer integral and we have:

$L_{\omega_o} = fL_i\int\limits_{\phi=-\pi}^{\pi}\frac{1}{2}\,d\phi = fL_i\pi$

Energy Conservation

So, in this scenario, we put in radiant intensity $L_i$ from all directions and get $L_{\omega_o} = fL_i\pi$ out. Did we conserve energy? Probably not, but let’s prove it.

Energy conservation means that the flux leaving the surface must not exceed the flux arriving at the surface. We’ll write this as:

$F_o \leq F_i$

Once again, we construct an imaginary test case that makes the math easier. If we have radiant intensity $L_i$ over the entire hemisphere, and intensity is flux per steradian (unit hemisphere area), then the incoming flux is given by:

$F_i = L_i\int\limits_\Omega\,d\omega = 2\pi L_i$

If you’re not sure what flux and intensity are, you may find my earlier post helpful.

Now lets calculate how much flux is leaving. We’ve just calculated $L_{\omega_o}$ , the flux per unit hemisphere area in any particular direction, so we can just integrate that over the outgoing hemisphere.

$F_o = \int\limits_\Omega L_{\omega_o}\,d\omega = \int\limits_\Omega fL_i\pi\,d\omega= fL_i\pi\int\limits_\Omega \,d\omega= 2\pi^2fL_i$

Now, let $f=1$ , meaning that we want $F_o=F_i$ , and plug back in:

$2\pi^2L_i = 2\pi L_i \rightarrow \leftarrow$

Ooops, we’ve got $\pi$ times as much flux as we started with, and that won’t do. In order to fix it, we need to divide out the extra $\pi$ in our BRDF. And that, dear friends, is where that mysterious $\pi$ comes from.

It’s There By Definition

You can also derive this directly from the conventional definition of the BRDF. The BRDF is the ratio:

$f(\omega_o,\omega_i) = \frac{L_{\omega_o}}{I_{\omega_i}}$

Where $L_{\omega_o}$ is the reflected intensity through solid angle ${\omega_o}$ and $I_{\omega_i}$ is the irradiance due to light entering from solid angle ${\omega_i}$ .