<p>There are multiple areas in path tracing that can be importance sampled. In addition, each of those areas can also use Multiple Importance Sampling, first proposed in Veach and Guibas's <a href="https://sites.fas.harvard.edu/~cs278/papers/veach.pdf" rel="noreferrer">1995 paper</a>. To better explain, let's look at a backwards path tracer:</p>

<pre><code>void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene-&gt;Camera-&gt;CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces &lt; maxBounces; ++bounces) {
        m_scene-&gt;Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene-&gt;BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene-&gt;GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene-&gt;GetLight(ray.GeomID);

        // If we hit a light, add the emission
        if (light != nullptr) {
            color += throughput * light-&gt;Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene-&gt;InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);


        // Get the new ray direction
        // Choose the direction based on the bsdf        
        material-&gt;bsdf-&gt;Sample(interaction, sampler);
        float pdf = material-&gt;bsdf-&gt;Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material-&gt;bsdf-&gt;Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;


        // Russian Roulette
        if (bounces &gt; 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler-&gt;NextFloat() &gt; p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene-&gt;Camera-&gt;FrameBufferData.SplatPixel(x, y, color);
}
</code></pre>

<p><strong>In English:</strong></p>

<ol>
<li>Shoot a ray through the scene</li>
<li>Check if we hit anything. If not we return the skybox color and break.</li>
<li>Check if we hit a light. If so, we add the light emission to our color accumulation</li>
<li>Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF</li>
<li>Evaluate the BRDF and accumulate it. Here we have to divide by the pdf of our chosen direction, in order to follow the Monte Carlo Algorithm.</li>
<li>Create a new ray based on our chosen direction and where we just came from</li>
<li>[Optional] Use Russian Roulette to choose if we should terminate the ray</li>
<li>Goto 1</li>
</ol>

<p>With this code, we only get color if the ray eventually hits a light. In addition, it doesn't support punctual light sources, since they have no area. </p>

<p>To fix this, we sample the lights directly at every bounce. We have to do a few small changes:</p>

<pre><code>void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene-&gt;Camera-&gt;CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces &lt; maxBounces; ++bounces) {
        m_scene-&gt;Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene-&gt;BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene-&gt;GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene-&gt;GetLight(ray.GeomID);

        // If this is the first bounce or if we just had a specular bounce,
        // we need to add the emmisive light
        if ((bounces == 0 || (interaction.SampledLobe &amp; BSDFLobe::Specular) != 0) &amp;&amp; light != nullptr) {
            color += throughput * light-&gt;Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene-&gt;InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);


        // Calculate the direct lighting
        color += throughput * SampleLights(sampler, interaction, material-&gt;bsdf, light);


        // Get the new ray direction
        // Choose the direction based on the bsdf        
        material-&gt;bsdf-&gt;Sample(interaction, sampler);
        float pdf = material-&gt;bsdf-&gt;Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material-&gt;bsdf-&gt;Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;


        // Russian Roulette
        if (bounces &gt; 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler-&gt;NextFloat() &gt; p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene-&gt;Camera-&gt;FrameBufferData.SplatPixel(x, y, color);
}
</code></pre>

<p>First, we add "color += throughput * SampleLights(...)". I'll go into detail about SampleLights() in a bit. But, essentially, it loops through all the lights, and returns their contribution to the color, attenuated by the BSDF.</p>

<p>This is great, but we need to make one more change in order to make it correct; specifically, what happens when we hit a light. In the old code, we added the light's emission to the color accumulation. But now we directly sample the light every bounce, so if we added the light's emission, we would "double dip". Therefore, the correct thing to do is... nothing; we skip accumulating the light's emission.</p>

<p>However, there are two corner cases: </p>

<ol>
<li>The first ray</li>
<li>Perfectly specular bounces (aka mirrors)</li>
</ol>

<p>If the first ray hits the light, you should see the light's emission directly. So if we skip it, all the lights will show up as black, even though the surfaces around them are lit.</p>

<p>When you hit a perfectly specular surfaces you can't directly sample a light, because an input ray has only one output. Well, technically, we <em>could</em> check if the input ray is going to hit a light, but there's no point; the main Path Tracing loop is going to do that anyway. Therefore, if we hit a light just after we hit a specular surface, we need to accumulate the color. If we don't, lights will be black in mirrors.</p>

<p>Now, let's delve into SampleLights(): </p>

<pre><code>float3 SampleLights(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene-&gt;NumLights();

    float3 L(0.0f);
    for (uint i = 0; i &lt; numLights; ++i) {
        Light *light = &amp;m_scene-&gt;Lights[i];

        // Don't let a light contribute light to itself
        if (light == hitLight) {
            continue;
        }

        L = L + EstimateDirect(light, sampler, interaction, bsdf);
    }

    return L;
}
</code></pre>

<p><strong>In English:</strong></p>

<ol>
<li>Loop through all the lights</li>
<li>Skip the light if we hit it

<ul>
<li>Don't double dip</li>
</ul></li>
<li>Accumulate the direct lighting from all the lights</li>
<li>Return the direct lighting</li>
</ol>

<p>Finally, EstimateDirect() is just evaluating $BSDF(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}})$</p>

<p>For punctual light sources, this is simple as:</p>

<pre><code>float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &amp;interaction, BSDF *bsdf) const {
    // Only sample if the BRDF is non-specular 
    if ((bsdf-&gt;SupportedLobes &amp; ~BSDFLobe::Specular) != 0) {
        return float3(0.0f);
    }

    interaction.InputDirection = normalize(light-&gt;Origin - interaction.Position);
    return bsdf-&gt;Eval(interaction) * light-&gt;Li;
}
</code></pre>

<p>However, if we want lights to have area, we first need to sample a point on the light. Therefore, the full definition is:</p>

<pre><code>float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &amp;interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);

    // Only sample if the BRDF is non-specular 
    if ((bsdf-&gt;SupportedLobes &amp; ~BSDFLobe::Specular) != 0) {
        float pdf;
        float3 Li = light-&gt;SampleLi(sampler, m_scene, interaction, &amp;pdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (pdf != 0.0f &amp;&amp; !all(Li)) {
            directLighting += bsdf-&gt;Eval(interaction) * Li / pdf;
        }
    }

    return directLighting;
}
</code></pre>

<p>We can implement light->SampleLi however we want; we can choose the point uniformly, or importance sample. In either case, we divide the radiosity by the pdf of choosing the point. Again, to satisfy the requirements of Monte Carlo.</p>

<p>If the BRDF is highly view dependent, it may be better to choose a point based on the BRDF, instead of a random point on the light. But how do we choose? Sample based on the light, or based on the BRDF?</p>

<p>Why not both? Enter Multiple Importance Sampling. In short, we evaluate $BSDF(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}})$ multiple times, using different sampling techniques, then average them together using weights based on their pdfs. In code this is:</p>

<pre><code>float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &amp;interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);
    float3 f;
    float lightPdf, scatteringPdf;


    // Sample lighting with multiple importance sampling
    // Only sample if the BRDF is non-specular 
    if ((bsdf-&gt;SupportedLobes &amp; ~BSDFLobe::Specular) != 0) {
        float3 Li = light-&gt;SampleLi(sampler, m_scene, interaction, &amp;lightPdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (lightPdf != 0.0f &amp;&amp; !all(Li)) {
            // Calculate the brdf value
            f = bsdf-&gt;Eval(interaction);
            scatteringPdf = bsdf-&gt;Pdf(interaction);

            if (scatteringPdf != 0.0f &amp;&amp; !all(f)) {
                float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf);
                directLighting += f * Li * weight / lightPdf;
            }
        }
    }


    // Sample brdf with multiple importance sampling
    bsdf-&gt;Sample(interaction, sampler);
    f = bsdf-&gt;Eval(interaction);
    scatteringPdf = bsdf-&gt;Pdf(interaction);
    if (scatteringPdf != 0.0f &amp;&amp; !all(f)) {
        lightPdf = light-&gt;PdfLi(m_scene, interaction);
        if (lightPdf == 0.0f) {
            // We didn't hit anything, so ignore the brdf sample
            return directLighting;
        }

        float weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf);
        float3 Li = light-&gt;Le();
        directLighting += f * Li * weight / scatteringPdf;
    }

    return directLighting;
}
</code></pre>

<p><strong>In English:</strong></p>

<ol>
<li>First, we sample the light

<ul>
<li>This updates interaction.InputDirection</li>
<li>Gives us the Li for the light</li>
<li>And the pdf of choosing that point on the light</li>
</ul></li>
<li>Check that the pdf is valid and the radiance is non-zero</li>
<li>Evaluate the BSDF using the sampled InputDirection</li>
<li>Calculate the pdf for the BSDF given the sampled InputDirection

<ul>
<li>Essentially, how likely is this sample, if we were to sample using the BSDF, instead of the light</li>
</ul></li>
<li>Calculate the weight, using the light pdf and the BSDF pdf

<ul>
<li>Veach and Guibas define a couple different ways to calculate the weight. Experimentally, they found the power heuristic with a power of 2 to work the best for most cases. I refer you to the paper for more details. The implementation is below</li>
</ul></li>
<li>Multiply the weight with the direct lighting calculation and divide by the light pdf. (For Monte Carlo) And add to the direct light accumulation.</li>
<li>Then, we sample the BRDF

<ul>
<li>This updates interaction.InputDirection</li>
</ul></li>
<li>Evaluate the BRDF</li>
<li>Get the pdf for choosing this direction based on the BRDF</li>
<li>Calculate the light pdf, given the sampled InputDirection

<ul>
<li>This is the mirror of before. How likely is this direction, if we were to sample the light</li>
</ul></li>
<li>If lightPdf == 0.0f, then the ray missed the light, so just return the direct lighting from the light sample.</li>
<li>Otherwise, calculate the weight, and add the BSDF direct lighting to the accumulation</li>
<li>Finally, return the accumulated direct lighting</li>
</ol>

<p>.</p>

<pre><code>inline float PowerHeuristic(uint numf, float fPdf, uint numg, float gPdf) {
    float f = numf * fPdf;
    float g = numg * gPdf;

    return (f * f) / (f * f + g * g);
}
</code></pre>

<p>There are a number of optimizations / improvements you can do in these functions, but I've pared them down to try to make them easier to comprehend. If you would like, I can share some of these improvements.</p>

<h2>Only Sampling One Light</h2>

<p>In SampleLights() we loop through all the lights, and get their contribution. For a small number of lights, this is fine, but for hundreds or thousands of lights, this gets expensive. Fortunately, we can exploit the fact that Monte Carlo Integration is a giant average. Example:</p>

<p>Let's define </p>

<p>$$h(x) = f(x) + g(x)$$</p>

<p>Currently, we're estimating $h(x)$ by:</p>

<p>$$h(x) = \frac{1}{N} \sum_{i=1}^N f(x_i) + g(x_i)$$</p>

<p>But, calculating both $f(x)$ and $g(x)$ is expensive, so instead we do:</p>

<p>$$h(x) = \frac{1}{N} \sum_{i=1}^N \frac{r(\zeta, x)}{pdf}$$</p>

<p>Where $\zeta$ is a uniform random variable, and $r(\zeta, x)$ is defined as:</p>

<p>$$r(\zeta, x) = 
\begin{cases}
  f(x), &amp; 0.0 \leq \zeta \lt 0.5 \\
  g(x), &amp; 0.5 \leq \zeta \lt 1.0
\end{cases}$$</p>

<p>In this case $pdf = \frac{1}{2}$ because the pdf must integrate to 1, and there are 2 functions to choose from.</p>

<p><strong>In English:</strong></p>

<ol>
<li>Randomly choose either $f(x)$ or $g(x)$ to evaluate. </li>
<li>Divide the result by $\frac{1}{2}$ (since there are two items)</li>
<li>Average</li>
</ol>

<p>As N gets large, the estimate will converge to the correct solution. </p>

<p>We can apply this same principle to light sampling. Instead of sampling every light, we randomly pick one, and multiply the result by the number of lights (This is the same as dividing by the fractional pdf):</p>

<pre><code>float3 SampleOneLight(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene-&gt;NumLights();

    // Return black if there are no lights
    // And don't let a light contribute light to itself
    // Aka, if we hit a light
    // This is the special case where there is only 1 light
    if (numLights == 0 || numLights == 1 &amp;&amp; hitLight != nullptr) {
        return float3(0.0f);
    }

    // Don't let a light contribute light to itself
    // Choose another one
    Light *light;
    do {
        light = m_scene-&gt;RandomOneLight(sampler);
    } while (light == hitLight);

    return numLights * EstimateDirect(light, sampler, interaction, bsdf);
}
</code></pre>

<p>In this code, all the lights have an equal chance of being picked. However, we can importance sample, if we like. For example, we can give larger lights a higher chance of being picked, or lights closer to the hit surface. You just have to divide the result by the pdf, which would no longer be $\frac{1}{\text{numLights}}$.</p>

<h2>Multiple Importance Sampling the "New Ray" Direction</h2>

<p>The current code only importance samples the "New Ray" direction based on the BSDF. What if we want to also importance sample based on the location of lights?</p>

<p>Taking from what we learned above, one method would be to shoot <strong>two</strong> "new" rays and weight each based on their pdfs. However, this is both computationally expensive, and hard to implement without recursion.</p>

<p>To overcome this, we can apply the same principles we learned by sampling only one light. That is, randomly choose one to sample, and divide by the pdf of choosing it.</p>

<pre><code>// Get the new ray direction

// Randomly (uniform) choose whether to sample based on the BSDF or the Lights
float p = sampler-&gt;NextFloat();

Light *light = m_scene-&gt;RandomLight();

if (p &lt; 0.5f) {
    // Choose the direction based on the bsdf 
    material-&gt;bsdf-&gt;Sample(interaction, sampler);
    float bsdfPdf = material-&gt;bsdf-&gt;Pdf(interaction);

    float lightPdf = light-&gt;PdfLi(m_scene, interaction);
    float weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material-&gt;bsdf-&gt;Eval(interaction) / bsdfPdf;

} else {
    // Choose the direction based on a light
    float lightPdf;
    light-&gt;SampleLi(sampler, m_scene, interaction, &amp;lightPdf);

    float bsdfPdf = material-&gt;bsdf-&gt;Pdf(interaction);
    float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material-&gt;bsdf-&gt;Eval(interaction) / lightPdf;
}
</code></pre>

<p>That all said, do we really <strong><em>want</em></strong> to importance sample the "New Ray" direction based on the light? For <em>direct</em> lighting, the radiosity is affected by both the BSDF of the surface, and the direction of the light. But for <em>indirect</em> lighting, the radiosity is almost exclusively defined by the BSDF of the surface hit before. So, adding light importance sampling doesn't give us anything. </p>

<p>Therefore, it is common to only importance sample the "New Direction" with the BSDF, but apply Multiple Importance Sampling to the direct lighting.</p>

CC BY-SA 3.0
2017-05-25T15:37:02.323Z
Thank you for the clarifying answer! I understand that if we were to use a path tracer without explicit light sampling, we would never hit a point light source. So, we can basically add its contribution. On the other hand, if we sample an area light source, we have to make sure we should not hit it again with the indirect lighting in order to avoid double dip
Mustafa Işık
CC BY-SA 3.0
2017-05-25T15:39:24.237Z
Exactly! Is there any part that you need clarification on? Or there isn't enough detail?
RichieSams
CC BY-SA 3.0
2017-05-25T15:43:01.690Z
Also, is multiple importance sampling used only for direct lighting calculation? Maybe I missed but I didn't see another example of it. If I shoot just one ray per bounce in my path tracer, it seems that I cannot do it for the indirect lighting calculation.
Mustafa Işık
CC BY-SA 3.0
2017-05-25T18:53:35.473Z
Multiple Importance Sampling can be applied anywhere you use importance sampling. The power of multiple importance sampling is that we can combine the benefits of multiple sampling techniques. For example, in some cases, light importance sampling will be better than BSDF sampling. In other cases, vice versa. MIS will combine the best of both worlds. However, if BSDF sampling will be better 100% of the time, there is no reason to add the complexity of MIS. I added some sections to the answer to expand upon this point
RichieSams
CC BY-SA 3.0
2017-05-25T19:50:22.677Z
It seems we separated incoming radiance sources into two parts as direct and indirect. We sample lights explicitly for the direct part and while sampling this part, it is reasonable to importance sample the lights as well as BSDFs. For the indirect part, however, we have no idea about which direction may potentially give us higher radiance values since it is the problem itself that we want to solve. However, we can say which direction can contribute more according to the cosine term and BSDF. This is what I understand. Correct me if I'm wrong and thank you for your awesome answer.
Mustafa Işık
CC BY-SA 3.0
2017-05-25T19:52:58.853Z
That's correct! Very good summary
RichieSams
CC BY-SA 3.0
2017-05-27T23:24:40.693Z
Hi Richie, I enjoyed reading your answer here. While ago I posted a question on MIS but didn't get answer. I wonder if I could have your thoughts as well, thanks. https://computergraphics.stackexchange.com/questions/4868/efficiently-sampling-specific-surfaces-using-mis-in-path-tracing
ali
CC BY-SA 4.0
2018-12-23T00:32:44.167Z
@RichieSams Am I missing something or is your code sample not casting any shadow rays?
FogleBird
CC BY-SA 4.0
2018-12-23T01:12:15.977Z
The first code block doesn't have shadow rays. However, the following code blocks add light sampling, which are a form of shadow rays
RichieSams
CC BY-SA 4.0
2021-06-22T12:23:31.737Z
Hi, how to treat transparent object when testing with direct light?
iaomw
CC BY-SA 4.0
2021-06-28T20:12:28.350Z
The "correct" (but will need lots more SPP) approach is to treat all "transparent" objects as solid, for purposes of visibility for NEE. There is active research in how to get around that. For example: https://www.ics.uci.edu/~yug10/projects/translucent/report.html
RichieSams
CC BY-SA 4.0
2024-06-17T18:58:39.360Z
Hi again @RichieSams, 3 years later. You mentioned "a number of optimizations / improvements", what are they? I'm curious. Also, the link you shared on handling transparent objects is dead. Do you have a "replacement" for it? I'm struggling to find resources on how to handle transmissive objects when evaluating the visibility of my shadow rays. Treating transmissive objects as opaque (shadowing the light) completely disables caustics in the scene and so my renders really lack something when glass materials are involved : (
Tom Clabault
CC BY-SA 3.0
2017-05-25T04:09:07.557Z
2017-05-25T18:56:47.003Z
2017-05-25T18:56:47.003Z
RichieSams
<p>I understood the logic behind the importance sampling for the BRDF part. However, when it comes to sampling light sources explicitly, all becomes confusing. For example, if I have one point light source in my scene and if I directly sample it at each frame constantly, should I count it as one more sample for the monte carlo integration? That is, I take one sample from cosine-weighted distribution and other from the point light. Is it two samples in total or just one? Also, should I divide the radiance coming from the direct sample to any term?</p>

<raytracing><pathtracing><global-illumination>
Progressive Path Tracing with Explicit Light Sampling
Mustafa Işık
# Progressive Path Tracing with Explicit Light Sampling

(Asked by: username_0 on 2017-05-24)

I understood the logic behind the importance sampling for the BRDF part. However, when it comes to sampling light sources explicitly, all becomes confusing. For example, if I have one point light source in my scene and if I directly sample it at each frame constantly, should I count it as one more sample for the monte carlo integration? That is, I take one sample from cosine-weighted distribution and other from the point light. Is it two samples in total or just one? Also, should I divide the radiance coming from the direct sample to any term?

**Tags:** raytracing, pathtracing, global-illumination


---

## 1 answers

### Answer by username_1 on 2017-05-25. Score: 30 | Accepted answer

There are multiple areas in path tracing that can be importance sampled. In addition, each of those areas can also use Multiple Importance Sampling, first proposed in Veach and Guibas's [1995 paper](https://sites.fas.harvard.edu/~cs278/papers/veach.pdf). To better explain, let's look at a backwards path tracer:

```
void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera->CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces < maxBounces; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.GeomID);

        // If we hit a light, add the emission
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Get the new ray direction
        // Choose the direction based on the bsdf        
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene->Camera->FrameBufferData.SplatPixel(x, y, color);
}

```

**In English:**

1. Shoot a ray through the scene
2. Check if we hit anything. If not we return the skybox color and break.
3. Check if we hit a light. If so, we add the light emission to our color accumulation
4. Choose a new direction for the next ray. We can do this uniformly, or importance sample based on the BRDF
5. Evaluate the BRDF and accumulate it. Here we have to divide by the pdf of our chosen direction, in order to follow the Monte Carlo Algorithm.
6. Create a new ray based on our chosen direction and where we just came from
7. [Optional] Use Russian Roulette to choose if we should terminate the ray
8. Goto 1

With this code, we only get color if the ray eventually hits a light. In addition, it doesn't support punctual light sources, since they have no area.

To fix this, we sample the lights directly at every bounce. We have to do a few small changes:

```
void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene->Camera->CalculateRayFromPixel(x, y, sampler);

    float3 color(0.0f);
    float3 throughput(1.0f);
    SurfaceInteraction interaction;

    // Bounce the ray around the scene
    const uint maxBounces = 15;
    for (uint bounces = 0; bounces < maxBounces; ++bounces) {
        m_scene->Intersect(ray);

        // The ray missed. Return the background color
        if (ray.GeomID == INVALID_GEOMETRY_ID) {
            color += throughput * m_scene->BackgroundColor;
            break;
        }

        // Fetch the material
        Material *material = m_scene->GetMaterial(ray.GeomID);
        // The object might be emissive. If so, it will have a corresponding light
        // Otherwise, GetLight will return nullptr
        Light *light = m_scene->GetLight(ray.GeomID);

        // If this is the first bounce or if we just had a specular bounce,
        // we need to add the emmisive light
        if ((bounces == 0 || (interaction.SampledLobe & BSDFLobe::Specular) != 0) && light != nullptr) {
            color += throughput * light->Le();
        }

        interaction.Position = ray.Origin + ray.Direction * ray.TFar;
        interaction.Normal = normalize(m_scene->InterpolateNormal(ray.GeomID, ray.PrimID, ray.U, ray.V));
        interaction.OutputDirection = normalize(-ray.Direction);

        // Calculate the direct lighting
        color += throughput * SampleLights(sampler, interaction, material->bsdf, light);

        // Get the new ray direction
        // Choose the direction based on the bsdf        
        material->bsdf->Sample(interaction, sampler);
        float pdf = material->bsdf->Pdf(interaction);

        // Accumulate the weight
        throughput = throughput * material->bsdf->Eval(interaction) / pdf;

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.Origin = interaction.Position;

        // Reset the other ray properties
        ray.Direction = interaction.InputDirection;
        ray.TNear = 0.001f;
        ray.TFar = infinity;

        // Russian Roulette
        if (bounces > 3) {
            float p = std::max(throughput.x, std::max(throughput.y, throughput.z));
            if (sampler->NextFloat() > p) {
                break;
            }

            throughput *= 1 / p;
        }
    }

    m_scene->Camera->FrameBufferData.SplatPixel(x, y, color);
}

```

First, we add "color += throughput \* SampleLights(...)". I'll go into detail about SampleLights() in a bit. But, essentially, it loops through all the lights, and returns their contribution to the color, attenuated by the BSDF.

This is great, but we need to make one more change in order to make it correct; specifically, what happens when we hit a light. In the old code, we added the light's emission to the color accumulation. But now we directly sample the light every bounce, so if we added the light's emission, we would "double dip". Therefore, the correct thing to do is... nothing; we skip accumulating the light's emission.

However, there are two corner cases:

1. The first ray
2. Perfectly specular bounces (aka mirrors)

If the first ray hits the light, you should see the light's emission directly. So if we skip it, all the lights will show up as black, even though the surfaces around them are lit.

When you hit a perfectly specular surfaces you can't directly sample a light, because an input ray has only one output. Well, technically, we *could* check if the input ray is going to hit a light, but there's no point; the main Path Tracing loop is going to do that anyway. Therefore, if we hit a light just after we hit a specular surface, we need to accumulate the color. If we don't, lights will be black in mirrors.

Now, let's delve into SampleLights():

```
float3 SampleLights(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene->NumLights();

    float3 L(0.0f);
    for (uint i = 0; i < numLights; ++i) {
        Light *light = &m_scene->Lights[i];

        // Don't let a light contribute light to itself
        if (light == hitLight) {
            continue;
        }

        L = L + EstimateDirect(light, sampler, interaction, bsdf);
    }

    return L;
}

```

**In English:**

1. Loop through all the lights
2. Skip the light if we hit it
   * Don't double dip
3. Accumulate the direct lighting from all the lights
4. Return the direct lighting

Finally, EstimateDirect() is just evaluating $BSDF(p, \omega\_{\text{i}}, \omega\_{\text{o}}) L\_{\text{i}}(p, \omega\_{\text{i}})$

For punctual light sources, this is simple as:

```
float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    // Only sample if the BRDF is non-specular 
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        return float3(0.0f);
    }

    interaction.InputDirection = normalize(light->Origin - interaction.Position);
    return bsdf->Eval(interaction) * light->Li;
}

```

However, if we want lights to have area, we first need to sample a point on the light. Therefore, the full definition is:

```
float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);

    // Only sample if the BRDF is non-specular 
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        float pdf;
        float3 Li = light->SampleLi(sampler, m_scene, interaction, &pdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (pdf != 0.0f && !all(Li)) {
            directLighting += bsdf->Eval(interaction) * Li / pdf;
        }
    }

    return directLighting;
}

```

We can implement light->SampleLi however we want; we can choose the point uniformly, or importance sample. In either case, we divide the radiosity by the pdf of choosing the point. Again, to satisfy the requirements of Monte Carlo.

If the BRDF is highly view dependent, it may be better to choose a point based on the BRDF, instead of a random point on the light. But how do we choose? Sample based on the light, or based on the BRDF?

Why not both? Enter Multiple Importance Sampling. In short, we evaluate $BSDF(p, \omega\_{\text{i}}, \omega\_{\text{o}}) L\_{\text{i}}(p, \omega\_{\text{i}})$ multiple times, using different sampling techniques, then average them together using weights based on their pdfs. In code this is:

```
float3 EstimateDirect(Light *light, UniformSampler *sampler, SurfaceInteraction &interaction, BSDF *bsdf) const {
    float3 directLighting = float3(0.0f);
    float3 f;
    float lightPdf, scatteringPdf;

    // Sample lighting with multiple importance sampling
    // Only sample if the BRDF is non-specular 
    if ((bsdf->SupportedLobes & ~BSDFLobe::Specular) != 0) {
        float3 Li = light->SampleLi(sampler, m_scene, interaction, &lightPdf);

        // Make sure the pdf isn't zero and the radiance isn't black
        if (lightPdf != 0.0f && !all(Li)) {
            // Calculate the brdf value
            f = bsdf->Eval(interaction);
            scatteringPdf = bsdf->Pdf(interaction);

            if (scatteringPdf != 0.0f && !all(f)) {
                float weight = PowerHeuristic(1, lightPdf, 1, scatteringPdf);
                directLighting += f * Li * weight / lightPdf;
            }
        }
    }

    // Sample brdf with multiple importance sampling
    bsdf->Sample(interaction, sampler);
    f = bsdf->Eval(interaction);
    scatteringPdf = bsdf->Pdf(interaction);
    if (scatteringPdf != 0.0f && !all(f)) {
        lightPdf = light->PdfLi(m_scene, interaction);
        if (lightPdf == 0.0f) {
            // We didn't hit anything, so ignore the brdf sample
            return directLighting;
        }

        float weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf);
        float3 Li = light->Le();
        directLighting += f * Li * weight / scatteringPdf;
    }

    return directLighting;
}

```

**In English:**

1. First, we sample the light
   * This updates interaction.InputDirection
   * Gives us the Li for the light
   * And the pdf of choosing that point on the light
2. Check that the pdf is valid and the radiance is non-zero
3. Evaluate the BSDF using the sampled InputDirection
4. Calculate the pdf for the BSDF given the sampled InputDirection
   * Essentially, how likely is this sample, if we were to sample using the BSDF, instead of the light
5. Calculate the weight, using the light pdf and the BSDF pdf
   * Veach and Guibas define a couple different ways to calculate the weight. Experimentally, they found the power heuristic with a power of 2 to work the best for most cases. I refer you to the paper for more details. The implementation is below
6. Multiply the weight with the direct lighting calculation and divide by the light pdf. (For Monte Carlo) And add to the direct light accumulation.
7. Then, we sample the BRDF
   * This updates interaction.InputDirection
8. Evaluate the BRDF
9. Get the pdf for choosing this direction based on the BRDF
10. Calculate the light pdf, given the sampled InputDirection
    * This is the mirror of before. How likely is this direction, if we were to sample the light
11. If lightPdf == 0.0f, then the ray missed the light, so just return the direct lighting from the light sample.
12. Otherwise, calculate the weight, and add the BSDF direct lighting to the accumulation
13. Finally, return the accumulated direct lighting

.

```
inline float PowerHeuristic(uint numf, float fPdf, uint numg, float gPdf) {
    float f = numf * fPdf;
    float g = numg * gPdf;

    return (f * f) / (f * f + g * g);
}

```

There are a number of optimizations / improvements you can do in these functions, but I've pared them down to try to make them easier to comprehend. If you would like, I can share some of these improvements.

## Only Sampling One Light

In SampleLights() we loop through all the lights, and get their contribution. For a small number of lights, this is fine, but for hundreds or thousands of lights, this gets expensive. Fortunately, we can exploit the fact that Monte Carlo Integration is a giant average. Example:

Let's define

$$h(x) = f(x) + g(x)$$

Currently, we're estimating $h(x)$ by:

$$h(x) = \frac{1}{N} \sum\_{i=1}^N f(x\_i) + g(x\_i)$$

But, calculating both $f(x)$ and $g(x)$ is expensive, so instead we do:

$$h(x) = \frac{1}{N} \sum\_{i=1}^N \frac{r(\zeta, x)}{pdf}$$

Where $\zeta$ is a uniform random variable, and $r(\zeta, x)$ is defined as:

$$r(\zeta, x) =
\begin{cases}
f(x), & 0.0 \leq \zeta \lt 0.5 \\
g(x), & 0.5 \leq \zeta \lt 1.0
\end{cases}$$

In this case $pdf = \frac{1}{2}$ because the pdf must integrate to 1, and there are 2 functions to choose from.

**In English:**

1. Randomly choose either $f(x)$ or $g(x)$ to evaluate.
2. Divide the result by $\frac{1}{2}$ (since there are two items)
3. Average

As N gets large, the estimate will converge to the correct solution.

We can apply this same principle to light sampling. Instead of sampling every light, we randomly pick one, and multiply the result by the number of lights (This is the same as dividing by the fractional pdf):

```
float3 SampleOneLight(UniformSampler *sampler, SurfaceInteraction interaction, BSDF *bsdf, Light *hitLight) const {
    std::size_t numLights = m_scene->NumLights();

    // Return black if there are no lights
    // And don't let a light contribute light to itself
    // Aka, if we hit a light
    // This is the special case where there is only 1 light
    if (numLights == 0 || numLights == 1 && hitLight != nullptr) {
        return float3(0.0f);
    }

    // Don't let a light contribute light to itself
    // Choose another one
    Light *light;
    do {
        light = m_scene->RandomOneLight(sampler);
    } while (light == hitLight);

    return numLights * EstimateDirect(light, sampler, interaction, bsdf);
}

```

In this code, all the lights have an equal chance of being picked. However, we can importance sample, if we like. For example, we can give larger lights a higher chance of being picked, or lights closer to the hit surface. You just have to divide the result by the pdf, which would no longer be $\frac{1}{\text{numLights}}$.

## Multiple Importance Sampling the "New Ray" Direction

The current code only importance samples the "New Ray" direction based on the BSDF. What if we want to also importance sample based on the location of lights?

Taking from what we learned above, one method would be to shoot **two** "new" rays and weight each based on their pdfs. However, this is both computationally expensive, and hard to implement without recursion.

To overcome this, we can apply the same principles we learned by sampling only one light. That is, randomly choose one to sample, and divide by the pdf of choosing it.

```
// Get the new ray direction

// Randomly (uniform) choose whether to sample based on the BSDF or the Lights
float p = sampler->NextFloat();

Light *light = m_scene->RandomLight();

if (p < 0.5f) {
    // Choose the direction based on the bsdf 
    material->bsdf->Sample(interaction, sampler);
    float bsdfPdf = material->bsdf->Pdf(interaction);

    float lightPdf = light->PdfLi(m_scene, interaction);
    float weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material->bsdf->Eval(interaction) / bsdfPdf;

} else {
    // Choose the direction based on a light
    float lightPdf;
    light->SampleLi(sampler, m_scene, interaction, &lightPdf);

    float bsdfPdf = material->bsdf->Pdf(interaction);
    float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);

    // Accumulate the throughput
    throughput = throughput * weight * material->bsdf->Eval(interaction) / lightPdf;
}

```

That all said, do we really ***want*** to importance sample the "New Ray" direction based on the light? For *direct* lighting, the radiosity is affected by both the BSDF of the surface, and the direction of the light. But for *indirect* lighting, the radiosity is almost exclusively defined by the BSDF of the surface hit before. So, adding light importance sampling doesn't give us anything.

Therefore, it is common to only importance sample the "New Direction" with the BSDF, but apply Multiple Importance Sampling to the direct lighting.

> **Comments:**
> **username_0**: Thank you for the clarifying answer! I understand that if we were to use a path tracer without explicit light sampling, we would never hit a point light source. So, we can basically add its contribution. On the other hand, if we sample an area light source, we have to make sure we should not hit it again with the indirect lighting in order to avoid double dip  
> **username_1**: Exactly! Is there any part that you need clarification on? Or there isn't enough detail?  
> **username_0**: Also, is multiple importance sampling used only for direct lighting calculation? Maybe I missed but I didn't see another example of it. If I shoot just one ray per bounce in my path tracer, it seems that I cannot do it for the indirect lighting calculation.  
> **username_1**: Multiple Importance Sampling can be applied anywhere you use importance sampling. The power of multiple importance sampling is that we can combine the benefits of multiple sampling techniques. For example, in some cases, light importance sampling will be better than BSDF sampling. In other cases, vice versa. MIS will combine the best of both worlds. However, if BSDF sampling will be better 100% of the time, there is no reason to add the complexity of MIS. I added some sections to the answer to expand upon this point  
> **username_0**: It seems we separated incoming radiance sources into two parts as direct and indirect. We sample lights explicitly for the direct part and while sampling this part, it is reasonable to importance sample the lights as well as BSDFs. For the indirect part, however, we have no idea about which direction may potentially give us higher radiance values since it is the problem itself that we want to solve. However, we can say which direction can contribute more according to the cosine term and BSDF. This is what I understand. Correct me if I'm wrong and thank you for your awesome answer.  
> **username_1**: That's correct! Very good summary  
> **username_2**: Hi Richie, I enjoyed reading your answer here. While ago I posted a question on MIS but didn't get answer. I wonder if I could have your thoughts as well, thanks. https://computergraphics.stackexchange.com/questions/4868/efficiently-sampling-specific-surfaces-using-mis-in-path-tracing  
> **username_3**: @username_1 Am I missing something or is your code sample not casting any shadow rays?  
> **username_1**: The first code block doesn't have shadow rays. However, the following code blocks add light sampling, which are a form of shadow rays  
> **username_4**: Hi, how to treat transparent object when testing with direct light?  
> **username_1**: The "correct" (but will need lots more SPP) approach is to treat all "transparent" objects as solid, for purposes of visibility for NEE. There is active research in how to get around that. For example: https://www.ics.uci.edu/~yug10/projects/translucent/report.html  
> **username_5**: Hi again @username_1, 3 years later. You mentioned "a number of optimizations / improvements", what are they? I'm curious. Also, the link you shared on handling transparent objects is dead. Do you have a "replacement" for it? I'm struggling to find resources on how to handle transmissive objects when evaluating the visibility of my shadow rays. Treating transmissive objects as opaque (shadowing the light) completely disables caustics in the scene and so my renders really lack something when glass materials are involved : (  

---

