# Introduction

For those familiar with smallpt by Kevin Beason:
In path tracing, when rays hit a surface, they are scattered and the theta angle (the angle between the normal to the surface and the scattered ray) depends on a probability density function. For cosine importance sampling, this probability depends on cos(theta). This is the case implemented in the basic smallpt. Multiple rays should be scattered for each hitting point, and their contribution needs to be calculated to return the color value of the specific pixel.

# Problem

Looking at how rays are scattered, this doesn't actually happen in the smallpt code. Multiple rays per pixel are shot at slightly different positions. That means, they don't always land in the same spot on the surface. For this reason, the image posted above does not represent what is happening. Indeed, it is not true that multiple rays are scattered FROM EVERY HITTING POINT, but multiple rays are shot for each pixel. Is my interpretation correct?

• this might help, Commented Apr 1, 2019 at 11:21
• @gallickgunner It helped indeed, thanks! It is exactly what I asked. So, the theory states that we need to shoot x number of rays for each intersection, but in practice, it would be too computationally expensive. Commented Apr 1, 2019 at 11:41
• Yes try reading Cook's paper on "Distributed Ray Tracing" 1984. There he clearly explains that instead of adding X number of rays for every dimension (shadows, light etc) we distribute them pre hand in space. That way we don't need to shoot multiple rays as they are already distributed in space. Commented Apr 1, 2019 at 11:45
• @gallickgunner I am reading it right now, thank you. I also noticed that, in the smallpt, the cos(theta) of the BRDF is not considered. The color is calculated as Vec f = obj.c, while according to the literature it should be multiplied by cos(theta). Is it something that you noticed as well? Commented Apr 1, 2019 at 12:09
• That might be because cosing weighted hemispherical sampling has a pdf of cos. The cos in the rendering equation and the pdf gets cancelled. It's micro optimization process. So don't take smallpt as a reference. Smaller isn't always better. Go check up on PBRT or some other site like scratchapixel. Commented Apr 1, 2019 at 12:26

## 1 Answer

Let me first address some misconceptions you have:

the theory states that we need to shoot x number of rays for each intersection

No, the "theory" doesn't state such a thing. Note also that the paper @gallickgunner is referring to is inapplicable in this case since smallpt is based on the rendering equation and not the limited variant present in Cook's paper. Thus you actually need to refer to Kajiya's paper or Veach's path formulation. What @gallickgunner is right about is the canceling of the cosine term (but there's a missing $$\pi$$ factor).

Now onto the bigger misconception you have:

Multiple rays should be scattered for each hitting point, and their contribution needs to be calculated to return the color value of the specific pixel.

There's no such requirement, in fact Kajiya proved that scattering multiple rays at each bounce is bad since you get an exponential explosion of the ray count later down the bounce tree where the actual contribution is actually very small. In MC terms you perform splitting where the contribution is small thus wasting computational power which would be a lot more useful in reducing variance in other places. To understand why this is the case, one has to dig deeper into the theory of how raytracing works and how you estimate the rendering equation.

You start with a Fredholm integral equation of the second kind:

$$L_o(x,\omega_o) = L_e(x, \omega_o) + \int_{\Omega}{f(\omega_o,x,\omega_i)L_i(x,\omega_i)\cos\theta\,d\omega_i}$$

As long as the integral converges (which happens when the operator is a contraction, or when it is not but enough energy can escape the scene for there to be a steady state), you can find the solution as the Liouville-Neumann expansion. To do this however I'll first rewrite the equation in operator form: $$L = L_e + TL$$ $$(I-T)L = L_e$$ $$L = (I-T)^{-1}L_e = \sum_{i=0}^{\infty}T^iL_e$$ Where the last equality holds as long as $$T$$ is a contraction ($$||T||<1$$). Note that $$T^0$$ represents $$0$$ bounces, $$T^1$$ represents $$1$$ bounce, and $$T^n$$ represents $$n$$ bounces. Thus you can look at the solution as accumulating contributions from all light paths of length $$0$$, $$1$$, $$2$$ etc. till infinity. Now to actually write this out explicitly it's more convenient to rewrite the rendering equation in terms of integration over all scene surface points $$S$$ rather than the upper hemisphere (this was actually how Kajiya derived it initially, not precisely the form below but in an equivalent form): $$L_o(x_1\rightarrow x_0) = L_e(x_1 \rightarrow x_0) + \int_{S}{f(x_2 \rightarrow x_1 \rightarrow x_0)L_o(x_2 \rightarrow x_1) G(x_1 \leftrightarrow x_2)\,dA(x_2)}$$ Where $$G$$ is defined as: $$G(x \leftrightarrow y) = \frac{|\cos\theta_x\cos\theta_y|}{||x-y||}V(x,y)$$ Where $$\theta_x$$ is the angle between the vector $$y-x$$ and the normal at $$x$$, and $$\theta_y$$ is the angle betwen the vector $$x-y$$ and the normal at $$y$$ (the cosines can be written as dot products provided you have the normals). $$V(x,y)$$ is the visibility function and it is 1 only if there's no surface point between $$x$$ and $$y$$ (if you shoot a ray from $$x$$ to $$y$$ the first thing you hit will be $$y$$), and $$0$$ otherwise. The visibility function is implicit in the solid angle formulation inside $$L_i$$ (through the ray-casting function). Now let us rewrite the operator expansion in terms of what we have above, we get: $$L_o(x_1\rightarrow x_0) = L_e(x_1 \rightarrow x_0) + \int_{S}{f(x_2 \rightarrow x_1 \rightarrow x_0)L_e(x_2 \rightarrow x_1) G(x_1 \leftrightarrow x_2)\,dA(x_2)} +$$ $$\int_{S}\int_{S}{{L_e(x_3 \rightarrow x_2)f(x_3 \rightarrow x_2 \rightarrow x_1) G(x_3 \leftrightarrow x_2)f(x_2 \rightarrow x_1 \rightarrow x_0) G(x_2 \leftrightarrow x_1)\,dA(x_3)\,dA(x_2)}}$$ $$+ ...$$ This sum continues till infinity with iterated integrals corresponding to $$T^i$$, the part expanded above is $$T^0L_e + T^1L_e + T^2L_e$$. Note that usually you have an extra integration over the area of the pixel (then there's a factor $$W$$ to account for the sensor sensitivity). When you do pathtracing (what smallpt does) you're actually estimating this infinite sum of integrals through Monte Carlo (this is also the reason why Monte Carlo is employed, since the convergence doesn't depend on the dimension, which in this case is infinite). You can estimate each term in the sum by Monte Carlo integration, that is for a $$n$$-iterated integral your samples are really paths starting from the camera with length $$n$$: $$x_0x_1x_2...x_{n+1}$$ for a path of length $$n$$. Note that you do not require that you split the path for each subsequent integral. What path tracing implementations do, is create this path $$x_0x_1...x_{n+1}$$ incrementally, and reuse the path for earlier terms in the sum. That means that you start by picking some point $$x_0$$ inside a pixel, pick some point $$x_1$$ in the scene compute the value of the integrand of the term $$T^0L_e$$ for the path $$x_0x_1$$, pick a point $$x_2$$ in the scene and compute the value of the integrand of the term $$T^1L_e$$ for the path $$x_0x_1x_2$$, and so on. In an implementation you usually use the solid angle formulation (at least for bsdf sampling), but there's not major difference, you still end up with paths of the form $$x_0x_1...x_{n+1}$$.

The takeaway from this is that pathtracing works with paths where, a path starts at the camera and ends at a lightsource. Using different paths that do not split (in general) to estimate the sum of integrals is usually more efficient than splitting each path into more paths at each bounce. I am not saying there's no place for splitting (actually smallpt uses it for the first few bounces for materials that refract), you just have to know when to apply it (when the contribution is large and is inadequately sampled by only 1 or a few samples). Note that the opposite of splitting is early path termination in the form of the Russian roullete - basically kill off paths with small contributions since it's not worth it wasting computational resources on those (it also lets you make your estimator unbiased, but in practice there's usually a ceiling to the number of bounces allowed, since for most scenes your estimator being "fully unbiased" is not worth it in terms of performance).

• Thank you so much for your amazing answer, it went even beyond my original question providing the mathematical framework. Commented Apr 1, 2019 at 14:35
• I don't think your question can really be answered (in a satisfactory manner) without mathematics, unless one just states that it's less efficient without explaining the basis for this. The basis is that you assume energy conservation of the brdf (otherwise the operator won't be a contraction), then with each bounce the contribution will be attenuated more and more, so you'd expected that bounces later down the bounce tree have a very small weight, thus you don't want to split there. The other obvious reason is that splitting results in exponential ray count growth. Commented Apr 1, 2019 at 14:40
• One cannot really talk about using Monte Carlo to solve the rendering equation without understanding how this is done. What my answer does is explain (albeit not in very many details) to what reformulation of the problem Monte Carlo is applied, which naturally address the misconception you had about the necessity of splitting. Note that both Whitted and Cook (and some other researchers), has the exact same idea you had, it was ultimately Kajiya that showed that this idea is not optimal (nor is it necessary for things to work out). My point is that your misunderstanding was in some way natural. Commented Apr 1, 2019 at 14:44
• It is perfectly clear now, thank you. Regarding the canceling of the cosine term in the BRDF, is there a source you know that provide the math for this micro-optimization? I have just calculated on paper the x,y and z coordinates both for uniform sampling, and for p(w) = cos(theta) / pi depending on two random uniform generators, and I obtained the same equations that smallpt already implements for the second case. Commented Apr 1, 2019 at 15:28
• There's not much math, you have the integrand piece for the current bounce: $f(x_{i+1}\rightarrow x_{i} \rightarrow x_{i-1})L_e(x_{i+1}\rightarrow x_i)\cos\theta$ which you must divide by the pdf as I have explained here computergraphics.stackexchange.com/questions/8676/… Since you have the cosine in the integrand and the cosine in the pdf they cancel out. In fact you should cancel them out otherwise you may get numerical imprecision or NaNs especially for the cosine equal or close to $0$ (remember computers have finite precision). Commented Apr 1, 2019 at 15:58