Beginning of 2006 I got a new graphics card at my work place. It had a pixel shaders 3.0 capable GPU, and I made few experiments on it. One of them was a shader based raytacer with some spheres and planes. I decided it would be great to add ambient occlussion to all those moving objects. For that I had to find an analytic way to compute the occlussion produced by planes and spheres in each other, since doing a montercarlo approximation was out the question. The occlussion produced by an infinite plane into an sphere is trivial: take the normal to the spheres surface point, and dot it with the plane normal, add one and divide by two, et voila. Now, the sphere to plane ambient occlussion is something else than a dot, and it would be very nice to know it since we could use if for triangles also... (hm, this sound interesting right?)
So this article is about analytically calculating the occlussion of a sphere to a point on an arbitrary suface. As we will see the result is surprisingly simple.
We start by assuming without any lose of generality that we have a suface point with normal pointing up in the y direction. We recall that the ambient occlussion means how much light from everywhere in the upper hemisphere centered in the point and with north pole in the normal direction does not arrive to the given point. It's the opposite of the "iradiance". Here we are going to use the blocking version, that is easier to work with. The problem is to find the blocking factor for a sphere of radius r that is at a distance d from the point. For now we don't need the angle between the normal and the position of the sphere since the blocking factor is independant of the angle. Later, we will take into account the cosine term of the lighting integral and we will need this angle. But, for now what we have is:
Now, the blocking factor is the area that the sphere projects into the hemisphere divided by the area of the hemisphere. If the projection (shadow) of the sphere was covering the complete hemisphere we would have both areas equal and thus the blocking factor would be one.
By assuming the hemisphere has radious one, the area of the hemisphere is two times pi. Working with radius one is nice beacause we can exchange the concept of area and solid angle (the areas of an object project on a sphere is just it's solid angle multiplied by the radius squared of the sphere). The solid angle is given in steradians (even if it is a unit-less measure). The image shows the setup of the problem. The blue sphere is projecting thru a yellow cone a shadow into the pink hemisphere. We have to calculate the solid angle w of the sphere (the small circle drawn by the intersection of the cone and the hemisphere).
|The solid angle of an object is calculated by the formula given below. The integral runs over all the surface S of the
object projecting "the shadow" (the sphere in our case). The idea behind the formula is not that complex: we have to take every surface element
(diferential) dA (note that it's a vector with the orientation of the surface and modulo the differential measure of area), and dot it
with the position vector. That way parts of the surface pointing to our point will occlude more than those being perpendicular, what makes much sense.
The contribution of each of these suface elements is divided by the square of the distance r (don't mess with the radious of our sphere)
to account for the projection of the object into the hemisphere: objects far away of the point have a smaller shadow on the hemisphere, so they occlude
Well, don't be afraid, we will very quickly convert this ugly integral into something that will make more sense. For that we are going to make a change on the setup. We realize that the shadow casted by the sphere is the same as the one casted by the disk resulting from the intersection of the sphere and the cone. Basically, we can change the integral to run on that disk, and that will make our life easier.
So, the first thing to do it to get the dimensions of that disk, based on the dimensions of the sphere. Following the next image bellow, we will call R (in uppercase) the radious of the disk, and D the distance from it's center to the point of interest. We can quickly find how R and D relate to r and d by applying Pitagora's theorem twice after noting that the radius of the disk line intersects the sphere position line in ninty degrees. The relations are the ones on the left of the image.
To integrate over the disk we are going to use polar (flat cylindrical) coordinates, that fit very well the geometry of the disk. So, the integral on the surface is going to be done by a double integral, one along a radial axis from the center of the disk until it's border, and another one over the complete circunference. Let's call these two coordinates for the angle and for the distance to the center of the disk. Then, given the following picture, we have to identify the symbols in the previous integral.
So first, r is just the hypotenuse of the triangle with sides and D, so
Next, we have to do the dot product of the numerator of the integral. Because the normal n has modulus one by definition, we have
For the cosinus, we just have to note that the angle beta between dA and n is the same as the angle formed in the triangle D--r. So,
Finally only |dA| remains. It is the differential area created when moving the integration point a small amount in the direction (d) and a small amount in a d rotation (this creates a ·d arc). It is the shaded area in the picture on the left. The result is
what is pretty simple after all, just a few operations per fragment and we are done. Note how the expression is just a function of r/d as a hole. This makes a lot of sense, since making a sphere bigger means we have to move it further to have the same projected shadow or occlussion factor. Basically it means that a sphere x times bigger will project the same shadow if it is x times futher away.
|We have to take care of correctly integrating the disk surface. First note that I called the angle between
the current integration point and the up direction. That is the angle for the lambert cosine term. Then, by reusing what we got in the previous analysis,
we have that:
The tricky part is to put in terms of and (the angle between the suface vector and the center of the sphere/disk direction). By drawing the thing on paper slowly, one arrives to:
so the complete integral is:
|The second integral evaluates to null, because of the||so we get:|
|meaning that for the disk,|
And making the final substitution
So the result is amazingly simple (and intuitive!) Just an inversesqrt() (RSQ), a dot product and few muls! The image below shows an example of the formula being applied in a fragment shader to get the AO produced by the sphere (note that I also included the plane occlussions, that are very simple to compute). The nice things are: 1] it's done per pixel, not per vertex, so no tesselation or clashing polygon artifacts 2] it's fully dynamic, so you can move the sphere. The technique is not restricted to planes of course, any 3D mesh can recieve the AO of the sphere ;)