Fix precision issues and a bug in vec_roll_to_mat3_normalized.
When the input vector gets close to -Y, y and theta becomes totally unreliable. It is thus necessary to compute the result in a different way based on x and z. The code already had a special case, but: - The threshold for using the special case was way too low. - The special case was not precise enough to extend the threshold. - The special case math had a sign error, resulting in a jump. This adds tests for the computation precision and fixes the issues by adjusting the threshold, and replacing the special case with one based on a quadratic Taylor expansion of sqrt instead of linear. Replacing the special case fixes the bug and results in a compatibility break, requiring versioning for the roll of affected bones. Differential Revision: https://developer.blender.org/D9551
This commit is contained in:
parent
df445cc571
commit
16eafdadf6
Notes:
blender-bot
2023-09-13 08:48:34 +02:00
Referenced by issue #82455, The conversion of roll to matrix breaks in some cases
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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. */
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -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. */
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue