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?
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
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.
