I spent my morning programming a small ambient occlusion generator for height maps and the quality of the output is surprisingly nice.  In the picture below, 2048 random rays are cast from the center of each pixel out towards the hemisphere.  Each pixel represents the ratio between rays that hit other locations of the height maps and those rays that are clear of all other terrain.

Several methods were used to improve the quality.  Again at all stages, I use sub pixel accuracy, using bi-cubic interpolation in all calculations.  Random rays are generated that maintain a strong random distribution over the entire hemisphere.  This is accomplished by setting a rays x and y components to a random value, between -1.0f and 1.0f.  The z component is set to a random value between 0.0f and 1.0f, to ensure it rises out of the terrain.  The statistical trick here is to discard the ray if it has a magnitude > 1.0f, and try to generate another random ray.  Once a ray is found with a magnitude of < 1.0f it can be normalized and then used.  This method was described in “Texturing and Modeling: A procedural Aproach”.