# 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, 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. 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. 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? 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. Apr 1, 2019 at 12:26

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.

$$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}$$.
• 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. 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). Apr 1, 2019 at 15:58