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

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?

• It accounts for the term that allows you to rewrite the integral over the hemisphere as an iterated integral. The iterated integral is over a rectangle $[0,2\pi)\times [0,\pi/2]$, this rectangle gets mapped to the hemisphere using spherical coordinates. Now if you take a small square at a point from this rectangle it gets mapped to a tangent parallelogram at the corresponding point on the sphere. The $|\sin\theta|$ term is the area of that parallelogram. It's basically a sirface integration term. You can derive it as $\sqrt{det[J^TJ]}$ where $J$ is the Jacobian of spherical coords wrt angles. Jun 29 at 17:17
• @lightxbulb I'm aware of that, what I'm asking about is why for the MC estimator we don't include the metric tensor of the spherical-to-cartesian mapping. Does this imply we don't use the mapping at all? Jun 29 at 18:30
• It's because you very sloppily failed to define the measure associated with your probability density function. Since the $|\sin\theta|$ is missing, for this to be correct the density must be defined w.r.t. the solid angle measure. For clarity and correctness I would have written instead: $I\approx \frac{1}{N}\sum_{k=1}^N \frac{albedo}{\pi} \frac{\cos\theta_k \sin\theta_k}{p(\phi_k, \theta_k)}$ where $(\phi_k, \theta_k)$ are distributed according to $p$ (wrt $d\phi d\theta$). In fact $p$ does not need to include a sine term if you pick other mappings, and would thus not cancel in general. Jun 29 at 21:21
• @lightxbulb I would really appreciate a fully-fledged answer with detailed explanations, which I'll upvote and accept. I cannot do so with a comment that mostly puts a blame on me. Jun 29 at 23:13

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.

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.

• I'm sorry, but could you please elaborate further and expand the answer? As this answer is currently, there is not enough information to address the question in its entirety, which in turn makes it a link-only answer, and this is against the rules of StackExchange network. Jun 29 at 18:35
• I think my answer is crystal clear. Instead of pointing out why I might be breaking some rules, why don't you just take 5 minutes to read my link and come back? Jun 30 at 9:36
• I also don't like Ocelot's entitled "approach", on the other hand pbrt is a pretty terrible resource when it comes to mathematical motivation - the above link is a good example of this (granted, there are much worse parts in the book). Jul 5 at 0:06