# Evaluating the Rendering Equation by hand (simple case)

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?

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$$.