Where does sin(theta) go in estimators of The Rendering Equation?

Question

Consider this equation, where we have a surface integral over hemisphere with Lambertian BRDF and cos(theta) from The Rendering Equation: $$ \int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}}\frac{albedo}{\pi}cos(\theta)sin(\theta)d\theta d\phi = 1 $$ However, in the code, we would usually do this, employing a Monte-Carlo estimator: $$ \frac{1}{N}\sum_{i=1}^{N}cos(\theta)\frac{albedo}{\pi}\frac{1}{p(x_i)} $$ where only cosine remains, but where does sin(theta) go?

Answer 1

Assume you want to estimate the following integral: $$I = \int_{\Omega} f d\mu$$ over the set $\Omega\subseteq \mathbb{R}^n$. Suppose you are given a bijective and continuously differentiable parametrization $\varphi^{-1}:D\to\Omega$ where $D=[a^1,b^1]\times\ldots\times [a^m,b^m]$ whose Jacobian has linearly independent columns almost everywhere (i.e. except on a set of measure zero). That is, we have an almost everywhere regular $m$-manifold embedded in $\mathbb{R}^n$: $$\Omega = \varphi^{-1}(D) = \{\varphi^{-1}(t)\in\mathbb{R}^n\,:\, t\in D\}.$$ You can now rewrite the integral using $m$ iterated integrals as follows: $$I = \int_{\Omega} f\,d\mu = \int_{a^1}^{b^1}\ldots \int_{a^m}^{b^m} f(\varphi^{-1}(t))\sqrt{det[J_{\varphi^{-1}}^T(t)J_{\varphi^{-1}}(t)]}\,dt^m\ldots dt^1.$$ The average value of the integrand: $$g(t) = f(\varphi^{-1}(t))\sqrt{det[J_{\varphi^{-1}}^T(t)J_{\varphi^{-1}}(t)]}$$ over $D$ can then be defined as $g_{avg} = I / |D|$, where $|D|$ is the measure/area/volume of $D$: $|D| =\prod_{j=1}^m(b^j-a^j)$. If we have a sequence of points $t_1,\ldots,t_N\in D$ that uniformly fills the box $D$ then we may estimate the integral using the following expression: $$I = |D|g_{avg} \approx |D|\frac{1}{N}\sum_{j=1}^N g(t_j).$$ Using random independent uniformly distributed points in $D$ you get the Monte Carlo estimator. Since the sampling is uniform over $D$, you can interpret the density as $p(t) = \frac{1}{|D|}$: $$I = |D|g_{avg} \approx |D|\frac{1}{N}\sum_{j=1}^N g(t_j)= \frac{1}{N}\sum_{j=1}^N \frac{g(t_j)}{p(t_j)}.$$ Generally, estimators resulting from integrands that are flatter (i.e. $g$ is close to $g_{avg}$ without much wiggling) have lower variance. This suggests that we may reduce the error for a fixed number of samples by applying a change of variables of our integral in order to flatten the integrand and get a "better" estimator. That is, assume that you have a bijective and continuously differentiable map $\psi:D\to D$ that has a non-singular Jacobian almost everywhere, such that $|det[J_{\psi}]|\propto |g|$ or at the very least it is close to proportional to $|g|$. Then by the inverse function theorem $J_{\psi^{-1}} = J^{-1}_{\psi}$ and from $det[A^{-1}] = (det[A])^{-1}$ we get $$|det[J_{\psi^{-1}}]| = |det[J^{-1}_{\psi}]| = \frac{1}{|det[J_{\psi}]|} \propto \frac{1}{|g|}.$$ Thus $|det[J_{\psi^{-1}}||g|$ is close to constant, and in the case where $g$ is non-negative we may also drop the absolute value around $g$. Let us look at the change of variables: $$\int_D g(t)\, dt = \int_D g(\psi^{-1}(s))|det[J_{\psi^{-1}}(s)]| \, ds.$$ We can now also compare the two estimators: $$|D|\frac{1}{N} \sum_{j=1}^N g(t_j), \quad |D|\frac{1}{N}\sum_{j=1}^N g(\psi^{-1}(s_j)) |det[J_{\psi^{-1}}(s_j)]|.$$ If $\psi$ was chosen as described then the latter estimator has a product close to constant in the sum, and in the ideal case it is constant - then only 1 sample is required to find the true solution. In probability and statistics the usage of such a map $\psi^{-1}$ is known as importance sampling, however note that this concept is more general and doesn't actually require any probabilistic considerations. Nevertheless, in the probabilistic setting we can interpret the probability density function (pdf) in the first case as being constant $p_1(t) = \frac{1}{|D|}$ and in the second case the pdf is $p_2(t) = \frac{1}{|D||det[J_{\psi^{-1}}(t)]|}$ where this density can be achieved by first generating uniform samples $s_1,\ldots,s_N\in D$ and then remapping them as $t_k = \psi^{-1}(s_k)$. The two estimators can then be rewritten as: $$\frac{1}{N} \sum_{j=1}^N \frac{g(t_j)}{p_1(t_j)}, \quad \frac{1}{N}\sum_{j=1}^N \frac{g(t_j)}{p_2(t_j)}.$$ Note that the $t_1,\ldots,t_N$ refer to points generated in a different manner for the 2 cases, in the first they are generated according to a density $p_1$, while in the second according to $p_2$.

Now consider the special case where $\Omega$ is the unit hemisphere around some normal vector $n$. Let $t,b,n$ be an arbitrary orthonormal basis ($n$ is fixed), then we may form the rotation matrix $R=\begin{bmatrix} t & b & n\end{bmatrix}$. We can parametrize the unit hemisphere around $n$ using $$\varphi^{-1}(\phi,\theta) = R\begin{bmatrix}\cos\phi\sin\theta\\\sin\phi\sin\theta\\\cos\theta\end{bmatrix}, \quad \phi\in[0,2\pi), \, \theta \in [0,\pi/2].$$ That is $\Omega = \varphi^{-1}(D)$ where $D=[0,2\pi)\times[0,\pi/2]$. We may compute the Jacobian of this map as follows: $$J_{\varphi^{-1}}(\phi,\theta) = R \begin{bmatrix} -\sin\phi\sin\theta & \cos\phi\cos\theta \\ \cos\phi \sin\theta & \sin\phi \cos\theta \\ 0 & -\sin\theta \end{bmatrix}.$$

It is clear that the Jacobian has linearly independent columns for $\theta\in(0,\pi/2]$ and those are linearly dependent only for $\theta = 0$ or $\theta = \pi$, so as desired the map is regular almost everywhere. You can now compute $J^TJ$ and the square root of its determinant as follows: $$J_{\varphi^{-1}}^T J_{\varphi^{-1}}= \begin{bmatrix} \sin^2\theta & 0 \\ 0 & 1 \end{bmatrix} \implies \sqrt{det[J^TJ]} = |\sin\theta|.$$ Then the integral over the hemisphere can be rewritten as: $$I = \int_{\Omega} f\,d\mu = \int_0^{2\pi}\int_0^{\pi/2} f(\varphi^{-1}(\phi, \theta))\sin\theta\,d\theta d\phi.$$ We are allowed to drop the absolute value around the sine since it is non-negative in $[0,\pi/2]$. If you were to directly construct an estimator from this you would get: $$\frac{\pi^2}{N}\sum_{j=1}^{N}f(\varphi^{-1}(\phi_j,\theta_j))\sin\theta_j = \frac{1}{N} \sum_{j=1}^N \frac{f(\varphi^{-1}(\phi_j,\theta_j))\sin\theta_j}{p_1(\phi_j, \theta_j)},\, p_1(\phi,\theta) = \frac{1}{\pi^2}.$$ This is the estimator for uniform sampling $D$. As discussed, ideally we want to flatten the integrand. One way to do this here (assuming $f$ is not inversely proportional to $\sin\theta$) would be to cancel out the $\sin\theta$ by using a change of variables with a map $\psi$ such $det[J_{\psi}] = \sin\theta$. Such a map $\psi$ can for example be constructed using inverse transform sampling, but there are many more ways to do so. The inverse transform sampling approach results in $\psi^{-1}(\phi,\theta) = (\phi, \arccos(1-2\theta/\pi))$, and you can verify that $det[J_{\psi^{-1}}(\phi,\theta)] = \frac{2}{\pi\sin\theta}$ as desired, then $p_2(\phi, \theta) = \frac{1}{|D||det[J_{\psi^{-1}}]|} = \frac{\sin\theta}{2\pi}$, and thus the $\sin\theta$ gets cancelled in the new estimator. When you have $\cos\theta\sin\theta$ people usually come up with a map $\psi$ that cancels out both. Similarly if there is a complex brdf people typically come up with a map that would cancel out some terms of it in order to get a better estimator.

A minor remark - in the literature one usually picks $D$ to be the unit box. Another remark is that in the literature they sometimes write the pdf w.r.t. different measures: $$p_{\omega}(\theta) d\omega = p_{\omega}(\theta)|\sin\theta|\,d\theta d\phi = p_{\phi,\theta}(\theta)\,d\theta d\phi.$$ From there you see that the relationship is $p_{\phi,\theta}(\theta) = p_{\omega}(\theta)|\sin\theta|$. So for example a constant pdf w.r.t. the solid angle measure is translated into a sine weighted pdf w.r.t. $\phi,\theta$. This corresponds to the $\arccos$ map I mentioned. A constant pdf w.r.t. $\phi, \theta$ corresponds to a uniform distribution in the rectangle $[0,2\pi)\times[0,\pi/2]$, but it is not uniform on the hemisphere. On the other hand a constant pdf wrt solid angle corresponds to a uniform distribution on the unit hemisphere but not in the rectangle $D$. In general if you have a parametrization $\varphi^{-1}$ of some manifold such that $d\mu = \sqrt{det[J^T_{\varphi^{-1}}(t)J_{\varphi^{-1}}(t)]}\,dt$, then a constant pdf wrt $d\mu$ will give you a uniform distribution on the manifold and will cancel out this square root of the determinant term. Of course it is not always trivial constructing a map $\psi^{-1}$ that produces the desired density. Things get even uglier when you consider the area formulation of the rendering equation and you have to compute weights for multiple importance sampling.

Answer 2

Your second equation is estimated with respect to solid angle $\Omega$, and the Jacobian of the transformation from spherical coordinates to solid angles is equal to $\sin\theta$. Since $\mathrm{d}\omega = \sin\theta\,\mathrm{d}\theta\mathrm{d}\phi$, you have that

$$ p(\theta,\phi) = \sin\theta\, p(\omega), $$

so you are only left with the cosine.

See PBRT for a more detailed explanation.