Book 3.4 / Listing 16: Confusion

[dimitry-ishenko]

Let's start with analytical solution. The paragraph just below listing 16 says:

The analytic answer is $\frac 4 3 \pi = 4.188790204786391$ — if you remember enough advanced calc, check me!

I remember very little advanced calc, but let me try...

NB: I've updated this post, which was initially using wrong formulas to compute surface integral.

General equation to compute an integral of function $f$ over surface $S$ is:

$$I_S = \iint_S f \cdot dS$$

Using spherical coordinates, surface $S$ can be described in terms of a parameterized vector function $\vec v(\theta, \phi)$, and it can be proven (I won't do that here) that:

$$dS = \bigg| \frac {\partial \vec v}{\partial \theta} \times \frac {\partial \vec v}{\partial \phi} \bigg| d\theta d\phi$$

For a sphere of radius $r$, we can write the following parametric equations:

$$\displaylines{ x = r \cdot \sin \theta \cdot \cos \phi \ y = r \cdot \sin \theta \cdot \sin \phi \ z = r \cdot \sin \theta }$$

$$\vec v(\theta, \phi) = (x, y, z) = (r \cdot \sin \theta \cdot \cos \phi, r \cdot \sin \theta \cdot \sin \phi, r \cdot \sin \theta)$$

If you take partial derivatives $\frac {\partial \vec v}{\partial \theta}$ and $\frac {\partial \vec v}{\partial \phi}$, and compute their cross-product, you will arrive at the following equation:

$$dS = r^2 \sin (\theta) d\theta d\phi$$

And, since we have a unit sphere, it becomes: $dS = \sin (\theta) d\theta d\phi$

So now, if we integrate our function $f(\theta, \phi) = \cos^2(\theta)$ over this unit sphere, we get:

$$I_S = \iint_S \cos^2(\theta) \sin(\theta) d\theta d\phi = \int_0^{2\pi} \bigg( \int_0^\pi \cos^2(\theta) \sin(\theta) d\theta \bigg) d\phi$$

Let's integrate over $\theta$ first:

$$I_\theta = \int_0^\pi \cos^2(\theta) \sin(\theta) d\theta$$

Using u-substitution with $u = \cos(\theta)$, we get $\frac {du}{d\theta} = -\sin(\theta)$ or $d\theta = -\frac {du}{\sin(\theta)}$. Substitute $u$ and $d\theta$ into the above formula, and we get:

$$I_\theta = \int_0^\pi -u^2 \sin(\theta) \frac {du}{\sin(\theta)} = \int_0^\pi -u^2 du$$

$$I_\theta = -\frac {u^3} 3 \bigg|_0^\pi = -\frac {\cos^3(\theta)} 3 \bigg|_0^\pi = \frac 2 3$$

Now we integrate over $\phi$:

$$I = \int0^{2\pi} I\theta d\phi = \int_0^{2\pi} \frac 2 3 d\phi = \frac 2 3 \phi \bigg|_0^{2\pi} = \frac 4 3 \pi$$

QED

[dimitry-ishenko]

Now, the confusing part...

As I've mentioned in #1534 we've made a very unfortunate decision of using letter $f$ to designate inverse of the CDF, as well as arbitrary functions that we want to integrate.

In Chapter 4 we are given function $f(\theta, \phi) = \cos^2(\theta)$ and "we want to integrate this function over the surface of the unit sphere".

Then, in Listing 16 we stick this function into:

double f(const vec3& d) {
    auto cosine_squared = d.z()*d.z();
    return cosine_squared;
}

which up until now (see Listings 13, 14 and 15) was used for inverse of the CDF...

Looking closer at the code in Listing 16, I've realized that while it looks almost identical to the code in Listings 13, 14 and 15, in this particular case double f(const vec3&) is the function that we want to integrate—not inverse of the CDF. And, since we have constant PDF, we don't need to map our random samples using inverse CDF. Mystery solved.

[dimitry-ishenko]

If I may suggest, the code in Listings 13, 14 and 15 should be adjusted as follows:

Listing 13:

-double f(double d) {
+double inv_cdf(double d) {
     return 2.0 * d;
 }

 double pdf(double x) {
     return 0.5;
 }
+
+double f(double x) {
+    return x * x;
+}

 int main() {
     ...
     for (int i = 0; i < N; i++) {
-        auto x = f(random_double());
-        sum += x*x / pdf(x);
+        auto x = inv_cdf(random_double());
+        sum += f(x) / pdf(x);
     }
     ...
 }

Listing 14:

-double f(double d) {
+double inv_cdfdouble d) {
     return std::sqrt(4.0 * d);
 }

 double pdf(double x) {
     return x / 2.0;
 }
+
+double f(double x) {
+    return x * x;
+}

 int main() {
     ...
     for (int i = 0; i < N; i++) {
         auto z = random_double();
         ...
-        auto x = f(z);
-        sum += x*x / pdf(x);
+        auto x = inv_cdf(z);
+        sum += f(x) / pdf(x);
     }
     ...
 }

Listing 15:

-double f(double d) {
+double inv_cdf(double d) {
     return 8.0 * std::pow(d, 1.0/3.0);
 }

 double pdf(double x) {
     return (3.0/8.0) * x*x;
 }
+
+double f(double x) {
+    return x * x;
+}

 int main() {
     ...
     for (int i = 0; i < N; i++) {
         auto z = random_double();
         ...
-        auto x = f(z);
-        sum += x*x / pdf(x);
+        auto x = inv_cdf(z);
+        sum += f(x) / pdf(x);
     }
     ...
 }

With the above changes, Listing 16 can actually remain the same.

The text surrounding Listings 13-15 might need to be tweaked as well to match up with the new function names.

[hollasch]

How do you feel about making the opposite change? That is, leave $f(x) = CDF^{-1}$, and instead introduce a different name for the function of directions? Obviously, $d$ is an appealing name, but can also be confused with direction vectors. Perhaps just go with $g(x)$ or some other single-letter function name?

[hollasch]

@dimitry-ishenko ↑

[hollasch]

@dimitry-ishenko — I thought about this a bit more yesterday, and now I'm planning on proceeding with a variant of your suggestion. I'll name the inverse cumulative probability distribution function $ICD(x)$, for Inverse Cumulative Distribution (the "F" is redundant, so I'll drop that). Thus we have $PDF(x)$, $ICD(x)$, pdf(x), and icd(x). Sound good?