<p>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.</p>

<p>(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.)</p>

CC BY-SA 3.0
2015-08-12T13:37:15.447Z
I'm not sure that's the correct takeaway. I would first double-check all the computations in the BRDF sampling branch. I think ShapeSphere::areaPdf looks suspect. It should depend on the distance squared and both the surface normal and the light normal.
John Calsbeek
CC BY-SA 3.0
2015-08-12T08:17:24.157Z
2015-08-12T08:17:24.157Z
yuriks
<p>(Note: This has been cross posted from my <a href="http://ompf2.com/viewtopic.php?f=15&amp;t=1977" rel="nofollow noreferrer">ompf2 post</a>.)</p>

<p>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.</p>

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

<p><img src="https://i.sstatic.net/kElUN.png" alt="BRDF sampling only"></p>

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

<p><img src="https://i.sstatic.net/H0JbB.png" alt="BRDF+Light sampling"></p>

<p>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:</p>

<p><a href="https://i.sstatic.net/zVdm7.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zVdm7.png" alt="BRDF+Light with 2 lights"></a></p>

<p>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:</p>

<pre class="lang-cpp prettyprint-override"><code>vec3 calc_light_incidence(const Scene&amp; scene, Rng&amp; rng, const Ray&amp; ray, int depth) {
    float ray_weight = 1.0f;

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

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

    vec3 color = vec3_0;

    const Optional&lt;Intersection&gt; 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 = &amp;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-&gt;normal);
            const vec3 reflectance = surface_hit-&gt;object-&gt;material.diffuse_brdf(*surface_hit, in_dir, out_dir);
            const float cos_term = dot(in_dir, surface_hit-&gt;normal);

            const float brdf_pdf = cos_term / pi;

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

            if (brdf_pdf != 0.f &amp;&amp; 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-&gt;shape-&gt;sampleArea(rng, surface_hit-&gt;position);
            const vec3 light_vec = light_sample.point - surface_hit-&gt;position;
            const vec3 in_dir = normalized(light_vec);
            const float cos_term = vmax(0.f, dot(in_dir, surface_hit-&gt;normal));

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

            const float brdf_pdf = cos_term / pi;

            const vec3 reflectance = surface_hit-&gt;object-&gt;material.diffuse_brdf(*surface_hit, in_dir, out_dir);
            if (!occluded &amp;&amp; light_pdf != 0.f &amp;&amp; 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-&gt;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-&gt;normal) &gt;= 0.f) {
            color += surface_hit-&gt;object-&gt;material.emmision-&gt;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;
}
</code></pre>

CC BY-SA 3.0
2015-09-03T21:01:10.693Z
Try randomly picking a light source *per sample* instead of per pixel.
mm201
CC BY-SA 3.0
2015-08-12T03:07:56.277Z
Do you get a similarly smooth image from only the second light? i.e. it's definitely an artifact of multiple lights, not just the geometry of the second light?
John Calsbeek
CC BY-SA 3.0
2015-08-12T03:36:00.983Z
Yes, using any kind of single light gives a smooth image. The noise only comes back when there are multiple lights to choose from. (And I suppose that the 3rd example image is somewhat misleading, but when there are 2 lights, the image looks almost as bad as the first one which uses the naive algorithm if I use only 1024 spp.)
yuriks
CC BY-SA 3.0
2015-08-12T03:39:30.193Z
To me, that suggests a bug. But it's rather difficult to debug code just by reading it. I'd try disabling pieces of code—for example, in this scene the importance weight of the BRDF sample should be extremely small since light sampling is far more likely to hit the light, so you should be able to set brdf_pdf to 0 and only use light sampling and see if the noise goes away.
John Calsbeek
CC BY-SA 3.0
2015-08-12T05:54:57.620Z
@JohnCalsbeek Looks like you gave me a good lead! This is how it looks like if I ignore BRDF sampling: i.imgur.com/4c7bRCo.png I'll try to further track down the issue later when I have more time. It seems like the issue is that the impact of the BRDF is overestimated, causing it to overwhelm the light contribution in some cases.
yuriks
CC BY-SA 3.0
2015-08-12T09:26:27.740Z
You should actually add the significant parts of your code into the question instead of linking them elsewhere.
Christian Rau
<raytracing><importance-sampling>
Multiple Importance Sampling with more than one light
yuriks
# Multiple Importance Sampling with more than one light

(Asked by: username_0 on 2015-08-12)

(Note: This has been cross posted from my [ompf2 post](http://ompf2.com/viewtopic.php?f=15&t=1977).)

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:

![BRDF sampling only](https://i.sstatic.net/kElUN.png)

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

![BRDF+Light sampling](https://i.sstatic.net/H0JbB.png)

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:

[![BRDF+Light with 2 lights](https://i.sstatic.net/zVdm7.png)](https://i.sstatic.net/zVdm7.png)

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:

```cpp
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;
}

```

**Tags:** raytracing, importance-sampling

> **Comments:**
> **username_2**: Do you get a similarly smooth image from only the second light? i.e. it's definitely an artifact of multiple lights, not just the geometry of the second light?  
> **username_0**: Yes, using any kind of single light gives a smooth image. The noise only comes back when there are multiple lights to choose from. (And I suppose that the 3rd example image is somewhat misleading, but when there are 2 lights, the image looks almost as bad as the first one which uses the naive algorithm if I use only 1024 spp.)  
> **username_2**: To me, that suggests a bug. But it's rather difficult to debug code just by reading it. I'd try disabling pieces of code—for example, in this scene the importance weight of the BRDF sample should be extremely small since light sampling is far more likely to hit the light, so you should be able to set brdf_pdf to 0 and only use light sampling and see if the noise goes away.  
> **username_0**: @username_2 Looks like you gave me a good lead! This is how it looks like if I ignore BRDF sampling: i.imgur.com/4c7bRCo.png I'll try to further track down the issue later when I have more time. It seems like the issue is that the impact of the BRDF is overestimated, causing it to overwhelm the light contribution in some cases.  
> **username_3**: You should actually add the significant parts of your code into the question instead of linking them elsewhere.  
> **username_1**: Try randomly picking a light source *per sample* instead of per pixel.  

---

## 1 answers

### Answer by username_0 on 2015-08-12. Score: 3

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

> **Comments:**
> **username_2**: I'm not sure that's the correct takeaway. I would first double-check all the computations in the BRDF sampling branch. I think ShapeSphere::areaPdf looks suspect. It should depend on the distance squared and both the surface normal and the light normal.  

---

