# Multiple Importance Sampling in Path tracer produces Dark Images

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

• You shouldn't include the cosine foreshortening term in your brdf_pdf. – Hubble Jan 7 at 6:26
• Im using cosine weighted hemispherical sampling so the BRDF PDF is supposed to be cos(theta)/pi – gallickgunner Jan 7 at 9:57
• May be a coincidence but are your diffuse surfaces exactly PI times darker ? – PaulHK Jan 11 at 7:38
• @PaulHK - What would that imply? I don't think it's actually PI times darker, its just mis_weight times darker/brighter. Because as I said as soon as i remove mis_wieght and just divide by light_pdf I get the brighter image. – gallickgunner Jan 11 at 10:33
• A shot in the dark really, I have seen other question on here were the solution was surprisingly that. – PaulHK Jan 14 at 2:07

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

• Thanks for the detailed write-up. So in short (1) I need to compute both light and brdf samples using integral over area not solid angle. And (2) I also need to compute the BRDF PDF w.r.t to the area. I was actually doing (1) but changed it to integral over solid angle and noticed images getting brighter. Well anyways, will report back shortly. – gallickgunner Feb 24 at 18:33
• You can decide on one of the two formulations: either area or solid angle. That means that all of your pdfs (even outside of the weights) should match this formulation, as well the rendering equation formulation that you use and consequently the estimator that you will compute. The good thing is that nothing except for the pdfs in your weights needs to change. Since I think that most of your stuff is already in solid angle formulation (your brdf pdf and the estimators), then just convert the light pdf wrt solid angle (mult by $r^2/\cos\theta_y$). – lightxbulb Feb 24 at 18:47
• Hey thanks for the answer. I think most of the problems have been resolved. Check the update images. The fireflies seem no less though or it's just that there are too many of them to make a difference. I'll try with different scene parameters. I also don't get what you said about the first ray should not use MIS. I think you meant for naive path tracing. I'm using NEE. – gallickgunner Feb 24 at 19:43
• The first ray (emanating from the camera) should not use NEE, since you will hit anything visible either way. The fact that the fireflies do not even diminish is strange (it's not impossible though), so I still have my doubts. – lightxbulb Feb 24 at 20:02
• About fireflies diminishing, I think the real reason for fireflies in my scene is because the reflective sphere shows an image of bright light source. This is very small as compared to the actual size of the light sources. Hence when the rays hitting the walls bounce, only a very few of them end up striking the light source image on the reflective sphere giving an incredibly high color in some samples. This just sounds like I am double dipping lights (1 from NEE, 2nd from reflective sphere) Sounds wrong but never seen such a case where people ignore the reflection of light in GI. – gallickgunner Feb 24 at 21:07

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.

• That answers my point mentioned in the edit portion. The main problem still remains though. – gallickgunner Jan 8 at 16:08