【转】sphere ambient occlusion
From：http://www.iquilezles.org/www/articles/sphereao/sphereao.htm
Intro
Beginning of 2006 I got a new graphics card at my workplace. 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 occlusion to all those moving objects. For that I had to find an analytic way to compute the occlusion produced by planes and spheres in each other, since doing a montercarlo approximation was out the question. The occlusion produced by an infinite plane into an sphere is trivial: take the normal to the sphere's surface point, and dot it with the plane normal, add one and divide by two, et voila. Now, the sphere to plane ambient occlusion 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 occlusion of a sphere to a point on an arbitrary surface. As we will see the result is surprisingly simple.
We start by assuming without any lose of generality that we have a surface point with normal pointing up in the y direction. We recall that the ambient occlusion 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 "irradiance". 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 unit radius is nice because we can exchange the concept of area and solid angle (the area of an object project on a sphere is just its 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 a shadow into the pink hemisphere through a yellow cone. 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 (differential) 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 surface 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 less. 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 its center to the point of interest. We can quickly find how R and D relate to r and d by applying Pythagoras theorem twice after noting that the radius of the disk line intersects the sphere position line in ninety 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 its border, and another one over the complete circumference. 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 cosine of beta, 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 |
We have all our ingredients now. The resulting integral is thus:
By substituting back r and d in R and D we get our final blocking factor for the sphere:
 |
or |  |
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 occlusion factor. Basically it means that a sphere x times bigger will project the same shadow if it is x times further away.
However this result is not physically correct, because we know that light coming from the sides will have less influence in out point that light coming directly from above. So, next step is to fix this. For that we have to include the cosine term into our integral. This will have the effect that if the disk is straight above the point, it will block more light even if the projected area is the same as if it was anywhere else. So now we need to know the position of the disk relative to the normal of the surface point.
 |
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 surface 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 realtime image below (drag with the mouse) 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 occlusions, that are very simple to compute). The nice things are: 1] it's done per pixel, not per vertex, so no tessellation 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 recive the AO of the sphere ;)
