Variance Shadow Mapping - theory and implementation

In my previous post I presented video from my implementation of VSM - Variance Shadow Mapping. Today I will share theory and code behind it. You can find original paper and sources by Donnelly and Lauritzen here.
The technique is very simple yet it produces very nice and quite natural shadows. It can be considered a modification to standard shadow mapping.

In shadow mapping it is normal to store depth value in a one channel floating-point texture like R32F. Then to produce soft shadows PCF (Percentage Close Filtering) can be used on a full size shadow map to compute influence on the scene by sampling shadow map around a given point n x n times (where n is a size of PCF kernel) and averaging the results. Although this technique has potential of producing very nice soft shadows it becomes very slow for larger n's.
Variance shadow mapping in turn requires two channel texture (16- or 32-bit floating point). One of the channels stores depth value while the other stores depth squared. This allows to consider shadow map as not storing depth value but rather depth distribution.
What we want is the average depth and average square depth. To obtain it we can blur the shadow map (using Gaussian blur or even simple box filter) or use filtering on it (eg. anisotropic). Having average values we can now easily derive (E[X]) ^ 2 and E[X^2] from them. With both values calculated we can compute variances as:

Var = E[X^2] - (E[X]) ^ 2.

Then using Chebyshev's inequality we calculate p as:

p = Var / (Var + (depth - E[x]) ^ 2)

Then we use this value as shadowing factor, i.e. how much shadow is visible.

Here's the code doing the last part:

...
float2 moments = tex2D(shadowMap, texCoord.xy).xy;
float mean = moments.x;
float meanSqr = moments.y;


float Ex_2 = mean * mean;
float E_x2 = meanSqr;
float variance = min(max(E_x2 - Ex_2, 0.0f) + shadowBias, 1.0f);
float m_d = (depth - mean);
float p = variance / (variance + m_d * m_d);


return max(p, depth <= mean);

where depth is value we compare against value stored in the shadow map. In my case it is depth stored in the G-Buffer as I use deferred shading.

So summing up there are just a few steps to be done to use Variance Shadow Mapping:

1. Calculate depth and depth * depth and store it in shadow map (2-channel 16- or 32-bit texture).
2. Blur shadow map using a filter of your choice. It is also possible to down-sample shadow map prior to blurring.
3. Apply shadowing to the scene using snippet above.

And voila!