Demystifying the Cosine Weighted Hemisphere Algorithm

So… ray tracing seems to be all the rage lately, and with that there have been a lot of talks on Probably Density Functions (PDF) and choosing good sampling rays rather than just casting rays in random directions. In fact, this isn’t the first time we’ve talked about this – there was a great paper published by Sony back in 2003 which covered these very topics, the goal was Spherical Harmonics, but the ray tracing concepts are the same. I just noticed that the author has a somewhat recent blog post which still links to the article, you can find it here.

If you’re reading this, you’ve probably heard the terms Hemisphere Sampling and Cosine Weighted Hemisphere Sampling, but if you’ve ever looked at the equation you might just take it as is with only a vague understand of what it’s doing. I started diving into it and surprisingly enough I couldn’t find a good description of the algorithm – I found a lot of reference which use it and explain the results, but not how the math actually works. I decided to step through and understand it and I figured I’d write about it here in case it helps someone else (or just to remind myself in the future). I’m not going to dive into the theory or more intense mathematical equations, if that interests you then you’ve probably already read them, but I’ll add them to the references at the bottom just in case.

To get a thorough understanding, let’s back up a step and try to understand what we want to achieve…

Imagine you have a ground plane and you want to cast rays from it to see how much sun is hitting it and how much it’s occluded from the sun. Obviously each ray cast is expensive so you have a limited amount you can cast, the first approach would be to randomly cast in all directions…

Such an equation would look something like this (I’m going to start in 2D, but I’ll move to 3D after the basic concepts are explained):

//assume rand.x is a random number from 0 to 2PI...
x = cos(rand.x);
y = sin(rand.x);

What does this do? Well, the full scope of cosine and sine are beyond this tutorial – but a quick summary:

Depending on the angle you give cosine (0 to 2PI) it will return the X coordinate of a line with the length of 1 which shoots out the direction of that angle (that is, where the line would intersect the ‘unit circle’; a circle with a radius of 1). Sine is similar to cosine, except it will return the Y coordinate of the same line.

Wikipedia has some nice visuals to help with this.

You can see from the Wikipedia article that when the angle is > PI the Y position [y = sin(rand.x)] is less than 0 (below the ground in our case), obviously these are wasted rays, which brings in the concept of hemisphere sampling. Hemisphere sampling is the concept of only shooting rays towards the hemisphere. In our 2D case we can accomplish this by clamping the angle from 0 to PI.

Now all of our rays at least have a chance of hitting the sun.

This is great for 2D land, but if we move to 3D, how do we generate our z (forward) ray, and what should it be? The approach I take is to flip our unit circle on its side and let cosine and sine determine the x and z coordinates the ray should shoot through. I’m not worried about up or down yet, so I want to use the full range of angles, 0 to 2PI…

//assume rand.x is a random number from 0 to 2PI...
x = cos(rand.x);
z = sin(rand.x);

For the image below, we’re in the sky looking straight down at the ground, much like if we looking straight down at a compass. Here, as the angle goes from 0 to 2PI, the rays will rotate similar to a compass needle:

Ok that helped our 3D case in terms of the right and forward vector, but what about up? Well as we talked about before, we don’t want our angle to point down and we know using sine from 0 to PI will keep us from pointing down, so we could introduce a second random number:

//assume rand.x is a random number from 0 to 2PI...
//assume rand.y is a random number from 0 to PI...
x = cos(rand.x);
y = sin(rand.y);
z = sin(rand.x);

So.. this would almost work, however the result will create an x, y, z which are not normalized – and our ray direction needs to be normalized. The x and z will be normalized with respect to each other because they use cosine, sine based off of the same number rand.x (which, remember cosine and sine will give the x,y intercept location of a unit circle so the length is guaranteed to be 1). However, z will be it’s own value from 0 to 1. This means we’d need to add this additional line to normalize it:

v = normalize(x, y, z);

This would further increase the cost (we’re already paying the cost of an additional sin call). But, with that said, let’s look at what the result would look like:

I want to show this in 3D so I used my crappy debug graphics, please ignore the horrid leafs on the arrows

Another non-ideal issue is that rand.x and rand.y have different ranges… it’s probably not best practice to assume whoever uses our methods will respect the correct ranges for each value. Let’s fix that part first…

A common practice is to send in two values, both in the range from 0 to 1, this changes the equation to something like this:

//assume rand.x is a random number from 0 to 1...
//assume rand.y is a random number from 0 to 1...

xt = rand.x * 2 * PI; //expand to 0 to 2PI
yt = rand.y * PI; //expand 0 to PI

x = cos(xt);
y = sin(yt);
z = sin(xt);

v = normalize(x, y, z);

This fixes the inconsistent input argument ranges, but it still requires a second call to sine and a normalization at the end which is troublesome (as described above). Another option is to simply use the y value passed in (remember, it’s already 0 to 1 which is what we want) and then shorten the x, z with respect to y so the values are normalized, this would remove the need for the second sine call along with the normalize call (this is based on Hammersley Points which is linked to in the references):

//assume rand.x is a random number from 0 to 1...
//assume rand.y is a random number from 0 to 1...

xt = rand.x * 2 * PI; //expand to 0 to 2PI
yt = sqrt(1.0 - rand.y * rand.y);

x = cos(xt) * yt;
y = rand.y;
z = sin(xt) * yt;

What the heck is going on here??? It’s not as strange as it might seem, remember for y we want a value from 0 to 1, originally we were calling sine for that. However, the random number passed in is already 0 to 1 so we simply choose to use that, thus saving a sine call. That’s what the y = rand.y is doing.

But what about yt? yt is being used to adjust (i.e. shorten) the x and z length of the ray depending on how long the y ray should be. Why do we need to do this? Let’s back up a step…

Remember our 2D ray at the beginning of this page?

x = cos(rand.x); y = sin(rand.x);

Remember those x and y are points on the unit circle, so taken together their distance from the origin is 1 unit. What’s the equation to calculate the distance? It’s based on the Pythagorean theorem. The equation is x^2 + y^2 = c^2, so for our 2D ray case (which we know equals 1) it’s…

x^2 + y^2 = 1^2; //which is the same as...
x^2 + y^2 = 1;

With our 3D ray we’re adding another dimension, but we still need the length to equal 1. With the raw result of x = cos, z = sin this would be impossible because their length equals 1 with no room for another variable. This means we need to modify their result by some number such that adding in the new y^2 will take us to a distance of 1. Let’s add a placeholder coefficient, Enter theta:

(x^2 + z^2) * theta + y^2 = 1;

Theta is going to adjust our x^2 + z^2 result so that adding y^2 will give the result of 1. So how do we determine what theta should be?

First let’s move theta in to directly multiply against x^2 and z^2 because, well, it’s easier for my brain…

(x^2)(theta) + (z^2)(theta) + y^2 = 1

Now apply some sqrt properties; we’re taking the square root of each value before multiplying them and then we’re squaring the result, this is an equivalent algorithm of above but a different way of looking at it…

(sqrt(x^2)sqrt(theta))^2 + (sqrt(z^2)sqrt(theta))^2 + y^2 = 1

Since the square root of x^2 is x and z^2 is z we can remove two square root calls…

((x)sqrt(theta))^2 + ((z)sqrt(theta))^2 + y^2 = 1

So now we see that if we multiply the x and z values by the square root of theta it will modify them such that including the additional dimension y will give us a length of 1.

But how do we find the sqrt(theta) in this slightly complicated equation?…
Well if we step back to the beginning we remember that

x^2 + z^2 = 1

Our initial algorithm…

((x^2) + (z^2))(theta) + y^2 = 1

can become…

(1)(theta) + y^2 = 1

which means…

theta + y^2 = 1

or…

theta = 1 - y^2

so going back to needing to multiply x/z by the sqrt of theta…

yt = sqrt(theta)

which is…

yt = sqrt(1 - y^2)

or, as our original code snippet puts it…

yt = sqrt(1 - y * y)

Neat huh???

The result will be a 3D direction vector (length of 1) which gives us 3D uniform hemisphere sampling that can be quite useful when casting rays:

However, we’re not done yet! While none of these rays will shoot down, many will still shoot off in adjacent angles which aren’t completely useful, especially when wanting to sample something in a certain direction like the sun. This is where Cosine Weighted Hemisphere Sampling comes into play.

Cosine Weighted Hemisphere Sampling biases the random rays towards the up vector which means more random rays in a direction we care about.

It takes random ray generation from something like this:

I increased the ray amount to help visualize the differences between the two images. The red arrow in denotes the up vector.

to this:

You can see how the rays are more biased to an up (but random) direction rather than uniformly scattered around the hemisphere. The genius of this equation is that it’s quite a simple modification on what we’ve already figured out:

//assume rand.x is a random number from 0 to 1...
//assume rand.y is a random number from 0 to 1...

xt = rand.x * 2 * PI; //expand to 0 to 2PI
yt = sqrt(1.0 - rand.y);

x = cos(xt) * yt;
y = sqrt(rand.y);
z = sin(xt) * yt;

Only two minor changes! First, we’re now using the square root of rand.y for our Y length. Since the value passed in, rand.y, is between 0 and 1 using its square root will increase the value rather than decrease it, thus biasing the ray more towards the up vector. The second change is to adjust the x/z length accordingly, since we’re taking the square root of rand.y for our y length we can remove a rand.y multiplication when we’re calculating yt.

Why can we do that? Well it’s a mathematical simplification, remember our uniform hemispherical equation had:

yt = sqrt(1 - rand.y * rand.y);

and…

y = rand.y;

so it could have been written as…

yt = sqrt(1 - y * y)

and remember this time y is defined as…

y = sqrt(rand.y);

so essentially yt is…

yt = sqrt(1 - sqrt(rand.y) * sqrt(rand.y));

which is the same as…

yt = sqrt(1 - rand.y);

which is what we used!

The last question you might be asking yourself is, well this is great if the light source is straight up; but what if it’s at an angle, like the sun in the morning? Fortunately this is a fairly trivial solution, simply take the resulting vector ray(x, y, z) and rotate it about the direction of your light source. For instance, if you have a rotation matrix defining your light source ‘light_rot_matrix’ you would just multiply:

Assuming row major…

rotated_ray = mul(ray, light_rot_matrix);

Assuming col major

rotated_ray = mul(light_rot_matrix, ray);

You can see the results here where I rotated my light source off to the side:

I hope that was helpful in demystifying how the Uniform Hemisphere Sampling and Cosine Weighted Hemisphere Sampling algorithms work. I’m not a math expert so please feel free to ping me if you notice anything which needs to be corrected!

References


Images