An Invertible Split Scheme for Cascading Shadow Maps

split-scheme.txt

           An Invertible Split Scheme for Cascading Shadow Maps

                           by Leonard Ritter
                              Duangle GbR
                              
                        Second revision: 2016/01/01


DISCLAIMER: This is original research and has not been verified by a
second or third party yet.


This is a proposal for a new function for the practical split 
scheme[1] for shadow map cascades by Zhang et al. The new 
function can be inverted in a shader to easily project 
distance z  back to layer index i. An intuitive replacement 
for the weighing factor sigma is proposed as well.


In the original article[1], a split scheme function was suggested to
find good boundaries for individual shadow map cascades, so that
given a layer index i, the optimal depth z for that index could
be derived; given a normal scalar x which is retrieved via

                                         i
                                x = -----------
                                    layer_count

and a sigma constant S in the range 0..1, which would allow the
user interpolate linearly between an equidistant depth distribution
and an ideal exponential one to bias the resolution distribution,
the original split scheme equation would be

                          f
                  z = pow(-, x) n S + ((f - n) x + n) (1 - S) = f(x)
                          n

n and f designate the distance of the near and far planes respectively.

Unfortunately this function can't be inverted. With this solution,
it is necessary to precompute and pass the layer boundaries as 
coefficients to the pixel shader. This also limits applicability,
as the number of addressable layers is bounded by storage. 

The New Split Scheme
====================

Through experimentation and some guesswork, I found an alternative
way to bias the distribution, which is purely exponential and
has no additive polynomial component, so inversion is trivial.

This bias effectively projects the near and far distances to a
slope further up the exponential curve, which dampens the initial
slow momentum and causes the function to approach a more linear
scaling.

The original sigma constant S is squared so that C = sqrt(S),
and the split scheme function is altered so that 

                         f                 f
            z = (f - n - -) (1 - pow(-------------, x)) + n = f(x)
                         C           (n - f) C + f

At the original proposed default of S = 0.5, the function starts
with near identical momentum, but converges marginally sooner
towards f; across all values of S, the first half of the function
distributes very similarly, so this function can easily be used 
as a drop-in replacement.

The only limitation with this modified function is that, while there
is a solution for C = 0, z and x can not be computed, which poses no
issue as S = 0 describes a simple linear mapping that is trivial to
invert, and undesirable in the context of efficient depth layer usage.

Analysis
========

http://i.imgur.com/9OeUMBz.png

The above picture shows the distribution for sigma = 0.5; the
green and yellow curves show the boundaries at sigma 0.0 and 1.0,
which are the same for both old and new functions.
The purple curve is the original split scheme, the blue curve
is the modified one.


Shader Usage
============

In the shader, the layer index can now be retrieved by resolving x from a
distance z and multiplying it with layer_count to retrieve the interpolated
index of the layer to be retrieved.

To calculate x from z, we invert the new split scheme function so that

                (z - far) C + far                  far
       x = log(--------------------) / log(--------------------) = f(z)
               (near - far) C + far        (near - far) C + far

For optimization, the inner factors and the second log term can be stored in
three constant variables U, V, W:

                           d = far * (1 - C) + near * C

                                         C
                                     U = -
                                         d

                                   far * (1 - C)
                               V = -------------
                                         d

                                                                far
                                 W = layer_count * log(2) / log(---)
                                                                 d

With these optimizations, the resolve reduces to

                              i = log2(z U + V) W


Future Work
===========

I did not find the original weighing parameter sigma to be particularly
easy to understand and apply. In practice, I found myself aiming to get
an optimal balance between resolution and the range of the first cascade
layer by guessing values for sigma; There was no clear relationship
between the parameter and its implications.

Ideally, I would prefer a new balancing variable z1 instead of sigma, which
defines the distance at which the first shadow map layer ends, and is
used to calculate the bias C.

Unfortunately, I could not find a way to solve the equation towards C,
so I leave this for further exploration. Until a simpler method can be found,
bisecting towards C would be the nearest choice; A routine would 
methodically attempt to solve the split scheme function by adjusting C until
z = z1 for x = 1 / layer_count.


References
==========

[1] Parallel-Split Shadow Maps on Programmable GPUs
    http://http.developer.nvidia.com/GPUGems3/gpugems3_ch10.html


Thanks to Sean Barrett and Fabian Giesen for help and suggestions.