fix shadow stability when eye is far from world origin

in stable mode the scale was ever so slightly varying with the 
camera position, because it was calculated from the camera frustum in
world-space, this variation was amplified when the camera is far from 
the origin, which eventually caused the modulo needed for snapping the
shadowmap projection to widely vary, leading to the instability.

We now calculate the camera frustum sphere in view space, which is
guaranteed to be constant. If "shadow caster mode" is chosen, we 
quantize the scale a little bit so it stays constant.

The snapping code itself has been cleaned.

filament/src/ShadowMap.cpp

@@ -110,7 +110,7 @@ math::mat4f ShadowMap::getPointLightViewMatrix(backend::TextureCubemapFace face,
 }
 
 ShadowMap::ShaderParameters ShadowMap::updateDirectional(FEngine& engine,
-        const FScene::LightSoa& lightData, size_t index,
+        FScene::LightSoa const& lightData, size_t index,
         filament::CameraInfo const& camera,
         ShadowMapInfo const& shadowMapInfo,
         SceneInfo const& sceneInfo) noexcept {
@@ -216,7 +216,7 @@ ShadowMap::ShaderParameters ShadowMap::updateDirectional(FEngine& engine,
     }
 
     // Now that we know the znear (-lsLightFrustumBounds.max.z), adjust the light's position such
-    // that znear = 0, this is only need for VSM, but doesn't hurt PCF.
+    // that znear = 0, this is only needed for VSM, but doesn't hurt PCF.
     const mat4f Mv = getDirectionalLightViewMatrix(direction, direction * -lsLightFrustumBounds.max.z);
 
     // near / far planes are specified relative to the direction the eye is looking at
@@ -240,8 +240,11 @@ ShadowMap::ShaderParameters ShadowMap::updateDirectional(FEngine& engine,
         const float4 shadowReceiverVolumeBoundingSphere = computeBoundingSphere(
                 wsShadowReceiversVolume.getCorners().data(), 8);
 
-        // in stable mode we simply take the view volume, bounding sphere
-        viewVolumeBoundingSphere = computeBoundingSphere(wsViewFrustumVertices, 8);
+        // in stable mode we simply take the view volume bounding sphere, but we calculate it
+        // in view space, so that it's perfectly stable.
+        float3 vertices[8];
+        computeFrustumCorners(vertices, inverse(cullingProjection), sceneInfo.csNearFar);
+        viewVolumeBoundingSphere = computeBoundingSphere(vertices, 8);
 
         if (shadowReceiverVolumeBoundingSphere.w < viewVolumeBoundingSphere.w) {
 
@@ -320,45 +323,65 @@ ShadowMap::ShaderParameters ShadowMap::updateDirectional(FEngine& engine,
     //
     //   In LiPSM mode, we're using the warped space here.
 
-    Aabb bounds;
-    if (params.options.stable && viewVolumeBoundingSphere.w > 0) {
-        bounds = compute2DBounds(Mv, viewVolumeBoundingSphere);
-    } else {
-        bounds = compute2DBounds(WLMpMv, wsClippedShadowReceiverVolume.data(), vertexCount);
-    }
-    lsLightFrustumBounds.min.xy = bounds.min.xy;
-    lsLightFrustumBounds.max.xy = bounds.max.xy;
-
+    float2 s, o;
     if (params.options.stable) {
-        // in stable mode we can't do anything that can change the scaling of the texture
+        if (viewVolumeBoundingSphere.w > 0) {
+            s = 1.0f / viewVolumeBoundingSphere.w;
+            o = mat4f::project(Mv * camera.model, viewVolumeBoundingSphere.xyz).xy;
+        } else {
+            Aabb const bounds = compute2DBounds(Mv,
+                    wsClippedShadowReceiverVolume.data(), vertexCount);
+            if (UTILS_UNLIKELY((bounds.min.x >= bounds.max.x) || (bounds.min.y >= bounds.max.y))) {
+                // this could happen if the only thing visible is a perfectly horizontal or
+                // vertical thin line
+                mHasVisibleShadows = false;
+                return {};
+            }
+            assert_invariant(bounds.min.x < bounds.max.x);
+            assert_invariant(bounds.min.y < bounds.max.y);
+
+            s = 2.0f / float2(bounds.max.xy - bounds.min.xy);
+            o = float2(bounds.max.xy + bounds.min.xy) * 0.5f;
+
+            // Quantize the scale in world-space units. This value can be very small because
+            // if it wasn't for floating-point imprecision, the scale would be a constant.
+            double2 const quantizer = 0.0625;
+            s = 1.0 / (ceil(1.0 / (s * quantizer)) * quantizer);
+        }
     } else {
+        Aabb const bounds = compute2DBounds(WLMpMv,
+                wsClippedShadowReceiverVolume.data(), vertexCount);
+        lsLightFrustumBounds.min.xy = bounds.min.xy;
+        lsLightFrustumBounds.max.xy = bounds.max.xy;
         // For directional lights, we further constraint the light frustum to the
         // intersection of the shadow casters & shadow receivers in light-space.
         // ** This relies on the 1-texel shadow map border **
         if (engine.debug.shadowmap.focus_shadowcasters) {
             intersectWithShadowCasters(lsLightFrustumBounds, WLMpMv, wsShadowCastersVolume);
         }
-    }
+        if (UTILS_UNLIKELY((lsLightFrustumBounds.min.x >= lsLightFrustumBounds.max.x) ||
+                           (lsLightFrustumBounds.min.y >= lsLightFrustumBounds.max.y))) {
+            // this could happen if the only thing visible is a perfectly horizontal or
+            // vertical thin line
+            mHasVisibleShadows = false;
+            return {};
+        }
+        assert_invariant(lsLightFrustumBounds.min.x < lsLightFrustumBounds.max.x);
+        assert_invariant(lsLightFrustumBounds.min.y < lsLightFrustumBounds.max.y);
 
-    if (UTILS_UNLIKELY((lsLightFrustumBounds.min.x >= lsLightFrustumBounds.max.x) ||
-                       (lsLightFrustumBounds.min.y >= lsLightFrustumBounds.max.y))) {
-        // this could happen if the only thing visible is a perfectly horizontal or
-        // vertical thin line
-        mHasVisibleShadows = false;
-        return {};
-    }
+        s = 2.0f / float2(bounds.max.xy - bounds.min.xy);
+        o = float2(bounds.max.xy + bounds.min.xy) * 0.5f;
 
-    assert_invariant(lsLightFrustumBounds.min.x < lsLightFrustumBounds.max.x);
-    assert_invariant(lsLightFrustumBounds.min.y < lsLightFrustumBounds.max.y);
+        // TODO: we could quantize `s` here to give some stability when lispsm is disabled,
+        //       however, the quantization paramater should probably be user settable.
+    }
 
-    // compute focus scale and offset
-    float2 s = 2.0f / float2(lsLightFrustumBounds.max.xy - lsLightFrustumBounds.min.xy);
-    float2 o =   -s * float2(lsLightFrustumBounds.max.xy + lsLightFrustumBounds.min.xy) * 0.5f;
+    // adjust offset for scale
+    o = -s * o;
 
-    if (params.options.stable) {
-        // Use the world origin as reference point, fixed w.r.t. the camera
-        snapLightFrustum(s, o, Mv, camera.worldOrigin[3].xyz,
-                1.0f / float(shadowMapInfo.shadowDimension));
+    if (!useLispsm) {
+        // stabilize the shadowmap in all modes, except lispsm which can never be stable
+        snapLightFrustum(s, o, Mv, camera.worldOrigin, shadowMapInfo.shadowDimension);
     }
 
     const mat4f F(mat4f::row_major_init {
@@ -825,21 +848,39 @@ void ShadowMap::computeFrustumCorners(float3* UTILS_RESTRICT out,
 }
 
 void ShadowMap::snapLightFrustum(float2& s, float2& o,
-        mat4f const& Mv, float3 worldOrigin, float2 shadowMapResolution) noexcept {
+        mat4f const& Mv, mat4 worldOrigin, int2 resolution) noexcept {
 
-    auto fmod = [](float2 x, float2 y) -> float2 {
-        auto mod = [](float x, float y) -> float { return std::fmod(x, y); };
-        return float2{ mod(x[0], y[0]), mod(x[1], y[1]) };
+    auto proj = [](mat4 m, double4 v) -> double3 {
+        // for directional light p.w == 1, exactly
+        auto p = m * v;
+        assert_invariant(p.w == 1.0);
+        return p.xyz;
     };
 
-    // This snaps the shadow map bounds to texels.
-    // The 2.0 comes from Mv having a NDC in the range -1,1 (so a range of 2).
-    const float2 r = 2.0f * shadowMapResolution;
-    o -= fmod(o, r);
+    auto fract = [](auto v) {
+        using namespace std;
+        using T = decltype(v);
+        return fmod(v, T{1});
+    };
+
+    const mat4 F(mat4::row_major_init {
+            s.x, 0.0, 0.0, o.x,
+            0.0, s.y, 0.0, o.y,
+            0.0, 0.0, 1.0, 0.0,
+            0.0, 0.0, 0.0, 1.0,
+    });
+
+    // The (resolution * 0.5) comes from Mv having a NDC in the range -1,1 (so a range of 2).
+
+    // focused light-space
+    mat4 const FMv{ F * Mv };
 
     // This offsets the texture coordinates, so it has a fixed offset w.r.t the world
-    const float2 lsOrigin = mat4f::project(Mv, worldOrigin).xy * s;
-    o -= fmod(lsOrigin, r);
+    double2 const lsOrigin = proj(FMv, worldOrigin[3]).xy;
+    double2 const d = (fract(lsOrigin * resolution * 0.5) * 2.0) / resolution;
+
+    // adjust offset
+    o -= d;
 }
 
 size_t ShadowMap::intersectFrustumWithBox(

filament/src/ShadowMap.h

@@ -222,7 +222,7 @@ class ShadowMap {
             const math::float3& dir);
 
     static inline void snapLightFrustum(math::float2& s, math::float2& o,
-            math::mat4f const& Mv, math::float3 worldOrigin, math::float2 shadowMapResolution) noexcept;
+            math::mat4f const& Mv, math::mat4 worldOrigin, math::int2 resolution) noexcept;
 
     static inline void computeFrustumCorners(math::float3* out,
             const math::mat4f& projectionViewInverse, math::float2 csNearFar = { -1.0f, 1.0f }) noexcept;

filament/src/ShadowMapManager.cpp

@@ -504,19 +504,13 @@ ShadowMapManager::ShadowTechnique ShadowMapManager::updateCascadeShadowMaps(FEng
             splitPercentages[i] = options.cascadeSplitPositions[i - 1];
         }
 
-        const CascadeSplits::Params p{
+        const CascadeSplits splits({
                 .proj = cameraInfo.cullingProjection,
                 .near = vsNear,
                 .far = vsFar,
                 .cascadeCount = cascadeCount,
                 .splitPositions = splitPercentages
-        };
-        if (p != mCascadeSplitParams) {
-            mCascadeSplits = CascadeSplits{ p };
-            mCascadeSplitParams = p;
-        }
-
-        const CascadeSplits& splits = mCascadeSplits;
+        });
 
         // The split positions uniform is a float4. To save space, we chop off the first split position
         // (which is the near plane, and doesn't need to be communicated to the shaders).
@@ -530,7 +524,7 @@ ShadowMapManager::ShadowTechnique ShadowMapManager::updateCascadeShadowMaps(FEng
 
         mShadowMappingUniforms.cascadeSplits = wsSplitPositionUniform;
 
-        // when computing the required bias we need a half-texel size, so we multiply by 0.5 here.
+        // When computing the required bias we need a half-texel size, so we multiply by 0.5 here.
         // note: normalBias is set to zero for VSM
         const float normalBias = shadowMapInfo.vsm ? 0.0f : 0.5f * lcm.getShadowNormalBias(0);
 

filament/src/ShadowMapManager.h

@@ -145,17 +145,8 @@ class ShadowMapManager {
             float far = 0.0f;
             size_t cascadeCount = 1;
             std::array<float, SPLIT_COUNT> splitPositions = { 0.0f };
-
-            bool operator!=(const Params& rhs) const {
-                return proj != rhs.proj ||
-                       near != rhs.near ||
-                       far != rhs.far ||
-                       cascadeCount != rhs.cascadeCount ||
-                       splitPositions != rhs.splitPositions;
-            }
         };
 
-        CascadeSplits() noexcept : CascadeSplits(Params{}) {}
         explicit CascadeSplits(Params const& params) noexcept;
 
         // Split positions in world-space.
@@ -186,9 +177,6 @@ class ShadowMapManager {
 
     SoftShadowOptions mSoftShadowOptions;
 
-    CascadeSplits::Params mCascadeSplitParams;
-    CascadeSplits mCascadeSplits;
-
     mutable TypedUniformBuffer<ShadowUib> mShadowUb;
     backend::Handle<backend::HwBufferObject> mShadowUbh;
 
@@ -204,8 +192,8 @@ class ShadowMapManager {
             utils::FixedCapacityVector<ShadowMap*>::with_capacity(
                     CONFIG_MAX_SHADOWMAPS - CONFIG_MAX_SHADOW_CASCADES) };
 
-    // inline storage for all our ShadowMap objects, we can't easily use a std::array<> directly.
-    // because ShadowMap doesn't have a default ctor, and we avoid out-of-line allocations.
+    // Inline storage for all our ShadowMap objects, we can't easily use a std::array<> directly.
+    // Because ShadowMap doesn't have a default ctor, and we avoid out-of-line allocations.
     // Each ShadowMap is currently 40 bytes (total of 2.5KB for 64 shadow maps)
     using ShadowMapStorage = std::aligned_storage<sizeof(ShadowMap), alignof(ShadowMap)>::type;
     std::array<ShadowMapStorage, CONFIG_MAX_SHADOWMAPS> mShadowMapCache;

libs/math/include/math/TVecHelpers.h

@@ -443,6 +443,37 @@ class TVecFunctions {
         return v;
     }
 
+    template<typename U>
+    friend inline
+    VECTOR<T> MATH_PURE fmod(VECTOR<T> const& x, VECTOR<U> const& y) {
+        VECTOR<T> r;
+        for (size_t i = 0; i < r.size(); i++) {
+            r[i] = std::fmod(x[i], y[i]);
+        }
+        return r;
+    }
+
+    template<typename U>
+    friend inline
+    VECTOR<T> MATH_PURE remainder(VECTOR<T> const& x, VECTOR<U> const& y) {
+        VECTOR<T> r;
+        for (size_t i = 0; i < r.size(); i++) {
+            r[i] = std::remainder(x[i], y[i]);
+        }
+        return r;
+    }
+
+    template<typename U>
+    friend inline
+    VECTOR<T> MATH_PURE remquo(VECTOR<T> const& x, VECTOR<U> const& y,
+            VECTOR<int>* q) {
+        VECTOR<T> r;
+        for (size_t i = 0; i < r.size(); i++) {
+            r[i] = std::remquo(x[i], y[i], &((*q)[i]));
+        }
+        return r;
+    }
+
     friend inline VECTOR<T> MATH_PURE inversesqrt(VECTOR<T> v) {
         for (size_t i = 0; i < v.size(); i++) {
             v[i] = T(1) / std::sqrt(v[i]);