<p>There are several misunderstandings:</p>
<ul>
<li><code>normalize(normal + uniformSphere())</code>: I don't think this is the correct way to implement cosine-weighted hemisphere sampling. Check this out: <a href="https://ameye.dev/notes/sampling-the-hemisphere/" rel="nofollow noreferrer">Sampling the hemisphere</a> and <a href="https://pbr-book.org/3ed-2018/Monte_Carlo_Integration/2D_Sampling_with_Multidimensional_Transformations" rel="nofollow noreferrer">13.6.3 Cosine-Weighted Hemisphere Sampling</a> (this is more involved). The method you implemented by normalizing the sum of normal and uniformly sampled vector is neither uniform hemisphere sampling nor cosine-weighted sampling.</li>
<li>PDF is not <code>cos(theta)</code>. PDF is <code>cos(theta) / pi</code>. If you integrate the PDF over the hemisphere then you will find that the latter one results in 1 (valid). Don't forget to add the term <span class="math-container">$\sin \theta$</span> to account for the measure conversion between solid angle and polar coords.</li>
</ul>
<p><strong>Getting brighter results is caused by not using unbiased sampling method and PDF pairs.</strong> Let's break it down (a little bit of math, but I will try to make it easy to understand):</p>
<p>MC integration is aimed at approximating the following integral <strong>in an unbiased way</strong>:</p>
<p><span class="math-container">$$
\int_{\Omega} f(x)d\mu(x)\tag{1}
$$</span></p>
<p>That is, for the following estimator (let's consider a one-sample case, N = 1):
<span class="math-container">$$
\hat{I} = \underbrace{f(x)}_{\text{evaluation}} / \underbrace{p(x)}_{\text{some PDF you choose to divide}}\tag{2}
$$</span>
If we take its expectation:
<span class="math-container">$$
\mathbb{E}(\hat{I}) = \int_{S}\frac{f(x)}{p(x)}p_{act}(x)d\mu(x)\tag{3}
$$</span></p>
<p>We want <span class="math-container">$p(x)$</span> (the PDF you choose as the denominator) and <span class="math-container">$p_{act}(x)$</span> (the actual PDF of the samples, determined by the sampling method you use) can cancel each other out, so that the result will be exactly <span class="math-container">$(1)$</span>. Now that we know this, for the specific case of yours:</p>
<p>If you choose <span class="math-container">$p(x)$</span> as <span class="math-container">$\cos$</span> or even <span class="math-container">$\cos / \pi$</span>: since the reflected direction given by your code is not drawn from the cosine-weighted PDF, <span class="math-container">$p_{act}$</span> won't cancel itself out with <span class="math-container">$p$</span> in the denominator, which results in bias (hence, brighter). Another extreme example is that if you choose <span class="math-container">$\cos /\pi$</span> as <span class="math-container">$p$</span> and use uniform hemisphere sampling, <span class="math-container">$p$</span> and <span class="math-container">$p_{act}$</span> won't cancel each other out, either. So, to get the desired output, always <strong>choose sampling method (<span class="math-container">$p_{act}$</span>) and the denominator <span class="math-container">$p(x)$</span> wisely.</strong></p>
<p>What you should do here, is to:</p>
<ol>
<li>Implement the local reflection direction sampling, according to <a href="https://ameye.dev/notes/sampling-the-hemisphere/#cosine-weighted-hemisphere" rel="nofollow noreferrer">Cosine-weighted hemisphere</a>.</li>
<li>Rotate the local vector to the global frame: since you know the normal, and your cosine-weighted sample is drawn w.r.t to vector (0, 0, 1): calculate the rotation from (0, 0, 1) to your normal and apply the rotation to your local vector.</li>
<li>Devide <span class="math-container">$f(x)$</span> by <span class="math-container">$\cos\theta / \pi$</span>, the <span class="math-container">$\cos\theta$</span> term should cancel the foreshortening term in <span class="math-container">$f(x)$</span>, so you will end up with a simple path throughput: <span class="math-container">$k_d$</span>.</li>
</ol>

CC BY-SA 4.0
2024-07-22T18:41:17.760Z
The thing that confuses me is that I don't see where the `P act` term comes from. In this case, `f(x)` just shoots a ray in that direction, so for my test scene `f(x) = 0.1`. So regardless of the probability of the ray in any given direction, the resulting samples aren't actually numerically changed, whereas dividing them by the PDF does change the result.
Puppy
CC BY-SA 4.0
2024-07-23T01:19:31.287Z
$p_{act}$ is given by the *act*ual distribution of the sample. The example is given in the post: if you sample the hemisphere by uniform sampling, then $p_{act}$ will be the PDF of uniform hemisphere distribution, regardless what you choose as the denominator.
Enigmatisms
CC BY-SA 4.0
2024-07-23T01:21:29.560Z
So, if you choose uniform hemisphere sampling to generate samples, while divide $f(x)$ by its PDF, you will have unbiased results. The same holds for cosine-weighted sampling. To understand why this is the case, always take the expectation of the stochastic estimator, as in Equation (3).
Enigmatisms
CC BY-SA 4.0
2024-07-23T01:23:02.893Z
You will get $\int f(x) p_{act}(x) d\mu(x)$ as expectation, if you don't account for the actual distribution of samples in your estimator. Yet what we desire, is $\int f(x) d\mu(x)$.
Enigmatisms
CC BY-SA 4.0
2024-07-23T18:38:50.863Z
I don't understand how that works. If I choose $f(x) = 0.1$ then all distributions produce the same result, so there's no guarantee that the actual distribution impacts the final result. Even in a more realistic case where I might get different results, it won't change e.g. the range. If $f(x)$ has the range 0-1 and $cos(\theta)$ has a maximum of 1, then $f(x)p_{act}(x)$ would have the range 0-$\frac{1}{\pi}$. But no matter how I distribute my samples, they still have the range 0-1. So changing the distribution doesn't seem to have the same properties as multiplying the sample value.
Puppy
CC BY-SA 4.0
2024-07-27T03:46:46.363Z
"all distributions produce the same result". This is not true. If you completely remove $p(x)$ (the denominator) and sample from "cosine-weighted hemisphere" and "uniform hemisphere" respectively, you will get very different results. Also, I've never stated that you should "multiplying the sample value", $p_{act}$ is a implicit PDF value, that is determined by your sampling distribution. It is NOT multiplied manually. So in the code, there is no where to find a $p_{act}$.
Enigmatisms
CC BY-SA 4.0
2024-07-27T03:54:55.947Z
See this discussion, I just created a room for this question. Comments are not supposed to be used for extensive discussion. Check this out: [Discussion for unbiased-brdf-sampling](https://chat.stackexchange.com/rooms/154331/unbiased-brdf-sampling)
Enigmatisms
CC BY-SA 4.0
2024-07-21T16:10:52.740Z
2024-07-22T01:20:13.173Z
2024-07-22T01:20:13.173Z
Enigmatisms
<p>It seems like the confusion comes from some assumptions that go unstated. The division by the PDF is given as a general rule for all Monte Carlo, but this doesn't actually seem to be true. In practice, PDFs which cancel out with part of <code>f</code> are used, not any arbitrary PDF, even for functions where this should work like <code>f(x) = 0.1</code>. PBRT is a little unhelpful on this front because first the division by PDF is given, then later on, the parts of <code>f</code> that it should cancel with. So it seems like the correct way of describing it should be that given that <code>f(x)</code> = <code>g(x)p(x)</code>, <em>then</em> sampling <code>g(x)</code> by <code>p(x)</code> and cancelling out <code>p(x)</code> is valid. Otherwise, it is not valid.</p>

CC BY-SA 4.0
2024-08-01T17:18:36.333Z
It is true, $\int_{\Omega}f(x) \,dx \approx \frac{1}{N}\sum_{k=1}^N \frac{f(x_k)}{p(x_k)}$ is valid as long as $x_1, \ldots, x_N$ are independently and identically distributed with density $p$. It doesn't matter whether $f$ was formed as $f(x) = p(x) g(x)$ originally. You can rewrite any function like this anyways.
lightxbulb
CC BY-SA 4.0
2024-08-01T16:08:44.733Z
2024-08-01T16:08:44.733Z
Puppy
<p>I've been looking at Monte Carlo with <a href="https://pbr-book.org/4ed/Monte_Carlo_Integration/Monte_Carlo_Basics#eq:MC-estimator" rel="nofollow noreferrer">PBRT</a>. Currently I'm aiming to finish up the basic MC so that I can move on to multiple importance sampling. For diffuse surfaces I'm calculating the random sample direction with cosine weighting i.e. <code>normalize(normal + uniformSphere())</code>.  PBRT states that I should divide each sample by the PDF, and as I'm using cosine weighting, the PDF should be <code>cos(theta)</code> or simply <code>dot(normal, direction)</code>.  The major issue I'm observing is that this value is always less than 1, so every sample ends up brighter than it should be. The test case I'm using has no actual light sources, only ambient light, so every sample produces the same constant result, however the image now has a large amount of noise where none existed before, because the ambient light can be sampled from many different directions which the PDF weights differently.</p>
<p>How do I make sure that weighting the samples doesn't impact their final brightness? I've seen various sources on this but none of them seem to have any kind of renormalisation which you would normally see for a weighted average.</p>

<pathtracing><monte-carlo>
Applying PDFs brightens all samples in monte carlo
Puppy
# Applying PDFs brightens all samples in monte carlo

(Asked by: username_0 on 2024-07-21)

I've been looking at Monte Carlo with [PBRT](https://pbr-book.org/4ed/Monte_Carlo_Integration/Monte_Carlo_Basics#eq:MC-estimator). Currently I'm aiming to finish up the basic MC so that I can move on to multiple importance sampling. For diffuse surfaces I'm calculating the random sample direction with cosine weighting i.e. `normalize(normal + uniformSphere())`. PBRT states that I should divide each sample by the PDF, and as I'm using cosine weighting, the PDF should be `cos(theta)` or simply `dot(normal, direction)`. The major issue I'm observing is that this value is always less than 1, so every sample ends up brighter than it should be. The test case I'm using has no actual light sources, only ambient light, so every sample produces the same constant result, however the image now has a large amount of noise where none existed before, because the ambient light can be sampled from many different directions which the PDF weights differently.

How do I make sure that weighting the samples doesn't impact their final brightness? I've seen various sources on this but none of them seem to have any kind of renormalisation which you would normally see for a weighted average.

**Tags:** pathtracing, monte-carlo


---

## 2 answers

### Answer by username_1 on 2024-07-21. Score: 3 | Accepted answer

There are several misunderstandings:

* `normalize(normal + uniformSphere())`: I don't think this is the correct way to implement cosine-weighted hemisphere sampling. Check this out: [Sampling the hemisphere](https://ameye.dev/notes/sampling-the-hemisphere/) and [13.6.3 Cosine-Weighted Hemisphere Sampling](https://pbr-book.org/3ed-2018/Monte_Carlo_Integration/2D_Sampling_with_Multidimensional_Transformations) (this is more involved). The method you implemented by normalizing the sum of normal and uniformly sampled vector is neither uniform hemisphere sampling nor cosine-weighted sampling.
* PDF is not `cos(theta)`. PDF is `cos(theta) / pi`. If you integrate the PDF over the hemisphere then you will find that the latter one results in 1 (valid). Don't forget to add the term $\sin \theta$ to account for the measure conversion between solid angle and polar coords.

**Getting brighter results is caused by not using unbiased sampling method and PDF pairs.** Let's break it down (a little bit of math, but I will try to make it easy to understand):

MC integration is aimed at approximating the following integral **in an unbiased way**:

$$
\int\_{\Omega} f(x)d\mu(x)\tag{1}
$$

That is, for the following estimator (let's consider a one-sample case, N = 1):
$$
\hat{I} = \underbrace{f(x)}\_{\text{evaluation}} / \underbrace{p(x)}\_{\text{some PDF you choose to divide}}\tag{2}
$$
If we take its expectation:
$$
\mathbb{E}(\hat{I}) = \int\_{S}\frac{f(x)}{p(x)}p\_{act}(x)d\mu(x)\tag{3}
$$

We want $p(x)$ (the PDF you choose as the denominator) and $p\_{act}(x)$ (the actual PDF of the samples, determined by the sampling method you use) can cancel each other out, so that the result will be exactly $(1)$. Now that we know this, for the specific case of yours:

If you choose $p(x)$ as $\cos$ or even $\cos / \pi$: since the reflected direction given by your code is not drawn from the cosine-weighted PDF, $p\_{act}$ won't cancel itself out with $p$ in the denominator, which results in bias (hence, brighter). Another extreme example is that if you choose $\cos /\pi$ as $p$ and use uniform hemisphere sampling, $p$ and $p\_{act}$ won't cancel each other out, either. So, to get the desired output, always **choose sampling method ($p\_{act}$) and the denominator $p(x)$ wisely.**

What you should do here, is to:

1. Implement the local reflection direction sampling, according to [Cosine-weighted hemisphere](https://ameye.dev/notes/sampling-the-hemisphere/#cosine-weighted-hemisphere).
2. Rotate the local vector to the global frame: since you know the normal, and your cosine-weighted sample is drawn w.r.t to vector (0, 0, 1): calculate the rotation from (0, 0, 1) to your normal and apply the rotation to your local vector.
3. Devide $f(x)$ by $\cos\theta / \pi$, the $\cos\theta$ term should cancel the foreshortening term in $f(x)$, so you will end up with a simple path throughput: $k\_d$.

> **Comments:**
> **username_0**: The thing that confuses me is that I don't see where the `P act` term comes from. In this case, `f(x)` just shoots a ray in that direction, so for my test scene `f(x) = 0.1`. So regardless of the probability of the ray in any given direction, the resulting samples aren't actually numerically changed, whereas dividing them by the PDF does change the result.  
> **username_1**: $p_{act}$ is given by the *act*ual distribution of the sample. The example is given in the post: if you sample the hemisphere by uniform sampling, then $p_{act}$ will be the PDF of uniform hemisphere distribution, regardless what you choose as the denominator.  
> **username_1**: So, if you choose uniform hemisphere sampling to generate samples, while divide $f(x)$ by its PDF, you will have unbiased results. The same holds for cosine-weighted sampling. To understand why this is the case, always take the expectation of the stochastic estimator, as in Equation (3).  
> **username_1**: You will get $\int f(x) p_{act}(x) d\mu(x)$ as expectation, if you don't account for the actual distribution of samples in your estimator. Yet what we desire, is $\int f(x) d\mu(x)$.  
> **username_0**: I don't understand how that works. If I choose $f(x) = 0.1$ then all distributions produce the same result, so there's no guarantee that the actual distribution impacts the final result. Even in a more realistic case where I might get different results, it won't change e.g. the range. If $f(x)$ has the range 0-1 and $cos(\theta)$ has a maximum of 1, then $f(x)p_{act}(x)$ would have the range 0-$\frac{1}{\pi}$. But no matter how I distribute my samples, they still have the range 0-1. So changing the distribution doesn't seem to have the same properties as multiplying the sample value.  
> **username_1**: "all distributions produce the same result". This is not true. If you completely remove $p(x)$ (the denominator) and sample from "cosine-weighted hemisphere" and "uniform hemisphere" respectively, you will get very different results. Also, I've never stated that you should "multiplying the sample value", $p_{act}$ is a implicit PDF value, that is determined by your sampling distribution. It is NOT multiplied manually. So in the code, there is no where to find a $p_{act}$.  
> **username_1**: See this discussion, I just created a room for this question. Comments are not supposed to be used for extensive discussion. Check this out: [Discussion for unbiased-brdf-sampling](https://chat.stackexchange.com/rooms/154331/unbiased-brdf-sampling)  

---

### Answer by username_0 on 2024-08-01. Score: 0

It seems like the confusion comes from some assumptions that go unstated. The division by the PDF is given as a general rule for all Monte Carlo, but this doesn't actually seem to be true. In practice, PDFs which cancel out with part of `f` are used, not any arbitrary PDF, even for functions where this should work like `f(x) = 0.1`. PBRT is a little unhelpful on this front because first the division by PDF is given, then later on, the parts of `f` that it should cancel with. So it seems like the correct way of describing it should be that given that `f(x)` = `g(x)p(x)`, *then* sampling `g(x)` by `p(x)` and cancelling out `p(x)` is valid. Otherwise, it is not valid.

> **Comments:**
> **username_2**: It is true, $\int_{\Omega}f(x) \,dx \approx \frac{1}{N}\sum_{k=1}^N \frac{f(x_k)}{p(x_k)}$ is valid as long as $x_1, \ldots, x_N$ are independently and identically distributed with density $p$. It doesn't matter whether $f$ was formed as $f(x) = p(x) g(x)$ originally. You can rewrite any function like this anyways.  

---

