diff options
author | Alexander Gavrilov <angavrilov@gmail.com> | 2020-11-21 21:45:14 +0300 |
---|---|---|
committer | Alexander Gavrilov <angavrilov@gmail.com> | 2021-10-20 12:58:19 +0300 |
commit | 16eafdadf6040fb84bacf657ac0bf16a78e1057e (patch) | |
tree | f4f89df1993cc18a60775f2a1a21bc6fd67878ea /source/blender/blenkernel/intern/armature.c | |
parent | df445cc571bd1cf7fab4c5c8474f5e185a757fe2 (diff) |
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
Diffstat (limited to 'source/blender/blenkernel/intern/armature.c')
-rw-r--r-- | source/blender/blenkernel/intern/armature.c | 45 |
1 files changed, 25 insertions, 20 deletions
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. */ |