Where that #!@&?$^*%! PI comes from

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.

sphere_rings

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}.

Again, assuming uniform incoming intensity L_i and an ideal diffuse reflector with albedo r, we have:

F_o = r F_i

There are 2\pi \omega_o‘s in a hemisphere, and input and output are both constant over the hemisphere, meaning:

L_{\omega_o} = \frac{F_o}{2\pi} = r \frac{F_i}{2\pi} = r L_i

And we’ve already derived irradiance earlier:

I_{\omega_i} = \int\limits_\Omega L_i cos(\theta_i)\,d\omega = ... = \pi L_i

So:

f = \frac{r L_i}{\pi L_i} = \frac{r}{\pi}