Horizon Based Ambient Occlusion

Introduction

In modern computer graphics there are many different ways to compute the lighting of a 3D scene. The most basic lighting models, which the reader encounters early in their studies, are the Lambert model, the Phong model, and models derived from them. The problem is that such models are local — computing the illumination of a given point in a 3D scene uses (typically) the position of that point, its normal, and the positions of the light source and the camera. Thus, when computing the illumination of one object, the presence of other objects in the scene is not taken into account at all.

Another way of computing illumination is shading models, in which the illumination of a point is determined by how much "it is occluded" from the light source by other objects in the scene. It is assumed that we have a global light source (for example, the sky), and then the more light that can reach a point, the more illuminated that point is. This way of computing illumination is called ambient occlusion.

The simplest way to implement such a lighting model is ray tracing: from the point for which the illumination is being computed, a bundle of rays is cast in all directions. The more rays that reach the light source, the more illuminated the point is. This method is not applicable to computer games, since it is computationally expensive and cannot be implemented in a way that provides the necessary frame rate.

Computer games use approximations that produce a roughly similar (more or less) result and can be implemented in "real time". One such approximation is SSAO (screen space ambient occlusion) technology. Proposed by Crytek in 2007, this technology uses the scene's depth map: for each point, the depth of surrounding points is estimated, and the resulting illumination is computed as the difference between the depth of the current point (for which the calculation is being performed) and the depth of the surrounding points of the 3D scene. This method of calculation is fairly crude and not physically correct, but it was used in computer games for quite a long time, since it gave (with certain settings) a decent result.

In this article we'll look at another, more correct way of approximating ambient occlusion, proposed by NVidia in 2008 and giving a more correct and pleasing result: horizon-based ambient occlusion (HBAO).

This article partially uses source code and terminology from the previous article on deferred lighting. If the concepts of the g-buffer, depth map, world and view-space coordinates don't mean anything to you, I recommend first reading the article on deferred lighting.

Theory

General case

In the general case, for a point P of a 3D scene the illumination can be computed using the following formula:

A = 1 - {1 \over 2\pi} \int_{\varOmega} V(\vec{\omega}) W(\vec{\omega}) \,d\omega

In this formula, the integration is performed over the unit hemisphere ω, oriented along the normal at point P.

The function V is the visibility function and returns 1 if a ray cast in a given direction from the hemisphere ω out of point P hits an obstacle, and 0 if it does not and reaches the light source.

The function W is a linear attenuation function: the farther the intersection is from point P, the less significant it is.

Horizon-Based Ambient Occlusion

The essence of the HBAO method is to find the maximum viewing angle that can be built from point P without intersecting geometry. The larger this angle, the more illuminated the point is. This is how the function V from the previous formula is approximated. This maximum angle is called the horizon, and the method is accordingly called horizon-based.

Let's look at an example in a one-dimensional illustration:

The algorithm consists of a step-by-step movement from point P, for which the illumination is being computed, along a given direction. At each step, the angle between the vector S_iP (where S_i is the scene point at step i) and the X axis is taken. It's important to note that if at a given step the angle is smaller than the one found previously, this step is skipped.

As can be seen in the illustration, at step S₁ the angle between the vector S₁P and the X axis is smaller than the angle between the vector S₀P and the X axis, so this step is skipped. The illustration shows that such a situation can arise when the scene geometry at step S₀ occludes the geometry at the step with vector S₁.

In actual calculations, the XY screen plane is used instead of the X axis, with all calculations carried out in camera coordinates (view space), where the origin of the coordinates is the camera position. In addition, the normal at point P needs to be accounted for in the calculations:

In the diagram above, a tangent vector (perpendicular to the normal) has appeared. Now the calculations will use not only the angle h (between the horizon vector and the XY plane), but also the angle t (between the tangent and the XY plane). Logically, this can be explained as follows: the larger the angle between the normal and the horizon vector, the more illuminated the point is. If the horizon vector coincides with the normal at the point, then that point is not illuminated at all.

Thus, the formula considered above for computing the general case of ambient occlusion can be rewritten as follows:

A = 1 - {1 \over 2\pi} \int_{\theta=-\pi}^\pi \bigg(sin(h(\theta)) - sin(t(\theta))\bigg) W(\theta) \,d\theta

In the formula above, the integration is performed over θ, i.e. over all directions (from -π to π). For each direction θ, the difference between the sine of the tangent and the sine of the horizon is taken. This difference is multiplied by the attenuation function W.

To simplify the calculation, a practical implementation of the formula above uses the Monte Carlo method: the number of directions is limited, the sum of the integrand over them is computed, and the result is divided by the number of directions.

Implementation

Data preparation

The input data for the HBAO algorithm is a g-buffer containing a color texture (albedo), a packed normal map, and a depth map:

Preparing the g-buffer and passing it to the shader is covered in detail in the article on deferred lighting.

Camera coordinates (view space)

For convenience, all calculations are performed in camera space (the so-called view space). To convert a vector from world space to view space, this vector must be multiplied by the camera matrix. Let's get this matrix in the script and pass it to the shader:

glsl

Matrix4x4 worldToCameraMatrix = camera.worldToCameraMatrix;
deffered_material.SetMatrix("_worldToView", worldToCameraMatrix);

In the shader, let's add a function to convert a coordinate from world space to camera space:

glsl

uniform float4x4 _worldToView;

float3 worldToView(float3 p) {
    float4 result = mul(_worldToView, float4(p, 1.0));
    return result.xyz /result.w;
}

Now, all that remains to be done before computing ambient occlusion itself is to convert the point's coordinate and normal into camera space:

glsl

float3 pointVS = worldToView(screenToWorld(uv));
float3 normalVS = normalize(mul((float3x3)_worldToView, normal).xyz);

Now all the data is prepared and we can move on to computing HBAO itself. Below are images that can be obtained if you output the coordinates and normals in view space. If in your case the color ranges differ significantly from the images below, you probably made a mistake somewhere in the calculations. I recommend repeatedly checking the correctness of converting coordinates from world space to view space — with an incorrect implementation, you can lose a lot of time trying to understand why nothing works, or why it works incorrectly.

A reliable way to check the correctness of the coordinate conversion is to implement Lambertian lighting in view space. As is known, Lambertian lighting is, in essence, the cosine between the normal at a point and the vector from the point to the light source. Accordingly, by converting the coordinate, the normal, and the light source position to view space and then computing the cosine, you can get a correctly computed Lambertian lighting. If the lighting looks correct and doesn't change as the camera moves/rotates (doesn't jump, doesn't drift, etc.), then your conversion to view space is working correctly.

Sampling along a single direction

Let's add variables responsible for configuring the algorithm:

radiusSS — the maximum distance at which the horizon is searched for. Specified as a fractional number, where 1.0 is the entire screen. At values greater than 0.25–0.3 the algorithm won't work sensibly.

stepsCount — the number of steps into which radiusSS is divided, in which the horizon angle is searched for. The larger this value, the more precisely the algorithm works and the higher-quality the result will be... and the more computationally expensive the algorithm will be.

deltaUV = float2(radiusSS / (stepsCount + 1.0), 0.0) — a two-dimensional vector that will be added to the texture coordinates in the loop, thereby providing movement from the source point to the point maximally distant (by radiusSS) from the source point. At each step of the loop, the horizon will be searched for.

Let's look at the source code that finds the horizon and computes the light attenuation occlusion along a single direction:

glsl

// the sensitivity of the algorithm to bumps depends on this value
float horizonAngle = 0.04;

// loop starts from 1, since at zero we'd get the source point
for (int j = 1; j <= stepsCount; j++) {
    // add the offset to the source point's texture coordinates
    float2 sampleUV = uv + j * deltaUV;
    // get the view-space coordinates for the sample point
    // screen-space -> world-space -> view-space
    float3 sampleVS = worldToView(screenToWorld(sampleUV));
    // build the vector from the sample point to the source point (the horizon vector)
    float3 sampleDirVS = sampleVS - pointVS;
    
    // compute the horizon for the sample point: if it's smaller than the one found earlier
    // do nothing, otherwise update the horizon angle
    float angle = (PI / 2.0) - acos(dot(normalVS, normalize(sampleDirVS)));  
    if (angle > horizonAngle) {
        float value = sin(angle) - sin(horizonAngle);
        float attenuation = clamp(1.0 - pow(length(sampleDirVS) / 2.0, 2.0), 0.0, 1.0);
        occlusion += value * attenuation;
        horizonAngle = angle;
    }
}

The following points are worth paying attention to:

If, at a given step of the loop, the horizon found is larger than the one found at the previous step, we add the difference between the sines of these angles to the resulting attenuation, multiplied by the attenuation coefficient attenuation.
The length of sampleDirVS equals the distance from the current sample point to the source point for which the attenuation is being computed, and can be used to construct the attenuation function.
The attenuation function is, as a rule, chosen manually. In this case it's an exponent of sampleDirVS with small coefficients for fine-tuning.
On the internet you can find implementations where the maximum horizon angle is found first, and then occlusion is computed based on it once and multiplied by the attenuation function. This approach is incorrect, since the attenuation function in that case would be computed only once, whereas in our implementation it is computed as many times as the maximum horizon angle is found during the loop pass. Thus, our implementation accounts for the fact that the closer the found horizon angle is to the source point, the more significant its contribution to the resulting light attenuation. This is very important.

Sampling over a sphere

We have one last step left: to move from sampling along a single direction to sampling over a sphere. To do this, let's add a variable directionsCount, denoting the number of directions along which sampling will occur. Let's also add a rotation matrix deltaRotationMatrix, by which the vector deltaUV will be multiplied, thereby rotating it by an angle of 2π/directionsCount:

glsl

const int directionsCount = 64;

float theta = 2.0 * PI / float(directionsCount);
float2x2 deltaRotationMatrix = float2x2(
    cos(theta), -sin(theta),
    sin(theta),  cos(theta)
);

Now the shader will have two loops: the outer one will pick the direction, and the inner, nested one will sample along the chosen direction:

glsl

float occlusion = 0.0;

for (int i = 0; i < directionsCount; i++) {
    float horizonAngle = 0.04;

    // rotate the deltaUV vector by an angle of 2PI/directionsCount each time,
    // thereby moving on to a new direction
    deltaUV = mul(deltaRotationMatrix, deltaUV);

    for (int j = 1; j <= stepsCount; j++) {
        float2 sampleUV = uv + j * deltaUV;
        float3 sampleVS = worldToView(screenToWorld(sampleUV));
        float3 sampleDirVS = sampleVS - pointVS;
        
        float angle = (PI / 2.0) - acos(dot(normalVS, normalize(sampleDirVS)));  
        if (angle > horizonAngle) {
            float value = sin(angle) - sin(horizonAngle);
            float attenuation = clamp(1.0 - pow(length(sampleDirVS) / 2.0, 2.0), 0.0, 1.0);
            occlusion += value * attenuation;
            horizonAngle = angle;
        }
    }
}

All that's left is to divide the resulting occlusion by the number of directions (in accordance with the Monte Carlo integration method), subtract it from 1 (the value corresponding to the brightness of a fully lit point), and fine-tune the result for your scene:

glsl

occlusion = 1.0 - occlusion / directionsCount;
occlusion = clamp(pow(occlusion, 2.7), 0.0, 1.0);

Final shader and result

Below is the full version of the shader implementing the horizon-based ambient occlusion lighting algorithm:

Result:

Implementation details

Parameter tuning

The HBAO algorithm is very sensitive to correct parameter tuning:

radiusSS — the sampling radius, specified as a fraction of the screen: the higher it is, the more distant objects will be taken into account when computing a point's occlusion. Typical values range from 0.05 to 0.15. At values greater than 0.25 the algorithm will work incorrectly.
directionsCount and stepsCount — the number of directions along which sampling is performed and the number of steps in which calculations are done for each direction. The higher these values, the higher-quality the result and the lower the FPS.
initial value of horizonAngle — chosen as 0.04 in this article. In other sources this parameter is sometimes called angleBias. The initial horizon value: the larger this value, the less the algorithm will take small details into account.
attenuation — the attenuation function. Usually chosen empirically and always depends on the distance between the source point and the sample point. Most often this distance is simply raised to a certain power.
Processing of the resulting value. The final occlusion value is also often raised to some empirically chosen power to increase ambient occlusion's contribution to the final result.

With incorrectly chosen parameters, it's very easy to get an unacceptable result even with a correct implementation. I'd like to once again caution readers and recommend checking the correctness of the coordinate and normal transformation from screen space to view space before implementing the HBAO part of the algorithm itself.

Speeding up the algorithm: noise and blur

The HBAO algorithm can be sped up by reducing the directionsCount and stepsCount parameters. Reducing these parameters leads to banding appearing in the final result. To eliminate this negative effect, a random value is added to the rotation angle of deltaUV and to the sample point. This typically eliminates the banding, but introduces noise, which is removed by ordinary blurring of the final lighting result.

In this article, noise and blur were not added so as not to overcomplicate the article, but in a real project these steps are mandatory — otherwise the algorithm will be too expensive.

Next steps

The HBAO algorithm belongs to the class of screen space algorithms: only the visible part of the frame is used to compute the result. This is a very, very serious limitation: if there is geometry outside the camera's view that occludes an object visible in the frame from the light source, the visible object will be rendered without that occlusion. Further study of the topic of ambient occlusion will lead the reader to algorithms that are more complex to implement, but also more correct: LSAO (Line-Sweep Ambient Obscurance) and GTAO (Ground Truth-based Ambient Occlusion).

Links and files

The author of this article bows and expresses sincere gratitude to Bogdan Mazurenko aka Che@ter for help in working through the theory and practical implementation of the HBAO algorithm. There are a great many incorrect implementations on the internet, and only a single document from NVidia describing the theoretical basis underlying the algorithm. Bogdan's expert help not only made it possible to implement the algorithm correctly, but also to understand its finer points, so they could later be shared with you.

Computer graphicsambient-occlusionhbaoshading