~/imallett (Ian Mallett)

Applying Monte Carlo Integration to Evaluate the Rendering Equation (Explicit Formulation)

Recall equation $(3)$ from Monte Carlo Integration (renumbered here):$\int_{domain}g(x) d x = E\left[\frac{g(x)}{pdf(x)}\right] ~~~~~~~~~~~~~~~(1)$. . . and the rendering equation (from my explicit formulation):$\vec{L}(\vec{x},\vec{\omega}_o) = \left[ \left[ \int_{\Omega_{\vec{x},\vec{N}_{\vec{x}}}} \vec{L}(\vec{y}_i,\,-\vec{\omega}_i) \left| \vec{N}_{\vec{x}} \cdot \vec{\omega}_i \right| \vec{f}_r(\vec{x},\vec{\omega}_i,\vec{\omega}_o) \, d\,\vec{\omega}_i \right] + \overrightarrow{Em}(\vec{x},\vec{\omega}_o) \right] \frac{1}{\left|\vec{N}_{\vec{x}}\cdot\vec{\omega}_o\right|} ~~~~~~~~~~~~~~~(2)$

We have a recursive integral equation. It's time to solve it numerically. The hard part is the integral, and we now have a Monte Carlo method to apply to find its average value!

Rename the integrand:$g(\vec{\omega}_i) = \vec{L}(\vec{y}_i,\,-\vec{\omega}_i) \left| \vec{N}_{\vec{x}} \cdot \vec{\omega}_i \right| \vec{f}_r(\vec{x},\vec{\omega}_i,\vec{\omega}_o) ~~~~~~~~~~~~~~~(3)$

Now, it should be easy to see that by these first three equations, the value $E\left[g(\vec{\omega}_i)~/~pdf(\vec{\omega}_i)\right]$ is the same as the integral. This produces our final equation:$\vec{L}(\vec{x},\vec{\omega}_o) = \left[E\left[ \frac{ \vec{L}(\vec{y}_i,\,-\vec{\omega}_i) \left| \vec{N}_{\vec{x}} \cdot \vec{\omega}_i \right| \vec{f}_r(\vec{x},\vec{\omega}_i,\vec{\omega}_o) } { pdf(\vec{\omega}_i) }\right] + \overrightarrow{Em}(\vec{x},\vec{\omega}_o) \right] \frac{1}{\left|\vec{N}_{\vec{x}}\cdot\vec{\omega}_o\right|} ~~~~~~~~~~~~~~~(4)$

BRDF Energy Conservation

It is important to note that all the light that enters the surface must be accounted for. If no absorption happens, all the light that enters must be reflected—so if the surface emits no light of its own, then the total radiance in will be the same as the total radiance out.

This can be represented by the following constraint (found by removing the emission term from $(4)$, setting the incoming radiance $\vec{L}(\vec{y}_i,\,-\vec{\omega}_i)$ and outgoing radiance $\vec{L}(\vec{x},\vec{\omega}_o)$ equal (the reason we can do that may be clear if you consider estimating the expected value with a single sample), and then rearranging):$\vec{k} \cdot \left|\vec{N}_{\vec{x}}\cdot\vec{\omega}_o\right| = E\left[ \frac{ \left| \vec{N}_{\vec{x}} \cdot \vec{\omega}_i \right| \vec{f}_r(\vec{x},\vec{\omega}_i,\vec{\omega}_o) } { pdf(\vec{\omega}_i) } \right] ~~~~~~~~~~~~~~~(5)$In the above, $\vec{k}$ is the absorption coefficient, which is (also) "contained" in the BRDF $\vec{f}_r$. Each channel of the coefficient ranges from $0$ through $1$.

Deriving Diffuse BRDF from Lambert's Cosine Law

Let's see an example of how $(5)$ works out for the real example of an importance-sampled diffuse distribution.

The BRDF comes from Lambert's Cosine Law:$\vec{f}_r(\vec{x},\vec{\omega}_i,\vec{\omega}_o) = \left|\vec{N}_{\vec{x}} \cdot \vec{\omega}_o\right| \cdot \vec{k}_d \cdot x ~~~~~~~~~~~~~~~(6)$Note the absorption coefficient $\vec{k}$ for diffuse has been named $\vec{k}_d$ here, and below. The $x$ is a constant normalization factor that will be found below (see equation $(9)$). The absolute values are given here for convenience; Lambert's Cosine Law is not defined for viewing angles below the ecliptic plane, so the choice doesn't matter.

Recall that we can choose any probability density function. Certain functions have some rays more likely than others. If these more-likely rays are (generally) where the BRDF is strongest, then this is called "importance sampling". The best function for diffuse is:$pdf(\vec{\omega}_i) = \frac{\left|\vec{N}_{\vec{x}} \cdot \vec{\omega}_i\right|}{\pi} ~~~~~~~~~~~~~~~(7)$Note that the $\pi$ here normalizes the PDF (since the integral of a PDF over its domain (here the unit hemisphere) must be $1$). Again, the absolute values do not affect the result since all dot products within the domain are positive.

By substituting the BRDF and the PDF into $(5)$ and simplifying a bit, one obtains:$\vec{k}_d \cdot \left|\vec{N}_{\vec{x}}\cdot\vec{\omega}_o\right| = \pi \cdot \vec{k}_d \cdot x \cdot E\left[ \left| \vec{N}_{\vec{x}} \cdot \vec{\omega}_o \right| \right] ~~~~~~~~~~~~~~~(8)$

Since $\vec{\omega}_o$ is not random, the expected value is just itself, and so almost all the rest cancels leaving $x=1/\pi$. Substituting this back in, for rays sampled according to the PDF above, the correct BRDF to use is:$\vec{f}_r(\vec{x},\vec{\omega}_i,\vec{\omega}_o) = \frac{ \left|\vec{N}_{\vec{x}} \cdot \vec{\omega}_o\right| \cdot \vec{k}_d }{\pi} ~~~~~~~~~~~~~~~(9)$

With $(4)$, $(6)$, and $(7)$ together, we have:$\vec{L}(\vec{x},\vec{\omega}_o) = \left[E\left[ \frac{ \vec{L}(\vec{y}_i,\,-\vec{\omega}_i) \left| \vec{N}_{\vec{x}} \cdot \vec{\omega}_i \right| \frac{\left|\vec{N}_{\vec{x}} \cdot \vec{\omega}_o\right| \cdot \vec{k}_d}{\pi} } { \frac{\left|\vec{N}_{\vec{x}} \cdot \vec{\omega}_i\right|}{\pi} }\right] + \overrightarrow{Em}(\vec{x},\vec{\omega}_o) \right] \frac{1}{\left|\vec{N}_{\vec{x}}\cdot\vec{\omega}_o\right|} ~~~~~~~~~~~~~~~(10)$A whole bunch of stuff cancels, and we have left:$\vec{L}(\vec{x},\vec{\omega}_o) = \vec{k}_d \cdot E\left[\vec{L}(\vec{y}_i,\,-\vec{\omega}_i)\right] + \overrightarrow{Em}(\vec{x},\vec{\omega}_o) \frac{1}{\left|\vec{N}_{\vec{x}}\cdot\vec{\omega}_o\right|} ~~~~~~~~~~~~~~~(11)$Equation $(11)$ is exactly what I reverse engineered from SmallPT. Like other sources, SmallPT doesn't have the reciprocal cosine scaling on the emission term explicitly. As I mentioned in my derivation of the rendering equation, linked above, it is already implicitly included: a diffuse emitter cancels out the scaling and becomes the constant term I always see everywhere.