|
Warning: This post is going to be pretty math heavy. If you suck at math, then go and read this first, then come back. If you still think that you suck, then I suggest going to Khan Academy and start watching videos. You won’t regret it!
|
I was reading through the AMD CubeMapGen source last week and came across the code for calculating the solid angle of a cube map texel. This code piqued my interest, since it seemed very terse for what I thought would be a horrific calculation.
static float32 AreaElement( float32 x, float32 y )
{
return atan2(x * y, sqrt(x * x + y * y + 1));
}
float32 TexelCoordSolidAngle(int32 a_FaceIdx, float32 a_U, float32 a_V, int32 a_Size)
{
//scale up to [-1, 1] range (inclusive), offset by 0.5 to point to texel center.
float32 U = (2.0f * ((float32)a_U + 0.5f) / (float32)a_Size ) - 1.0f;
float32 V = (2.0f * ((float32)a_V + 0.5f) / (float32)a_Size ) - 1.0f;
float32 InvResolution = 1.0f / a_Size;
// U and V are the -1..1 texture coordinate on the current face.
// Get projected area for this texel
float32 x0 = U - InvResolution;
float32 y0 = V - InvResolution;
float32 x1 = U + InvResolution;
float32 y1 = V + InvResolution;
float32 SolidAngle = AreaElement(x0, y0) - AreaElement(x0, y1) - AreaElement(x1, y0) + AreaElement(x1, y1);
return SolidAngle;
}
The source code for this particular part is well documented, and points you towards this thesis by Manne Öhrström (@manneohrstrom) where he gives a high level overview of the derivation. I was interested in finding out some more of the details, so I had a go myself, and this post is the result.
Why is it Useful?
When processing cube maps (for example, generating a diffuse irradiance map, or spherical harmonic approximation), you need to be able to integrate the texel values over a sphere.
One way of approximating this integral is to use a Monte Carlo estimator. This is a statistical technique that may oversample some texels and undersample other. This seems a bit wasteful considering that we have a finite number of input values. Ideally we’d like to use each texel value just once.
A naive approach to analytical integration where each texel has the same weight would result in overly bright values in the corner areas. This is because the texels in the corners project to smaller and smaller areas on the sphere. The correct approach is to factor in the solid angle during the integral, and this is what CubeMapGen does.
The Plan
Imagine a single cube map face placed at (0,0,1) and scaled such that the texel locations are all in [-1,1]. For any texel in this cube map, we want to project it onto a unit sphere sitting at the origin, then work out the area on the sphere. This area corresponds to the solid angle because the sphere is a unit sphere.
We can repeat this same calculation for any of the other cube map faces by first transforming them into the same range.
This is the high-level game plan for calculating out the solid angle:
- Determine a formula for projecting a position from texture-space onto the sphere.
- Work out how this projected position changes as the texture-space coordinates change in x and y.
- Imagining that these position change vectors define two sides of a microscopic quadrilateral, then calculate the microscopic area of this quad using the magnitude of the cross product.
- Integrate the microscopic area using the corner coordinates of a texel to calculate its area on the sphere, and solid angle.
The Details
We start off with the formula for projecting a point from its location on the texture face (x, y, 1) onto the unit sphere. This is just a standard vector normalization.
Note: I’ll be switching back and forth between negative and fractional exponents as I see fit. This makes things easier. Remember,  , and  .
|
We want to calculate how this projected point changes as the texture-space x and y coordinates change. We can do this separately for each axis using partial derivatives. First we’ll start by calculating how the projected z component changes along to the texture-space x axis.
Projected Z Change According to X
The z-component of p is simply:
We need to differentiate this equation with respect to x. Because of the exponent, we need to use the chain rule to do this. The chain rule is a method for finding the derivative of the composition of two functions. First, we can reformulate the equation a little bit to make the two functions a bit clearer:
In our case our first function is a function of 
and our second function is a function of 
and 
. Given this, the chain rule says:
We can apply this rule very easily to our reformulated functions:
![\begin{align*} \frac{\partial p_z}{\partial u} &=-\frac{u^{-\frac{3}{2}}}{2}\\ &=-\frac{1}{2(x^2+y^2+1)^{\frac{3}{2}}} \\[10] \frac{\partial u}{\partial x}&=2x\\[10] \frac{\partial p_z}{\partial x} &= -\frac{x}{(x^2+y^2+1)^{\frac{3}{2}}} \end{align*}](http://www.rorydriscoll.com/wp-content/ql-cache/quicklatex.com-ae0b2428e50b1171daec7bafca098554_l3.png "Rendered by QuickLaTeX.com")
This equation tells us exactly how the z component of the projected point changes as the texture-space position moves along the x axis.
Projected X Change According to X
Now we’ve found the projected z component derivative, it’s going to make finding the x and y components a little easier. Why? Because we can express the x and y components in terms of the z component.
We don’t have the same ‘composition of functions’ setup that we did last time, so we can’t use the chain rule to differentiate this. Instead, we can use the product rule. The product rule in our case says:
Applying this to equation for for the projected x component:
Projected Y Change According to X
We have a very similar derivation for the projected y derivative:
We can use the product rule again:
Projected Position Change According to X and Y
Putting this all together, we have our equation showing how the projected position changes as the texture-space position changes in the x direction. We can use the exact same process to work out how it moves in the y direction.
Differential Area
The next step is to calculate the differential (microscopic) area of the projected point using the partial derivatives we just calculated. Clearly at a normal scale, we wouldn’t be able to take the cross product of two projected vectors on a sphere and expect the magnitude to be the area they define on the sphere. But at this differential scale, we can treat the surface as if it is flat, so this works.
The first thing we need to do is to calculate the cross product of the partial derivatives.
Calculating the cross product for each of the components relatively straightforward.
If you’re reading carefully, you’ll notice that each component has a factor of 
on the top and the bottom, so we can divide through. Combining all the components back together again, we arrive at the final equation for the perpendicular vector.
Now we simply need to take the length of the result of the cross product to find the differential area on the sphere.
![\begin{align*} \partial A &= \sqrt{\vec{r} \cdot \vec{r}}\\[5] &= \sqrt{\frac{x^2+y^2+1}{(x^2 + y^2+1)^4}}\\[5] &= \frac{1}{(x^2 + y^2+1)^{\frac{3}{2}}} \end{align*}](http://www.rorydriscoll.com/wp-content/ql-cache/quicklatex.com-d01ca9cd76b91c518393ab4a9ad617bf_l3.png "Rendered by QuickLaTeX.com")
Solid Angle
The final step is to integrate the differential area over our range of texture-space values to get the solid angle of the texel. We can start by calculating the integral between 
and some point 
on the cube map face.
From this formula, we can calculate the area of any texel in the cube map face by adding together the two right-diagonal corners, A and C, and subtracting the left-diagonal corners, B and D.
You can see on the image below that the added areas in green are canceled out by the subtracted areas in red.
That should look familiar, since that’s exactly what the CubeMapGen code does. If you look at the surrounding source code to TexelCoordSolidAngle, then you’ll notice that there’s another method mentioned for calculating the solid angle of a texel. This method is based on Girard’s theorem, which describes how to calculate the area of a spherical triangle based on the excess of the sum of its interior angles. This method was also suggested to me on Twitter by Ignacio Castaño (@castano). I haven’t actually tried it, but it looks fascinating!
Is it Correct?
It’s always a bit daunting to get to the end of a derivation like this, and not know if the answer is correct or not. In this case, it’s pretty easy to verify if this result is correct.
Remember that the texture-space coordinates are in range [-1,1]. If we set our 
values to 1, that corresponds to the top right quarter of the cube map face. We know that there are 
steradians in a sphere, so that means that each face gets 
steradians. Since we’re only calculating for a quarter of a face we expect our result to be 
.
And it does. Thanks to the various people on twitter (@SebLagarde, @mattpharr, @manneohrstrom, @castano and @ChristerEricson) for engaging me in conversation over this.
Please let me know in the comments if you spot an error in this post, or if anything needs to be explained better or more easily.
