<p>Sampling the entire surface would produce a correct result, even if some contributions are zero, as long as you properly account for those being zero. Consider a similar example in 1D, where you are given a function <span class="math-container">$g$</span> which is nonzero in <span class="math-container">$[a,b]$</span> (i.e. <span class="math-container">$\operatorname{supp}(g)=[a,b]$</span> - the support of <span class="math-container">$g$</span> is <span class="math-container">$[a,b]$</span>). Say you want to numerically estimate the integral <span class="math-container">$\int_a^bg(t)\,dt$</span>. Assume that for some reason you aren't able to generate samples in <span class="math-container">$[a,b]$</span> but can instead generate those in <span class="math-container">$[a,b+\delta],\, \delta&gt;0$</span>. Since <span class="math-container">$g\big((b,b+\delta]\big) = 0$</span> then it follows that:</p>
<p><span class="math-container">$$\int_a^{b+\delta}g(t)\,dt = \int_a^bg(t)\,dt + \int_b^{b+\delta}g(t)\,dt = \int_a^bg(t)\,dt + \int_b^{b+\delta}0\,dt = \int_a^bg(t)\,dt.$$</span></p>
<p>Granted, when you numerically estimate <span class="math-container">$\int_a^{b+\delta}g(t)\,dt$</span> you will have higher variance since samples in <span class="math-container">$(b,b+\delta]$</span> would be used for estimating the zero function.</p>
<p>While the above example is not very realistic, the case in the rendering equation is a realistic one. There the visibility function <span class="math-container">$V(x,y)$</span> is zero over some of the domain of integration, but the non-zero part is hard to explicitly integrate over, since the evaluation of <span class="math-container">$V(x,y)$</span> generally requires a ray-tracing operation.</p>

CC BY-SA 4.0
2021-11-27T22:43:23.477Z
2021-11-27T22:43:23.477Z
lightxbulb
<p>I'm writing a Monte Carlo path tracer, and I'm trying to allow any mesh to be an emitter, but I'm not entirely sure about the probabilities to use when I sample them.</p>
<p>Right now, my algorithm for light sampling does the following:</p>
<ol>
<li>Pick a light <span class="math-container">$L$</span> in the scene with some PDF <span class="math-container">$p_{\mathcal{L}}$</span>.</li>
<li>Pick a random point <span class="math-container">$y$</span> on the mesh with uniform probability <span class="math-container">$p_{\mathcal{S}}(y) = \frac1{S_L}$</span>, where <span class="math-container">$S_L$</span> is the total area of the mesh (computed summing over all the areas of all the triangles composing the mesh).</li>
<li>Cast a shadow ray from the current position <span class="math-container">$x$</span> to the sampled point <span class="math-container">$y$</span>.</li>
<li>Check visibility.</li>
<li>Compute the PDF wrt solid angle (I need it wrt to solid angle for my integrator); the PDF wrt area mesure should now be just <span class="math-container">$p_A(y) = p_{\mathcal{L}}(L) \cdot p_{\mathcal{S}}(y)$</span> (the choice of light and of a point on the light are independent events), hence the actual PDF I need is EDIT: not <span class="math-container">$p_\omega(\omega) = \frac{\cos \theta_y}{\lvert y-x \rvert ^2} p_A(y)$</span> but <span class="math-container">$p_\omega(\omega) = \frac{\lvert y-x \rvert ^2}{\cos \theta_y} p_A(y)$</span> :).</li>
<li>Compute the light contribution to the current point as <span class="math-container">$\frac{E_L \cdot f_r \cdot \cos(\theta_x)}{p_\omega(y)}\cdot \mathrm{thr}$</span>, where <span class="math-container">$\mathrm{thr}$</span> is the ray throughput, <span class="math-container">$E_L$</span> is the emissive value of the light <span class="math-container">$L$</span>, <span class="math-container">$f_r = f_r(-\omega_{\mathrm{incoming}},\omega_{\mathrm{shadow}})$</span> is the BRDF of the surface containing <span class="math-container">$x$</span>, <span class="math-container">$\omega_{\mathrm{incoming}}$</span> is the unit vector of the current ray direction, and <span class="math-container">$\omega_{\mathrm{shadow}} = \frac{y-x}{\lvert y-x \rvert ^2}$</span> is the direction of the shadow ray.</li>
<li>Use that light contribution for next event estimation; use <span class="math-container">$p_\omega$</span> as my PDF for the NEE when computing the weights for multiple importance sampling.</li>
</ol>
<p>I think this should account correctly for the fact that only the projected area of my mesh on the plane of the triangle <span class="math-container">$x$</span> lies on should contribute to the illumination of <span class="math-container">$x$</span>, as the visibility check rules out all the triangles &quot;hidden&quot; by the geometry of the mesh, but I'd like to be sure that I'm not forgetting anything.</p>
<p>EDIT: In particular, the thing that makes me a bit unsure is using the entire surface area of the mesh, even though the part facing the point <span class="math-container">$x$</span> might be much smaller.</p>

CC BY-SA 4.0
2021-11-26T11:32:54.793Z
How do you pick a point uniformly on the triangular mesh? Also $p_\omega(\omega)\,d\omega = p_\omega(\omega)\frac{\cos\theta_y}{\|y-x\|^2}\,dA = p_A(y)\,dA$ and then $p_\omega(\omega) = \frac{\|y-x\|^2}{\cos\theta_y}p_A(y)$, so your pdf was wrong unless I'm trippin.
lightxbulb
CC BY-SA 4.0
2021-11-26T11:48:09.160Z
Entire surface is correct since you filter out 0 contribution through the visibility function. It is like having a function $g$ that is nonzero in $[a,b]$ (i.e. $\operatorname{supp}(f) = [a,b]$) and instead of numerically approximating $\int_{a}^{b}g(t)\,dt$ you approximate $\int_{a}^{b+\delta}g(t)\,dt= \int_{a}^{b}g(t)\,dt + \int_{b}^{b+\delta}g(t)\,dt =\int_a^b g(t)\,dt$, with $\delta>0$. But as you can see those are equal (since $g((b,b+\delta]) = 0$) Granted you'll have higher variance when approximating $\int_a^{b+\delta}g(t)\,dt$ because some samples will end up with zero contribution.
lightxbulb
CC BY-SA 4.0
2021-11-27T21:13:36.753Z
@lightxbulb 2) to randomly sample the mesh, I select a random number $r \in [0,S_L)$, then select a triangle with the inversion method (so that the PDF for picking a triangle is weighted by the surface of each triangle). I have two arrays pre-computed at mesh instantiation: one with the surfaces of each triangle, and one with their cumulative sums (that's the CDF I'm inverting). I do a binary search for $r$ to get the biggest value $\leq r$ in the array with the cumulative sums and pick the corresponding triangle. Then I uniformly sample the triangle.
uhwo
CC BY-SA 4.0
2021-11-27T21:35:12.533Z
`so that the PDF for picking a triangle is weighted by the surface of each triangle` - I was worried that may not have been the case. It's all good if you're doing this (it's called inverse transform sampling btw). Note that if you non-uniformly rescale your mesh, the CDF becomes invalid. Also I hope you properly uniformly sample the triangle - there's some math in the Global Illumination Compendium, in PBRT, in ray-tracing gems, and in some other places.
lightxbulb
<pathtracing><monte-carlo><importance-sampling><pdf>
PDF for sampling emissive meshes for NEE/MIS
uhwo
# PDF for sampling emissive meshes for NEE/MIS

(Asked by: username_0 on 2021-11-23)

I'm writing a Monte Carlo path tracer, and I'm trying to allow any mesh to be an emitter, but I'm not entirely sure about the probabilities to use when I sample them.

Right now, my algorithm for light sampling does the following:

1. Pick a light $L$ in the scene with some PDF $p\_{\mathcal{L}}$.
2. Pick a random point $y$ on the mesh with uniform probability $p\_{\mathcal{S}}(y) = \frac1{S\_L}$, where $S\_L$ is the total area of the mesh (computed summing over all the areas of all the triangles composing the mesh).
3. Cast a shadow ray from the current position $x$ to the sampled point $y$.
4. Check visibility.
5. Compute the PDF wrt solid angle (I need it wrt to solid angle for my integrator); the PDF wrt area mesure should now be just $p\_A(y) = p\_{\mathcal{L}}(L) \cdot p\_{\mathcal{S}}(y)$ (the choice of light and of a point on the light are independent events), hence the actual PDF I need is EDIT: not $p\_\omega(\omega) = \frac{\cos \theta\_y}{\lvert y-x \rvert ^2} p\_A(y)$ but $p\_\omega(\omega) = \frac{\lvert y-x \rvert ^2}{\cos \theta\_y} p\_A(y)$ :).
6. Compute the light contribution to the current point as $\frac{E\_L \cdot f\_r \cdot \cos(\theta\_x)}{p\_\omega(y)}\cdot \mathrm{thr}$, where $\mathrm{thr}$ is the ray throughput, $E\_L$ is the emissive value of the light $L$, $f\_r = f\_r(-\omega\_{\mathrm{incoming}},\omega\_{\mathrm{shadow}})$ is the BRDF of the surface containing $x$, $\omega\_{\mathrm{incoming}}$ is the unit vector of the current ray direction, and $\omega\_{\mathrm{shadow}} = \frac{y-x}{\lvert y-x \rvert ^2}$ is the direction of the shadow ray.
7. Use that light contribution for next event estimation; use $p\_\omega$ as my PDF for the NEE when computing the weights for multiple importance sampling.

I think this should account correctly for the fact that only the projected area of my mesh on the plane of the triangle $x$ lies on should contribute to the illumination of $x$, as the visibility check rules out all the triangles "hidden" by the geometry of the mesh, but I'd like to be sure that I'm not forgetting anything.

EDIT: In particular, the thing that makes me a bit unsure is using the entire surface area of the mesh, even though the part facing the point $x$ might be much smaller.

**Tags:** pathtracing, monte-carlo, importance-sampling, pdf

> **Comments:**
> **username_1**: How do you pick a point uniformly on the triangular mesh? Also $p_\omega(\omega)\,d\omega = p_\omega(\omega)\frac{\cos\theta_y}{\|y-x\|^2}\,dA = p_A(y)\,dA$ and then $p_\omega(\omega) = \frac{\|y-x\|^2}{\cos\theta_y}p_A(y)$, so your pdf was wrong unless I'm trippin.  
> **username_1**: Entire surface is correct since you filter out 0 contribution through the visibility function. It is like having a function $g$ that is nonzero in $[a,b]$ (i.e. $\operatorname{supp}(f) = [a,b]$) and instead of numerically approximating $\int_{a}^{b}g(t)\,dt$ you approximate $\int_{a}^{b+\delta}g(t)\,dt= \int_{a}^{b}g(t)\,dt + \int_{b}^{b+\delta}g(t)\,dt =\int_a^b g(t)\,dt$, with $\delta>0$. But as you can see those are equal (since $g((b,b+\delta]) = 0$) Granted you'll have higher variance when approximating $\int_a^{b+\delta}g(t)\,dt$ because some samples will end up with zero contribution.  
> **username_0**: @username_1 2) to randomly sample the mesh, I select a random number $r \in [0,S_L)$, then select a triangle with the inversion method (so that the PDF for picking a triangle is weighted by the surface of each triangle). I have two arrays pre-computed at mesh instantiation: one with the surfaces of each triangle, and one with their cumulative sums (that's the CDF I'm inverting). I do a binary search for $r$ to get the biggest value $\leq r$ in the array with the cumulative sums and pick the corresponding triangle. Then I uniformly sample the triangle.  
> **username_1**: `so that the PDF for picking a triangle is weighted by the surface of each triangle` - I was worried that may not have been the case. It's all good if you're doing this (it's called inverse transform sampling btw). Note that if you non-uniformly rescale your mesh, the CDF becomes invalid. Also I hope you properly uniformly sample the triangle - there's some math in the Global Illumination Compendium, in PBRT, in ray-tracing gems, and in some other places.  

---

## 1 answers

### Answer by username_1 on 2021-11-27. Score: 0 | Accepted answer

Sampling the entire surface would produce a correct result, even if some contributions are zero, as long as you properly account for those being zero. Consider a similar example in 1D, where you are given a function $g$ which is nonzero in $[a,b]$ (i.e. $\operatorname{supp}(g)=[a,b]$ - the support of $g$ is $[a,b]$). Say you want to numerically estimate the integral $\int\_a^bg(t)\,dt$. Assume that for some reason you aren't able to generate samples in $[a,b]$ but can instead generate those in $[a,b+\delta],\, \delta>0$. Since $g\big((b,b+\delta]\big) = 0$ then it follows that:

$$\int\_a^{b+\delta}g(t)\,dt = \int\_a^bg(t)\,dt + \int\_b^{b+\delta}g(t)\,dt = \int\_a^bg(t)\,dt + \int\_b^{b+\delta}0\,dt = \int\_a^bg(t)\,dt.$$

Granted, when you numerically estimate $\int\_a^{b+\delta}g(t)\,dt$ you will have higher variance since samples in $(b,b+\delta]$ would be used for estimating the zero function.

While the above example is not very realistic, the case in the rendering equation is a realistic one. There the visibility function $V(x,y)$ is zero over some of the domain of integration, but the non-zero part is hard to explicitly integrate over, since the evaluation of $V(x,y)$ generally requires a ray-tracing operation.


---

