Having major problems transforming frustum from camera space to light space for shadow mapping

Question

I'm trying to implement cascaded shadow mapping in my OpenGL/C++ renderer. I've successfully managed to implement directional shadow mapping when the ortho matrix is built from arbitrary numbers and stays at the origin of the scene. From everything I can see, the problem I'm having is deciding the bounds of the ortho matrices based on the slice of the view frustum.

I've found two ways of calculating it, both are related; The first way creates a projection matrix using the near/far planes for the current slice, then inverts it. Then I take the corners of a cube in NDC space (each axis ranging from [-1:+1]), multiply it through the inverted projection matrix, the inverted view matrix, then the light space matrix; then divide the whole thing by the corner's 'w'.

vec4 corners[8] =
    {
        // Near plane
        { 1, 1, 1, 1 },
        { -1, 1, 1, 1 },
        { 1, -1, 1, 1 },
        { -1, -1, 1, 1 },
        // Far plane
        { 1, 1, -1, 1 },
        { -1, 1, -1, 1 },
        { 1, -1, -1, 1 },
        { -1, -1, -1, 1 }
    };

    // Create a projection matrix for this cascade.
    // This will transform the corner from camera space to view space.
    mat4 proj = mat4::perspective(fov, aspect, near_plane, far_plane);

    // B.inverted()*A.inverted() == C.inverted() == (A*B).inverted()
    // This is now an inverted viewproj matrix, which transforms the corner from camera space to world space.
    mat4 inv_viewproj = (proj * view_matrix).inverted();

    // This is the matrix that transforms each corner from camera space to light space.
    // (this_frame.view is the light-space view matrix)
    mat4 matrix = this_frame.view * matrix;
    
    float minx = FLT_MAX;
    float maxx = -FLT_MAX;
    float miny = FLT_MAX;
    float maxy = -FLT_MAX;
    float minz = FLT_MAX;
    float maxz = -FLT_MAX;

    // For each frustum corner...
    for (uint32_t i = 0; i < 8; ++i)
    {
        corners[i] = matrix * corners[i];
        corners[i] /= corners[i].w;

        minx = min(minx, corners[i].x);
        maxx = max(maxx, corners[i].x);
        miny = min(miny, corners[i].y);
        maxy = max(maxy, corners[i].y);
        minz = min(minz, corners[i].z);
        maxz = max(maxz, corners[i].z);
    }

    this_frame.proj = mat4::ortho(minx, maxx, miny, maxy, minz, maxz);

The second way is similar, it starts with a set of corners in view-space, alleviating the need to divide by w and saving a few matrix inversions.

mat4 inv_view = view_matrix.inverted();
float tan_half_hfov = tanf(to_radians(fov * 0.5f));
float tan_half_vfov = tanf(to_radians((fov * 0.5f) / aspect));
//...

    float xn = near_plane * tan_half_hfov;
    float yn = near_plane * tan_half_vfov;
    float xf = far_plane * tan_half_hfov;
    float yf = far_plane * tan_half_vfov;

    // These are the corners of the screen in view space.
    // We'll be transforming these in order to get the dimensions of the view frustum.
    vec4 corners[8] =
    {
        // Near plane
        { xn, yn, -near_plane, 1 },
        { -xn, yn, -near_plane, 1 },
        { xn, -yn, -near_plane, 1 },
        { -xn, -yn, -near_plane, 1 },
        // Far plane
        { xf, yf, -far_plane, 1 },
        { -xf, yf, -far_plane, 1 },
        { xf, -yf, -far_plane, 1 },
        { -xf, -yf, -far_plane, 1 }
    };

    // This is the matrix that transforms each corner from camera space to light space.
    mat4 matrix = this_frame.view * inv_view;
    
    float minx = FLT_MAX;
    float maxx = -FLT_MAX;
    float miny = FLT_MAX;
    float maxy = -FLT_MAX;
    float minz = FLT_MAX;
    float maxz = -FLT_MAX;

    // For each frustum corner...
    for (uint32_t i = 0; i < 8; ++i)
    {
        corners[i] = matrix * corners[i];

        minx = min(minx, corners[i].x);
        maxx = max(maxx, corners[i].x);
        miny = min(miny, corners[i].y);
        maxy = max(maxy, corners[i].y);
        minz = min(minz, corners[i].z);
        maxz = max(maxz, corners[i].z);
    }

    this_frame.proj = mat4::ortho(minx, maxx, miny, maxy, minz, maxz);

Both produce incorrect, but kinda-sorta-close-to-correct results. The second method appears to produce a tighter box.

Interestingly, both appear to mostly work when the camera moves along the x-axis in light-space. Whenever the camera moves along the y/z axes in light space, you get some very wrong results.

Here's some images of what I'm talking about. The first of each pair is the scene according to the normal camera perspective, the second is from the view of the light (first cascade, I'm not even going to try the second cascade until the first works correctly). The red region is what should be included in the first cascade, the green region is later cascades, and the unshaded region is outside of the view frustum.

(These screenshots were taken using method #1)

1: http://imgur.com/FCxvDy1, http://imgur.com/WUFoEV4

2: http://imgur.com/wUTzkwO, http://imgur.com/fVvTTlG

So far so good, right? well...

3: http://imgur.com/Mxunbwg, http://imgur.com/yeFYxVp

4: http://imgur.com/J1FAuKv, http://imgur.com/bzrorwe

For a while I thought that the problem might be in my mat4::inverted() function, but I've tried 2 different implementations and they both give the same results. The one I'm currently using is:

mat4 mat4::inverted() const
{
float invdet = 1.0f / determinant();

// Note: This constructor for mat4 accepts a series of 16 floats as a row-major matrix,
// as that's simply more natural to write out,
// then converts it automatically to the correct column-major storage.
return {
    (m12*m23*m31 - m13*m22*m31 + m13*m21*m32 - m11*m23*m32 - m12*m21*m33 + m11*m22*m33) * invdet,
    (m03*m22*m31 - m02*m23*m31 - m03*m21*m32 + m01*m23*m32 + m02*m21*m33 - m01*m22*m33) * invdet,
    (m02*m13*m31 - m03*m12*m31 + m03*m11*m32 - m01*m13*m32 - m02*m11*m33 + m01*m12*m33) * invdet,
    (m03*m12*m21 - m02*m13*m21 - m03*m11*m22 + m01*m13*m22 + m02*m11*m23 - m01*m12*m23) * invdet,
    (m13*m22*m30 - m12*m23*m30 - m13*m20*m32 + m10*m23*m32 + m12*m20*m33 - m10*m22*m33) * invdet,
    (m02*m23*m30 - m03*m22*m30 + m03*m20*m32 - m00*m23*m32 - m02*m20*m33 + m00*m22*m33) * invdet,
    (m03*m12*m30 - m02*m13*m30 - m03*m10*m32 + m00*m13*m32 + m02*m10*m33 - m00*m12*m33) * invdet,
    (m02*m13*m20 - m03*m12*m20 + m03*m10*m22 - m00*m13*m22 - m02*m10*m23 + m00*m12*m23) * invdet,
    (m11*m23*m30 - m13*m21*m30 + m13*m20*m31 - m10*m23*m31 - m11*m20*m33 + m10*m21*m33) * invdet,
    (m03*m21*m30 - m01*m23*m30 - m03*m20*m31 + m00*m23*m31 + m01*m20*m33 - m00*m21*m33) * invdet,
    (m01*m13*m30 - m03*m11*m30 + m03*m10*m31 - m00*m13*m31 - m01*m10*m33 + m00*m11*m33) * invdet,
    (m03*m11*m20 - m01*m13*m20 - m03*m10*m21 + m00*m13*m21 + m01*m10*m23 - m00*m11*m23) * invdet,
    (m12*m21*m30 - m11*m22*m30 - m12*m20*m31 + m10*m22*m31 + m11*m20*m32 - m10*m21*m32) * invdet,
    (m01*m22*m30 - m02*m21*m30 + m02*m20*m31 - m00*m22*m31 - m01*m20*m32 + m00*m21*m32) * invdet,
    (m02*m11*m30 - m01*m12*m30 - m02*m10*m31 + m00*m12*m31 + m01*m10*m32 - m00*m11*m32) * invdet,
    (m01*m12*m20 - m02*m11*m20 + m02*m10*m21 - m00*m12*m21 - m01*m10*m22 + m00*m11*m22) * invdet };
}


float mat4::determinant() const
{
    return 
        m03*m12*m21*m30 - m02*m13*m21*m30 - m03*m11*m22*m30 + m01*m13*m22*m30 +
        m02*m11*m23*m30 - m01*m12*m23*m30 - m03*m12*m20*m31 + m02*m13*m20*m31 +
        m03*m10*m22*m31 - m00*m13*m22*m31 - m02*m10*m23*m31 + m00*m12*m23*m31 +
        m03*m11*m20*m32 - m01*m13*m20*m32 - m03*m10*m21*m32 + m00*m13*m21*m32 +
        m01*m10*m23*m32 - m00*m11*m23*m32 - m02*m11*m20*m33 + m01*m12*m20*m33 +
        m02*m10*m21*m33 - m00*m12*m21*m33 - m01*m10*m22*m33 + m00*m11*m22*m33;
}

What am I doing wrong?

Answer 1

vec4 corners[8] =
    {
        // Near plane
        { 1, 1, 1, 1 },
        { -1, 1, 1, 1 },
        { 1, -1, 1, 1 },
        { -1, -1, 1, 1 },
        // Far plane
        { 1, 1, -1, 1 },
        { -1, 1, -1, 1 },
        { 1, -1, -1, 1 },
        { -1, -1, -1, 1 }
    };

Let me tell you why here is wrong. you are looking to a screen. The near is depth 0. and the far is depth 1.

vec4 corners[8] =
    {
        // Near plane
        { 1, 1, 0, 1 },
        { -1, 1, 0, 1 },
        { 1, -1, 0, 1 },
        { -1, -1, 0, 1 },
        // Far plane
        { 1, 1, 1, 1 },
        { -1, 1, 1, 1 },
        { 1, -1, 1, 1 },
        { -1, -1, 1, 1 }
    };

.....// somestuff

    // For each frustum corner...
    for (uint32_t i = 0; i < 8; ++i)
    {
        corners[i] = invVP * corners[i];
        corners[i] /= corners[i].w;

        ... go on....

    }

Let me know if works.