From dfa1c7e554fcc3ddd40780fb8555cdd6e90eaba3 Mon Sep 17 00:00:00 2001 From: Alexander Gavrilov Date: Sat, 21 Nov 2020 16:06:02 +0300 Subject: 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 --- source/blender/blenkernel/intern/armature_test.cc | 246 +++++++++++++++++----- source/blender/blenlib/BLI_math_matrix.h | 2 + source/blender/blenlib/intern/math_matrix.c | 22 ++ 3 files changed, 219 insertions(+), 51 deletions(-) (limited to 'source') 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 { diff --git a/source/blender/blenlib/BLI_math_matrix.h b/source/blender/blenlib/BLI_math_matrix.h index 2b0c3db21ee..241acebffa3 100644 --- a/source/blender/blenlib/BLI_math_matrix.h +++ b/source/blender/blenlib/BLI_math_matrix.h @@ -57,6 +57,7 @@ void copy_m4_m4_db(double m1[4][4], const double m2[4][4]); void copy_m3_m3d(float m1[3][3], const double m2[3][3]); /* float->double */ +void copy_m3d_m3(double m1[3][3], const float m2[3][3]); void copy_m4d_m4(double m1[4][4], const float m2[4][4]); void swap_m3m3(float m1[3][3], float m2[3][3]); @@ -291,6 +292,7 @@ float determinant_m3( float a1, float a2, float a3, float b1, float b2, float b3, float c1, float c2, float c3); float determinant_m3_array(const float m[3][3]); float determinant_m4_mat3_array(const float m[4][4]); +double determinant_m3_array_db(const double m[3][3]); float determinant_m4(const float m[4][4]); #define PSEUDOINVERSE_EPSILON 1e-8f diff --git a/source/blender/blenlib/intern/math_matrix.c b/source/blender/blenlib/intern/math_matrix.c index 554506e90e7..b6d80d76be1 100644 --- a/source/blender/blenlib/intern/math_matrix.c +++ b/source/blender/blenlib/intern/math_matrix.c @@ -180,6 +180,21 @@ void copy_m4_m2(float m1[4][4], const float m2[2][2]) m1[3][3] = 1.0f; } +void copy_m3d_m3(double m1[3][3], const float m2[3][3]) +{ + m1[0][0] = m2[0][0]; + m1[0][1] = m2[0][1]; + m1[0][2] = m2[0][2]; + + m1[1][0] = m2[1][0]; + m1[1][1] = m2[1][1]; + m1[1][2] = m2[1][2]; + + m1[2][0] = m2[2][0]; + m1[2][1] = m2[2][1]; + m1[2][2] = m2[2][2]; +} + void copy_m4d_m4(double m1[4][4], const float m2[4][4]) { m1[0][0] = m2[0][0]; @@ -1113,6 +1128,13 @@ float determinant_m4_mat3_array(const float m[4][4]) m[2][0] * (m[0][1] * m[1][2] - m[0][2] * m[1][1])); } +double determinant_m3_array_db(const double m[3][3]) +{ + return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - + m[1][0] * (m[0][1] * m[2][2] - m[0][2] * m[2][1]) + + m[2][0] * (m[0][1] * m[1][2] - m[0][2] * m[1][1])); +} + bool invert_m3_ex(float m[3][3], const float epsilon) { float tmp[3][3]; -- cgit v1.2.3