<p>Assume you want to estimate the following integral:
<span class="math-container">$$I = \int_{\Omega} f d\mu$$</span>
over the set <span class="math-container">$\Omega\subseteq \mathbb{R}^n$</span>. Suppose you are given a bijective and continuously differentiable parametrization <span class="math-container">$\varphi^{-1}:D\to\Omega$</span> where <span class="math-container">$D=[a^1,b^1]\times\ldots\times [a^m,b^m]$</span> 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 <span class="math-container">$m$</span>-manifold embedded in <span class="math-container">$\mathbb{R}^n$</span>:
<span class="math-container">$$\Omega = \varphi^{-1}(D) = \{\varphi^{-1}(t)\in\mathbb{R}^n\,:\, t\in D\}.$$</span>
You can now rewrite the integral using <span class="math-container">$m$</span> iterated integrals as follows:
<span class="math-container">$$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.$$</span>
The average value of the integrand:
<span class="math-container">$$g(t) = f(\varphi^{-1}(t))\sqrt{det[J_{\varphi^{-1}}^T(t)J_{\varphi^{-1}}(t)]}$$</span>
over <span class="math-container">$D$</span> can then be defined as <span class="math-container">$g_{avg} = I / |D|$</span>, where <span class="math-container">$|D|$</span> is the measure/area/volume of <span class="math-container">$D$</span>: <span class="math-container">$|D| =\prod_{j=1}^m(b^j-a^j)$</span>. If we have a sequence of points <span class="math-container">$t_1,\ldots,t_N\in D$</span> that uniformly fills the box <span class="math-container">$D$</span> then we may estimate the integral using the following expression:
<span class="math-container">$$I = |D|g_{avg} \approx |D|\frac{1}{N}\sum_{j=1}^N g(t_j).$$</span>
Using random independent uniformly distributed points in <span class="math-container">$D$</span> you get the Monte Carlo estimator. Since the sampling is uniform over <span class="math-container">$D$</span>, you can interpret the density as <span class="math-container">$p(t) = \frac{1}{|D|}$</span>:
<span class="math-container">$$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)}.$$</span>
Generally, estimators resulting from integrands that are flatter (i.e. <span class="math-container">$g$</span> is close to <span class="math-container">$g_{avg}$</span> 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 &quot;better&quot; estimator. That is, assume that you have a bijective and continuously differentiable map <span class="math-container">$\psi:D\to D$</span> that has a non-singular Jacobian almost everywhere, such that <span class="math-container">$|det[J_{\psi}]|\propto |g|$</span> or at the very least it is close to proportional to <span class="math-container">$|g|$</span>. Then by the inverse function theorem <span class="math-container">$J_{\psi^{-1}} = J^{-1}_{\psi}$</span> and from <span class="math-container">$det[A^{-1}] = (det[A])^{-1}$</span> we get
<span class="math-container">$$|det[J_{\psi^{-1}}]| = |det[J^{-1}_{\psi}]| = \frac{1}{|det[J_{\psi}]|} \propto \frac{1}{|g|}.$$</span>
Thus <span class="math-container">$|det[J_{\psi^{-1}}||g|$</span> is close to constant, and in the case where <span class="math-container">$g$</span> is non-negative we may also drop the absolute value around <span class="math-container">$g$</span>. Let us look at the change of variables:
<span class="math-container">$$\int_D g(t)\, dt = \int_D g(\psi^{-1}(s))|det[J_{\psi^{-1}}(s)]| \, ds.$$</span>
We can now also compare the two estimators:
<span class="math-container">$$|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)]|.$$</span>
If <span class="math-container">$\psi$</span> 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 <span class="math-container">$\psi^{-1}$</span> 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 <span class="math-container">$p_1(t) = \frac{1}{|D|}$</span> and in the second case the pdf is <span class="math-container">$p_2(t) = \frac{1}{|D||det[J_{\psi^{-1}}(t)]|}$</span> where this density can be achieved by first generating uniform samples <span class="math-container">$s_1,\ldots,s_N\in D$</span> and then remapping them as <span class="math-container">$t_k = \psi^{-1}(s_k)$</span>. The two estimators can then be rewritten as:
<span class="math-container">$$\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)}.$$</span>
Note that the <span class="math-container">$t_1,\ldots,t_N$</span> refer to points generated in a different manner for the 2 cases, in the first they are generated according to a density <span class="math-container">$p_1$</span>, while in the second according to <span class="math-container">$p_2$</span>.</p>
<p>Now consider the special case where <span class="math-container">$\Omega$</span> is the unit hemisphere around some normal vector <span class="math-container">$n$</span>. Let <span class="math-container">$t,b,n$</span> be an arbitrary orthonormal basis (<span class="math-container">$n$</span> is fixed), then we may form the rotation matrix <span class="math-container">$R=\begin{bmatrix} t &amp; b &amp; n\end{bmatrix}$</span>. We can parametrize the unit hemisphere around <span class="math-container">$n$</span> using
<span class="math-container">$$\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].$$</span>
That is <span class="math-container">$\Omega = \varphi^{-1}(D)$</span> where <span class="math-container">$D=[0,2\pi)\times[0,\pi/2]$</span>. We may compute the Jacobian of this map as follows:
<span class="math-container">$$J_{\varphi^{-1}}(\phi,\theta) = R \begin{bmatrix} -\sin\phi\sin\theta &amp; \cos\phi\cos\theta \\ \cos\phi \sin\theta &amp; \sin\phi \cos\theta \\ 0 &amp; -\sin\theta \end{bmatrix}.$$</span></p>
<p>It is clear that the Jacobian has linearly independent columns for <span class="math-container">$\theta\in(0,\pi/2]$</span> and those are linearly dependent only for <span class="math-container">$\theta = 0$</span> or <span class="math-container">$\theta = \pi$</span>, so as desired the map is regular almost everywhere.
You can now compute <span class="math-container">$J^TJ$</span> and the square root of its determinant as follows:
<span class="math-container">$$J_{\varphi^{-1}}^T J_{\varphi^{-1}}=
\begin{bmatrix} 
\sin^2\theta &amp; 0 \\ 0 &amp; 1
\end{bmatrix} \implies \sqrt{det[J^TJ]} = |\sin\theta|.$$</span>
Then the integral over the hemisphere can be rewritten as:
<span class="math-container">$$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.$$</span>
We are allowed to drop the absolute value around the sine since it is non-negative in <span class="math-container">$[0,\pi/2]$</span>.
If you were to directly construct an estimator from this you would get:
<span class="math-container">$$\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}.$$</span>
This is the estimator for uniform sampling <span class="math-container">$D$</span>. As discussed, ideally we want to flatten the integrand. One way to do this here (assuming <span class="math-container">$f$</span> is not inversely proportional to <span class="math-container">$\sin\theta$</span>) would be to cancel out the <span class="math-container">$\sin\theta$</span> by using a change of variables with a map <span class="math-container">$\psi$</span> such <span class="math-container">$det[J_{\psi}] = \sin\theta$</span>. Such a map <span class="math-container">$\psi$</span> 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 <span class="math-container">$\psi^{-1}(\phi,\theta) = (\phi, \arccos(1-2\theta/\pi))$</span>, and you can verify that <span class="math-container">$det[J_{\psi^{-1}}(\phi,\theta)] = \frac{2}{\pi\sin\theta}$</span> as desired, then <span class="math-container">$p_2(\phi, \theta) = \frac{1}{|D||det[J_{\psi^{-1}}]|} = \frac{\sin\theta}{2\pi}$</span>, and thus the <span class="math-container">$\sin\theta$</span> gets cancelled in the new estimator. When you have <span class="math-container">$\cos\theta\sin\theta$</span> people usually come up with a map <span class="math-container">$\psi$</span> 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.</p>
<p>A minor remark - in the literature one usually picks <span class="math-container">$D$</span> to be the unit box. Another remark is that in the literature they sometimes write the pdf w.r.t. different measures:
<span class="math-container">$$p_{\omega}(\theta) d\omega = p_{\omega}(\theta)|\sin\theta|\,d\theta d\phi = p_{\phi,\theta}(\theta)\,d\theta d\phi.$$</span>
From there you see that the relationship is <span class="math-container">$p_{\phi,\theta}(\theta) = p_{\omega}(\theta)|\sin\theta|$</span>. So for example a constant pdf w.r.t. the solid angle measure is translated into a sine weighted pdf w.r.t. <span class="math-container">$\phi,\theta$</span>. This corresponds to the <span class="math-container">$\arccos$</span> map I mentioned. A constant pdf w.r.t. <span class="math-container">$\phi, \theta$</span> corresponds to a uniform distribution in the rectangle <span class="math-container">$[0,2\pi)\times[0,\pi/2]$</span>, 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 <span class="math-container">$D$</span>. In general if you have a parametrization <span class="math-container">$\varphi^{-1}$</span> of some manifold such that <span class="math-container">$d\mu = \sqrt{det[J^T_{\varphi^{-1}}(t)J_{\varphi^{-1}}(t)]}\,dt$</span>, then a constant pdf wrt <span class="math-container">$d\mu$</span> 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 <span class="math-container">$\psi^{-1}$</span> 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.</p>

CC BY-SA 4.0
2023-06-30T01:55:22.867Z
2023-06-30T14:14:55.830Z
2023-06-30T14:14:55.830Z
lightxbulb
<p>Your second equation is estimated with respect to solid angle <span class="math-container">$\Omega$</span>, and the Jacobian of the transformation from spherical coordinates to solid angles is equal to <span class="math-container">$\sin\theta$</span>. Since <span class="math-container">$\mathrm{d}\omega = \sin\theta\,\mathrm{d}\theta\mathrm{d}\phi$</span>, you have that</p>
<p><span class="math-container">$$
p(\theta,\phi) = \sin\theta\, p(\omega),
$$</span></p>
<p>so you are only left with the cosine.</p>
<p>See <a href="https://pbr-book.org/3ed-2018/Monte_Carlo_Integration/Transforming_between_Distributions" rel="nofollow noreferrer">PBRT</a> for a more detailed explanation.</p>

CC BY-SA 4.0
2023-06-29T18:35:41.790Z
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.
Ocelot
CC BY-SA 4.0
2023-06-30T09:36:55.200Z
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?
Hubble
CC BY-SA 4.0
2023-07-05T00:06:18.530Z
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).
lightxbulb
CC BY-SA 4.0
2023-06-29T15:50:44.797Z
2023-06-29T15:50:44.797Z
Hubble
<p>Consider this equation, where we have a surface integral over hemisphere with Lambertian BRDF and cos(theta) from The Rendering Equation:
<span class="math-container">$$
\int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}}\frac{albedo}{\pi}cos(\theta)sin(\theta)d\theta d\phi = 1
$$</span>
However, in the code, we would usually do this, employing a Monte-Carlo estimator:
<span class="math-container">$$
\frac{1}{N}\sum_{i=1}^{N}cos(\theta)\frac{albedo}{\pi}\frac{1}{p(x_i)}
$$</span>
where only cosine remains, but where does sin(theta) go?</p>

CC BY-SA 4.0
2023-06-29T17:17:51.660Z
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.
lightxbulb
CC BY-SA 4.0
2023-06-29T18:30:15.970Z
@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?
Ocelot
CC BY-SA 4.0
2023-06-29T21:21:11.677Z
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.
lightxbulb
CC BY-SA 4.0
2023-06-29T23:13:08.450Z
@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.
Ocelot
<raytracing><mathematics><integral>
Where does sin(theta) go in estimators of The Rendering Equation?
Ocelot
# Where does sin(theta) go in estimators of The Rendering Equation?

(Asked by: username_0 on 2023-06-29)

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?

**Tags:** raytracing, mathematics, integral

> **Comments:**
> **username_1**: 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.  
> **username_0**: @username_1 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?  
> **username_1**: 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.  
> **username_0**: @username_1 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.  

---

## 2 answers

### Answer by username_1 on 2023-06-30. Score: 1 | Accepted answer

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 by username_2 on 2023-06-29. Score: 1

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](https://pbr-book.org/3ed-2018/Monte_Carlo_Integration/Transforming_between_Distributions) for a more detailed explanation.

> **Comments:**
> **username_0**: 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.  
> **username_2**: 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?  
> **username_1**: I also don't like username_0'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).  

---

