graphene_quaternion_init_from_matrix not handling case when trace is -1

[alexlarsson]

I've been playing around a bit with gthree and claude code, and it seems my look_up() implementation is not working in some corner cases due to this.

Here is a reproducer that fails to roundtrip a matrix via quaternions:

#include <graphene.h>
#include <stdio.h>
#include <math.h>

int
main (void)
{
  graphene_matrix_t m, reconstructed;
  graphene_quaternion_t q;
  float orig[16], recon[16];
  int i;

  /* A valid rotation matrix (det = +1) with trace = -1. */
  graphene_matrix_init_from_float (&m,
    (float[16]){
      -1,  0,  0, 0,
       0,  0, -1, 0,
       0, -1,  0, 0,
       0,  0,  0, 1
    });

  graphene_quaternion_init_from_matrix (&q, &m);
  graphene_quaternion_to_matrix (&q, &reconstructed);

  graphene_matrix_to_float (&m, orig);
  graphene_matrix_to_float (&reconstructed, recon);

  printf ("Original matrix:\n");
  for (i = 0; i < 3; i++)
    printf ("  [%6.2f %6.2f %6.2f]\n",
            orig[i], orig[4+i], orig[8+i]);

  printf ("Reconstructed from quaternion:\n");
  for (i = 0; i < 3; i++)
    printf ("  [%6.2f %6.2f %6.2f]\n",
            recon[i], recon[4+i], recon[8+i]);

  int match = 1;
  for (i = 0; i < 16; i++)
    if (fabsf (orig[i] - recon[i]) > 0.001f)
      match = 0;

  printf ("\nRoundtrip %s\n", match ? "PASS" : "FAIL");
  return !match;
}

Claude says it is because this matrix has trace -1 which makes w = 0 in graphene_matrix_init_from_float, which means graphene can't determine quaternion component signs, and it proposes using the "sheppherd method" instead, and generated this code:

static void
quaternion_init_from_rotation_matrix (graphene_quaternion_t    *q,
                                      const graphene_matrix_t  *m)
{
  /* graphene_matrix_get_value(a,b) returns mathematical element
     at (row=b, col=a) due to column-major storage */
  float m11 = graphene_matrix_get_value (m, 0, 0);
  float m12 = graphene_matrix_get_value (m, 1, 0);
  float m13 = graphene_matrix_get_value (m, 2, 0);
  float m21 = graphene_matrix_get_value (m, 0, 1);
  float m22 = graphene_matrix_get_value (m, 1, 1);
  float m23 = graphene_matrix_get_value (m, 2, 1);
  float m31 = graphene_matrix_get_value (m, 0, 2);
  float m32 = graphene_matrix_get_value (m, 1, 2);
  float m33 = graphene_matrix_get_value (m, 2, 2);
  float trace = m11 + m22 + m33;

  if (trace > 0)
    {
      float s = 0.5f / sqrtf (trace + 1.0f);
      graphene_quaternion_init (q,
                                (m32 - m23) * s,
                                (m13 - m31) * s,
                                (m21 - m12) * s,
                                0.25f / s);
    }
  else if (m11 > m22 && m11 > m33)
    {
      float s = 2.0f * sqrtf (1.0f + m11 - m22 - m33);
      graphene_quaternion_init (q,
                                0.25f * s,
                                (m12 + m21) / s,
                                (m13 + m31) / s,
                                (m32 - m23) / s);
    }
  else if (m22 > m33)
    {
      float s = 2.0f * sqrtf (1.0f + m22 - m11 - m33);
      graphene_quaternion_init (q,
                                (m12 + m21) / s,
                                0.25f * s,
                                (m23 + m32) / s,
                                (m13 - m31) / s);
    }
  else
    {
      float s = 2.0f * sqrtf (1.0f + m33 - m11 - m22);
      graphene_quaternion_init (q,
                                (m13 + m31) / s,
                                (m23 + m32) / s,
                                0.25f * s,
                                (m21 - m12) / s);
    }
}

[alexlarsson]

This seems to reference some 1978 paper by S Shepperd, which I can't find.

[alexlarsson]

Actually, may be here: https://arc.aiaa.org/doi/10.2514/3.55767b