Multiple Importance Sampling with more than one light

Question

(Note: This has been cross posted from my ompf2 post.)

Recently I've implemented Multiple Importance Sampling for the sampling of surfaces in my ray tracer. I do this by, on each intersection, sampling both a random direction from the BRDF and also sampling a random light, then combining both using the power heuristic.

When I have only a single light, the results are beautiful. Whereas previously I had this:

Now I have this, with the same amount (1024 per pixel) of samples:

So the algorithm seems to work well. The problem is when I add a second light source to the scene, the noise comes back! This is a render with a second light source added:

I am probably using the technique incorrectly. Perhaps I need to sample all light sources at once, instead of randomly picking one like I'm doing now? That doesn't seem like it would scale well as the number of lights go up, however. The implementation was loosely based on pbrt's DirectLightingIntegrator, and I suspect the reason why it doesn't work properly is that I also have multi-bounce indirect lighting. This is my core integration function:

vec3 calc_light_incidence(const Scene& scene, Rng& rng, const Ray& ray, int depth) {
    float ray_weight = 1.0f;

    depth += 1;
    if (depth > 2) {
        const float live_probability = 0.75f;

        if (rng.canonical() > live_probability) {
            return vec3_0;
        } else {
            ray_weight = 1.0f / live_probability;
        }
    }

    vec3 color = vec3_0;

    const Optional<Intersection> surface_hit = find_nearest_intersection(scene, ray);
    if (surface_hit) {
        const vec3 out_dir = -normalized(ray.direction);

        const size_t light_index = (size_t)(rng.canonical() * scene.lights.size());
        const SceneObject* light = &scene.objects[scene.lights[light_index]];
        const float light_weight = (float)scene.lights.size();

        // Sample BRDF
        {
            const vec3 in_dir = cosine_weighted_point_on_hemisphere(rng.canonical(), rng.canonical(), surface_hit->normal);
            const vec3 reflectance = surface_hit->object->material.diffuse_brdf(*surface_hit, in_dir, out_dir);
            const float cos_term = dot(in_dir, surface_hit->normal);

            const float brdf_pdf = cos_term / pi;

            const Optional<Intersection> light_hit = light->shape->intersect(ray_from_surface(*surface_hit, in_dir));
            const float light_pdf = (light_hit ? light->shape->areaPdf(surface_hit->position, in_dir) : 0.f);

            if (brdf_pdf != 0.f && reflectance != vec3_0) {
                const vec3 illuminance = calc_light_incidence(scene, rng, ray_from_surface(*surface_hit, in_dir), depth);
                color += (1.0f / brdf_pdf) * cos_term * reflectance * illuminance * power_heuristic(brdf_pdf, light_pdf);
            }
        }

        // Sample lights
        {
            const ShapeSample light_sample = light->shape->sampleArea(rng, surface_hit->position);
            const vec3 light_vec = light_sample.point - surface_hit->position;
            const vec3 in_dir = normalized(light_vec);
            const float cos_term = vmax(0.f, dot(in_dir, surface_hit->normal));

            const Optional<Intersection> light_hit = find_nearest_intersection(scene, ray_from_surface(*surface_hit, light_vec));
            bool occluded = !light_hit || light_hit->object != light;
            const float light_pdf = light_weight * light_sample.pdf;

            const float brdf_pdf = cos_term / pi;

            const vec3 reflectance = surface_hit->object->material.diffuse_brdf(*surface_hit, in_dir, out_dir);
            if (!occluded && light_pdf != 0.f && reflectance != vec3_0) {
                const vec3 illuminance = calc_light_incidence(scene, rng, ray_from_surface(*surface_hit, in_dir), depth);
                const float differential_area = -dot(light_hit->normal, in_dir) / length_sqr(light_vec);
                color += light_weight * (1.0f / light_pdf) * differential_area * cos_term * reflectance * illuminance * power_heuristic(light_pdf, brdf_pdf);
            }
        }

        if (dot(out_dir, surface_hit->normal) >= 0.f) {
            color += surface_hit->object->material.emmision->getValue(*surface_hit);
        }
    } else {
        //color = lerp(mvec3(1.0f, 0.2f, 0.0f), mvec3(0.35f, 0.9f, 1.0f), 1.0f - std::pow(1.0f - vmax(0.0f, dot(ray.direction, vec3_y)), 2)) * 0.5f;
        color = lerp(mvec3(0.02f, 0.06f, 0.36f), mvec3(0.0f, 0.0f, 0.0f), 1.0f - std::pow(1.0f - vmax(0.0f, dot(ray.direction, vec3_y)), 2)) * 0.5f;
    }

    return color * ray_weight;
}

Answer 1

It seems that the answer to my question is that my approach inherently can't work. After doing thinking about it some more and researching existing renderers, none seem to implement what I'm doing, and I think the noise comes from contributions from lights other than the one randomly picked will not be estimated with the PDF. To do it correctly I would need to loop through all the lights, which gives me essentially path tracing.

(Note however, that at Mitsuba and PBRT implement a Direct Illumination integrator which uses the approach I do, but only works when excluding indirect illumination, which is also true of my approach.)