[mutter/gbsneto/graphene-matrix: 13/47] cogl/matrix: Workaround precision loss

[Georges Basile Stavracas Neto <gbsneto src gnome org>]

commit 9bb280096cb66b5105036c0d13335a1038f36d4f
Author: Georges Basile Stavracas Neto <georges stavracas gmail com>
Date:   Thu Sep 10 18:28:27 2020 -0300

    cogl/matrix: Workaround precision loss
    
    This is a complicated issue, to please bear with me.
    
    The way to detect whether a matrix is invertible or not (i.e.
    whether it's not a "singular" matrix, or not) is by checking
    if the determinant equals 0. So far, so good.
    
    Both graphene_matrix_t and CoglMatrix use single-precision
    floats to store their 4x4 matrices. Graphene uses vectorized
    operations to optimize determinant calculation, while Cogl
    tries to keep track of the matrix type and has special-purpose
    determinant functions for different matrix types (the most
    common one being a 3D matrix).
    
    Cogl, however, has a fundamentally flawed check for whether
    the matrix is invertible or not. Have a look:
    
    ```
    float det;
    
    …
    
    if (det*det < 1e-25)
       return FALSE;
    ```
    
    Notice that 1e-25 is *way* smaller than FLT_EPSILON. This
    check is fundamentally flawed.
    
    "In practice, what does it break?", the reader might ask.
    Well, in this case, the answer is opposite of that: Cogl
    inverts matrices that should not be invertible. Let's see
    an example: the model-view-projection of a 4K monitor. It
    looks like this:
    
    ```
    | +0,002693 +0,000000 +0,000000 +0,000000 |
    | +0,000000 -0,002693 +0,000000 +0,000000 |
    | +0,000000 +0,000000 +0,002693 +0,000000 |
    | -5,169809 +2,908017 -5,036834 +1,000000 |
    ```
    
    The determinant of this matrix is -0.000000019530306557.
    It evidently is smaller than FLT_EPSILON. In this situation,
    Cogl would happily calculate the inverse matrix, whereas
    Graphene (correctly) bails out and thinks it's a singular
    matrix.
    
    This commit works around that by exploiting the maths around
    it. The basis of it is:
    
      inverse(scalar * M) = (1/scalar) * M'
    
    which can be extrapolated to:
    
      inverse(M) = scalar * inverse(scalar * M) = M'
    
    In other words, scaling the to-be-inversed matrix, then
    scaling the inverse matrix by the same factor, gives us
    the desired inverse.
    
    I'm sorry for everyone that has to read through this :(
    
    https://gitlab.gnome.org/GNOME/mutter/-/merge_requests/1439

 cogl/cogl/cogl-matrix.c | 32 +++++++++++++++++++++++++++-----
 1 file changed, 27 insertions(+), 5 deletions(-)
---
diff --git a/cogl/cogl/cogl-matrix.c b/cogl/cogl/cogl-matrix.c
index 43c535cfaa..4ae1c07bf9 100644
--- a/cogl/cogl/cogl-matrix.c
+++ b/cogl/cogl/cogl-matrix.c
@@ -135,16 +135,38 @@ cogl_debug_matrix_print (const CoglMatrix *matrix)
  * given matrix type.  In case of failure, updates the MAT_FLAG_SINGULAR flag,
  * and copies the identity matrix into CoglMatrix::inv.
  */
+
+#define SCALE 1000.F
+
+static inline gboolean
+calculate_inverse (CoglMatrix *matrix)
+{
+  graphene_matrix_t scaled;
+  graphene_matrix_t m;
+  gboolean invertible;
+
+  graphene_matrix_init_scale (&scaled, SCALE, SCALE, SCALE);
+
+  /* Float precision is a limiting factor */
+  graphene_matrix_init_from_matrix (&m, &matrix->m);
+  graphene_matrix_multiply (&m, &scaled, &m);
+
+  invertible = graphene_matrix_inverse (&m, &matrix->inv);
+
+  if (invertible)
+    graphene_matrix_multiply (&scaled, &matrix->inv, &matrix->inv);
+  else
+    graphene_matrix_init_identity (&matrix->inv);
+
+  return invertible;
+}
+
 static gboolean
 _cogl_matrix_update_inverse (CoglMatrix *matrix)
 {
   if (matrix->flags & COGL_MATRIX_FLAG_DIRTY_INVERSE)
     {
-      gboolean invertible;
-
-      invertible = graphene_matrix_inverse (&matrix->m, &matrix->inv);
-
-      if (invertible)
+      if (calculate_inverse (matrix))
         matrix->flags &= ~COGL_MATRIX_FLAG_SINGULAR;
       else
         matrix->flags |= COGL_MATRIX_FLAG_SINGULAR;