<p>Let me first address some misconceptions you have:</p>

<blockquote>
  <p>the theory states that we need to shoot x number of rays for each
  intersection</p>
</blockquote>

<p>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 <span class="math-container">$\pi$</span> factor).</p>

<p>Now onto the bigger misconception you have:</p>

<blockquote>
  <p>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. </p>
</blockquote>

<p>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.</p>

<p>You start with a Fredholm integral equation of the second kind:</p>

<p><span class="math-container">$$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}$$</span></p>

<p>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:
<span class="math-container">$$L = L_e + TL$$</span>
<span class="math-container">$$(I-T)L = L_e$$</span>
<span class="math-container">$$L = (I-T)^{-1}L_e = \sum_{i=0}^{\infty}T^iL_e$$</span>
Where the last equality holds as long as <span class="math-container">$T$</span> is a contraction (<span class="math-container">$||T||&lt;1$</span>). Note that <span class="math-container">$T^0$</span> represents <span class="math-container">$0$</span> bounces, <span class="math-container">$T^1$</span> represents <span class="math-container">$1$</span> bounce, and <span class="math-container">$T^n$</span> represents <span class="math-container">$n$</span> bounces. Thus you can look at the solution as accumulating contributions from all light paths of length <span class="math-container">$0$</span>, <span class="math-container">$1$</span>, <span class="math-container">$2$</span> 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 <span class="math-container">$S$</span> rather than the upper hemisphere (this was actually how Kajiya derived it initially, not precisely the form below but in an equivalent form):
<span class="math-container">$$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)}$$</span>
Where <span class="math-container">$G$</span> is defined as:
<span class="math-container">$$G(x \leftrightarrow y) = \frac{|\cos\theta_x\cos\theta_y|}{||x-y||}V(x,y)$$</span>
Where <span class="math-container">$\theta_x$</span> is the angle between the vector <span class="math-container">$y-x$</span> and the normal at <span class="math-container">$x$</span>, and <span class="math-container">$\theta_y$</span> is the angle betwen the vector <span class="math-container">$x-y$</span> and the normal at <span class="math-container">$y$</span> (the cosines can be written as dot products provided you have the normals). <span class="math-container">$V(x,y)$</span> is the visibility function and it is 1 only if there's no surface point between <span class="math-container">$x$</span> and <span class="math-container">$y$</span> (if you shoot a ray from <span class="math-container">$x$</span> to <span class="math-container">$y$</span> the first thing you hit will be <span class="math-container">$y$</span>), and <span class="math-container">$0$</span> otherwise. The visibility function is implicit in the solid angle formulation inside <span class="math-container">$L_i$</span> (through the ray-casting function). Now let us rewrite the operator expansion in terms of what we have above, we get:
<span class="math-container">$$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)} + $$</span>
<span class="math-container">$$\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)}}$$</span>
<span class="math-container">$$+ ...$$</span>
This sum continues till infinity with iterated integrals corresponding to <span class="math-container">$T^i$</span>, the part expanded above is <span class="math-container">$T^0L_e + T^1L_e + T^2L_e$</span>. Note that usually you have an extra integration over the area of the pixel (then there's a factor <span class="math-container">$W$</span> 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 <span class="math-container">$n$</span>-iterated integral your samples are really paths starting from the camera with length <span class="math-container">$n$</span>: <span class="math-container">$x_0x_1x_2...x_{n+1}$</span> for a path of length <span class="math-container">$n$</span>. Note that you do not require that you split the path for each subsequent integral.
What path tracing implementations do, is create this path <span class="math-container">$x_0x_1...x_{n+1}$</span> incrementally, and reuse the path for earlier terms in the sum. That means that you start by picking some point <span class="math-container">$x_0$</span> inside a pixel, pick some point <span class="math-container">$x_1$</span> in the scene compute the value of the integrand of the term <span class="math-container">$T^0L_e$</span> for the path <span class="math-container">$x_0x_1$</span>, pick a point <span class="math-container">$x_2$</span> in the scene and compute the value of the integrand of the term <span class="math-container">$T^1L_e$</span> for the path <span class="math-container">$x_0x_1x_2$</span>, 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 <span class="math-container">$x_0x_1...x_{n+1}$</span>.</p>

<p>The takeaway from this is that <strong>path</strong>tracing works with <strong>paths</strong> 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).</p>

CC BY-SA 4.0
2019-04-01T14:35:45.377Z
Thank you so much for your amazing answer, it went even beyond my original question providing the mathematical framework.
maurock
CC BY-SA 4.0
2019-04-01T14:40:27.470Z
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.
lightxbulb
CC BY-SA 4.0
2019-04-01T14:44:11.920Z
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.
lightxbulb
CC BY-SA 4.0
2019-04-01T15:28:50.140Z
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.
maurock
CC BY-SA 4.0
2019-04-01T15:58:17.920Z
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 https://computergraphics.stackexchange.com/questions/8676/monte-carlo-importance-sampling/8679#8679
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).
lightxbulb
CC BY-SA 4.0
2019-04-01T16:21:10.643Z
Great, so if I understand correctly, the returned reflectance should be multiplied in each bounce for `pi`, since we are dividing the integrand you have just written by `cos(theta)\pi`. I mean, in the Diffusion case: `return obj.e + f.mult(radiance(Ray(x,d),depth,Xi)) * pi`;
maurock
CC BY-SA 4.0
2019-04-01T16:35:34.487Z
The $\pi$ (it is a constant after all) can be absorbed into the albedo depending on how you define your albedo, which I think is what they did in smallpt. But yes, your idea is correct if you consider $\pi$ to not be absorbed. You should also try to rewrite your raytracer with a while loop rather than a recursion to gain further insight into what actually happens (I would recommend also checking out the path formulation in vech's thesis, or the same thing in Iliyan's thesis: http://www.iliyan.com/publications/ThesisPhD).
lightxbulb
CC BY-SA 4.0
2019-04-01T17:34:09.827Z
I multiplied for π as I wrote above, but the image returned is almost completely white. As you said, π is probably already absorbed in the albedo. Thank you for those sources!
maurock
CC BY-SA 4.0
2019-04-01T18:16:16.217Z
If you multiplied it in smallpt, then as I explained there it is accounted for, so another multiplication simply yields a non-energy conserving brdf most likely.
lightxbulb
CC BY-SA 4.0
2019-04-02T09:04:13.990Z
Just for the sake of completeness: I don't think the `π` is absorbed in the albedo (at least, in the diffusion case). I think that, since the BRDF for the Lambertian surface is `reflectivity  / π `, and `pdf = cos(theta) / π`, `π ` cancels out along with `cos(theta)`.
maurock
CC BY-SA 4.0
2019-04-02T10:17:32.853Z
The point is that the brdf is $f = r/\pi$, so when you write $f * L$ you're really writing $r/\pi * L$ thus it hasn't really been canceled. If one is to be entirely precise if he considers the $\pi$ to be absorbed, then he must call it $f_{\pi}$ or something. This becomes especially problematic when you use other sampling strategies where the constant factor is not $\pi$, and then people often make mistakes because of this. I recommend things being explicit, to make it explicit in smallpt you'll just have to divide the albedo by $\pi$ where the material is defined, then you can mult by $\pi$.
lightxbulb
CC BY-SA 4.0
2019-04-01T13:58:41.347Z
2019-04-01T14:36:11.170Z
2019-04-01T14:36:11.170Z
lightxbulb
<h1>Introduction</h1>

<p>For those familiar with <b>smallpt</b> by Kevin Beason: <br>
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. </p>

<p><a href="https://i.sstatic.net/Sv5IS.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Sv5IS.png" alt="enter image description here"></a></p>

<h1>Problem</h1>

<p>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?</p>

CC BY-SA 4.0
2019-04-01T11:21:50.093Z
[this](https://computergraphics.stackexchange.com/questions/7933/path-tracer-sampling-dimensions-confusion) might help,
gallickgunner
CC BY-SA 4.0
2019-04-01T11:41:17.983Z
@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.
maurock
CC BY-SA 4.0
2019-04-01T11:45:56.457Z
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.
gallickgunner
CC BY-SA 4.0
2019-04-01T12:09:31.793Z
@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?
maurock
CC BY-SA 4.0
2019-04-01T12:26:51.210Z
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.
gallickgunner
<raytracing><pathtracing><c++>
smallpt: ray scattering and Importance Sampling
maurock
# smallpt: ray scattering and Importance Sampling

(Asked by: username_0 on 2019-04-01)

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

[![enter image description here](https://i.sstatic.net/Sv5IS.png)](https://i.sstatic.net/Sv5IS.png)

# 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?

**Tags:** raytracing, pathtracing, c++

> **Comments:**
> **username_1**: [this](https://computergraphics.stackexchange.com/questions/7933/path-tracer-sampling-dimensions-confusion) might help,  
> **username_0**: @username_1 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.  
> **username_1**: 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.  
> **username_0**: @username_1 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?  
> **username_1**: 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.  

---

## 1 answers

### Answer by username_2 on 2019-04-01. Score: 2 | Accepted 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 @username_1 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 @username_1 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 **path**tracing 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).

> **Comments:**
> **username_0**: Thank you so much for your amazing answer, it went even beyond my original question providing the mathematical framework.  
> **username_2**: 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.  
> **username_2**: 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.  
> **username_0**: 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.  
> **username_2**: 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 https://computergraphics.stackexchange.com/questions/8676/monte-carlo-importance-sampling/8679#8679 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).  
> **username_0**: Great, so if I understand correctly, the returned reflectance should be multiplied in each bounce for `pi`, since we are dividing the integrand you have just written by `cos(theta)\pi`. I mean, in the Diffusion case: `return obj.e + f.mult(radiance(Ray(x,d),depth,Xi)) * pi`;  
> **username_2**: The $\pi$ (it is a constant after all) can be absorbed into the albedo depending on how you define your albedo, which I think is what they did in smallpt. But yes, your idea is correct if you consider $\pi$ to not be absorbed. You should also try to rewrite your raytracer with a while loop rather than a recursion to gain further insight into what actually happens (I would recommend also checking out the path formulation in vech's thesis, or the same thing in Iliyan's thesis: http://www.iliyan.com/publications/ThesisPhD).  
> **username_0**: I multiplied for π as I wrote above, but the image returned is almost completely white. As you said, π is probably already absorbed in the albedo. Thank you for those sources!  
> **username_2**: If you multiplied it in smallpt, then as I explained there it is accounted for, so another multiplication simply yields a non-energy conserving brdf most likely.  
> **username_0**: Just for the sake of completeness: I don't think the `π` is absorbed in the albedo (at least, in the diffusion case). I think that, since the BRDF for the Lambertian surface is `reflectivity  / π `, and `pdf = cos(theta) / π`, `π ` cancels out along with `cos(theta)`.  
> **username_2**: The point is that the brdf is $f = r/\pi$, so when you write $f * L$ you're really writing $r/\pi * L$ thus it hasn't really been canceled. If one is to be entirely precise if he considers the $\pi$ to be absorbed, then he must call it $f_{\pi}$ or something. This becomes especially problematic when you use other sampling strategies where the constant factor is not $\pi$, and then people often make mistakes because of this. I recommend things being explicit, to make it explicit in smallpt you'll just have to divide the albedo by $\pi$ where the material is defined, then you can mult by $\pi$.  

---

