Multiple Importance Sampling in Path tracer produces Dark Images

Question

So I recently implemented Multiple Importance Sampling in my path tracer which was based on next event estimation.

The problem is without MIS I get images like,

This is obtained by setting
light_sample *= 1/light_pdf; // Note no mis_weight
and just returning this sample alone. No BRDF sampling.

Where as with MIS as defined above I get darker images like,

The reason is surely the mis_weight factor. From what I gathered over the internet my MIS code is ok but theoretically it doesn't seem right. For example suppose we performed Light Sampling first and obtained a mis_weight. After that when we tried BRDF sampling, the ray didn't intersect any light source resulting in light_pdf = 0. We neglected this estimate. Since we neglected this estimate isn't it wrong to weight our light sample by mis_weight when we didn't even use multiple estimators, plus how is the sum of weights equal to 1 when we didnt even estimate anything using the BRDF pdf?

Imo, mis_weight factor in light_sample should only be used when BRDF sampling also results in a light ray that intersects the light source.

Can anybody explain if these results are ok or there is something wrong with the code?

EDIT:- There is another case which is a little confusing. I am currently using a heuristic (basically power/distance) to choose 1 light out of multiple lights. What if Light sampling and BRDF sampling choose different light sources. Is MIS still valid then? Since the PDF would change resulting the weights not being able to sum to 1 anymore

EDIT2:- So Stefan pointed out a mistake here, that I was clamping radiance. And this solved that issue and changed the results I'm getting for MIS. I have double checked PBRT's implementation and it is similar to mine. I'm still getting darker images using MIS however now I'm getting more fireflies and I think I read that MIS reduces fireflies? Updated the images. It seems fireflies are less in MIS however, the reflection seems harder to converge.

Note the reflection is harder to converge. Both images are 1000 spp. Top is without MIS.

I'm removing the snippet and adding a link to github for the whole kernel. The core functions are evaluateDirectLighting, shading, sampleLights. Link to code is "here"

--Found a mistake. I was using the direct lighting equation when calculating the brdf sample. (multiplying by $\cos(\theta^{\prime})$ and dividing by $r^2$). Removed it since we initially sampled through BSDF which is integral over solid angle not area. The image got a little brighter. Still don't know if MIS is working as intended. See answer.

EDIT3:- Added a release as well. So if anybody wants, they can try changing the code in "*.cl" file and run the program to see the results. (You must have an OpenCL 1.1 supported GPU or CPU)

EDIT4:- Here's an overview of what I'm doing now.

First I choose a single light source out of multiple using simple heuristic scheme like the distance, intensity, area and the cosine falloff angles. I appropriately set the weights for the light like this

$light\_pdf = weight/area$

where weight is in range 0-1.

Next I trace the ray to see if light source is visible. If it is I calculate the $light\_sample$ using the direct lighting equation (integral over Area).

Then I calculate the BRDF PDF for this given ray. However I use Lafortune's algorithm for it. If a random number falls under the specular color I sample through the modfied Phong PDF else through Cosine.

The weights are computed using power heuristic, $ weight = light\_pdf^2/(light\_pdf^2+brdf\_pdf^2)$

The MIS estimator is then calculated as

$light\_sample = light\_sample * weight/light\_pdf$

After that I come to BRDF sampling. I again sample through either Phong/Cosine based on what I did earlier during light sample calculation. If the sampled ray doesn't hit any light source or not the same one. I set the $brdf\_sample$ to zero. If it hits, I set the $light\_pdf$ to same value as before. Calculate the weights like mentioned above and calculate brdf sample using the original equation (integral over solid angle).

EDIT5:- After lightxbulb suggestions, I think the problem has resolved.

Note the images might look a whole lot different but that's cause I implemented Tonemapping + Gamma Correction in the meantime :) With MIS 700spp

Without MIS 700spp

Answer 1

Throughout my answer I'll sometimes refer to some results in https://sites.fas.harvard.edu/~cs278/papers/veach.pdf by using [MIS,section_number].

You can skip the following derivation if you don't care about the mathematical explanation of why using MIS to combine estimators is valid. I'll have to start with what the purpose of MIS is. The general idea is that you want to estimate some integral $I=\int_{\Omega}{f(x)\,d\mu(x)}$ through Monte Carlo, and you have various sampling pdfs (in our case a pdf for sampling points on the light, and a pdf for sampling directions from the brdf). Additionally you want to sample using more than one technique and then combine the contributions optimally (in some sense). For simplicity I will restrict myself to the case where you have only two pdfs: $p_L(x), p_B(x)$ with respect to the same measure $\mu(x)$.

Assume that you also have the weighting functions $w_1(x), w_2(x) : \Omega \rightarrow [0,1], w_1(x) + w_2(x) = 1$ and also $n$ independently and identically distributed (i.i.d.) samples $x_{p_L,i}: i=1,...,n$ according to $p_L$, and $m$ i.i.d. samples $x_{p_B,j} : j=1,...,m$ according to $p_B$, we can use Monte Carlo with MIS to estimate the integral $I$: $$ I = \int_{\Omega}{f(x)\,d\mu(x)} = \int_{\Omega}{f(x)(w_1(x)+w_2(x))\,d\mu(x)} \\ = \int_{\Omega}{\frac{w_1(x)f(x)}{p_L(x)}p_L(x)\,d\mu(x)} + \int_{\Omega}{\frac{w_2(x)f(x)}{p_B(x)}p_B(x)\,d\mu(x)} \\ = E\left[\frac{w_1(X_{p_L})f(X_{p_L})}{p_L(X_{p_L})}\right] + E\left[\frac{w_2(X_{p_B})f(X_{p_B})} {p_B(X_{p_B})}\right] \\ = \frac{1}{n}\sum_{i=1}^{n}{E\left[\frac{w_1(X_{p_L,i})f(X_{p_L,i})}{p_L(X_{p_L,i})}\right]} + \frac{1}{m}\sum_{j=1}^{m}{E\left[\frac{w_2(X_{p_B,j})f(X_{p_B,j})} {p_B(X_{p_B,j})}\right]} \\ = E\left[\frac{1}{n}\sum_{i=1}^{n}{\frac{w_1(X_{p_L,i})f(X_{p_L,i})}{p_L(X_{p_L,i})}}\right] + E\left[\frac{1}{m}\sum_{j=1}^{m}{\frac{w_2(X_{p_B,j})f(X_{p_B,j})} {p_B(X_{p_B,j})}}\right] \\ \approx \frac{1}{n}\sum_{i=1}^{n}{\frac{w_1(x_{p_L,i})f(x_{p_L,i})}{p_L(x_{p_L,i})}} + \frac{1}{m}\sum_{j=1}^{m}{\frac{w_2(x_{p_B,j})f(x_{p_B,j})} {p_B(x_{p_B,j})}} $$ Where the second equality holds because $w_1(x)+w_2(x)=1$, the fourth follows from the definition of expected value, the fifth holds since $X_{p_L,i}$ are independently and identically distributed (i.i.d.) according to $p_L$ and similarly $X_{p_B,i}$ are i.i.d. according to $p_B$. The sixth holds because of the properties of the expectation (which follows from the fact that integration is a linear operator). The last approximate equality follows from the strong law of large numbers (the average of the samples converges almost surely to the expected value with the number of samples going to infinity).

Having this result the first thing to note is that both pdfs are defined with respect to the same measure $\mu(x)$. However I believe that in your code and overview you pdf for sampling the light $p_L$ is defined with respect to the area measure, whereas the pdf for sampling the brdf $p_B$ is defined with respect to the solid angle measure. The relationship between those being $d\sigma(x\rightarrow y) = \frac{\cos\theta_y}{||x-y||^2}dA(y)$, where $x \in \mathbb{R}^3$ is the current point for which you are computing the illumination (that is your intersection point), and $y \in \mathbb{R}^3$ is a point on the surface of some light source. Additionally $(x\rightarrow y) = \frac{y-x}{||y-x||}$ is simply the unit direction vector from $x$ to $y$, $\cos\theta_y = n_y \cdot (y\rightarrow x)$ is the cosine of the angle between the normal $n_y$ (also unit length) of the light surface at $y$ and $-(x\rightarrow y)$. More specifically what you have produced as light pdf is really the probability of picking the light $weights[i]/sum$ multiplied by the probability of picking uniformly a point on the light source that you have chosen $1/area$. The solid angle to area (or vice versa) conversion seems to be missing in your code and your overview, that is your pdf is with respect to the area measure, and you combine it with a pdf (the brdf's) that is with respect to the solid angle measure, which is clearly wrong as you can see from my derivation above. To get the brdf pdf with respect to the area measure and not the solid angle measure you can use: $p_{\omega}(\omega)\,d\sigma(\omega) = p_A(y)\,dA(y)$ then using the relationship between the measures $p_A(y) = p_{\omega}(\omega)\,d\sigma(\omega) / dA(y) = p_{\omega}(\omega) \frac{\cos\theta_y}{||x-y||^2}$. You can refer to [MIS, 2.3] equation (9). For a formal derivation of the relationship between the two measures you can refer to http://www.dgp.toronto.edu/~lessig/dissertation/files/area_formulation.pdf , another possible derivation of the same fact can be done through the divergence theorem.

Now as far as the practical part goes - you need to know when to apply the mis weighting. There are a few cases which occur. Initially for the very first ray, you should not use MIS since you shouldn't sample any light directly. Additionally if the last bounce was fully deterministic (ideal reflection/refraction) then you also should not use MIS, since sampling the light is useless in that case. Finally if you hit a light through sampling the bsdf, you should use MIS to add this contribution, and convert the bsdf pdf wrt the area measure, and then combine it with the probability of sampling this point on this light by using the light pdf and use the area formulation estimator of the rendering equation. When you sample a light, you should also use MIS, once again computing the area formulation estimator. Finally if you have point lights, those can be sampled only through NEE, so you should not use MIS.

Note that there may be more mistakes that you have. Also I didn't understand anything of your explanation past "After that I come to BRDF sampling.".

Answer 2

One of the first things that many people get wrong with MIS of direct lighting is that you have to always consider the same light source for both light sampling and BSDF sampling. For example, if you sampled Light $L_i$ during light sampling and the ray spawned from BSDF sampling hit $L_j$, you cannot mix them together since the probabilities cannot be used together for weights.