Split and extend unit tests for vec_roll_to_mat3_normalized.

Separate the huge test into huge logical parts and add more cases
to check. Also add a utility to check that the matrix is orthogonal,
with arbitrary epsilon values and calculations in double.

A couple of tests deliberately fail, to be fixed in following commits.

Ref D9551

1 files changed, 195 insertions, 51 deletions

diff --git a/source/blender/blenkernel/intern/armature_test.cc b/source/blender/blenkernel/intern/armature_test.cc
index 2994563175f..8ebb91ffc74 100644
--- a/source/blender/blenkernel/intern/armature_test.cc
+++ b/source/blender/blenkernel/intern/armature_test.cc
@@ -30,6 +30,36 @@ namespace blender::bke::tests {
 
 static const float FLOAT_EPSILON = 1.2e-7;
 
+static const float SCALE_EPSILON = 3.71e-5;
+static const float ORTHO_EPSILON = 5e-5;
+
+/** Test that the matrix is orthogonal, i.e. has no scale or shear within acceptable precision. */
+static double EXPECT_M3_ORTHOGONAL(const float mat[3][3],
+                                   double epsilon_scale,
+                                   double epsilon_ortho)
+{
+  /* Do the checks in double precision to avoid precision issues in the checks themselves. */
+  double dmat[3][3];
+  copy_m3d_m3(dmat, mat);
+
+  /* Check individual axis scaling. */
+  EXPECT_NEAR(len_v3_db(dmat[0]), 1.0, epsilon_scale);
+  EXPECT_NEAR(len_v3_db(dmat[1]), 1.0, epsilon_scale);
+  EXPECT_NEAR(len_v3_db(dmat[2]), 1.0, epsilon_scale);
+
+  /* Check orthogonality. */
+  EXPECT_NEAR(dot_v3v3_db(dmat[0], dmat[1]), 0.0, epsilon_ortho);
+  EXPECT_NEAR(dot_v3v3_db(dmat[0], dmat[2]), 0.0, epsilon_ortho);
+  EXPECT_NEAR(dot_v3v3_db(dmat[1], dmat[2]), 0.0, epsilon_ortho);
+
+  /* Check determinant to detect flipping and as a secondary volume change check. */
+  double determinant = determinant_m3_array_db(dmat);
+
+  EXPECT_NEAR(determinant, 1.0, epsilon_ortho);
+
+  return determinant;
+}
+
 TEST(mat3_vec_to_roll, UnitMatrix)
 {
   float unit_matrix[3][3];
@@ -93,71 +123,185 @@ TEST(mat3_vec_to_roll, Rotationmatrix)
   }
 }
 
-TEST(vec_roll_to_mat3_normalized, Rotationmatrix)
+/** Generic function to test vec_roll_to_mat3_normalized. */
+static double test_vec_roll_to_mat3_normalized(const float input[3],
+                                               float roll,
+                                               const float expected_roll_mat[3][3],
+                                               bool normalize = true)
 {
-  float negative_y_axis[3][3];
-  unit_m3(negative_y_axis);
-  negative_y_axis[0][0] = negative_y_axis[1][1] = -1.0f;
-
-  const float roll = 0.0f;
+  float input_normalized[3];
   float roll_mat[3][3];
 
-  /* If normalized_vector is -Y, simple symmetry by Z axis. */
-  {
-    const float normalized_vector[3] = {0.0f, -1.0f, 0.0f};
-    vec_roll_to_mat3_normalized(normalized_vector, roll, roll_mat);
-    EXPECT_M3_NEAR(roll_mat, negative_y_axis, FLT_EPSILON);
+  if (normalize) {
+    /* The vector is renormalized to replicate the actual usage. */
+    normalize_v3_v3(input_normalized, input);
+  }
+  else {
+    copy_v3_v3(input_normalized, input);
   }
 
-  /* If normalized_vector is far enough from -Y, apply the general case. */
-  {
-    const float expected_roll_mat[3][3] = {{1.000000f, 0.000000f, 0.000000f},
-                                           {0.000000f, -0.999989986f, -0.000000f},
-                                           {0.000000f, 0.000000f, 1.000000f}};
+  vec_roll_to_mat3_normalized(input_normalized, roll, roll_mat);
 
-    const float normalized_vector[3] = {0.0f, -1.0f + 1e-5f, 0.0f};
-    vec_roll_to_mat3_normalized(normalized_vector, roll, roll_mat);
+  EXPECT_V3_NEAR(roll_mat[1], input_normalized, FLT_EPSILON);
+
+  if (expected_roll_mat) {
     EXPECT_M3_NEAR(roll_mat, expected_roll_mat, FLT_EPSILON);
   }
 
-#if 0
-  /* TODO: This test will pass after fixing T82455) */
-  /* If normalized_vector is close to -Y and
-   * it has X and Z values above a threshold,
-   * apply the special case. */
-  {
-    const float expected_roll_mat[3][3] = {{0.000000f, -9.99999975e-06f, 1.000000f},
-                                           {9.99999975e-06f, -0.999999881f, 9.99999975e-06f},
-                                           {1.000000f, -9.99999975e-06, 0.000000f}};
-    const float normalized_vector[3] = {1e-24, -0.999999881, 0};
-    vec_roll_to_mat3_normalized(normalized_vector, roll, roll_mat);
-    EXPECT_M3_NEAR(roll_mat, expected_roll_mat, FLT_EPSILON);
+  return EXPECT_M3_ORTHOGONAL(roll_mat, SCALE_EPSILON, ORTHO_EPSILON);
+}
+
+/** Binary search to test where the code switches to the most degenerate special case. */
+static double find_flip_boundary(double x, double z)
+{
+  /* Irrational scale factor to ensure values aren't 'nice', have a lot of rounding errors,
+   * and can't accidentally produce the exact result returned by the special case. */
+  const double scale = M_1_PI / 10;
+  double theta = x * x + z * z;
+  double minv = 0, maxv = 1e-2;
+
+  while (maxv - minv > FLT_EPSILON * 1e-3) {
+    double mid = (minv + maxv) / 2;
+
+    float roll_mat[3][3];
+    float input[3] = {float(x * mid * scale),
+                      -float(sqrt(1 - theta * mid * mid) * scale),
+                      float(z * mid * scale)};
+
+    normalize_v3(input);
+    vec_roll_to_mat3_normalized(input, 0, roll_mat);
+
+    /* The special case assigns exact constants rather than computing. */
+    if (roll_mat[0][0] == -1 && roll_mat[0][1] == 0 && roll_mat[2][1] == 0) {
+      minv = mid;
+    }
+    else {
+      maxv = mid;
+    }
   }
-#endif
 
+  return sqrt(theta) * (minv + maxv) * 0.5;
+}
+
+TEST(vec_roll_to_mat3_normalized, FlippedBoundary1)
+{
+  EXPECT_NEAR(find_flip_boundary(0, 1), 2.40e-4, 0.01e-4);
+}
+
+TEST(vec_roll_to_mat3_normalized, FlippedBoundary2)
+{
+  EXPECT_NEAR(find_flip_boundary(1, 1), 3.39e-4, 0.01e-4);
+}
+
+/* Test cases close to the -Y axis. */
+TEST(vec_roll_to_mat3_normalized, Flipped1)
+{
+  /* If normalized_vector is -Y, simple symmetry by Z axis. */
+  const float input[3] = {0.0f, -1.0f, 0.0f};
+  const float expected_roll_mat[3][3] = {
+      {-1.0f, 0.0f, 0.0f}, {0.0f, -1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat, false);
+}
+
+TEST(vec_roll_to_mat3_normalized, Flipped2)
+{
+  /* If normalized_vector is close to -Y and
+   * it has X and Z values below a threshold,
+   * simple symmetry by Z axis. */
+  const float input[3] = {1e-24, -0.999999881, 0};
+  const float expected_roll_mat[3][3] = {
+      {-1.0f, 0.0f, 0.0f}, {0.0f, -1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat, false);
+}
+
+TEST(vec_roll_to_mat3_normalized, Flipped3)
+{
   /* If normalized_vector is in a critical range close to -Y, apply the special case. */
-  {
-    const float expected_roll_mat[3][3] = {{0.000000f, -9.99999975e-06f, 1.000000f},
-                                           {9.99999975e-06f, -0.999999881f, 9.99999975e-06f},
-                                           {1.000000f, -9.99999975e-06f, 0.000000f}};
+  const float input[3] = {2.5e-4f, -0.999999881f, 2.5e-4f}; /* Corner Case. */
+  const float expected_roll_mat[3][3] = {{0.000000f, -2.5e-4f, 1.000000f},
+                                         {2.5e-4f, -0.999999881f, 2.5e-4f},
+                                         {1.000000f, -2.5e-4f, 0.000000f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat, false);
+}
 
-    const float normalized_vector[3] = {1e-5f, -0.999999881f, 1e-5f}; /* Corner Case. */
-    vec_roll_to_mat3_normalized(normalized_vector, roll, roll_mat);
-    EXPECT_M3_NEAR(roll_mat, expected_roll_mat, FLT_EPSILON);
-  }
+/* Test 90 degree rotations. */
+TEST(vec_roll_to_mat3_normalized, Rotate90_Z_CW)
+{
+  /* Rotate 90 around Z. */
+  const float input[3] = {1, 0, 0};
+  const float expected_roll_mat[3][3] = {{0, -1, 0}, {1, 0, 0}, {0, 0, 1}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
 
-  /* If normalized_vector is far enough from -Y, apply the general case. */
-  {
-    const float expected_roll_mat[3][3] = {{0.788675129f, -0.577350259f, -0.211324856f},
-                                           {0.577350259f, 0.577350259f, 0.577350259f},
-                                           {-0.211324856f, -0.577350259f, 0.788675129f}};
-
-    const float vector[3] = {1.0f, 1.0f, 1.0f}; /* Arbitrary Value. */
-    float normalized_vector[3];
-    normalize_v3_v3(normalized_vector, vector);
-    vec_roll_to_mat3_normalized(normalized_vector, roll, roll_mat);
-    EXPECT_M3_NEAR(roll_mat, expected_roll_mat, FLT_EPSILON);
-  }
+TEST(vec_roll_to_mat3_normalized, Rotate90_Z_CCW)
+{
+  /* Rotate 90 around Z. */
+  const float input[3] = {-1, 0, 0};
+  const float expected_roll_mat[3][3] = {{0, 1, 0}, {-1, 0, 0}, {0, 0, 1}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+TEST(vec_roll_to_mat3_normalized, Rotate90_X_CW)
+{
+  /* Rotate 90 around X. */
+  const float input[3] = {0, 0, -1};
+  const float expected_roll_mat[3][3] = {{1, 0, 0}, {0, 0, -1}, {0, 1, 0}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+TEST(vec_roll_to_mat3_normalized, Rotate90_X_CCW)
+{
+  /* Rotate 90 around X. */
+  const float input[3] = {0, 0, 1};
+  const float expected_roll_mat[3][3] = {{1, 0, 0}, {0, 0, 1}, {0, -1, 0}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+/* Test the general case when the vector is far enough from -Y. */
+TEST(vec_roll_to_mat3_normalized, Generic1)
+{
+  const float input[3] = {1.0f, 1.0f, 1.0f}; /* Arbitrary Value. */
+  const float expected_roll_mat[3][3] = {{0.788675129f, -0.577350259f, -0.211324856f},
+                                         {0.577350259f, 0.577350259f, 0.577350259f},
+                                         {-0.211324856f, -0.577350259f, 0.788675129f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+TEST(vec_roll_to_mat3_normalized, Generic2)
+{
+  const float input[3] = {1.0f, -1.0f, 1.0f}; /* Arbitrary Value. */
+  const float expected_roll_mat[3][3] = {{0.211324856f, -0.577350259f, -0.788675129f},
+                                         {0.577350259f, -0.577350259f, 0.577350259f},
+                                         {-0.788675129f, -0.577350259f, 0.211324856f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+TEST(vec_roll_to_mat3_normalized, Generic3)
+{
+  const float input[3] = {-1.0f, -1.0f, 1.0f}; /* Arbitrary Value. */
+  const float expected_roll_mat[3][3] = {{0.211324856f, 0.577350259f, 0.788675129f},
+                                         {-0.577350259f, -0.577350259f, 0.577350259f},
+                                         {0.788675129f, -0.577350259f, 0.211324856f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+TEST(vec_roll_to_mat3_normalized, Generic4)
+{
+  const float input[3] = {-1.0f, -1.0f, -1.0f}; /* Arbitrary Value. */
+  const float expected_roll_mat[3][3] = {{0.211324856f, 0.577350259f, -0.788675129f},
+                                         {-0.577350259f, -0.577350259f, -0.577350259f},
+                                         {-0.788675129f, 0.577350259f, 0.211324856f}};
+  test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat);
+}
+
+/* Test roll. */
+TEST(vec_roll_to_mat3_normalized, Roll1)
+{
+  const float input[3] = {1.0f, 1.0f, 1.0f}; /* Arbitrary Value. */
+  const float expected_roll_mat[3][3] = {{0.211324856f, 0.577350259f, -0.788675129f},
+                                         {0.577350259f, 0.577350259f, 0.577350259f},
+                                         {0.788675129f, -0.577350259f, -0.211324856f}};
+  test_vec_roll_to_mat3_normalized(input, float(M_PI * 0.5), expected_roll_mat);
 }
 
 class BKE_armature_find_selected_bones_test : public testing::Test {