Okay, I’ll state this up front: I’m probably not going to use this approach in my own engine because of many issues inherent with Variance Shadow Mapping. However, I think I did end up with some interesting results to play with, so if VSM with fixed penumbra sizes (or just for filtering) is working well for you, the article may still be useful anyway.
Further worth noting is that most soft shadow articles discuss point lights, so I’ve done things with directional lights.
Introduction
Penumbra widening with distance
If you need an introduction to Variance Shadow Maps, I recommend Andrew Lauritzen’s classic article in GPU Gems3: Summed-Area Variance Shadow Maps. It also contains a technique for very nice precise soft shadows. So, yeah, VSM uses probability theory to estimate whether or not a point is in shadow. Groovy! Unlike standard shadow mapping, this allows for texture filtering the same way regular texture sampling does (bilinear/anisotropic sampling, mip-mapping, etc), and you can use anti-aliasing while rendering the shadow maps. What’s more, the shadow maps can be pre-filtered with a separable blur. This way we can eliminate jaggies using a small filter kernel, or create very soft shadows (with a fixed penumbra size) using a larger kernel.
Real shadows, however, do not have a fixed size penumbra size; they get “softer” further away from the occluding object. No matter the size of the penumbra, we will need to filter the shadow map to get a blurred version of the original. The two general approaches are:
- Using mip-maps: sample the mip-levels since they already contain further filtered data. By default, this starts looking very boxy with larger penumbrae. To alleviate this, you’d have to generate the mip levels with more expensive filtering (more or less like blurring the mip-maps as you’re generating them).
- Using Summed-Area Tables as in the GPU Gems article, resulting in very high-quality results. You can generate SATs like this: “Hensley [2005]: Fast Summed-Area Table Generation and its Applications”. Armed with this, any convolution using a SAT just takes 4 texture samples.
Either way, you end up pre-filtering the shadow maps, the cost of which is dependent on the amount and resolution of the shadow maps (cubic shadow maps, different cascades for directional lights, etc…); not something I wanted to spend too much frame time on. So instead of pre-filtering, I wanted to try and combine it with “post-filtering” in the lighting shader in a way similar to PCSS but without the crazy amount of samples. However, a standard mip-map chain still needs to be generated.
Overview
Calculating soft shadows with shadow maps is not exactly physically correct, but it does result in a visually pleasing approximation: penumbrae get wider with distance. Conceptually, we’ll be using the same method described in NVidia’s Percentage-Closer Soft Shadows paper (so review it ;) ) but of course the implementation will be quite different. As a recap, the steps involved are:
- Find the search area where potential occluders could be
- Find the average occluder depth
- Calculate the penumbra size from the average depth
- Test the shadow map to find the percentage (or in our case: probability) of occluders in the penumbra region.
The main building block
The main building block for our approach, unsurprisingly, uses the Chebyshev’s inequality theorem to find an upper bound for the probability that a sample is in the light. This is the default VSM fare:
// reference contains the depth value of the point to be compared
float UpperBoundShadow(float2 moments, float referenceDepth)
{
float variance = moments.y - moments.x * moments.x;
// clamp to some minimum small variance value for numerical stability
variance = max(variance, MIN_VARIANCE);
float diff = referenceDepth - moments.x;
// Chebyshev's inequality theorem
float upperBound = variance / (variance + diff*diff);
// The upper bound is only correct when referenceDepth < moments.x (if not, return 1.0, ie: fully lit)
return max(upperBound, referenceDepth < moments.x);
}
Finding the occluder search area
This one is exactly the same as for PCSS with an exception for directional lights. If actual directional lights would exist, there would be no penumbra. After all, all light rays are parallel! Also, the traditional PCSS way of back-projecting makes little sense either because the light doesn’t have an actual position in space. To get some handle on it, we’ll settle for a fixed search area instead.
Find average occluder depth
The original PCSS approach tests shadow map samples in the search area to figure out whether they’re occluders. Then, the average depth for the occluders is calculated. Our approach will use Chebyshev’s inequality theory to again get an upper probability bound of occlusion for the entire search area. From this probability, we can calculate the average depth (see [Yang 2010] Variance Soft Shadow Mapping). 
is the total average depth, 
and 
are the average depths for occluders and non-occluders, respectively. 
is the probability of a point being lit (ie: a non-occluder). Then, we can make the following observation:
![\[ d_{total} = d_{occluders} * (1 - p) + d_{nonOccluders} * p \]](https://www.derschmale.com/blog/wp-content/ql-cache/quicklatex.com-2a22ee4ceacb17b6b6788b88d865b0cc_l3.png "Rendered by QuickLaTeX.com")
![\[ \iff \]](https://www.derschmale.com/blog/wp-content/ql-cache/quicklatex.com-7a937a1e383c7f6082223acc1cd2f969_l3.png "Rendered by QuickLaTeX.com")
![\[ d_{occluders} = \frac{d_{total} - d_{nonOccluders} * p}{(1 - p)} \]](https://www.derschmale.com/blog/wp-content/ql-cache/quicklatex.com-cbefad9e25127a55a51ce425ff6a4552_l3.png "Rendered by QuickLaTeX.com")
Since the area search approach is based on the simplification that receiving and casting planes are parallel to the shadow map, is the reference depth. The only thing left to do is calculate the probability and the average depth for the entire search area. We’ll do this again very coarsely: we’ll use a single sample in a coarser mip level. It’s not exactly precise, but seems to work well enough. The average depth is already right there in the red channel of the shadow map, and for the probability we’ll use the upper bound again. The code below is for illustration, don’t expect any optimizations:
// shadowMapSize is the size of the shadow map in texels
// shadowMapCoord is the shadow map coord projected into directional light space (so z contains its depth)
float GetAverageOccluderDepth(float searchAreaSize, int shadowMapSize, float4 shadowMapCoord)
{
// calculate the mip level corresponding to the search area
// Really, mipLevel would be a passed in as a constant.
float mipLevel = log2(searchAreaSize * shadowMapSize);
// retrieve the distribution's moments for the entire area
// shadowMapSampler is a trilinear sampler, not a comparison sampler
float2 moments = shadowMap.SampleLevel(shadowMapSampler, shadowMapCoord.xy, mipLevel);
float averageTotalDepth = moments.x; // assign for semantic clarity
float probability = UpperBoundShadow(moments, shadowMapCoord.z);
// prevent numerical issues
if (probability > .99) return 0.0;
// calculate the average occluder depth
return (averageTotalDepth - probability * shadowMapCoord.z) / (1.0 - probability);
}
Calculate penumbra size from average depth
We calculate the penumbra size in the same exact way as in the PCSS. For directional lights, this again doesn’t hold up very well (ah you missing light position!) Instead, we can simply use the distance to the average occluder as a scale factor instead. It’s fun when things get simpler!
// shadowMapSize is the size of the shadow map in texels
// shadowMapCoord is the shadow map coord projected into directional light space (so z contains its depth)
// penumbraScale is a value describing how fast the penumbra should go soft. It can also be used to control the world space fall-off (by projecting world space distances to depth values)
float EstimatePenumbraSize(float lightSize, int shadowMapSize, float4 shadowMapCoord, float penumbraScale)
{
// the search area covers twice the light size
float averageOccluderDepth = GetAverageOccluderDepth(lightSize, shadowMapSize, shadowMapCoord);
float penumbraSize = lightSize * (shadowMapCoord.z - averageOccluderDepth) * penumbraScale;
// clamp to the maximum softness, which matches the search area
return min(penumbraSize, lightSize);
}
Calculate occluder probability
Instead of using a pre-blurred mip-map chain or a SAT table, we’ll perform the filtering on the fly. We’ll start by sampling a fixed number of points in a Poisson disk distribution to get the (approximate) moments of the entire filter region (ie: the penumbra size). We’ll rotate the sample points randomly to reduce banding in favour of noise. This is essentially the same as percentage closer filtering, but using probabilities instead. So, a first draft:
float penumbraSize = EstimatePenumbraSize(lightSize, shadowMapSize, shadowMapCoord, penumbraScale);
float2 moments = 0.0;
// ditherTexture contains 2d rotation matrix (cos, -sin, sin, cos), this will tile the texture across the screen
float4 rotation = ditherTexture.SampleLevel(nearestWrapSampler, screenUV * screenSize / ditherTextureSize, 0) * 2.0 - 1.0;
for (int i = 0; i < numShadowSamples; ++i) {
// poissonDiskValues contain the sampling offsets in the unit circle
// scale by penumbraSize / 2 to get samples within the penumbra radius (penumbraSize is diameter)
float2 sampleOffset = poissonDiskValues[i] * penumbraSize / 2;
float4 coord = shadowMapCoord;
// add rotated sample offset using dithered sample
coord.x += sampleOffset.x * rotation.x + sampleOffset.y * rotation.y;
coord.y += sampleOffset.x * rotation.z + sampleOffset.y * rotation.w;
// shadowMapSampler is a trilinear sampler, not a comparison sampler
moments += shadowMap.Sample(shadowMapSampler, shadowMapCoord.xy);
}
moments /= numShadowSamples;
float lightContribution = UpperBoundShadow(moments, shadowMapCoord.z);
But we can do better, observing that when sampling the disk distribution, we’d get a better approximation if we could get the average over every disk’s area instead of only at the sample point. Again, we can use the mip levels to get an approximation. A Poisson disk distribution has a minimum distance between any two points, so we can use this to calculate the mip level to sample from. Let’s replace some of the shader code:
float penumbraSize = EstimatePenumbraSize(lightSize, shadowMapSize, shadowMapCoord);
float2 moments = 0.0;
// ditherTexture contains 2d rotation matrix (cos, -sin, sin, cos), this will tile the texture across the screen
float4 rotation = ditherTexture.SampleLevel(nearestWrapSampler, screenUV * screenSize / ditherTextureSize, 0) * 2.0 - 1.0;
// calculate the mip level for the disk sample's area
// Sample points are expected to be penumbraSize * poissonRadius * shadowMapSize texels apart
// poissonRadius is half the minimum distance in the disk distribution
float mipLevel = log2(penumbraSize * poissonRadius * shadowMapSize);
for (int i = 0; i < numShadowSamples; ++i) {
// poissonDiskValues contain the sampling offsets in the unit circle
// scale by penumbraSize / 2 to get samples within the penumbra radius (penumbraSize is diameter)
float2 sampleOffset = poissonDiskValues[i] * penumbraSize / 2;
float4 coord = shadowMapCoord;
// add rotated sample offset using dithered sample
coord.x += sampleOffset.x * rotation.x + sampleOffset.y * rotation.y;
coord.y += sampleOffset.x * rotation.z + sampleOffset.y * rotation.w;
// shadowMapSampler is a trilinear sampler, not a comparison sampler
moments += shadowMap.SampleLevel(shadowMapSampler, shadowMapCoord.xy, mipLevel);
}
moments /= numShadowSamples;
float lightContribution = UpperBoundShadow(moments, shadowMapCoord.z);
Not only does this give us a better estimate, it also reduces the noise from the random rotations because samples are expected to differ less. And what’s more, we don’t have to take all that many samples, even for quite large filter sizes! Another way of looking at this approach is as blurring the mip levels in the lighting shader – just enough to remove the jaggies – instead of doing so on the shadow map directly.
Shadow map bounds
As usual with soft shadows, there’s an issue with sampling outside the shadow map boundaries. For this reason, it may be required to extend the shadow maps to accommodate the largest penumbra size (our “lightSize” value). You might also want to keep the light size within certain limits so that most of the shadow map usage isn’t just there to provide the area not on the screen.
Conclusion
Note the difference between the vase’s shadow and that of the distant flag pole next to it.
A problem with contact shadows for thin casters
As I’ve said before, chances are slim that I’ll actually use this approach in my own code – unless if it’s for something very specific and manageable. Variance shadow maps have too many issues for me:
- Light leaking: while this can be ameliorated easily, but not completely avoided, solutions have a strong impact on the softness of the penumbra, destroying some of our hard work.
- Thin caster leaks: the closer a point gets to the occluder, the smaller the upper bound becomes (as it’s less and less likely for it to be in the shadow). This creates severe light leaking close to thin casters.
But, again, VSMs have been used with success in the past, so who knows this article still may be of use to someone. You might run into other problems too, if you’re up for pursuing this.
And perhaps VSMs could be used only to perform the area search, and PCF sampling for the occlusion tests, which should remove any light leaking problems. Anyway, I’m up to receive any ideas, comments or poisonous arrows!
