PDF for sampling emissive meshes for NEE/MIS

Question

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.

Answer 1

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.