smallpt: ray scattering and Importance Sampling

Question

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?

Answer 1

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).