<p><strong>Probability Density Function</strong></p>
<p>Say you want to uniformly generate directions only within <span class="math-container">$\theta \in [0,\delta)$</span>. We know that the probability density w.r.t. the solid angle measure must be constant and integrate to one:
<span class="math-container">$$1 = \int_{\Omega_{\delta}} p(\omega)\,d\omega = \int_{\Omega_{\delta}} C\,d\omega =\int_0^{2\pi} \int_0^{\delta} C \sin\theta\,d\theta\,d\phi = 2\pi (1-\cos\delta) C\implies C = \frac{1}{2\pi (1-\cos\delta)}.$$</span></p>
<p>Using <span class="math-container">$\delta = \arcsin(R/d)$</span> and
<span class="math-container">$$\cos\arcsin x= \sqrt{1-(\sin\arcsin x)^2} = \sqrt{1-x^2},$$</span>
you get <span class="math-container">$\cos\arcsin \delta = \sqrt{1-R^2/d^2}$</span>. Then your probability density function is:
<span class="math-container">$$p(\omega) = \frac{1}{2\pi \left(1-\frac{\sqrt{d^2-R^2}}{d}\right)}.$$</span>
There is no <span class="math-container">$2\pi$</span> divided by it. For some reason you decided that the pdf is the ratio of the areas, but as you can see that is not the case. It is the reciprocal of the area in <span class="math-container">$[0,2\pi)\times [0,\delta)$</span>. Note that in <span class="math-container">$2\pi (1-\cos\delta)$</span> if you set <span class="math-container">$\delta = \pi/2$</span> you will recover the area of the whole hemisphere.</p>
<p><strong>Non-constant BRDF</strong></p>
<p>Let's take some strange diffuse brdf that is easy to integrate: <span class="math-container">$f(\theta) = \frac{\alpha+2}{2\pi}\cos^{\alpha}\theta$</span> (for simplicitly <span class="math-container">$\alpha&gt;0$</span>). You can now compute the integral explicitly:
<span class="math-container">\begin{align}
L_r(x,\omega_o) &amp;= \int_{\Omega} f(\theta) L_i(x,\theta) \cos\theta_i \,d\omega_i \\
&amp;= L_e\int_0^{2\pi} \int_0^{\delta} \frac{\alpha+2}{2\pi}\cos^{\alpha+1}\theta_i \sin\theta_i \,d\theta_i \,d\phi_i \\
&amp;= -(\alpha+2)L_e \int_0^{\delta}\cos^{\alpha+1}\theta_i \,d\cos\theta_i \\
&amp;= -(\alpha+2)L_e \left[\frac{\cos^{\alpha+2}\theta_i}{\alpha+2}\right]^{\delta}_0 = L_e (\cos^{\alpha+2} 0 - \cos^{\alpha+2}\delta) \\
&amp;= L_e(1-(1-R^2/d^2)^{\frac{\alpha}{2}+1})
\end{align}</span></p>
<p><strong>Numerical Integration</strong></p>
<p>Numerically you can be as general as you like:
<span class="math-container">$$L_r(x,\omega_o) = \int_{\Omega} f(\omega_o, x, \omega_i) L_i(x,\omega_i) \cos\theta_i\,d\omega_i \approx \frac{1}{N}\sum_{k=1}^N \frac{f(\omega_o, x, \omega_i^k) L_i(x,\omega^k_i)\cos\theta^k_i}{p(\omega_i^k)}.$$</span>
Here <span class="math-container">$\omega^k_i$</span> is the <span class="math-container">$k$</span>-th sampled direction according to the pdf <span class="math-container">$p$</span>, and <span class="math-container">$\cos\theta^k_i = N_x \cdot \omega_i^k$</span>.</p>

CC BY-SA 4.0
2024-02-17T09:07:50.143Z
2024-02-17T22:04:12.930Z
2024-02-17T22:04:12.930Z
lightxbulb
<p>I'm trying to make intuitive sense of the rendering equation and so I'm trying to work out a trivial example by hand, though ultimately I'm just trying to understand how monte carlo importance sampling is applied.</p>
<p>Lets ignore for emittance, and so the rendering equation is simply:</p>
<p><span class="math-container">$$
L_o(x, \omega_o) = \int_{\Omega} f(x,\omega_o,\omega_i)L_i(x,\omega_i) \left(\omega_i \cdot \hat{n}\right) \ d\omega_i
$$</span></p>
<p>Now for my purposes, I'm considering a spherical light with a radius <span class="math-container">$R$</span> and a radiant power of <span class="math-container">$\Phi$</span>, that is <span class="math-container">$d$</span> away from point <span class="math-container">$x$</span> directly along the surface normal vector <span class="math-container">$\hat{n}$</span>.  I'm also assuming the BRDF is lambertian:</p>
<p><span class="math-container">$$
f_r = \frac{\rho}{\pi}
$$</span></p>
<p>where <span class="math-container">$\rho$</span> is the surface's albedo.  In this case, the integral simply becomes:</p>
<p><span class="math-container">$$
L_r = \frac{\rho}{\pi}\int_{\Omega}L_i(x,\omega_i) \left(\omega_i \cdot \hat{n}\right) \ d\omega_i
$$</span></p>
<p>where <span class="math-container">$L_r$</span> is the total reflected radiance.  If I assume that the radiant power of the light is constant across its entire surface (and thus across the entire solid angle it subtends) I can pull that out of the integral as well, leaving me with:</p>
<p><span class="math-container">$$
L_r = \frac{\rho}{\pi} L_i \int_{\Omega} \left(\omega_i \cdot \hat{n}\right) \ d\omega_i
$$</span></p>
<p>Given that the light is located directly above the point, along <span class="math-container">$\hat{n}$</span>, then all of the light incident on <span class="math-container">$x$</span> is within <span class="math-container">$\delta$</span> of <span class="math-container">$\hat{n}$</span>, where <span class="math-container">$\delta$</span> is the half the angular diameter of the light.  The angular diameter is defined as:</p>
<p><span class="math-container">$$
\delta = \sin^{-1}\left(\frac{R}{d}\right)
$$</span></p>
<p>Therefore, our integral can be rewritten as:
<span class="math-container">$$
L_r = \frac{\rho}{\pi} L_i \int_0^{2\pi} \int_0^{\delta} \cos{\theta_i}\sin{\theta_i} \ d\theta d\phi
$$</span></p>
<p><span class="math-container">$$
= \frac{\rho}{\pi} L_i \int_0^{2\pi} \left[\frac{-\cos^2{\theta}}{2}\right]_0^\delta  \ d\phi
$$</span></p>
<p><span class="math-container">$$
= \frac{\rho}{\pi} L_i \left(\pi \sin^2\left(\sin^{-1}\frac{R}{d}\right)\right)
$$</span></p>
<p><span class="math-container">$$
\boxed{L_r = \rho L_i \frac{R^2}{d^2}}
$$</span></p>
<p>This looks as though it could be right, but I'm not sure.  Now lets say I want to setup a monte carlo estimator:</p>
<p><span class="math-container">$$
F = \int f(x) dx \Rightarrow F \approx \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{p(x_i)}
$$</span></p>
<p>For our particular problem here, that would look like this:</p>
<p><span class="math-container">$$
L_r \approx \frac{\rho}{\pi} \frac{1}{N} \sum_{i=1}^N \frac{L_i(\omega_i)\cos{\theta_i}}{p(\omega_i)} 
$$</span></p>
<p>If I uniformly sample <span class="math-container">$\omega_i$</span> from the hemisphere, then <span class="math-container">$p(\omega_i) = \frac{1}{2\pi}$</span> and so the estimator simplifies to:</p>
<p><span class="math-container">$$
L_r \approx 2\rho \frac{1}{N} \sum_{i=1}^N L_i(\omega_i)\cos{\theta_i} 
$$</span></p>
<p>If instead, I were to use importance sampling and directly sample the light, I would need a more complex <span class="math-container">$p(\omega_i)$</span>.  My intuition tells me it should be 1 over the solid angle subtended by the light divided by <span class="math-container">$2\pi$</span> (the solid angle of the entire hemisphere):</p>
<p><span class="math-container">$$
p(\omega_i) = \frac{1}{\Omega_{light} / 2\pi} = \frac{2\pi}{\Omega_{light}}
$$</span></p>
<p>where <span class="math-container">$\Omega_{light}$</span> can be calculated using:</p>
<p><span class="math-container">$$
\Omega_{light} = 2\pi\left(1 - \frac{\sqrt{d^2 - R^2}}{d}\right)
$$</span></p>
<p>My rationale being, if every direction in the hemisphere is equally likely, then the probability of a ray hitting the light should be equal to the proportion of the hemisphere the light occupies.</p>
<p>Introducing this to our monte carlo estimator gives us an importance sampled version written as:</p>
<p><span class="math-container">$$
L_r \approx \frac{\rho}{\pi} \frac{\Omega_{light}}{2\pi} \frac{1}{N} \sum_{i=1}^N L_i(\omega_i)\cos{\theta_i} 
$$</span></p>
<p>Coding this up quickly, if I assume <span class="math-container">$\Phi = 100W$</span>, <span class="math-container">$d = 100m$</span>, <span class="math-container">$R = 10m$</span>, <span class="math-container">$\rho = 0.4$</span>, my analytical solution yields:</p>
<p><span class="math-container">$$
L = 1.0132\times10^{-4} \frac{W}{sr\cdot m^2}
$$</span></p>
<p>If I use the uniformly sampled monte carlo method with a million samples, I get:</p>
<p><span class="math-container">$$
L_{uniform} = 1.0141\times10^{-4} \frac{W}{sr\cdot m^2}
$$</span></p>
<p>And if I use the importance sampled version, I get:</p>
<p><span class="math-container">$$
L_{importance} = 1.6126\times10^{-5} \frac{W}{sr\cdot m^2}
$$</span></p>
<p>So <span class="math-container">$L_{uniform}$</span> matches the analytical solution to within only .09%.  So that gives me some confidence in both solutions.  However the estimate from importance sampling is off by a clean factor of <span class="math-container">$2\pi$</span>:</p>
<p><span class="math-container">$$
\frac{1.0132\times10^{-4}}{1.6126\times10^{-5}} = 6.2832 \approx 2\pi
$$</span></p>
<p>I'm not entirely sure where I went wrong with my methods.  This discrepancy, if not an issue in my derivation/integration, would imply that the probability is actually:</p>
<p><span class="math-container">$$
p(\omega_i) = \frac{1}{\Omega}
$$</span></p>
<p>But I'm not sure how that makes sense.</p>
<p>Further, I'm confused of what happens when the BRDF is not constant.  How is it possible to properly integrate this (either analytically or numerically) when your integral also depends on a BRDF that depends on <span class="math-container">$\omega_o$</span> and <span class="math-container">$\omega_i$</span> as well?</p>

<raytracing><pathtracing><sampling><importance-sampling><monte-carlo>
Evaluating the Rendering Equation by hand (simple case)
Chris Gnam
# Evaluating the Rendering Equation by hand (simple case)

(Asked by: username_0 on 2024-02-17)

I'm trying to make intuitive sense of the rendering equation and so I'm trying to work out a trivial example by hand, though ultimately I'm just trying to understand how monte carlo importance sampling is applied.

Lets ignore for emittance, and so the rendering equation is simply:

$$
L\_o(x, \omega\_o) = \int\_{\Omega} f(x,\omega\_o,\omega\_i)L\_i(x,\omega\_i) \left(\omega\_i \cdot \hat{n}\right) \ d\omega\_i
$$

Now for my purposes, I'm considering a spherical light with a radius $R$ and a radiant power of $\Phi$, that is $d$ away from point $x$ directly along the surface normal vector $\hat{n}$. I'm also assuming the BRDF is lambertian:

$$
f\_r = \frac{\rho}{\pi}
$$

where $\rho$ is the surface's albedo. In this case, the integral simply becomes:

$$
L\_r = \frac{\rho}{\pi}\int\_{\Omega}L\_i(x,\omega\_i) \left(\omega\_i \cdot \hat{n}\right) \ d\omega\_i
$$

where $L\_r$ is the total reflected radiance. If I assume that the radiant power of the light is constant across its entire surface (and thus across the entire solid angle it subtends) I can pull that out of the integral as well, leaving me with:

$$
L\_r = \frac{\rho}{\pi} L\_i \int\_{\Omega} \left(\omega\_i \cdot \hat{n}\right) \ d\omega\_i
$$

Given that the light is located directly above the point, along $\hat{n}$, then all of the light incident on $x$ is within $\delta$ of $\hat{n}$, where $\delta$ is the half the angular diameter of the light. The angular diameter is defined as:

$$
\delta = \sin^{-1}\left(\frac{R}{d}\right)
$$

Therefore, our integral can be rewritten as:
$$
L\_r = \frac{\rho}{\pi} L\_i \int\_0^{2\pi} \int\_0^{\delta} \cos{\theta\_i}\sin{\theta\_i} \ d\theta d\phi
$$

$$
= \frac{\rho}{\pi} L\_i \int\_0^{2\pi} \left[\frac{-\cos^2{\theta}}{2}\right]\_0^\delta \ d\phi
$$

$$
= \frac{\rho}{\pi} L\_i \left(\pi \sin^2\left(\sin^{-1}\frac{R}{d}\right)\right)
$$

$$
\boxed{L\_r = \rho L\_i \frac{R^2}{d^2}}
$$

This looks as though it could be right, but I'm not sure. Now lets say I want to setup a monte carlo estimator:

$$
F = \int f(x) dx \Rightarrow F \approx \frac{1}{N} \sum\_{i=1}^N \frac{f(x\_i)}{p(x\_i)}
$$

For our particular problem here, that would look like this:

$$
L\_r \approx \frac{\rho}{\pi} \frac{1}{N} \sum\_{i=1}^N \frac{L\_i(\omega\_i)\cos{\theta\_i}}{p(\omega\_i)}
$$

If I uniformly sample $\omega\_i$ from the hemisphere, then $p(\omega\_i) = \frac{1}{2\pi}$ and so the estimator simplifies to:

$$
L\_r \approx 2\rho \frac{1}{N} \sum\_{i=1}^N L\_i(\omega\_i)\cos{\theta\_i}
$$

If instead, I were to use importance sampling and directly sample the light, I would need a more complex $p(\omega\_i)$. My intuition tells me it should be 1 over the solid angle subtended by the light divided by $2\pi$ (the solid angle of the entire hemisphere):

$$
p(\omega\_i) = \frac{1}{\Omega\_{light} / 2\pi} = \frac{2\pi}{\Omega\_{light}}
$$

where $\Omega\_{light}$ can be calculated using:

$$
\Omega\_{light} = 2\pi\left(1 - \frac{\sqrt{d^2 - R^2}}{d}\right)
$$

My rationale being, if every direction in the hemisphere is equally likely, then the probability of a ray hitting the light should be equal to the proportion of the hemisphere the light occupies.

Introducing this to our monte carlo estimator gives us an importance sampled version written as:

$$
L\_r \approx \frac{\rho}{\pi} \frac{\Omega\_{light}}{2\pi} \frac{1}{N} \sum\_{i=1}^N L\_i(\omega\_i)\cos{\theta\_i}
$$

Coding this up quickly, if I assume $\Phi = 100W$, $d = 100m$, $R = 10m$, $\rho = 0.4$, my analytical solution yields:

$$
L = 1.0132\times10^{-4} \frac{W}{sr\cdot m^2}
$$

If I use the uniformly sampled monte carlo method with a million samples, I get:

$$
L\_{uniform} = 1.0141\times10^{-4} \frac{W}{sr\cdot m^2}
$$

And if I use the importance sampled version, I get:

$$
L\_{importance} = 1.6126\times10^{-5} \frac{W}{sr\cdot m^2}
$$

So $L\_{uniform}$ matches the analytical solution to within only .09%. So that gives me some confidence in both solutions. However the estimate from importance sampling is off by a clean factor of $2\pi$:

$$
\frac{1.0132\times10^{-4}}{1.6126\times10^{-5}} = 6.2832 \approx 2\pi
$$

I'm not entirely sure where I went wrong with my methods. This discrepancy, if not an issue in my derivation/integration, would imply that the probability is actually:

$$
p(\omega\_i) = \frac{1}{\Omega}
$$

But I'm not sure how that makes sense.

Further, I'm confused of what happens when the BRDF is not constant. How is it possible to properly integrate this (either analytically or numerically) when your integral also depends on a BRDF that depends on $\omega\_o$ and $\omega\_i$ as well?

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


---

## 1 answers

### Answer by username_1 on 2024-02-17. Score: 1 | Accepted answer

**Probability Density Function**

Say you want to uniformly generate directions only within $\theta \in [0,\delta)$. We know that the probability density w.r.t. the solid angle measure must be constant and integrate to one:
$$1 = \int\_{\Omega\_{\delta}} p(\omega)\,d\omega = \int\_{\Omega\_{\delta}} C\,d\omega =\int\_0^{2\pi} \int\_0^{\delta} C \sin\theta\,d\theta\,d\phi = 2\pi (1-\cos\delta) C\implies C = \frac{1}{2\pi (1-\cos\delta)}.$$

Using $\delta = \arcsin(R/d)$ and
$$\cos\arcsin x= \sqrt{1-(\sin\arcsin x)^2} = \sqrt{1-x^2},$$
you get $\cos\arcsin \delta = \sqrt{1-R^2/d^2}$. Then your probability density function is:
$$p(\omega) = \frac{1}{2\pi \left(1-\frac{\sqrt{d^2-R^2}}{d}\right)}.$$
There is no $2\pi$ divided by it. For some reason you decided that the pdf is the ratio of the areas, but as you can see that is not the case. It is the reciprocal of the area in $[0,2\pi)\times [0,\delta)$. Note that in $2\pi (1-\cos\delta)$ if you set $\delta = \pi/2$ you will recover the area of the whole hemisphere.

**Non-constant BRDF**

Let's take some strange diffuse brdf that is easy to integrate: $f(\theta) = \frac{\alpha+2}{2\pi}\cos^{\alpha}\theta$ (for simplicitly $\alpha>0$). You can now compute the integral explicitly:
\begin{align}
L\_r(x,\omega\_o) &= \int\_{\Omega} f(\theta) L\_i(x,\theta) \cos\theta\_i \,d\omega\_i \\
&= L\_e\int\_0^{2\pi} \int\_0^{\delta} \frac{\alpha+2}{2\pi}\cos^{\alpha+1}\theta\_i \sin\theta\_i \,d\theta\_i \,d\phi\_i \\
&= -(\alpha+2)L\_e \int\_0^{\delta}\cos^{\alpha+1}\theta\_i \,d\cos\theta\_i \\
&= -(\alpha+2)L\_e \left[\frac{\cos^{\alpha+2}\theta\_i}{\alpha+2}\right]^{\delta}\_0 = L\_e (\cos^{\alpha+2} 0 - \cos^{\alpha+2}\delta) \\
&= L\_e(1-(1-R^2/d^2)^{\frac{\alpha}{2}+1})
\end{align}

**Numerical Integration**

Numerically you can be as general as you like:
$$L\_r(x,\omega\_o) = \int\_{\Omega} f(\omega\_o, x, \omega\_i) L\_i(x,\omega\_i) \cos\theta\_i\,d\omega\_i \approx \frac{1}{N}\sum\_{k=1}^N \frac{f(\omega\_o, x, \omega\_i^k) L\_i(x,\omega^k\_i)\cos\theta^k\_i}{p(\omega\_i^k)}.$$
Here $\omega^k\_i$ is the $k$-th sampled direction according to the pdf $p$, and $\cos\theta^k\_i = N\_x \cdot \omega\_i^k$.


---

