<p>If you have a 1D function <span class="math-container">$f(x)$</span> and you want to integrate this function from say 0 to 1, one way to perform this integration is by taking N random samples in range [0, 1], evaluate <span class="math-container">$f(x)$</span> for each sample and calculate the average of the samples. However, this &quot;naive&quot; Monte Carlo integration is said to &quot;converge slowly&quot;, i.e. you need a large number of samples to get close to the ground truth, particularly if the function has high frequencies.</p>
<p>With importance sampling, instead of taking N random samples in [0, 1] range, you take more samples in the &quot;important&quot; regions of <span class="math-container">$f(x)$</span> that contribute most to the final result. However, because you bias sampling towards the important regions of the function, these samples must be weighted less to counter the bias, which is where the PDF (probability density function) comes along. PDF tells the probability of a sample at a given position and is used to calculate the weighted average of the samples by dividing each sample with the PDF value at each sample position.</p>
<p>With Cook-Torrance importance sampling the common practice is to distribute samples based on the normal distribution function NDF. If NDF is already normalized, it can serve directly as PDF, which is convenient since it cancels the term out from the BRDF evaluation. Only thing you need to do then is to distribute sample positions based on PDF and evaluate BRDF without the NDF term, i.e.
<span class="math-container">$$f=\frac{FG}{\pi(n\cdot\omega_i)(n\cdot\omega_o)}$$</span>
And calculate the average of the sample results multiplied by the solid angle of the domain you integrate over (e.g. <span class="math-container">$2\pi$</span> for hemisphere).</p>
<p>For NDF you need to calculate <a href="https://en.wikipedia.org/wiki/Cumulative_distribution_function" rel="nofollow noreferrer">Cumulative Distribution Function</a> of the PDF to convert uniformly distributed sample position to PDF weighted sample position. For isotropic NDF this simplifies to 1D function because of the symmetry of the function. For more details about the CDF derivation you can check <a href="http://http.developer.nvidia.com/GPUGems3/gpugems3_ch20.html" rel="nofollow noreferrer">this old GPU Gems article</a>.</p>

CC BY-SA 4.0
2017-04-12T17:08:08.050Z
2022-10-28T10:20:25.967Z
2022-10-28T10:20:25.967Z
JarkkoL
<h1>Short answer:</h1>
<p><strong>Importance sampling</strong> is a method to reduce variance in Monte Carlo Integration by choosing an estimator close to the shape of the actual function.</p>
<p><strong>PDF</strong> is an abbreviation for <a href="https://en.wikipedia.org/wiki/Probability_density_function" rel="noreferrer">Probability Density Function</a>. A <span class="math-container">$pdf(x)$</span> gives the probability of a random sample generated being <span class="math-container">$x$</span>.</p>
<h1>Long Answer:</h1>
<p>To start, let's review what Monte Carlo Integration is, and what it looks like mathematically.</p>
<p>Monte Carlo Integration is an technique to estimate the value of an integral. It's typically used when there isn't a closed form solution to the integral. It looks like this:</p>
<p><span class="math-container">$$\int f(x) \, dx \approx \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{pdf(x_{i})}$$</span></p>
<p>In english, this says that you can approximate an integral by averaging successive random samples from the function. As <span class="math-container">$N$</span> gets large, the approximation gets closer and closer to the solution. <span class="math-container">$pdf(x_{i})$</span> represents the probability density function of each random sample.</p>
<p>Let's do an example: Calculate the value of the integral <span class="math-container">$I$</span>.</p>
<p><span class="math-container">$$I = \int_{0}^{2\pi} e^{-x} \sin(x) dx$$</span></p>
<p>Let's use Monte Carlo Integration:</p>
<p><span class="math-container">$$I \approx \frac{1}{N} \sum_{i=1}^N \frac{e^{-x} \sin(x_i)}{pdf(x_{i})}$$</span></p>
<p>A simple python program to calculate this is:</p>
<pre><code>import random
import math

N = 200000
TwoPi = 2.0 * math.pi

sum = 0.0

for i in range(N):
    x = random.uniform(0, TwoPi)
    
    fx = math.exp(-x) * math.sin(x)
    pdf = 1 / (TwoPi - 0.0)
    
    sum += fx / pdf

I = (1 / N) * sum
print(I)
</code></pre>
<p>If we run the program we get <span class="math-container">$I = 0.4986941$</span></p>
<p>Using separation by parts, we can get the exact solution:</p>
<p><span class="math-container">$$I = \frac{1}{2} (1 − e−2π) = 0.4990663$$</span></p>
<p>You'll notice that the Monte Carlo Solution is not quite correct. This is because it is an estimate. That said, as <span class="math-container">$N$</span> goes to infinity, the estimate should get closer and closer to the correct answer. Already at <span class="math-container">$N = 2000$</span> some runs are almost identical to the correct answer.</p>
<p>A note about the PDF: In this simple example, we always take a uniform random sample. A uniform random sample means every sample has the exact same probability of being chosen. We sample in the range <span class="math-container">$[0, 2\pi]$</span> so, <span class="math-container">$pdf(x) = 1 / (2\pi - 0)$</span></p>
<p>Importance sampling works by <em>not</em> uniformly sampling. Instead we try to choose more samples that contribute a lot to the result (important), and less samples that only contribute a little to the result (less important). Hence the name, importance sampling.</p>
<p>If you choose a sampling function whose pdf very closely matches the shape of <span class="math-container">$f$</span>, you can greatly reduce the variance, which means you can take less samples. However, if you choose a sampling function whose value is very different from <span class="math-container">$f$</span>, you can <em><strong>increase</strong></em> the variance. See the picture below:
<a href="https://i.sstatic.net/XGyXh.png" rel="noreferrer"><img src="https://i.sstatic.net/XGyXh.png" alt="Comparison of good sampling vs bad sampling" /></a>
Image from Wojciech Jarosz's <a href="https://www.cs.dartmouth.edu/%7Ewjarosz/publications/dissertation/appendixA.pdf" rel="noreferrer">Dissertation Appendix A</a></p>
<p>One example of importance sampling in Path Tracing is how to choose the direction of a ray after it hits a surface. If the surface is not perfectly specular (ie. a mirror or glass), the outgoing ray can be anywhere in the hemisphere.</p>
<p><a href="https://i.sstatic.net/6crnH.png" rel="noreferrer"><img src="https://i.sstatic.net/6crnH.png" alt="Outgoing rays can go anywhere in the hemisphere" /></a></p>
<p>We <em><strong>could</strong></em> uniformly sample the hemisphere to generate the new ray. However, we can exploit the fact that the rendering equation has a cosine factor in it:</p>
<p><span class="math-container">$$L_{\text{o}}(p, \omega_{\text{o}}) =  L_{e}(p, \omega_{\text{o}}) + \int_{\Omega} f(p, \omega_{\text{i}}, \omega_{\text{o}}) L_{\text{i}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}}$$</span></p>
<p>Specifically, we know that any rays at the horizon will be heavily attenuated (specifically, <span class="math-container">$\cos(x)$</span> ). So, rays generated near the horizon will not contribute very much to the final value.</p>
<p>To combat this, we use importance sampling. If we generate rays according to a cosine weighted hemisphere, we ensure that more rays are generated well above the horizon, and less near the horizon. This will lower variance and reduce noise.</p>
<p>In your case, you specified that you will be using a Cook-Torrance, microfacet-based BRDF. The common form being:</p>
<p><span class="math-container">$$f(p, \omega_{\text{i}}, \omega_{\text{o}}) = \frac{F(\omega_{\text{i}}, h) G(\omega_{\text{i}}, \omega_{\text{o}}, h) D(h)}{4 \cos(\theta_{i}) \cos(\theta_{o})}$$</span></p>
<p>where</p>
<p><span class="math-container">$$F(\omega_{\text{i}}, h) = \text{Fresnel function} \\
G(\omega_{\text{i}}, \omega_{\text{o}}, h) = \text{Geometry Masking and Shadowing function} \\
D(h) = \text{Normal Distribution Function}$$</span></p>
<p>The blog &quot;A Graphic's Guy's Note&quot; has an excellent write up on how to sample Cook-Torrance BRDFs. I will refer you to <a href="https://agraphicsguy.wordpress.com/2015/11/01/sampling-microfacet-brdf/" rel="noreferrer">his blog post</a>. That said, I will try to create a brief overview below:</p>
<p>The NDF is generally the dominant portion of the Cook-Torrance BRDF, so if we are going to importance sample, then we should sample based on the NDF.</p>
<p>Cook-Torrance doesn't specify a specific NDF to use; we are free to choose whichever one suits our fancy. That said, there are a few popular NDFs:</p>
<ul>
<li>GGX</li>
<li>Beckmann</li>
<li>Blinn</li>
</ul>
<p>Each NDF has it's own formula, thus each must be sampled differently. I am only going to show the final sampling function for each. If you would like to see how the formula is derived, see the blog post.</p>
<p><strong>GGX</strong> is defined as:</p>
<p><span class="math-container">$$D_{GGX}(m) = \frac{\alpha^2}{\pi((\alpha^2-1) \cos^2(\theta) + 1)^2}$$</span></p>
<p>To sample the spherical coordinates angle <span class="math-container">$\theta$</span>, we can use the formula:</p>
<p><span class="math-container">$$\theta = \arccos \left( \sqrt{\frac{\alpha^2}{\xi_1 (\alpha^2 - 1) + 1}} \right)$$</span></p>
<p>where <span class="math-container">$\xi$</span> is a uniform random variable.</p>
<p>We assume that the NDF is isotropic, so we can sample <span class="math-container">$\phi$</span> uniformly:</p>
<p><span class="math-container">$$\phi = \xi_{2}$$</span></p>
<p><strong>Beckmann</strong> is defined as:</p>
<p><span class="math-container">$$D_{Beckmann}(m) = \frac{1}{\pi \alpha^2\cos^4(\theta)} e^{-\frac{\tan^2(\theta)}{\alpha^2}}$$</span></p>
<p>Which can be sampled with:</p>
<p><span class="math-container">$$\theta = \arccos \left(\sqrt{\frac{1}{1 = \alpha^2 \ln(1 - \xi_1)}} \right) \\
\phi = \xi_2$$</span></p>
<p>Lastly, <strong>Blinn</strong> is defined as:</p>
<p><span class="math-container">$$D_{Blinn}(m) = \frac{\alpha + 2}{2 \pi} (\cos(\theta))^{\alpha}$$</span></p>
<p>Which can be sampled with:</p>
<p><span class="math-container">$$\theta = \arccos \left(\frac{1}{\xi_{1}^{\alpha + 1}} \right) \\
\phi = \xi_2$$</span></p>
<h2>Putting it in Practice</h2>
<p>Let's look at a basic backwards path tracer:</p>
<pre><code>void RenderPixel(uint x, uint y, UniformSampler *sampler) {
    Ray ray = m_scene-&gt;Camera.CalculateRayFromPixel(x, y, sampler);

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

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

        // The ray missed. Return the background color
        if (ray.geomID == RTC_INVALID_GEOMETRY_ID) {
            color += throughput * float3(0.846f, 0.933f, 0.949f);
            break;
        }

        // We hit an object
        
        // 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 emmisive light
        if (light != nullptr) {
            color += throughput * light-&gt;Le();
        }

        float3 normal = normalize(ray.Ng);
        float3 wo = normalize(-ray.dir);
        float3 surfacePos = ray.org + ray.dir * ray.tfar;

        // Get the new ray direction
        // Choose the direction based on the material
        float3 wi = material-&gt;Sample(wo, normal, sampler);
        float pdf = material-&gt;Pdf(wi, normal);

        // Accumulate the brdf attenuation
        throughput = throughput * material-&gt;Eval(wi, wo, normal) / pdf;
    

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.org = surfacePos;

        // Reset the other ray properties
        ray.dir = wi;
        ray.tnear = 0.001f;
        ray.tfar = embree::inf;
        ray.geomID = RTC_INVALID_GEOMETRY_ID;
        ray.primID = RTC_INVALID_GEOMETRY_ID;
        ray.instID = RTC_INVALID_GEOMETRY_ID;
        ray.mask = 0xFFFFFFFF;
        ray.time = 0.0f;
    }

    m_scene-&gt;Camera.FrameBuffer.SplatPixel(x, y, color);
}
</code></pre>
<p>IE. we bounce around the scene, accumulating color and light attenuation as we go. At each bounce, we have to choose a new direction for the ray. As mentioned above, we <em><strong>could</strong></em> uniformly sample the hemisphere to generate the new ray. However, the code is smarter; it importance samples the new direction based on the BRDF. (Note: This is the input direction, because we are a backwards path tracer)</p>
<pre><code>// Get the new ray direction
// Choose the direction based on the material
float3 wi = material-&gt;Sample(wo, normal, sampler);
float pdf = material-&gt;Pdf(wi, normal);
</code></pre>
<p>Which could be implemented as:</p>
<pre><code>void LambertBRDF::Sample(float3 outputDirection, float3 normal, UniformSampler *sampler) {
    float rand = sampler-&gt;NextFloat();
    float r = std::sqrtf(rand);
    float theta = sampler-&gt;NextFloat() * 2.0f * M_PI;

    float x = r * std::cosf(theta);
    float y = r * std::sinf(theta);

    // Project z up to the unit hemisphere
    float z = std::sqrtf(1.0f - x * x - y * y);

    return normalize(TransformToWorld(x, y, z, normal));
}

float3a TransformToWorld(float x, float y, float z, float3a &amp;normal) {
    // Find an axis that is not parallel to normal
    float3a majorAxis;
    if (abs(normal.x) &lt; 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(1, 0, 0);
    } else if (abs(normal.y) &lt; 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(0, 1, 0);
    } else {
        majorAxis = float3a(0, 0, 1);
    }

    // Use majorAxis to create a coordinate system relative to world space
    float3a u = normalize(cross(normal, majorAxis));
    float3a v = cross(normal, u);
    float3a w = normal;


    // Transform from local coordinates to world coordinates
    return u * x +
           v * y +
           w * z;
}

float LambertBRDF::Pdf(float3 inputDirection, float3 normal) {
    return dot(inputDirection, normal) * M_1_PI;
}
</code></pre>
<p>After we sample the inputDirection ('wi' in the code), we use that to calculate the value of the BRDF. And then we divide by the pdf as per the Monte Carlo formula:</p>
<pre><code>// Accumulate the brdf attenuation
throughput = throughput * material-&gt;Eval(wi, wo, normal) / pdf;
</code></pre>
<p>Where <em>Eval()</em> is just the BRDF function itself (Lambert, Blinn-Phong, Cook-Torrance, etc.):</p>
<pre><code>float3 LambertBRDF::Eval(float3 inputDirection, float3 outputDirection, float3 normal) const override {
    return m_albedo * M_1_PI * dot(inputDirection, normal);
}
</code></pre>

CC BY-SA 3.0
2017-04-13T08:22:53.040Z
Nice answer. The OP also asked about Cook-Torrance importance sampling which this answer doesn't touch upon.
PeteUK
CC BY-SA 3.0
2017-04-13T14:25:05.987Z
I updated the answer to add a section about Cook-Torrance
RichieSams
CC BY-SA 3.0
2017-04-19T15:22:30.663Z
For example GGX, to sample spherical coordinates angle cos(θ) we use the importance sampled formula to calculate the angle and use that in GGX as usual right? Or does the formula replace GGX entirely?
Arjan Singh
CC BY-SA 3.0
2017-04-19T21:26:38.637Z
I added a section to help answer your questions. But, in short, your first method is correct. You use the sampling formula to generate a direction, then you use that new direction in the normal GGX formula *and* to get the pdf for the Monte Carlo formula.
RichieSams
CC BY-SA 3.0
2018-03-13T08:00:17.400Z
>throughput = throughput * material->Eval(wi, wo, normal) / pdf;

Should this line be `throughput + material->...` ?
Cu2S
CC BY-SA 3.0
2018-03-13T15:31:49.263Z
No. Multiplying is correct. The throughput tracks what kinds of surfaces you've already bounced off of. IE, the allowable color of the bounce lighting. For example, if we have pure white light, and bounce off a pure red surface, any remaining bounce lighting will only be red. Because all green and blue light is absorbed by the red surface. Put another way, if we shine pure blue light at a pure red surface, it will be completely black. (This assumes the simplified model of RGB light, where each channel is independent)
RichieSams
CC BY-SA 4.0
2019-01-31T08:33:22.720Z
Wow, wonderful and thorough answer! And I could even understand it!
Benny Jobigan
CC BY-SA 4.0
2019-11-13T19:56:45.297Z
Where does this `radius` variable is coming from in the `LambertBRDF::Sample` function?
dblouis
CC BY-SA 4.0
2019-11-13T20:00:12.933Z
That's an error. Removing
RichieSams
CC BY-SA 3.0
2017-06-13T14:59:53.260Z
For GGX how would I calculate/sample `wi`? I understand how to sample the spherical coordinates angle θ but for the actual direction vector how is that done?
Arjan Singh
CC BY-SA 3.0
2017-06-14T15:05:38.347Z
GGX is assumed to be isotropic so $\phi = \xi_{2}$. With $\phi$ and $\theta$, you can transform to cartesian coordinates to get $w_i$ http://tutorial.math.lamar.edu/Classes/CalcIII/SphericalCoords.aspx
RichieSams
CC BY-SA 4.0
2021-02-03T04:55:01.470Z
@RichieSams I was looking for a solution for another problem and somehow landed here. I have a problem with the way I sample new direction vectors. Could you please take a look at my question [here](https://computergraphics.stackexchange.com/questions/10624/path-tracing-how-to-ensure-the-new-direction-vector-is-a-valid-direction-vector)? I'd really appreciate if you can help me figure out the correct way of sampling new direction vectors ?
Amir
CC BY-SA 4.0
2021-03-26T01:19:53.320Z
This is an excellent answer. Very interesting and informative.
Julien Guertault
CC BY-SA 4.0
2023-09-25T11:10:56.353Z
The best possible explanation is found on the internet.
bim
CC BY-SA 4.0
2017-04-12T22:48:41.963Z
2022-12-06T16:14:20.853Z
2022-12-06T16:14:20.853Z
RichieSams
<p>What is importance sampling? Every article I read about it mentions 'PDF' what is that as well?</p>

<p>From what I gather, importance sampling is a technique to only sample areas on a hemisphere that matter more than others. So, ideally, I should sample rays towards light sources to reduce noise and increase speed. Also, some BRDF's at grazing angles have little no difference in the calculation so using importance sampling to avoid that is good? </p>

<p>If I were to implement importance sampling for a Cook-Torrance BRDF how could I do this? </p>

CC BY-SA 3.0
2017-12-05T21:25:07.703Z
This is link a good read that explains what a PDF is. TL;DR a PDF is a function that describes the probability of (continuous aka floating point) random numbers.  Generating random numbers from a specific PDF can be challenging and there are a few techniques for doing so.  This talks about one of them.  The article after this one talks about another way.
https://blog.demofox.org/2017/08/05/generating-random-numbers-from-a-specific-distribution-by-inverting-the-cdf/
Alan Wolfe
CC BY-SA 4.0
2022-11-26T10:12:05.580Z
To my best knowledge, this blog @AlanWolfe shared is really excellent to understand the whole path tracing concept in easier language
bim
<brdf><importance-sampling>
What is Importance Sampling?
Arjan Singh
# What is Importance Sampling?

(Asked by: username_0 on 2017-04-10)

What is importance sampling? Every article I read about it mentions 'PDF' what is that as well?

From what I gather, importance sampling is a technique to only sample areas on a hemisphere that matter more than others. So, ideally, I should sample rays towards light sources to reduce noise and increase speed. Also, some BRDF's at grazing angles have little no difference in the calculation so using importance sampling to avoid that is good?

If I were to implement importance sampling for a Cook-Torrance BRDF how could I do this?

**Tags:** brdf, importance-sampling

> **Comments:**
> **username_1**: This is link a good read that explains what a PDF is. TL;DR a PDF is a function that describes the probability of (continuous aka floating point) random numbers.  Generating random numbers from a specific PDF can be challenging and there are a few techniques for doing so.  This talks about one of them.  The article after this one talks about another way. https://blog.demofox.org/2017/08/05/generating-random-numbers-from-a-specific-distribution-by-inverting-the-cdf/  
> **username_2**: To my best knowledge, this blog @username_1 shared is really excellent to understand the whole path tracing concept in easier language  

---

## 2 answers

### Answer by username_4 on 2017-04-12. Score: 70 | Accepted answer

# Short answer:

**Importance sampling** is a method to reduce variance in Monte Carlo Integration by choosing an estimator close to the shape of the actual function.

**PDF** is an abbreviation for [Probability Density Function](https://en.wikipedia.org/wiki/Probability_density_function). A $pdf(x)$ gives the probability of a random sample generated being $x$.

# Long Answer:

To start, let's review what Monte Carlo Integration is, and what it looks like mathematically.

Monte Carlo Integration is an technique to estimate the value of an integral. It's typically used when there isn't a closed form solution to the integral. It looks like this:

$$\int f(x) \, dx \approx \frac{1}{N} \sum\_{i=1}^N \frac{f(x\_i)}{pdf(x\_{i})}$$

In english, this says that you can approximate an integral by averaging successive random samples from the function. As $N$ gets large, the approximation gets closer and closer to the solution. $pdf(x\_{i})$ represents the probability density function of each random sample.

Let's do an example: Calculate the value of the integral $I$.

$$I = \int\_{0}^{2\pi} e^{-x} \sin(x) dx$$

Let's use Monte Carlo Integration:

$$I \approx \frac{1}{N} \sum\_{i=1}^N \frac{e^{-x} \sin(x\_i)}{pdf(x\_{i})}$$

A simple python program to calculate this is:

```
import random
import math

N = 200000
TwoPi = 2.0 * math.pi

sum = 0.0

for i in range(N):
    x = random.uniform(0, TwoPi)
    
    fx = math.exp(-x) * math.sin(x)
    pdf = 1 / (TwoPi - 0.0)
    
    sum += fx / pdf

I = (1 / N) * sum
print(I)

```

If we run the program we get $I = 0.4986941$

Using separation by parts, we can get the exact solution:

$$I = \frac{1}{2} (1 − e−2π) = 0.4990663$$

You'll notice that the Monte Carlo Solution is not quite correct. This is because it is an estimate. That said, as $N$ goes to infinity, the estimate should get closer and closer to the correct answer. Already at $N = 2000$ some runs are almost identical to the correct answer.

A note about the PDF: In this simple example, we always take a uniform random sample. A uniform random sample means every sample has the exact same probability of being chosen. We sample in the range $[0, 2\pi]$ so, $pdf(x) = 1 / (2\pi - 0)$

Importance sampling works by *not* uniformly sampling. Instead we try to choose more samples that contribute a lot to the result (important), and less samples that only contribute a little to the result (less important). Hence the name, importance sampling.

If you choose a sampling function whose pdf very closely matches the shape of $f$, you can greatly reduce the variance, which means you can take less samples. However, if you choose a sampling function whose value is very different from $f$, you can ***increase*** the variance. See the picture below:
[![Comparison of good sampling vs bad sampling](https://i.sstatic.net/XGyXh.png)](https://i.sstatic.net/XGyXh.png)
Image from Wojciech Jarosz's [Dissertation Appendix A](https://www.cs.dartmouth.edu/%7Ewjarosz/publications/dissertation/appendixA.pdf)

One example of importance sampling in Path Tracing is how to choose the direction of a ray after it hits a surface. If the surface is not perfectly specular (ie. a mirror or glass), the outgoing ray can be anywhere in the hemisphere.

[![Outgoing rays can go anywhere in the hemisphere](https://i.sstatic.net/6crnH.png)](https://i.sstatic.net/6crnH.png)

We ***could*** uniformly sample the hemisphere to generate the new ray. However, we can exploit the fact that the rendering equation has a cosine factor in it:

$$L\_{\text{o}}(p, \omega\_{\text{o}}) = L\_{e}(p, \omega\_{\text{o}}) + \int\_{\Omega} f(p, \omega\_{\text{i}}, \omega\_{\text{o}}) L\_{\text{i}}(p, \omega\_{\text{i}}) \left | \cos \theta\_{\text{i}} \right | d\omega\_{\text{i}}$$

Specifically, we know that any rays at the horizon will be heavily attenuated (specifically, $\cos(x)$ ). So, rays generated near the horizon will not contribute very much to the final value.

To combat this, we use importance sampling. If we generate rays according to a cosine weighted hemisphere, we ensure that more rays are generated well above the horizon, and less near the horizon. This will lower variance and reduce noise.

In your case, you specified that you will be using a Cook-Torrance, microfacet-based BRDF. The common form being:

$$f(p, \omega\_{\text{i}}, \omega\_{\text{o}}) = \frac{F(\omega\_{\text{i}}, h) G(\omega\_{\text{i}}, \omega\_{\text{o}}, h) D(h)}{4 \cos(\theta\_{i}) \cos(\theta\_{o})}$$

where

$$F(\omega\_{\text{i}}, h) = \text{Fresnel function} \\
G(\omega\_{\text{i}}, \omega\_{\text{o}}, h) = \text{Geometry Masking and Shadowing function} \\
D(h) = \text{Normal Distribution Function}$$

The blog "A Graphic's Guy's Note" has an excellent write up on how to sample Cook-Torrance BRDFs. I will refer you to [his blog post](https://agraphicsguy.wordpress.com/2015/11/01/sampling-microfacet-brdf/). That said, I will try to create a brief overview below:

The NDF is generally the dominant portion of the Cook-Torrance BRDF, so if we are going to importance sample, then we should sample based on the NDF.

Cook-Torrance doesn't specify a specific NDF to use; we are free to choose whichever one suits our fancy. That said, there are a few popular NDFs:

* GGX
* Beckmann
* Blinn

Each NDF has it's own formula, thus each must be sampled differently. I am only going to show the final sampling function for each. If you would like to see how the formula is derived, see the blog post.

**GGX** is defined as:

$$D\_{GGX}(m) = \frac{\alpha^2}{\pi((\alpha^2-1) \cos^2(\theta) + 1)^2}$$

To sample the spherical coordinates angle $\theta$, we can use the formula:

$$\theta = \arccos \left( \sqrt{\frac{\alpha^2}{\xi\_1 (\alpha^2 - 1) + 1}} \right)$$

where $\xi$ is a uniform random variable.

We assume that the NDF is isotropic, so we can sample $\phi$ uniformly:

$$\phi = \xi\_{2}$$

**Beckmann** is defined as:

$$D\_{Beckmann}(m) = \frac{1}{\pi \alpha^2\cos^4(\theta)} e^{-\frac{\tan^2(\theta)}{\alpha^2}}$$

Which can be sampled with:

$$\theta = \arccos \left(\sqrt{\frac{1}{1 = \alpha^2 \ln(1 - \xi\_1)}} \right) \\
\phi = \xi\_2$$

Lastly, **Blinn** is defined as:

$$D\_{Blinn}(m) = \frac{\alpha + 2}{2 \pi} (\cos(\theta))^{\alpha}$$

Which can be sampled with:

$$\theta = \arccos \left(\frac{1}{\xi\_{1}^{\alpha + 1}} \right) \\
\phi = \xi\_2$$

## Putting it in Practice

Let's look at a basic 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);

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

        // The ray missed. Return the background color
        if (ray.geomID == RTC_INVALID_GEOMETRY_ID) {
            color += throughput * float3(0.846f, 0.933f, 0.949f);
            break;
        }

        // We hit an object
        
        // 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 emmisive light
        if (light != nullptr) {
            color += throughput * light->Le();
        }

        float3 normal = normalize(ray.Ng);
        float3 wo = normalize(-ray.dir);
        float3 surfacePos = ray.org + ray.dir * ray.tfar;

        // Get the new ray direction
        // Choose the direction based on the material
        float3 wi = material->Sample(wo, normal, sampler);
        float pdf = material->Pdf(wi, normal);

        // Accumulate the brdf attenuation
        throughput = throughput * material->Eval(wi, wo, normal) / pdf;
    

        // Shoot a new ray

        // Set the origin at the intersection point
        ray.org = surfacePos;

        // Reset the other ray properties
        ray.dir = wi;
        ray.tnear = 0.001f;
        ray.tfar = embree::inf;
        ray.geomID = RTC_INVALID_GEOMETRY_ID;
        ray.primID = RTC_INVALID_GEOMETRY_ID;
        ray.instID = RTC_INVALID_GEOMETRY_ID;
        ray.mask = 0xFFFFFFFF;
        ray.time = 0.0f;
    }

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

```

IE. we bounce around the scene, accumulating color and light attenuation as we go. At each bounce, we have to choose a new direction for the ray. As mentioned above, we ***could*** uniformly sample the hemisphere to generate the new ray. However, the code is smarter; it importance samples the new direction based on the BRDF. (Note: This is the input direction, because we are a backwards path tracer)

```
// Get the new ray direction
// Choose the direction based on the material
float3 wi = material->Sample(wo, normal, sampler);
float pdf = material->Pdf(wi, normal);

```

Which could be implemented as:

```
void LambertBRDF::Sample(float3 outputDirection, float3 normal, UniformSampler *sampler) {
    float rand = sampler->NextFloat();
    float r = std::sqrtf(rand);
    float theta = sampler->NextFloat() * 2.0f * M_PI;

    float x = r * std::cosf(theta);
    float y = r * std::sinf(theta);

    // Project z up to the unit hemisphere
    float z = std::sqrtf(1.0f - x * x - y * y);

    return normalize(TransformToWorld(x, y, z, normal));
}

float3a TransformToWorld(float x, float y, float z, float3a &normal) {
    // Find an axis that is not parallel to normal
    float3a majorAxis;
    if (abs(normal.x) < 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(1, 0, 0);
    } else if (abs(normal.y) < 0.57735026919f /* 1 / sqrt(3) */) {
        majorAxis = float3a(0, 1, 0);
    } else {
        majorAxis = float3a(0, 0, 1);
    }

    // Use majorAxis to create a coordinate system relative to world space
    float3a u = normalize(cross(normal, majorAxis));
    float3a v = cross(normal, u);
    float3a w = normal;

    // Transform from local coordinates to world coordinates
    return u * x +
           v * y +
           w * z;
}

float LambertBRDF::Pdf(float3 inputDirection, float3 normal) {
    return dot(inputDirection, normal) * M_1_PI;
}

```

After we sample the inputDirection ('wi' in the code), we use that to calculate the value of the BRDF. And then we divide by the pdf as per the Monte Carlo formula:

```
// Accumulate the brdf attenuation
throughput = throughput * material->Eval(wi, wo, normal) / pdf;

```

Where *Eval()* is just the BRDF function itself (Lambert, Blinn-Phong, Cook-Torrance, etc.):

```
float3 LambertBRDF::Eval(float3 inputDirection, float3 outputDirection, float3 normal) const override {
    return m_albedo * M_1_PI * dot(inputDirection, normal);
}

```

> **Comments:**
> **username_5**: Nice answer. The OP also asked about Cook-Torrance importance sampling which this answer doesn't touch upon.  
> **username_4**: I updated the answer to add a section about Cook-Torrance  
> **username_0**: For example GGX, to sample spherical coordinates angle cos(θ) we use the importance sampled formula to calculate the angle and use that in GGX as usual right? Or does the formula replace GGX entirely?  
> **username_4**: I added a section to help answer your questions. But, in short, your first method is correct. You use the sampling formula to generate a direction, then you use that new direction in the normal GGX formula *and* to get the pdf for the Monte Carlo formula.  
> **username_0**: For GGX how would I calculate/sample `wi`? I understand how to sample the spherical coordinates angle θ but for the actual direction vector how is that done?  
> **username_4**: GGX is assumed to be isotropic so $\phi = \xi_{2}$. With $\phi$ and $\theta$, you can transform to cartesian coordinates to get $w_i$ http://tutorial.math.lamar.edu/Classes/CalcIII/SphericalCoords.aspx  
> **username_6**: >throughput = throughput * material->Eval(wi, wo, normal) / pdf;  Should this line be `throughput + material->...` ?  
> **username_4**: No. Multiplying is correct. The throughput tracks what kinds of surfaces you've already bounced off of. IE, the allowable color of the bounce lighting. For example, if we have pure white light, and bounce off a pure red surface, any remaining bounce lighting will only be red. Because all green and blue light is absorbed by the red surface. Put another way, if we shine pure blue light at a pure red surface, it will be completely black. (This assumes the simplified model of RGB light, where each channel is independent)  
> **username_7**: Wow, wonderful and thorough answer! And I could even understand it!  
> **username_8**: Where does this `radius` variable is coming from in the `LambertBRDF::Sample` function?  
> **username_4**: That's an error. Removing  
> **username_9**: @username_4 I was looking for a solution for another problem and somehow landed here. I have a problem with the way I sample new direction vectors. Could you please take a look at my question [here](https://computergraphics.stackexchange.com/questions/10624/path-tracing-how-to-ensure-the-new-direction-vector-is-a-valid-direction-vector)? I'd really appreciate if you can help me figure out the correct way of sampling new direction vectors ?  
> **username_10**: This is an excellent answer. Very interesting and informative.  
> **username_2**: The best possible explanation is found on the internet.  

---

### Answer by username_3 on 2017-04-12. Score: 12

If you have a 1D function $f(x)$ and you want to integrate this function from say 0 to 1, one way to perform this integration is by taking N random samples in range [0, 1], evaluate $f(x)$ for each sample and calculate the average of the samples. However, this "naive" Monte Carlo integration is said to "converge slowly", i.e. you need a large number of samples to get close to the ground truth, particularly if the function has high frequencies.

With importance sampling, instead of taking N random samples in [0, 1] range, you take more samples in the "important" regions of $f(x)$ that contribute most to the final result. However, because you bias sampling towards the important regions of the function, these samples must be weighted less to counter the bias, which is where the PDF (probability density function) comes along. PDF tells the probability of a sample at a given position and is used to calculate the weighted average of the samples by dividing each sample with the PDF value at each sample position.

With Cook-Torrance importance sampling the common practice is to distribute samples based on the normal distribution function NDF. If NDF is already normalized, it can serve directly as PDF, which is convenient since it cancels the term out from the BRDF evaluation. Only thing you need to do then is to distribute sample positions based on PDF and evaluate BRDF without the NDF term, i.e.
$$f=\frac{FG}{\pi(n\cdot\omega\_i)(n\cdot\omega\_o)}$$
And calculate the average of the sample results multiplied by the solid angle of the domain you integrate over (e.g. $2\pi$ for hemisphere).

For NDF you need to calculate [Cumulative Distribution Function](https://en.wikipedia.org/wiki/Cumulative_distribution_function) of the PDF to convert uniformly distributed sample position to PDF weighted sample position. For isotropic NDF this simplifies to 1D function because of the symmetry of the function. For more details about the CDF derivation you can check [this old GPU Gems article](http://http.developer.nvidia.com/GPUGems3/gpugems3_ch20.html).


---

