<p>The summation doesn't include the BxDF that was picked for sampling. Look again at this line:</p>

<pre><code>if (bxdfs[i] != bxdf &amp;&amp; bxdfs[i]-&gt;MatchesFlags(type))
</code></pre>

<p>Here <code>bxdf</code> is the one that was sampled earlier, so when it iterates to that one it skips it.</p>

CC BY-SA 4.0
2019-11-13T21:20:43.117Z
2019-11-13T21:20:43.117Z
Nathan Reed
<p>In PBRT, a BSDF is represented as an array of BxDF lobes. Each lobe sets a bitflag representing what kind of lobe it is. (Specular, Diffuse, Reflection, Transmission, etc.) I'm investigating how PBRT implemented sampling a BSDF for an input direction, wi.</p>

<p>Following standard Monte Carlo sampling, I was expecting the code to randomly sample an appropriate lobe, then set the <span class="math-container">$pdf = pdf_{picked}/n$</span> where <span class="math-container">$n$</span> is the number of lobes of the appropriate type. However this isn't what PBRT does. My question is why? To go in futher, this is the PBRT function:</p>

<pre><code>Spectrum BSDF::Sample_f(const Vector3f &amp;woWorld, Vector3f *wiWorld,
                        const Point2f &amp;u, Float *pdf, BxDFType type,
                        BxDFType *sampledType) const {
    // Choose which _BxDF_ to sample
    int matchingComps = NumComponents(type);
    if (matchingComps == 0) {
        *pdf = 0;
        if (sampledType)
            *sampledType = BxDFType(0);
        return Spectrum(0);
    }
    int comp = std::min((int)std::floor(u[0] * matchingComps), matchingComps - 1);

    // Get _BxDF_ pointer for chosen component
    BxDF *bxdf = nullptr;
    int count = comp;
    for (int i = 0; i &lt; nBxDFs; ++i)
        if (bxdfs[i]-&gt;MatchesFlags(type) &amp;&amp; count-- == 0) {
            bxdf = bxdfs[i];
            break;
        }

    // Remap _BxDF_ sample _u_ to <span class="math-container">$[0,1)^2$</span>
    Point2f uRemapped(std::min(u[0] * matchingComps - comp, OneMinusEpsilon),
                      u[1]);

    // Sample chosen _BxDF_
    Vector3f wi, wo = WorldToLocal(woWorld);
    if (wo.z == 0)
        return 0.;
    *pdf = 0;
    if (sampledType) *sampledType = bxdf-&gt;type;
    Spectrum f = bxdf-&gt;Sample_f(wo, &amp;wi, uRemapped, pdf, sampledType);

    if (*pdf == 0) {
        if (sampledType) *sampledType = BxDFType(0);
        return 0;
    }
    *wiWorld = LocalToWorld(wi);

    // Compute overall PDF with all matching _BxDF_s
    if (!(bxdf-&gt;type &amp; BSDF_SPECULAR) &amp;&amp; matchingComps &gt; 1)
        for (int i = 0; i &lt; nBxDFs; ++i)
            if (bxdfs[i] != bxdf &amp;&amp; bxdfs[i]-&gt;MatchesFlags(type))
                *pdf += bxdfs[i]-&gt;Pdf(wo, wi);
    if (matchingComps &gt; 1)
        *pdf /= matchingComps;

    // Compute value of BSDF for sampled direction
    if (!(bxdf-&gt;type &amp; BSDF_SPECULAR)) {
        bool reflect = Dot(*wiWorld, ng) * Dot(woWorld, ng) &gt; 0;
        f = 0.;
        for (int i = 0; i &lt; nBxDFs; ++i)
            if (bxdfs[i]-&gt;MatchesFlags(type) &amp;&amp;
                ((reflect &amp;&amp; (bxdfs[i]-&gt;type &amp; BSDF_REFLECTION)) ||
                 (!reflect &amp;&amp; (bxdfs[i]-&gt;type &amp; BSDF_TRANSMISSION))))
                f += bxdfs[i]-&gt;f(wo, wi);
    }

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

<p>In pseudo-code, it looks to be doing the following:</p>

<pre><code>// Get the number of lobes that match the ones we asked for
numMatchingLobes = NumComponents(type);

// Early out check
if (numMatchingLobes == 0):
    return

// Randomly pick one of the matching lobes
pickedBxDF = Rand()

// Transform wo from world space to local space
// Then sample the chosen BxDF to get f, pdf, and wi
f, pdf, wi = pickedBxDF.sample_f(woWorld.ToLocal())

// Early out check
if (pdf == 0):
    return

// Transform wi to world space
wiWorld = wi.ToWorld()

// Add all of the pdfs of the matching lobes if the picked BxDF wasn't specular
if pickedBxDF.type not Specular &amp;&amp; numMatchingLobes &gt; 1:
    for BxDF in matchingLobes:
        pdf += BxDF.pdf(wi, wo)

// Divide the resulting pdf by the number of matching lobes
pdf /= numMatchingLobes;


// Calculate the combined f if the picked BxDF isn't specular
if pickedBxDF.type not Specular:
    // Reset f to zero
    f = 0

    // Add each BxDF's f, iff it matches reflection / transmission
    reflect = wo and wi on same hemisphere as geoNormal
    for BxDF in matchingLobes:
        // Add the 
        if reflect and BxDF.type == Reflection:
            f += BxDF.f(wo, wi)
        elif !reflect and BxDF.type == Transmission:
            f += BxDF.f(wo, wi)
</code></pre>

<p>The part I'm not understanding is why it adds the pdf's of all the matching lobes on top of the already calculated pdf. This seems like it would double count the pdf of the picked lobe.</p>

<p>For example, if there are two matching lobes, called A and B. A is picked.</p>

<p>We calculate f, pdf, and wi with A. Then with a for loop, we add the pdf for A and B and divide by 2. Aka:</p>

<p><span class="math-container">$$
\frac{A + A + B}{2}
$$</span></p>

<p>This doesn't seem right. Is PBRT code wrong? Should it be resetting the pdf, similar to the f calculation code below it?</p>

<sampling><monte-carlo><bsdf><pbrt>
How to calculate the PDF of a multi-lobe BSDF
RichieSams
# How to calculate the PDF of a multi-lobe BSDF

(Asked by: username_0 on 2019-11-13)

In PBRT, a BSDF is represented as an array of BxDF lobes. Each lobe sets a bitflag representing what kind of lobe it is. (Specular, Diffuse, Reflection, Transmission, etc.) I'm investigating how PBRT implemented sampling a BSDF for an input direction, wi.

Following standard Monte Carlo sampling, I was expecting the code to randomly sample an appropriate lobe, then set the $pdf = pdf\_{picked}/n$ where $n$ is the number of lobes of the appropriate type. However this isn't what PBRT does. My question is why? To go in futher, this is the PBRT function:

```
Spectrum BSDF::Sample_f(const Vector3f &woWorld, Vector3f *wiWorld,
                        const Point2f &u, Float *pdf, BxDFType type,
                        BxDFType *sampledType) const {
    // Choose which _BxDF_ to sample
    int matchingComps = NumComponents(type);
    if (matchingComps == 0) {
        *pdf = 0;
        if (sampledType)
            *sampledType = BxDFType(0);
        return Spectrum(0);
    }
    int comp = std::min((int)std::floor(u[0] * matchingComps), matchingComps - 1);

    // Get _BxDF_ pointer for chosen component
    BxDF *bxdf = nullptr;
    int count = comp;
    for (int i = 0; i < nBxDFs; ++i)
        if (bxdfs[i]->MatchesFlags(type) && count-- == 0) {
            bxdf = bxdfs[i];
            break;
        }

    // Remap _BxDF_ sample _u_ to $[0,1)^2$
    Point2f uRemapped(std::min(u[0] * matchingComps - comp, OneMinusEpsilon),
                      u[1]);

    // Sample chosen _BxDF_
    Vector3f wi, wo = WorldToLocal(woWorld);
    if (wo.z == 0)
        return 0.;
    *pdf = 0;
    if (sampledType) *sampledType = bxdf->type;
    Spectrum f = bxdf->Sample_f(wo, &wi, uRemapped, pdf, sampledType);

    if (*pdf == 0) {
        if (sampledType) *sampledType = BxDFType(0);
        return 0;
    }
    *wiWorld = LocalToWorld(wi);

    // Compute overall PDF with all matching _BxDF_s
    if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1)
        for (int i = 0; i < nBxDFs; ++i)
            if (bxdfs[i] != bxdf && bxdfs[i]->MatchesFlags(type))
                *pdf += bxdfs[i]->Pdf(wo, wi);
    if (matchingComps > 1)
        *pdf /= matchingComps;

    // Compute value of BSDF for sampled direction
    if (!(bxdf->type & BSDF_SPECULAR)) {
        bool reflect = Dot(*wiWorld, ng) * Dot(woWorld, ng) > 0;
        f = 0.;
        for (int i = 0; i < nBxDFs; ++i)
            if (bxdfs[i]->MatchesFlags(type) &&
                ((reflect && (bxdfs[i]->type & BSDF_REFLECTION)) ||
                 (!reflect && (bxdfs[i]->type & BSDF_TRANSMISSION))))
                f += bxdfs[i]->f(wo, wi);
    }

    return f;
}

```

In pseudo-code, it looks to be doing the following:

```
// Get the number of lobes that match the ones we asked for
numMatchingLobes = NumComponents(type);

// Early out check
if (numMatchingLobes == 0):
    return

// Randomly pick one of the matching lobes
pickedBxDF = Rand()

// Transform wo from world space to local space
// Then sample the chosen BxDF to get f, pdf, and wi
f, pdf, wi = pickedBxDF.sample_f(woWorld.ToLocal())

// Early out check
if (pdf == 0):
    return

// Transform wi to world space
wiWorld = wi.ToWorld()

// Add all of the pdfs of the matching lobes if the picked BxDF wasn't specular
if pickedBxDF.type not Specular && numMatchingLobes > 1:
    for BxDF in matchingLobes:
        pdf += BxDF.pdf(wi, wo)

// Divide the resulting pdf by the number of matching lobes
pdf /= numMatchingLobes;

// Calculate the combined f if the picked BxDF isn't specular
if pickedBxDF.type not Specular:
    // Reset f to zero
    f = 0

    // Add each BxDF's f, iff it matches reflection / transmission
    reflect = wo and wi on same hemisphere as geoNormal
    for BxDF in matchingLobes:
        // Add the 
        if reflect and BxDF.type == Reflection:
            f += BxDF.f(wo, wi)
        elif !reflect and BxDF.type == Transmission:
            f += BxDF.f(wo, wi)

```

The part I'm not understanding is why it adds the pdf's of all the matching lobes on top of the already calculated pdf. This seems like it would double count the pdf of the picked lobe.

For example, if there are two matching lobes, called A and B. A is picked.

We calculate f, pdf, and wi with A. Then with a for loop, we add the pdf for A and B and divide by 2. Aka:

$$
\frac{A + A + B}{2}
$$

This doesn't seem right. Is PBRT code wrong? Should it be resetting the pdf, similar to the f calculation code below it?

**Tags:** sampling, monte-carlo, bsdf, pbrt


---

## 1 answers

### Answer by username_1 on 2019-11-13. Score: 1 | Accepted answer

The summation doesn't include the BxDF that was picked for sampling. Look again at this line:

```
if (bxdfs[i] != bxdf && bxdfs[i]->MatchesFlags(type))

```

Here `bxdf` is the one that was sampled earlier, so when it iterates to that one it skips it.


---

