I have started on getting a multiple importance sampling system going in my ray-tracer project for the light sampling. Rather than having a special "light" object type, I've decided to simply make all lights emissive meshes. There is a pre-processing step before BVH construction during which an internal "lightlist" is built - the renderer checks every triangle in the scene, and if any assigned material has the "emits light" flag enabled a pointer to the triangle is added to the lightlist.
During shader evaluation, lights are looped through from the lightlist and sampled directly. I had 2 possible ways to implement this in an unbiased and energy conserving way:
1. Two separate loops - a "bxdf sample" loop and a "light sample" loop. During bxdf sampling all rays that strike the light will return black (the bxdf sampling loop only returns indirect light contribution). During light sampling all rays will be exclusively directed at the lights and only return a value if striking a light. The math works out to be incredibly simple - the light loop and bxdf loop can happen independently and each one can have an arbitrary amount of samples. It looks like this:
Total_Bxdf_Irradiance += Bxdf_Irradiance / number_of_bxdf_samples
for each light:
Total_Light_Irradiance += (Light_Irradiance / distance_to_light_squared) / number_of_light_samples
Total_Irradiance += Total_Bxdf_Irradiance + Total_Light_Irradiance
2. One loop which randomly chooses either a bxdf sample or a light sample, shoots that out, then downweighs it probabilistically based on the probability density function (pdf) of that sampling strategy. This is called multiple importance sampling and the math is more tricky than the previous method. Assuming "ratio" is a zero-to-one ratio value between how many bxdf samples vs light samples are sent out, it looks like this:
if random_number(0..1) < ratio:
Total_Irradiance += Bxdf_Irradiance / ( PDF_Bxdf * ratio + PDF_Light * (1 - ratio) )
Total_Irradiance += Light_Irradiance / ( PDF_Bxdf * ratio + PDF_Light * (1 - ratio) )
Total_Irradiance /= number_of_samples
I have decided to try and implement the second method as it is more challenging, and avoids fireflies where light sources get near objects. For the moment I have a hacked solution and have not nailed the PDF calculation of the light sampling yet. I also wish to add a "probabilistic light selection" rather than always sampling every single light in the scene, and I am not sure how this will effect the PDF of the light sampling. Hopefully it is simple and the PDF remains unchanged, but I need to find and read more papers about this or flip through PBRT again and see if they cover this.