diff options
| -rw-r--r-- | source/blender/blenkernel/BKE_blender_version.h | 4 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/armature.c | 45 | ||||
| -rw-r--r-- | source/blender/blenkernel/intern/armature_test.cc | 67 | ||||
| -rw-r--r-- | source/blender/blenloader/intern/versioning_300.c | 138 |
4 files changed, 216 insertions, 38 deletions
diff --git a/source/blender/blenkernel/BKE_blender_version.h b/source/blender/blenkernel/BKE_blender_version.h index 5ee86f9446e..899d21683e4 100644 --- a/source/blender/blenkernel/BKE_blender_version.h +++ b/source/blender/blenkernel/BKE_blender_version.h @@ -39,13 +39,13 @@ extern "C" { /* Blender file format version. */ #define BLENDER_FILE_VERSION BLENDER_VERSION -#define BLENDER_FILE_SUBVERSION 35 +#define BLENDER_FILE_SUBVERSION 36 /* Minimum Blender version that supports reading file written with the current * version. Older Blender versions will test this and show a warning if the file * was written with too new a version. */ #define BLENDER_FILE_MIN_VERSION 300 -#define BLENDER_FILE_MIN_SUBVERSION 26 +#define BLENDER_FILE_MIN_SUBVERSION 36 /** User readable version string. */ const char *BKE_blender_version_string(void); diff --git a/source/blender/blenkernel/intern/armature.c b/source/blender/blenkernel/intern/armature.c index a266718dcfc..0fa4c6e47e8 100644 --- a/source/blender/blenkernel/intern/armature.c +++ b/source/blender/blenkernel/intern/armature.c @@ -2227,39 +2227,47 @@ void mat3_vec_to_roll(const float mat[3][3], const float vec[3], float *r_roll) * </pre> * * When y is close to -1, computing 1 / (1 + y) will cause severe numerical instability, - * so we ignore it and normalize M instead. + * so we use a different approach based on x and z as inputs. * We know `y^2 = 1 - (x^2 + z^2)`, and `y < 0`, hence `y = -sqrt(1 - (x^2 + z^2))`. * - * Since x and z are both close to 0, we apply the binomial expansion to the first order: - * `y = -sqrt(1 - (x^2 + z^2)) = -1 + (x^2 + z^2) / 2`. Which gives: + * Since x and z are both close to 0, we apply the binomial expansion to the second order: + * `y = -sqrt(1 - (x^2 + z^2)) = -1 + (x^2 + z^2) / 2 + (x^2 + z^2)^2 / 8`, which allows + * eliminating the problematic `1` constant. + * + * A first order expansion allows simplifying to this, but second order is more precise: * <pre> * ┌ z^2 - x^2, -2 * x * z ┐ * M* = 1 / (x^2 + z^2) * │ │ * └ -2 * x * z, x^2 - z^2 ┘ * </pre> + * + * P.S. In the end, this basically is a heavily optimized version of Damped Track +Y. */ void vec_roll_to_mat3_normalized(const float nor[3], const float roll, float r_mat[3][3]) { - const float SAFE_THRESHOLD = 1.0e-5f; /* theta above this value has good enough precision. */ - const float CRITICAL_THRESHOLD = 1.0e-9f; /* above this is safe under certain conditions. */ + const float SAFE_THRESHOLD = 6.1e-3f; /* theta above this value has good enough precision. */ + const float CRITICAL_THRESHOLD = 2.5e-4f; /* true singularity if xz distance is below this. */ const float THRESHOLD_SQUARED = CRITICAL_THRESHOLD * CRITICAL_THRESHOLD; const float x = nor[0]; const float y = nor[1]; const float z = nor[2]; - const float theta = 1.0f + y; /* remapping Y from [-1,+1] to [0,2]. */ - const float theta_alt = x * x + z * z; /* Helper value for matrix calculations.*/ + float theta = 1.0f + y; /* remapping Y from [-1,+1] to [0,2]. */ + const float theta_alt = x * x + z * z; /* squared distance from origin in x,z plane. */ float rMatrix[3][3], bMatrix[3][3]; BLI_ASSERT_UNIT_V3(nor); - /* When theta is close to zero (nor is aligned close to negative Y Axis), + /* Determine if the input is far enough from the true singularity of this type of + * transformation at (0,-1,0), where roll becomes 0/0 undefined without a limit. + * + * When theta is close to zero (nor is aligned close to negative Y Axis), * we have to check we do have non-null X/Z components as well. * Also, due to float precision errors, nor can be (0.0, -0.99999994, 0.0) which results * in theta being close to zero. This will cause problems when theta is used as divisor. */ - if (theta > SAFE_THRESHOLD || (theta > CRITICAL_THRESHOLD && theta_alt > THRESHOLD_SQUARED)) { + if (theta > SAFE_THRESHOLD || theta_alt > THRESHOLD_SQUARED) { /* nor is *not* aligned to negative Y-axis (0,-1,0). */ bMatrix[0][1] = -x; @@ -2268,18 +2276,15 @@ void vec_roll_to_mat3_normalized(const float nor[3], const float roll, float r_m bMatrix[1][2] = z; bMatrix[2][1] = -z; - if (theta > SAFE_THRESHOLD) { - /* nor differs significantly from negative Y axis (0,-1,0): apply the general case. */ - bMatrix[0][0] = 1 - x * x / theta; - bMatrix[2][2] = 1 - z * z / theta; - bMatrix[2][0] = bMatrix[0][2] = -x * z / theta; - } - else { - /* nor is close to negative Y axis (0,-1,0): apply the special case. */ - bMatrix[0][0] = (x + z) * (x - z) / -theta_alt; - bMatrix[2][2] = -bMatrix[0][0]; - bMatrix[2][0] = bMatrix[0][2] = 2.0f * x * z / theta_alt; + if (theta <= SAFE_THRESHOLD) { + /* When nor is close to negative Y axis (0,-1,0) the theta precision is very bad, + * so recompute it from x and z instead, using the series expansion for sqrt. */ + theta = theta_alt * 0.5f + theta_alt * theta_alt * 0.125f; } + + bMatrix[0][0] = 1 - x * x / theta; + bMatrix[2][2] = 1 - z * z / theta; + bMatrix[2][0] = bMatrix[0][2] = -x * z / theta; } else { /* nor is very close to negative Y axis (0,-1,0): use simple symmetry by Z axis. */ diff --git a/source/blender/blenkernel/intern/armature_test.cc b/source/blender/blenkernel/intern/armature_test.cc index 8ebb91ffc74..3d22351e9a6 100644 --- a/source/blender/blenkernel/intern/armature_test.cc +++ b/source/blender/blenkernel/intern/armature_test.cc @@ -185,12 +185,12 @@ static double find_flip_boundary(double x, double z) TEST(vec_roll_to_mat3_normalized, FlippedBoundary1) { - EXPECT_NEAR(find_flip_boundary(0, 1), 2.40e-4, 0.01e-4); + EXPECT_NEAR(find_flip_boundary(0, 1), 2.50e-4, 0.01e-4); } TEST(vec_roll_to_mat3_normalized, FlippedBoundary2) { - EXPECT_NEAR(find_flip_boundary(1, 1), 3.39e-4, 0.01e-4); + EXPECT_NEAR(find_flip_boundary(1, 1), 2.50e-4, 0.01e-4); } /* Test cases close to the -Y axis. */ @@ -218,9 +218,9 @@ TEST(vec_roll_to_mat3_normalized, Flipped3) { /* If normalized_vector is in a critical range close to -Y, apply the special case. */ 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}, + 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}}; + {-1.000000f, -2.5e-4f, 0.000000f}}; test_vec_roll_to_mat3_normalized(input, 0.0f, expected_roll_mat, false); } @@ -304,6 +304,65 @@ TEST(vec_roll_to_mat3_normalized, Roll1) test_vec_roll_to_mat3_normalized(input, float(M_PI * 0.5), expected_roll_mat); } +/** Test that the matrix is orthogonal for an input close to -Y. */ +static double test_vec_roll_to_mat3_orthogonal(double s, double x, double z) +{ + const float input[3] = {float(x), float(s * sqrt(1 - x * x - z * z)), float(z)}; + + return test_vec_roll_to_mat3_normalized(input, 0.0f, NULL); +} + +/** Test that the matrix is orthogonal for a range of inputs close to -Y. */ +static void test_vec_roll_to_mat3_orthogonal(double s, double x1, double x2, double y1, double y2) +{ + const int count = 5000; + double delta = 0; + double tmax = 0; + + for (int i = 0; i <= count; i++) { + double t = double(i) / count; + double det = test_vec_roll_to_mat3_orthogonal(s, interpd(x2, x1, t), interpd(y2, y1, t)); + + /* Find and report maximum error in the matrix determinant. */ + double curdelta = abs(det - 1); + if (curdelta > delta) { + delta = curdelta; + tmax = t; + } + } + + printf(" Max determinant error %.10f at %f.\n", delta, tmax); +} + +#define TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(name, s, x1, x2, y1, y2) \ + TEST(vec_roll_to_mat3_normalized, name) \ + { \ + test_vec_roll_to_mat3_orthogonal(s, x1, x2, y1, y2); \ + } + +/* Moving from -Y towards X. */ +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_000_005, -1, 0, 0, 3e-4, 0.005) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_000_010, -1, 0, 0, 0.005, 0.010) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_000_050, -1, 0, 0, 0.010, 0.050) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_000_100, -1, 0, 0, 0.050, 0.100) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_000_200, -1, 0, 0, 0.100, 0.200) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_000_300, -1, 0, 0, 0.200, 0.300) + +/* Moving from -Y towards X and Y. */ +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_005_005, -1, 3e-4, 0.005, 3e-4, 0.005) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_010_010, -1, 0.005, 0.010, 0.005, 0.010) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_050_050, -1, 0.010, 0.050, 0.010, 0.050) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_100_100, -1, 0.050, 0.100, 0.050, 0.100) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoN_200_200, -1, 0.100, 0.200, 0.100, 0.200) + +/* Moving from +Y towards X. */ +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoP_000_005, 1, 0, 0, 0, 0.005) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoP_000_100, 1, 0, 0, 0.005, 0.100) + +/* Moving from +Y towards X and Y. */ +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoP_005_005, 1, 0, 0.005, 0, 0.005) +TEST_VEC_ROLL_TO_MAT3_ORTHOGONAL(OrthoP_100_100, 1, 0.005, 0.100, 0.005, 0.100) + class BKE_armature_find_selected_bones_test : public testing::Test { protected: bArmature arm; diff --git a/source/blender/blenloader/intern/versioning_300.c b/source/blender/blenloader/intern/versioning_300.c index 4c6b70982a4..6b57bdf9a9c 100644 --- a/source/blender/blenloader/intern/versioning_300.c +++ b/source/blender/blenloader/intern/versioning_300.c @@ -47,6 +47,7 @@ #include "BKE_action.h" #include "BKE_animsys.h" +#include "BKE_armature.h" #include "BKE_asset.h" #include "BKE_collection.h" #include "BKE_deform.h" @@ -1097,6 +1098,112 @@ static void version_geometry_nodes_add_attribute_input_settings(NodesModifierDat } } +/* Copy of the function before the fixes. */ +static void legacy_vec_roll_to_mat3_normalized(const float nor[3], + const float roll, + float r_mat[3][3]) +{ + const float SAFE_THRESHOLD = 1.0e-5f; /* theta above this value has good enough precision. */ + const float CRITICAL_THRESHOLD = 1.0e-9f; /* above this is safe under certain conditions. */ + const float THRESHOLD_SQUARED = CRITICAL_THRESHOLD * CRITICAL_THRESHOLD; + + const float x = nor[0]; + const float y = nor[1]; + const float z = nor[2]; + + const float theta = 1.0f + y; /* remapping Y from [-1,+1] to [0,2]. */ + const float theta_alt = x * x + z * z; /* Helper value for matrix calculations.*/ + float rMatrix[3][3], bMatrix[3][3]; + + BLI_ASSERT_UNIT_V3(nor); + + /* When theta is close to zero (nor is aligned close to negative Y Axis), + * we have to check we do have non-null X/Z components as well. + * Also, due to float precision errors, nor can be (0.0, -0.99999994, 0.0) which results + * in theta being close to zero. This will cause problems when theta is used as divisor. + */ + if (theta > SAFE_THRESHOLD || (theta > CRITICAL_THRESHOLD && theta_alt > THRESHOLD_SQUARED)) { + /* nor is *not* aligned to negative Y-axis (0,-1,0). */ + + bMatrix[0][1] = -x; + bMatrix[1][0] = x; + bMatrix[1][1] = y; + bMatrix[1][2] = z; + bMatrix[2][1] = -z; + + if (theta > SAFE_THRESHOLD) { + /* nor differs significantly from negative Y axis (0,-1,0): apply the general case. */ + bMatrix[0][0] = 1 - x * x / theta; + bMatrix[2][2] = 1 - z * z / theta; + bMatrix[2][0] = bMatrix[0][2] = -x * z / theta; + } + else { + /* nor is close to negative Y axis (0,-1,0): apply the special case. */ + bMatrix[0][0] = (x + z) * (x - z) / -theta_alt; + bMatrix[2][2] = -bMatrix[0][0]; + bMatrix[2][0] = bMatrix[0][2] = 2.0f * x * z / theta_alt; + } + } + else { + /* nor is very close to negative Y axis (0,-1,0): use simple symmetry by Z axis. */ + unit_m3(bMatrix); + bMatrix[0][0] = bMatrix[1][1] = -1.0; + } + + /* Make Roll matrix */ + axis_angle_normalized_to_mat3(rMatrix, nor, roll); + + /* Combine and output result */ + mul_m3_m3m3(r_mat, rMatrix, bMatrix); +} + +static void correct_bone_roll_value(const float head[3], + const float tail[3], + const float check_x_axis[3], + const float check_y_axis[3], + float *r_roll) +{ + const float SAFE_THRESHOLD = 1.0e-5f; + float vec[3], bone_mat[3][3], vec2[3]; + + /* Compute the Y axis vector. */ + sub_v3_v3v3(vec, tail, head); + normalize_v3(vec); + + /* Only correct when in the danger zone. */ + if (1.0f + vec[1] < SAFE_THRESHOLD * 2 && (vec[0] || vec[2])) { + /* Use the armature matrix to double-check if adjustment is needed. + * This should minimize issues if the file is bounced back and forth between + * 2.92 and 2.91, provided Edit Mode isn't entered on the armature in 2.91. */ + vec_roll_to_mat3(vec, *r_roll, bone_mat); + + BLI_assert(dot_v3v3(bone_mat[1], check_y_axis) > 0.999f); + + if (dot_v3v3(bone_mat[0], check_x_axis) < 0.999f) { + /* Recompute roll using legacy code to interpret the old value. */ + legacy_vec_roll_to_mat3_normalized(vec, *r_roll, bone_mat); + mat3_to_vec_roll(bone_mat, vec2, r_roll); + BLI_assert(compare_v3v3(vec, vec2, FLT_EPSILON)); + } + } +} + +/* Update the armature Bone roll fields for bones very close to -Y direction. */ +static void do_version_bones_roll(ListBase *lb) +{ + LISTBASE_FOREACH (Bone *, bone, lb) { + /* Parent-relative orientation (used for posing). */ + correct_bone_roll_value( + bone->head, bone->tail, bone->bone_mat[0], bone->bone_mat[1], &bone->roll); + + /* Absolute orientation (used for Edit mode). */ + correct_bone_roll_value( + bone->arm_head, bone->arm_tail, bone->arm_mat[0], bone->arm_mat[1], &bone->arm_roll); + + do_version_bones_roll(&bone->childbase); + } +} + /* NOLINTNEXTLINE: readability-function-size */ void blo_do_versions_300(FileData *fd, Library *UNUSED(lib), Main *bmain) { @@ -1839,18 +1946,7 @@ void blo_do_versions_300(FileData *fd, Library *UNUSED(lib), Main *bmain) } } - /** - * Versioning code until next subversion bump goes here. - * - * \note Be sure to check when bumping the version: - * - "versioning_userdef.c", #blo_do_versions_userdef - * - "versioning_userdef.c", #do_versions_theme - * - * \note Keep this message at the bottom of the function. - */ - { - /* Keep this block, even when empty. */ - + if (!MAIN_VERSION_ATLEAST(bmain, 300, 36)) { /* Update the `idnames` for renamed geometry and function nodes. */ LISTBASE_FOREACH (bNodeTree *, ntree, &bmain->nodetrees) { if (ntree->type != NTREE_GEOMETRY) { @@ -1871,5 +1967,23 @@ void blo_do_versions_300(FileData *fd, Library *UNUSED(lib), Main *bmain) version_node_id(ntree, GEO_NODE_SET_MATERIAL, "GeometryNodeSetMaterial"); version_node_id(ntree, GEO_NODE_SPLIT_EDGES, "GeometryNodeSplitEdges"); } + + /* Update bone roll after a fix to vec_roll_to_mat3_normalized. */ + LISTBASE_FOREACH (bArmature *, arm, &bmain->armatures) { + do_version_bones_roll(&arm->bonebase); + } + } + + /** + * Versioning code until next subversion bump goes here. + * + * \note Be sure to check when bumping the version: + * - "versioning_userdef.c", #blo_do_versions_userdef + * - "versioning_userdef.c", #do_versions_theme + * + * \note Keep this message at the bottom of the function. + */ + { + /* Keep this block, even when empty. */ } } |