<p>This is a very good question. There is a common misconception that Monte Carlo, or integration is applied "recursively" on the rendering equation. That is not what's happening. Numerical integration methods are tailored to problems of the form:</p>

<p><span class="math-container">$$I = \int_{\Omega}{f(x)d\mu(x)} \approx \sum_{k=0}^{N-1}w(x_k)f(x_k)$$</span></p>

<p>Note that this is not the case for the rendering equation. Namely the fact that it is an equation, and the unknown radiance function <span class="math-container">$L$</span> shows up on both sides. In order to apply numerical integration techniques (such as Monte Carlo), one has to bring the rendering equation to the above form. Fortunately this is a well studied problem. The rendering equation is a specific case of a Fredholm integral equation of the second kind. If one can show that the scattering operator in that equation is a contraction <span class="math-container">$||T||&lt;1$</span> (which holds if brdfs are energy conserving, there are cases where this is not required but that goes outside the scope of my answer), then we may perform a Neumann expansion (<a href="https://en.wikipedia.org/wiki/Liouville%E2%80%93Neumann_series" rel="noreferrer">https://en.wikipedia.org/wiki/Liouville%E2%80%93Neumann_series</a>).</p>

<p>The commonly used form of the rendering equation is the solid angle formulation (here <span class="math-container">$\sigma$</span> is the solid angle measure, <span class="math-container">$\sigma(\omega) = \sin\theta\,d\theta d\phi$</span>):</p>

<p><span class="math-container">$$L(x,\omega_o) = L_e(x,\omega_o) + \int_{\Omega}{f(\omega_o,x,\omega_i)L_i(x,\omega_i)\cos\theta_i\,d\sigma(\omega_i)}$$</span></p>

<p>In order to perform the expansion, I will rewrite this in the area formulation (here <span class="math-container">$\mu(x)$</span> is the Lebesgue area measure):</p>

<p><span class="math-container">$$L(x_1 \rightarrow x_0) = L_e(x1 \rightarrow x_0) + \int_{M}{f(x_2 \rightarrow x_1 \rightarrow x_0)L(x_2 \rightarrow x_1)\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)}$$</span></p>

<p>Where <span class="math-container">$M$</span> is the set of all surface points in the scene, <span class="math-container">$V(x,y)$</span> is the visibility function which is <span class="math-container">$1$</span> if there is nothing between <span class="math-container">$x$</span> and <span class="math-container">$y$</span> and <span class="math-container">$0$</span> otherwise. If the normals of the surfaces at point <span class="math-container">$x_1$</span> and <span class="math-container">$x_2$</span> are respectively <span class="math-container">$N_{x_1}$</span> and <span class="math-container">$N_{x_2}$</span>, then <span class="math-container">$\cos\theta_{x_1} = N_{x_1} \cdot \frac{x_2-x1}{||x_2-x_1||}$</span> and <span class="math-container">$\cos\theta_{x_2} = N_{x_2} \cdot \frac{x_1-x2}{||x_2-x_1||}$</span>. As for the radiance, <span class="math-container">$L(x_2 \rightarrow x_1)$</span> gives the radiance arriving at <span class="math-container">$x_1$</span> from the direction of <span class="math-container">$x_2$</span>. And the brdf notation relationship is: <span class="math-container">$f(x_2 \rightarrow x_1 \rightarrow x_0) = f(x_1 \rightarrow x_0, x_1, x_1 \rightarrow x_2)$</span> (that means <span class="math-container">$\omega_o = x_1 \rightarrow x_0$</span> and <span class="math-container">$\omega_i = x_1 \rightarrow x_2$</span>).</p>

<p>I will rewrite the above for simplicity and conciseness as:</p>

<p><span class="math-container">$$L(x_1 \rightarrow x_0) = L_e(x1 \rightarrow x_0) + L_r(x_1 \rightarrow x_0)$$</span></p>

<p>Let us now split the incoming radiance into direct illumination arriving at <span class="math-container">$x_1$</span> and indirect (at least one bounce) illumination. The direct illumination is obviously due to direct light source rays arriving at <span class="math-container">$x_1$</span>:</p>

<p><span class="math-container">$$L(x_1 \rightarrow x_0) = L_e(x1 \rightarrow x_0) + \int_{M}{f(x_2 \rightarrow x_1 \rightarrow x_0)(L_e(x_2 \rightarrow x_1) + L_r(x_2 \rightarrow x_1))\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)} = 
L_e(x1 \rightarrow x_0) + \int_{M}{f(x_2 \rightarrow x_1 \rightarrow x_0)L_e(x_2 \rightarrow x_1)\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)} + $$</span></p>

<p><span class="math-container">$$\int_{M}f(x_2 \rightarrow x_1 \rightarrow x_0)\Bigg[\int_{M}f(x_3 \rightarrow x_2 \rightarrow x_1)L(x_3 \rightarrow x_2)\cos\theta_{x_2}$$</span></p>

<p><span class="math-container">$$\frac{\cos\theta_{x_3}}{||x_3-x_2||^2}V(x_3,x_2)\,d\mu(x_3)\Bigg]\cos\theta_{x_1}\frac{\cos\theta_{x_2}}{||x_2-x_1||^2}V(x_2,x_1)\,d\mu(x_2)$$</span></p>

<p>All I did here was separate the sum <span class="math-container">$L_e(x_2 \rightarrow x_1) + L_r(x_2 \rightarrow x_1)$</span> into two integrals, and then I additionally expanded <span class="math-container">$L_r$</span> into the corresponding integral using the recursive definition. You can expand this until infinity. For conciseness, I will rewrite the integration as an operator: <span class="math-container">$L_r = TL$</span>. Now we can see the relationship to the Neumann expansion. Since the rendering equation can be written as:</p>

<p><span class="math-container">$$L = L_e + TL$$</span></p>

<p>The solution is formally given as: </p>

<p><span class="math-container">$$(I-T)L = L_e$$</span>
<span class="math-container">$$L = (I-T)^{-1}L_e$$</span></p>

<p>Applying the Neumann expansion yields:</p>

<p><span class="math-container">$$L = \sum_{i=0}^{\infty}T^iL_e$$</span></p>

<p>Note that each term of the sum is an increasingly dimensional integral. The first term is obviously just <span class="math-container">$L_e(x_1 \rightarrow x_0)$</span>, the second term <span class="math-container">$TL_e$</span> is the integral for <span class="math-container">$L_e(x_2 \rightarrow x_2)$</span> that I wrote when I split <span class="math-container">$L_e$</span> and <span class="math-container">$L_r$</span>, and so on. What each term gives you is the energy coming from a light source after <span class="math-container">$i$</span> bounces. The first term gives you the radiance emitted from <span class="math-container">$x_1$</span> towards <span class="math-container">$x_0$</span> (paths of length 0), the second term gives you the radiance that is due to direct illumination of <span class="math-container">$x_1$</span> scattered towards <span class="math-container">$x_0$</span> (paths of length 1), then you have the radiance due to lights two bounces away (paths of length 2) etc. With each bounce each new integral in the sum gains two dimensions (if we are just integrating over incoming directions). I would also like to emphasize that the integration variable <span class="math-container">$x_2$</span> in <span class="math-container">$TL_e$</span> is not the same <span class="math-container">$x_2$</span> as the one in <span class="math-container">$T^iL_e$</span> (this means that in the example with integration, <span class="math-container">$x_2$</span> in the first integral is not the same as in the second one - they are just integration variables). Another important fact is that within each sum each path starts at the film <span class="math-container">$x_0$</span> and ends at a light source (<span class="math-container">$L_e(x_2\rightarrow x_1)$</span> is non-zero only for points <span class="math-container">$x_2$</span> lying on light sources).</p>

<p>Now that we have a sum of integrals, we can apply our numerical techniques to each integral to estimate the sum. An obvious optimization, is reusing the samples used to compute <span class="math-container">$T^iL_e$</span> to compute <span class="math-container">$T^{i+1}L_e = T(T^iL_e)$</span>. This yields a formulation where we can formally write out the sum of integrals as one integral with integration over all paths that start at the camera film and end at a light source. This also illustrates why Monte Carlo is a preferred technique: since we want to estimate an infinite dimensional integral, and Monte Carlo's convergence doesn't depend on the dimensionality unlike standard quadrature rules (there's also the point that it doesn't care about the smoothness of the integrand, but that's a double edged sword). </p>

<p>As you may have noticed, with the path integration formulation an exponential path count growth is not required. Splitting a path into multiple new ones is a technique known as splitting. While it is beneficial for high energy paths, on average with each bounce the energy gets lower due to attenuation, so this is in most cases counterproductive. On the other hand early path termination (for low energy paths) is more often than not beneficial in terms of efficiency (that's simply cutting off the sum at some point and not computing the remaining infinitely many terms, because you deem they have a contribution that is too low), and that's where Russian roulette comes in. </p>

<p>All of this is actually contained in Kajiya's paper on the rendering equation, even though he doesn't go into the details I went. He just refers to Rubinstein's book: <a href="https://dl.acm.org/citation.cfm?id=539488" rel="noreferrer">https://dl.acm.org/citation.cfm?id=539488</a> (note - the first edition, later editions do not have the part he is referring to). In this book how to solve Fredholm integrals of the second kind is described (what I explained above, but more formally).</p>

<p>I hope my explanation was useful, and will continue being useful for future readers.</p>

CC BY-SA 4.0
2019-07-19T09:14:10.740Z
2019-07-19T19:14:59.497Z
2019-07-19T19:14:59.497Z
lightxbulb
<p>Path tracing seems to perform Monte Carlo integration where the space being integrated over is the space of paths (e.g. chains of rays) but the rendering equation is stated in terms of outgoing rays. Naively you would have to perform the integration at each bounce to model the rendering equation but this would take exponential time. It seems that there should be a proof that either a) integration over the space of paths is equivalent to recursive integration over the space of rays or b) sampling paths is a valid numeric technique for evaluating the rendering equation. While somewhat intuitive neither of those facts is obvious to me. Is there a proof of some such claim?</p>

<pathtracing><monte-carlo><importance-sampling>
Rendering equation in terms of paths rather than directions
Jake
# Rendering equation in terms of paths rather than directions

(Asked by: username_0 on 2019-07-19)

Path tracing seems to perform Monte Carlo integration where the space being integrated over is the space of paths (e.g. chains of rays) but the rendering equation is stated in terms of outgoing rays. Naively you would have to perform the integration at each bounce to model the rendering equation but this would take exponential time. It seems that there should be a proof that either a) integration over the space of paths is equivalent to recursive integration over the space of rays or b) sampling paths is a valid numeric technique for evaluating the rendering equation. While somewhat intuitive neither of those facts is obvious to me. Is there a proof of some such claim?

**Tags:** pathtracing, monte-carlo, importance-sampling


---

## 1 answers

### Answer by username_1 on 2019-07-19. Score: 8 | Accepted answer

This is a very good question. There is a common misconception that Monte Carlo, or integration is applied "recursively" on the rendering equation. That is not what's happening. Numerical integration methods are tailored to problems of the form:

$$I = \int\_{\Omega}{f(x)d\mu(x)} \approx \sum\_{k=0}^{N-1}w(x\_k)f(x\_k)$$

Note that this is not the case for the rendering equation. Namely the fact that it is an equation, and the unknown radiance function $L$ shows up on both sides. In order to apply numerical integration techniques (such as Monte Carlo), one has to bring the rendering equation to the above form. Fortunately this is a well studied problem. The rendering equation is a specific case of a Fredholm integral equation of the second kind. If one can show that the scattering operator in that equation is a contraction $||T||<1$ (which holds if brdfs are energy conserving, there are cases where this is not required but that goes outside the scope of my answer), then we may perform a Neumann expansion (<https://en.wikipedia.org/wiki/Liouville%E2%80%93Neumann_series>).

The commonly used form of the rendering equation is the solid angle formulation (here $\sigma$ is the solid angle measure, $\sigma(\omega) = \sin\theta\,d\theta d\phi$):

$$L(x,\omega\_o) = L\_e(x,\omega\_o) + \int\_{\Omega}{f(\omega\_o,x,\omega\_i)L\_i(x,\omega\_i)\cos\theta\_i\,d\sigma(\omega\_i)}$$

In order to perform the expansion, I will rewrite this in the area formulation (here $\mu(x)$ is the Lebesgue area measure):

$$L(x\_1 \rightarrow x\_0) = L\_e(x1 \rightarrow x\_0) + \int\_{M}{f(x\_2 \rightarrow x\_1 \rightarrow x\_0)L(x\_2 \rightarrow x\_1)\cos\theta\_{x\_1}\frac{\cos\theta\_{x\_2}}{||x\_2-x\_1||^2}V(x\_2,x\_1)\,d\mu(x\_2)}$$

Where $M$ is the set of all surface points in the scene, $V(x,y)$ is the visibility function which is $1$ if there is nothing between $x$ and $y$ and $0$ otherwise. If the normals of the surfaces at point $x\_1$ and $x\_2$ are respectively $N\_{x\_1}$ and $N\_{x\_2}$, then $\cos\theta\_{x\_1} = N\_{x\_1} \cdot \frac{x\_2-x1}{||x\_2-x\_1||}$ and $\cos\theta\_{x\_2} = N\_{x\_2} \cdot \frac{x\_1-x2}{||x\_2-x\_1||}$. As for the radiance, $L(x\_2 \rightarrow x\_1)$ gives the radiance arriving at $x\_1$ from the direction of $x\_2$. And the brdf notation relationship is: $f(x\_2 \rightarrow x\_1 \rightarrow x\_0) = f(x\_1 \rightarrow x\_0, x\_1, x\_1 \rightarrow x\_2)$ (that means $\omega\_o = x\_1 \rightarrow x\_0$ and $\omega\_i = x\_1 \rightarrow x\_2$).

I will rewrite the above for simplicity and conciseness as:

$$L(x\_1 \rightarrow x\_0) = L\_e(x1 \rightarrow x\_0) + L\_r(x\_1 \rightarrow x\_0)$$

Let us now split the incoming radiance into direct illumination arriving at $x\_1$ and indirect (at least one bounce) illumination. The direct illumination is obviously due to direct light source rays arriving at $x\_1$:

$$L(x\_1 \rightarrow x\_0) = L\_e(x1 \rightarrow x\_0) + \int\_{M}{f(x\_2 \rightarrow x\_1 \rightarrow x\_0)(L\_e(x\_2 \rightarrow x\_1) + L\_r(x\_2 \rightarrow x\_1))\cos\theta\_{x\_1}\frac{\cos\theta\_{x\_2}}{||x\_2-x\_1||^2}V(x\_2,x\_1)\,d\mu(x\_2)} =
L\_e(x1 \rightarrow x\_0) + \int\_{M}{f(x\_2 \rightarrow x\_1 \rightarrow x\_0)L\_e(x\_2 \rightarrow x\_1)\cos\theta\_{x\_1}\frac{\cos\theta\_{x\_2}}{||x\_2-x\_1||^2}V(x\_2,x\_1)\,d\mu(x\_2)} + $$

$$\int\_{M}f(x\_2 \rightarrow x\_1 \rightarrow x\_0)\Bigg[\int\_{M}f(x\_3 \rightarrow x\_2 \rightarrow x\_1)L(x\_3 \rightarrow x\_2)\cos\theta\_{x\_2}$$

$$\frac{\cos\theta\_{x\_3}}{||x\_3-x\_2||^2}V(x\_3,x\_2)\,d\mu(x\_3)\Bigg]\cos\theta\_{x\_1}\frac{\cos\theta\_{x\_2}}{||x\_2-x\_1||^2}V(x\_2,x\_1)\,d\mu(x\_2)$$

All I did here was separate the sum $L\_e(x\_2 \rightarrow x\_1) + L\_r(x\_2 \rightarrow x\_1)$ into two integrals, and then I additionally expanded $L\_r$ into the corresponding integral using the recursive definition. You can expand this until infinity. For conciseness, I will rewrite the integration as an operator: $L\_r = TL$. Now we can see the relationship to the Neumann expansion. Since the rendering equation can be written as:

$$L = L\_e + TL$$

The solution is formally given as:

$$(I-T)L = L\_e$$
$$L = (I-T)^{-1}L\_e$$

Applying the Neumann expansion yields:

$$L = \sum\_{i=0}^{\infty}T^iL\_e$$

Note that each term of the sum is an increasingly dimensional integral. The first term is obviously just $L\_e(x\_1 \rightarrow x\_0)$, the second term $TL\_e$ is the integral for $L\_e(x\_2 \rightarrow x\_2)$ that I wrote when I split $L\_e$ and $L\_r$, and so on. What each term gives you is the energy coming from a light source after $i$ bounces. The first term gives you the radiance emitted from $x\_1$ towards $x\_0$ (paths of length 0), the second term gives you the radiance that is due to direct illumination of $x\_1$ scattered towards $x\_0$ (paths of length 1), then you have the radiance due to lights two bounces away (paths of length 2) etc. With each bounce each new integral in the sum gains two dimensions (if we are just integrating over incoming directions). I would also like to emphasize that the integration variable $x\_2$ in $TL\_e$ is not the same $x\_2$ as the one in $T^iL\_e$ (this means that in the example with integration, $x\_2$ in the first integral is not the same as in the second one - they are just integration variables). Another important fact is that within each sum each path starts at the film $x\_0$ and ends at a light source ($L\_e(x\_2\rightarrow x\_1)$ is non-zero only for points $x\_2$ lying on light sources).

Now that we have a sum of integrals, we can apply our numerical techniques to each integral to estimate the sum. An obvious optimization, is reusing the samples used to compute $T^iL\_e$ to compute $T^{i+1}L\_e = T(T^iL\_e)$. This yields a formulation where we can formally write out the sum of integrals as one integral with integration over all paths that start at the camera film and end at a light source. This also illustrates why Monte Carlo is a preferred technique: since we want to estimate an infinite dimensional integral, and Monte Carlo's convergence doesn't depend on the dimensionality unlike standard quadrature rules (there's also the point that it doesn't care about the smoothness of the integrand, but that's a double edged sword).

As you may have noticed, with the path integration formulation an exponential path count growth is not required. Splitting a path into multiple new ones is a technique known as splitting. While it is beneficial for high energy paths, on average with each bounce the energy gets lower due to attenuation, so this is in most cases counterproductive. On the other hand early path termination (for low energy paths) is more often than not beneficial in terms of efficiency (that's simply cutting off the sum at some point and not computing the remaining infinitely many terms, because you deem they have a contribution that is too low), and that's where Russian roulette comes in.

All of this is actually contained in Kajiya's paper on the rendering equation, even though he doesn't go into the details I went. He just refers to Rubinstein's book: <https://dl.acm.org/citation.cfm?id=539488> (note - the first edition, later editions do not have the part he is referring to). In this book how to solve Fredholm integrals of the second kind is described (what I explained above, but more formally).

I hope my explanation was useful, and will continue being useful for future readers.


---

