Turns out my ray tracer uses a biased estimator, so I unbiased it.
The issue
When using Monte Carlo integration, the samples need to be weighted by the reciprocal of the probability of choosing that sample. This ensures that the expected value of the sample is the integrand, i.e. that it is unbiased:
$ \begin{aligned} I = \frac{g(x)}{p(x)} \end{aligned} $
Here, $g$ contains all terms contributing to the terms inside the integral in the rendering equation at the current vertex and $p$ is the probability of choosing $x$ as the sample.
Expanding $g$:
$ \begin{aligned} I = \frac{f(x,\omega_o,\omega_i)L_i(x,\omega_i)cos(\theta_i)}{p(x)} \end{aligned} $
where $f$ is the BRDF transport function, $L_i$ the incoming light at point $x$ from direction $\omega_i$, $\omega_o$ is the outgoing direction and $\theta_i$ is the angle between the surface normal and the incoming light direction.
Since I am only dealing with diffuse surfaces for now, $f$ is independent of $\omega_i$ and $\omega_o$.
My error was to simply assume that $\frac{f(x)}{p(x)}$ simply cancel out to give the material’s albedo. This however is not correct, as $f$ is in reality:
$ \begin{aligned} f(x) = \frac{\rho}{\pi} \end{aligned} $
Where $\rho$ is the albedo. The $\pi$ term is a result of energy conservation - it ensures that the total exiting light is $\rho$ times the total incoming light. For a more in depth dive on this check out the 9.2 Diffuse Reflection in the PBRT book.
Therefore the estimator should be
$ \begin{aligned} I = \frac{\frac{\rho}{\pi}L_i(x,\omega_i)cos(\theta_i)}{\frac{1}{2\pi}} = L_i(x,\omega_i)cos(\theta_i) * 2 \end{aligned} $
Which means I have been missing a factor of $2$ at each bounce, which can be understood as a result of the lambertian surface scattering light in a cosine weighted semi-spherical distribution, while my sampling is uniform on the semi-sphere.
The solution
Missing a factor of 2 at each bounce resulted in longer paths being underweighted; they each have a lower chance of being randomly chosen, so their values are more “important”.
Adding this factor to my code is exceedingly easy:
const auto weakening_factor = dot(normal, new_v);
const auto path_sampling_probability = 1 / (2*pi);
const auto new_v_radiance = 1 / pi;
const auto beta = weakening_factor * new_v_radiance / path_sampling_probability;
transmission = transmission * beta * diffuse_material->diffuse_reflectance;
Simply adding the above lines to the code path that handles DiffuseMaterial is all I needed.
Results
After rendering the same scene with the unbiased integrator the result is visibly different:
There is less contrast and the light seems to reach corners more.
Closing thoughts
This was an important step, since later more complex materials will require getting this formula right even more.
Next up I want to look at ways of reducing variance to be able to use less samples and support smaller light sources.