New algorithm on computing bone matrix with improved accuracy #39470

New Issue

Bastien Montagne · 2014-03-27T20:40:49+01:00

tippisum commented

2014-03-27 20:40:49 +01:00

(Migrated from localhost:3001)

The bone matrices calculations in Blender hardly produce accurate results.
Sometimes it is annoying finding numbers like "7.549789415861596e-08" or "0.9999999403953552" spread out while the params assigned to that bone is indeed trival (e.g. a (1.0, 0.0, 0.0) with no roll or something).

And when I started to write some importers, things become worse.
Because Blender does not support assigning matrices to bones directly (while I cannot figure out what the reason is), I have to do some dirty math myself, and frustrated realize that I just have no way to import them and export later without messing up bone matrices.

Just a short time ago, I digged into Blender's source code and wanted to see what actually happens during the bone matrix calculating.
The key functions seems to be "vec_roll_to_mat3()" and its counterpart "mat3_to_vec_roll()".

There seems be quite a few problems with these functions, as the code comments said. In fact, they use extremely inaccurate method for solving rotation matrices, let alone the fact that all internal calculations are conducted in Single precision (I used to think it was Double).

The purpose of vec_roll_to_mat3() is to work out a rotation matrix for a certain bone vector and roll. This is started by trying to find an axis and angle that rotates Y-axis to the bone vector, then multiplying it with the roll. The later part is some sort of classic axis-angle conversion and not likely the major cause. But the algorithm applied to the former part does seem prone to numerical deviation.

At the first place, it uses traditional cross product to determine the rotating axis, with a threshold check (which changes over time, due to variant bugs related to it).
Well, cross product is less stable, and in fact it is unnecessary in this case. Let's begin with some solid geometry:

Assuming we want a rotation matrix corresponding to a bone vector v = (v.x, v.y, v.z), which is already normalized.
The rotation axis must be perpendicular to both the bone vector and the target.
Actually, we already know our target: It is the Y-axis. And it effectively means that the axis just lays on XZ-plane.
We recall that perpendicular to a vector is equivalent to perpendicular to its projection. So we project the bone vector to the XZ-plane.
Then it turns out to be a trivial plane geometry problem, and we can immediately write down the axis: a = (v.z, 0, -v.x). No actual calculations are really needed in order to get the rotation axis.

Unfortunately, the above strategy only solves a tiny part of the whole problem. The majority of numerical deviation comes from the highly inaccurate angle_normalized_v3v3() and axis_angle_normalized_to_mat3() call. They are inherently inaccurate, with more bad news that they both utilizes deviation-prone algorithm, and all operation are performed in Single precision.

Indeed, what we need is just a rotation matrix. We do not need the angle as it is useless to both Blender and the user. It is just an intermediate result, and converting it costs heavily on the numerical precision.

Is it possible to generate rotation matrix without calculating and converting rotation angle? The answer is YES.
with some time spent on the geometry and linear algebra of this specific rotation matrix, I finally got a effective algorithm computing it directly from bone vectors, and thus eliminates the precision cost of inverse trigonometric functions and radian measurement.

Still, assume we have a bone vector (v.x, v.y, v.z), which is already normalized. And we want to calculate a 3 * 3 rotation matrix, call it M, that will rotate Y-axis to this specific vector, in a pure linear algebra way.
A 3 * 3 rotation matrix contains 9 elements. However, we already have a not-bad start point. As stated before, this matrix will rotate Y-axis to the bone vector, means M * y = v. This immediately gives us 3 elements of the matrix for free. And the implementation of mat3_to_vec_roll() already made use of this characteristic.
Another 3 constraints come from rotation geometry, from which we know that a rotation axis is an eigenvector of its matrix, with eigenvalue 1. That gives us M * a = a, with a = (v.z, 0, -v.y) covered before. The representation of a is not normalized, but it does not matter since the equation still holds.
The final 3 is the most hard and tricky part. We know that this matrix will transform Y-axis to bone vector, and will keep axis invariant, and we still need one more to finally give out the result.
Now, let's have a imagine about the rotation. Y-axis is transformed to bone vector, then what vector (call it w) will now point to Y-axis?
The answer of this is not unique in general cases. But in this special situation, it can be resolved. The rotation axis is perpendicular to the plane determined by Y-axis and the bone vector, so w must resident in this plane too. And since all vectors are rotated by a same angle, v and w are symmetric to Y-axis, given w = (-v.x, v.y, -v.z) and the final equation M * u = y.

With all of these equation collected, it is now possible to finally solve the matrix M. Surprisingly, solving these equation is quite easy and the result looks reasonably nice:

    ┌ (x²y+z²)/(x²+z²), x,  xz(y-1)/(x²+z²) ┐
M = │  x(y²-1)/(x²+z²), y,  z(y²-1)/(x²+z²) │
    └  xz(y-1)/(x²+z²), z, (x²+z²y)/(x²+z²) ┘

It is easy to prove the correctness of this result, and such matrix is guaranteed to exist except when x² + z² = 0, which indicates the bone is parallel to Y-axis and need special treatment. And since all operations are internal, they can be done in a higher intermediate precision (e.g. Double). As long as x² + z² is not too small, this algorithm is likely to yield accuracy faith result.

A experimental patch was attached. But I'm not familiar with Blender's build system and thus unable to build and test it now. Also I am unable to use git so this is actually a diff from 2.69 release source. I'm sorry for that...
armature.diff

**The bone matrices calculations in Blender hardly produce accurate results.** Sometimes it is annoying finding numbers like "`7.549789415861596e-08`" or "`0.9999999403953552`" spread out while the params assigned to that bone is indeed trival (e.g. a `(1.0, 0.0, 0.0)` with no roll or something). And when I started to write some importers, things become worse. Because Blender does not support assigning matrices to bones directly (while I cannot figure out what the reason is), I have to do some dirty math myself, and frustrated realize that I just have no way to import them and export later without messing up bone matrices. Just a short time ago, I digged into Blender's source code and wanted to see what actually happens during the bone matrix calculating. The key functions seems to be "`vec_roll_to_mat3()`" and its counterpart "`mat3_to_vec_roll()`". There seems be quite a few problems with these functions, as the code comments said. In fact, they use extremely inaccurate method for solving rotation matrices, let alone the fact that all internal calculations are conducted in Single precision (I used to think it was Double). The purpose of `vec_roll_to_mat3()` is to work out a rotation matrix for a certain bone vector and roll. This is started by trying to find an axis and angle that rotates Y-axis to the bone vector, then multiplying it with the roll. The later part is some sort of classic axis-angle conversion and not likely the major cause. **But the algorithm applied to the former part does seem prone to numerical deviation.** At the first place, it uses traditional cross product to determine the rotating axis, with a threshold check (which changes over time, due to variant bugs related to it). Well, cross product is less stable, and in fact it is unnecessary in this case. Let's begin with some solid geometry: * Assuming we want a rotation matrix corresponding to a bone vector `v = (v.x, v.y, v.z)`, which is already normalized. * The rotation axis must be perpendicular to both the bone vector and the target. * Actually, we already know our target: **It is the Y-axis**. And it effectively means that the axis just lays on XZ-plane. * We recall that perpendicular to a vector is equivalent to perpendicular to its projection. So we project the bone vector to the XZ-plane. * Then it turns out to be a trivial plane geometry problem, and we can immediately write down the axis: `a = (v.z, 0, -v.x)`. **No actual calculations are really needed in order to get the rotation axis.** Unfortunately, the above strategy only solves a tiny part of the whole problem. **The majority of numerical deviation comes from the highly inaccurate `angle_normalized_v3v3()` and `axis_angle_normalized_to_mat3()` call.** They are inherently inaccurate, with more bad news that they both utilizes deviation-prone algorithm, and all operation are performed in Single precision. Indeed, what we need is just a rotation matrix. We do not need the angle as it is useless to both Blender and the user. It is just an intermediate result, and converting it costs heavily on the numerical precision. **Is it possible to generate rotation matrix without calculating and converting rotation angle? The answer is YES.** with some time spent on the geometry and linear algebra of this specific rotation matrix, I finally got a effective algorithm computing it directly from bone vectors, and thus eliminates the precision cost of inverse trigonometric functions and radian measurement. * Still, assume we have a bone vector `(v.x, v.y, v.z)`, which is already normalized. And we want to calculate a 3 * 3 rotation matrix, call it `M`, that will rotate Y-axis to this specific vector, in a pure linear algebra way. * A 3 * 3 rotation matrix contains 9 elements. However, we already have a not-bad start point. As stated before, this matrix will rotate Y-axis to the bone vector, means `M * y = v`. This immediately gives us 3 elements of the matrix for free. And the implementation of `mat3_to_vec_roll()` already made use of this characteristic. * Another 3 constraints come from rotation geometry, from which we know that a rotation axis is an eigenvector of its matrix, with eigenvalue 1. That gives us `M * a = a`, with `a = (v.z, 0, -v.y)` covered before. The representation of `a` is not normalized, but it does not matter since the equation still holds. * The final 3 is the most hard and tricky part. We know that this matrix will transform Y-axis to bone vector, and will keep axis invariant, and we still need one more to finally give out the result. * Now, let's have a imagine about the rotation. Y-axis is transformed to bone vector, then what vector (call it `w`) will now point to Y-axis? * The answer of this is not unique in general cases. But in this special situation, it can be resolved. The rotation axis is perpendicular to the plane determined by Y-axis and the bone vector, so `w` must resident in this plane too. And since all vectors are rotated by a same angle, `v` and `w` are symmetric to Y-axis, given `w = (-v.x, v.y, -v.z)` and the final equation `M * u = y`. With all of these equation collected, it is now possible to finally solve the matrix `M`. Surprisingly, solving these equation is quite easy and the result looks reasonably nice: ``` ┌ (x²y+z²)/(x²+z²), x, xz(y-1)/(x²+z²) ┐ M = │ x(y²-1)/(x²+z²), y, z(y²-1)/(x²+z²) │ └ xz(y-1)/(x²+z²), z, (x²+z²y)/(x²+z²) ┘ ``` It is easy to prove the correctness of this result, and such matrix is guaranteed to exist except when `x² + z² = 0`, which indicates the bone is parallel to Y-axis and need special treatment. And since all operations are internal, they can be done in a higher intermediate precision (e.g. Double). As long as `x² + z²` is not too small, this algorithm is likely to yield accuracy faith result. A experimental patch was attached. But I'm not familiar with Blender's build system and thus unable to build and test it now. Also I am unable to use git so this is actually a diff from 2.69 release source. I'm sorry for that... [armature.diff](https://archive.blender.org/developer/F83036/armature.diff)

tippisum commented

2014-03-27 20:40:49 +01:00

(Migrated from localhost:3001)

Changed status to: 'Open'

tippisum commented

2014-03-27 20:40:49 +01:00

(Migrated from localhost:3001)

Added subscriber: @tippisum

Bastien Montagne self-assigned this 2014-03-27 22:56:27 +01:00

Bastien Montagne commented

2014-03-27 22:59:25 +01:00

Wow… This sounds very promising!

Too late here to follow 100% the reasoning and test it, but if it does give us more accuracy in this mess, would be a great thing!

Will try it in the coming days, in any case many thanks for the work and patch! :D

Wow… This sounds very promising! Too late here to follow 100% the reasoning and test it, but if it does give us more accuracy in this mess, would be a great thing! Will try it in the coming days, in any case many thanks for the work and patch! :D

tippisum commented

2014-03-28 10:51:38 +01:00

(Migrated from localhost:3001)

When I reviewed it myself, I have found something more interesting and the algorithm is now further simplified.

As the bone vector is assumed to be normalized, the condition x² + y² + z² = 1 always holds, which is equivalent to x² + z² = 1 - y² = (1 - y)(1 + y).
By take this into account, we can simplify the above result even further.

(x²y + z²)/(x² + z²) = (x² + z² + x²(y - 1))/(x² + z²) = 1 + x²(y - 1)/((1 - y)(1 + y)) = 1 - x²/(1 + y)
xz(y - 1)/(x² + z²) = -xz(1 - y)/((1 - y)(1 + y)) = -xz/(1 + y)
x(y² - 1)/(x² + z²) = x(y² - 1)/(1 - y²) = -x
z(y² - 1)/(x² + z²) = z(y² - 1)/(1 - y²) = -z
(x² + z²y)/(x² + z²) = (x² + z² + z²(y - 1))/(x² + z²) = 1 + z²(y - 1)/((1 - y)(1 + y)) = 1 - z²/(1 + y)

Then we rewrite M as

    ┌ 1 - x²/(1 + y), x,    -xz/(1 + y) ┐
    │                                   │
M = │             -x, y,             -z │
    │                                   │
    └    -xz/(1 + y), z, 1 - z²/(1 + y) ┘

The final result is amazingly beautiful in mathematics, and eventually we only have 4 elements that needs an actual calculation, and their evaluations is indeed trivial and possible to be managed with both good performance and precision.
Yes, at least they must be FAR BETTER than the present one which invokes trigonometric functions that we all know cannot be accurate in its very nature.
And it even ruled out one of the singular point y = 1 which yields a unit matrix. by reducing (1 - y) factor from every elements, it is no longer singular. Now this algorithm is numerically stable when bone is pointing to +Y. The only case still needs to be treated separately is when y = -1, which is a intrinsic one from the arbitrariness of selecting rotation axis under such circumstance.

However, special caution is still needed when we are dealing with vectors near the singular point. Naive implementation may raise plenty of bugs because of its numerical instability.

Let's take a brief look about the asymptotic behavior when bone vector is reaching the limit of y = -1.

It is easy to check that each of the four corner elements is possible to vary from -1 to 1, depends on the axis chosen for doing the rotation.
And it can be figured out that, the "rotation" here is in fact established by mirroring XZ-plane by that certain axis, then inverse the Y-axis.
For sufficiently small x, z and with y approaching -1, all elements but the 4 corner ones of M will degenerate. So let's now focus on these corner elements.

We rewrite M so that it only contains its 4 corner elements, and combine the 1/(1 + y) factor.

               ┌                        ┐
        1      │ 1 + y - x²,        -xz │
M* = ------- * │                        │
      1 + y    │        -xz, 1 + y - z² │
               └                        ┘

In the circumstance that y is close to -1, mindlessly calculating the factor 1/(1 + y) will cause severe numerical instability. As a solution, we can just ignore it and normalize M instead.
Since bone vector is already normalized, y² = 1 - (x² + z²). And y is approaching -1 so we use the minus branch of square root: y = -√(1 - (x² + z²)).
As x and z are both small quantities, we apply the binomial expansion to the first order, which gives
y = -√(1 - (x² + z²)) = -1 + 1/2(x² + z²). With this substitution we then have:

                 ┌                  ┐
         1       │ z² - x²,    -2xz │
M* = --------- * │                  │
      x² + z²    │    -2xz, x² - z² │
                 └                  ┘

We may then "normalize" the rotation axis vector (which is previously represented as (z, 0, -x)), and calculate matrix this way. It will get a much more stable behavior when bone vectors are approaching -Y.

In conclusion, the new algorithm will take two threshold. As long as y is not close to -1, the default way is used, which should give a accurate enough result(and should even be faster that the present one). When y is approaching -1, the vector is re-normalized in XZ-plane(or the axis is normalized, which have same effect), and the modified method is used in order to maintain numerical stability. And when y is really close to -1, then we just treated as it is indeed anti-parallel. The actual value of these threshold may need some test and tweaks, but it does seem less error-prone than it currently does.

I am working on a patch of this new algorithm now, but since I really can't figure out how to compile Blender myself, I just can't do any real test. And the previous patch is indeed write in blind. (Yes, I even don't know whether it actually compiles or not, since I cannot try it on my own)
Can anybody give me some idea about how to compile Blender on Windows? It does not seem happy with the Visual Studio 2012 compiler.

Nonetheless, the algorithm itself should be okay. And I tested it under Blender's Python Console, where Matrix are using Blender's current implementation and computations issued in Python environment are conducted by Python using Double precision. My algorithm always give much more accurate and numerically stable result.
I think lots of people who want to handle import / export correctly will benefit from this improved algorithm (including myself). By the way, I wonder why it is now impossible to assign matrices to bones directly, since it seems possible back to 2.4x? Things can be a lot easier if this is possible...

When I reviewed it myself, I have found something more interesting and the algorithm is now further simplified. As the bone vector is assumed to be normalized, the condition `x² + y² + z² = 1` always holds, which is equivalent to `x² + z² = 1 - y² = (1 - y)(1 + y)`. By take this into account, we can simplify the above result even further. * `(x²y + z²)/(x² + z²) = (x² + z² + x²(y - 1))/(x² + z²) = 1 + x²(y - 1)/((1 - y)(1 + y)) = 1 - x²/(1 + y)` * `xz(y - 1)/(x² + z²) = -xz(1 - y)/((1 - y)(1 + y)) = -xz/(1 + y)` * `x(y² - 1)/(x² + z²) = x(y² - 1)/(1 - y²) = -x` * `z(y² - 1)/(x² + z²) = z(y² - 1)/(1 - y²) = -z` * `(x² + z²y)/(x² + z²) = (x² + z² + z²(y - 1))/(x² + z²) = 1 + z²(y - 1)/((1 - y)(1 + y)) = 1 - z²/(1 + y)` Then we rewrite `M` as ``` ┌ 1 - x²/(1 + y), x, -xz/(1 + y) ┐ │ │ M = │ -x, y, -z │ │ │ └ -xz/(1 + y), z, 1 - z²/(1 + y) ┘ ``` The final result is amazingly beautiful in mathematics, and eventually we only have 4 elements that needs an actual calculation, and their evaluations is indeed trivial and possible to be managed with both good performance and precision. **Yes, at least they must be FAR BETTER than the present one which invokes trigonometric functions that we all know cannot be accurate in its very nature.** And it even ruled out one of the singular point `y = 1` which yields a unit matrix. by reducing `(1 - y)` factor from every elements, it is no longer singular. **Now this algorithm is numerically stable when bone is pointing to +Y.** The only case still needs to be treated separately is when `y = -1`, which is a intrinsic one from the arbitrariness of selecting rotation axis under such circumstance. However, special caution is still needed when we are dealing with vectors near the singular point. Naive implementation may raise plenty of bugs because of its numerical instability. Let's take a brief look about the asymptotic behavior when bone vector is reaching the limit of `y = -1`. * It is easy to check that each of the four corner elements is possible to vary from `-1` to `1`, depends on the axis chosen for doing the rotation. * And it can be figured out that, the "rotation" here is in fact established by mirroring XZ-plane by that certain axis, then inverse the Y-axis. * For sufficiently small `x`, `z` and with `y` approaching `-1`, all elements but the 4 corner ones of `M` will degenerate. So let's now focus on these corner elements. We rewrite `M` so that it only contains its 4 corner elements, and combine the 1/(1 + y) factor. ``` ┌ ┐ 1 │ 1 + y - x², -xz │ M* = ------- * │ │ 1 + y │ -xz, 1 + y - z² │ └ ┘ ``` In the circumstance that `y` is close to `-1`, mindlessly calculating the factor `1/(1 + y)` will cause severe numerical instability. As a solution, we can just ignore it and normalize `M` instead. Since bone vector is already normalized, `y² = 1 - (x² + z²)`. And `y` is approaching `-1` so we use the minus branch of square root: `y = -√(1 - (x² + z²))`. As `x` and `z` are both small quantities, we apply the binomial expansion to the first order, which gives `y = -√(1 - (x² + z²)) = -1 + 1/2(x² + z²)`. With this substitution we then have: ``` ┌ ┐ 1 │ z² - x², -2xz │ M* = --------- * │ │ x² + z² │ -2xz, x² - z² │ └ ┘ ``` We may then "normalize" the rotation axis vector (which is previously represented as `(z, 0, -x)`), and calculate matrix this way. It will get a much more stable behavior when bone vectors are approaching -Y. In conclusion, the new algorithm will take two threshold. **As long as `y` is not close to `-1`, the default way is used, which should give a accurate enough result(and should even be faster that the present one).** When `y` is approaching `-1`, the vector is re-normalized in XZ-plane(or the axis is normalized, which have same effect), and the modified method is used in order to maintain numerical stability. And when `y` is really close to `-1`, then we just treated as it is indeed anti-parallel. The actual value of these threshold may need some test and tweaks, but it does seem less error-prone than it currently does. I am working on a patch of this new algorithm now, but since I really can't figure out how to compile Blender myself, I just can't do any real test. And the previous patch is indeed write in blind. (Yes, I even don't know whether it actually compiles or not, since I cannot try it on my own) Can anybody give me some idea about how to compile Blender on Windows? It does not seem happy with the Visual Studio 2012 compiler. Nonetheless, the algorithm itself should be okay. And I tested it under Blender's Python Console, where `Matrix` are using Blender's current implementation and computations issued in Python environment are conducted by Python using Double precision. My algorithm always give much more accurate and numerically stable result. I think lots of people who want to handle import / export correctly will benefit from this improved algorithm (including myself). By the way, **I wonder why it is now impossible to assign matrices to bones directly, since it seems possible back to 2.4x?** Things can be a lot easier if this is possible...

Bastien Montagne commented

2014-03-28 11:12:27 +01:00

Wow^2 :P

Ok, just tested previous patch this morning, had to do a few minors changes to get it working (mostly, fallback to 100% float for now, else we'd need to add more double/float conversions). From a very quick test, it seems to work quite well, though having different behavior than previous code regarding roll results, when reaching some positions… Will have to test this further, but will wait for your next patch now. :)

About building Blender under windows, http://wiki.blender.org/index.php/Dev:Doc/Building_Blender/Windows is the official ref. Note officially supported compiler remains VC2008, though we plan to switch to VC2013 in coming weeks. Not sure about VC2012 (since we have to compile a whole set of libs for each versions, we probably won't support it ever :/ ).

Wow^2 :P Ok, just tested previous patch this morning, had to do a few minors changes to get it working (mostly, fallback to 100% float for now, else we'd need to add more double/float conversions). From a very quick test, it seems to work quite well, though having different behavior than previous code regarding roll results, when reaching some positions… Will have to test this further, but will wait for your next patch now. :) About building Blender under windows, http://wiki.blender.org/index.php/Dev:Doc/Building_Blender/Windows is the official ref. Note officially supported compiler remains VC2008, though we plan to switch to VC2013 in coming weeks. Not sure about VC2012 (since we have to compile a whole set of libs for each versions, we probably won't support it ever :/ ).

tippisum commented

2014-03-28 12:53:42 +01:00

(Migrated from localhost:3001)

Here comes the new patch...

Now it seems that coding this algorithm is actually pretty easy...
but people will probably never know what I was doing without first seeing this ... LOL

Chance are big that it will deliver adequate accuracy even with Single precision.
(I just personally prefer double ...)

Any hope for the Python API to support assign bone matrices directly?
armature.diff

Here comes the new patch... Now it seems that coding this algorithm is actually pretty easy... but people will probably never know what I was doing without first seeing this ... LOL Chance are big that it will deliver adequate accuracy even with Single precision. (I just personally prefer double ...) Any hope for the Python API to support assign bone matrices directly? [armature.diff](https://archive.blender.org/developer/F83190/armature.diff)

Bastien Montagne commented

2014-03-30 17:33:16 +02:00

Ok, tested the code, seems to work rather well. Also made a nice diff for it. :)

@tippisum I have a question, though. When we come close to -Y, is it really better for stability to compute 1 / (x² + z²), rather than 1 / (1 + y)? Does not seems obvious to me…

Ok, tested the code, seems to work rather well. Also made a nice diff for it. :) @tippisum I have a question, though. When we come close to -Y, is it really better for stability to compute `1 / (x² + z²)`, rather than `1 / (1 + y)`? Does not seems obvious to me…

Karja Krähwald commented

2014-03-30 19:51:05 +02:00

Added subscriber: @karja

Karja Krähwald commented

2014-03-30 19:51:05 +02:00

Is this algorithm also applicable on #39326 somehow, or bone-only? It looks related.

tippisum commented

2014-03-30 19:55:39 +02:00

(Migrated from localhost:3001)

Yes, it should be.
Because in the previous discussion, we just assumed that the bone vector is already normalized, but in actual use case, We will need to normalize it.

When vector is approaching -Y, calculating 1/(1 + y) after normalizing will suffer severely from round-off error, just as the same reason that we should never try to subtract two entity that we already know they are close.

Imagine this situation: The vector is very close to -Y, so that x and z are both extremely small compared to y, and y will become -0.99999994 or something after normalizing. (or you will simply get a -1 because we are talking about Single precision). Evaluating 1 + y will obviously destroy every significant digit we have in the original vector or even give us unacceptable behavior. By evaluating x² + z² (which is sufficiently close to 1/2(1 + y) when y -> -1) instead, we can avoid this problem. Since the deviation of binomial expansion is in a higher order of (1 + y), using 1.0e-5 can be a reasonable threshold because the deviation will get rounded under Single precision (which hols 9 significant digits at most).

However, it should be clarified that this patch is not intended to fix #23954 or #27675. It should have some effect on them, I suppose (benefit from the higher accuracy vector->matrix converting we have now). But the root cause remains the same so do not put too much expect on it.

Nonetheless, I certainly have some models whose analyzing results suggested similar pattern of flaw in bone matrices as Blender used to. And they likely come from 3ds max (I'm not 100% sure). So it is possible that the bone matrix calculation in 3ds max suffers from numerical instability too. And with this patch equipped, we'll get the leading accuracy in bone matrix calculation (since I can't see what algorithm can outperform mine regarding precision). Enjoy it. :)

Yes, it should be. Because in the previous discussion, we just **assumed** that the bone vector is already normalized, but in actual use case, We will need to normalize it. When vector is approaching -Y, calculating `1/(1 + y)` after normalizing will suffer severely from round-off error, just as the same reason that we should never try to subtract two entity that we already know they are close. Imagine this situation: The vector is very close to -Y, so that `x` and `z` are both extremely small compared to `y`, and `y` will become `-0.99999994` or something after normalizing. (or you will simply get a `-1` because we are talking about **Single precision**). Evaluating `1 + y` will obviously destroy every significant digit we have in the original vector or even give us unacceptable behavior. By evaluating `x² + z²` (which is sufficiently close to `1/2(1 + y)` when `y -> -1`) instead, we can avoid this problem. Since the deviation of binomial expansion is in a higher order of `(1 + y)`, using `1.0e-5` can be a reasonable threshold because the deviation will get rounded under Single precision (which hols 9 significant digits at most). However, it should be clarified that this patch is not intended to fix #23954 or #27675. It should have some effect on them, I suppose (benefit from the higher accuracy vector->matrix converting we have now). But the root cause remains the same so do not put too much expect on it. Nonetheless, I certainly have some models whose analyzing results suggested similar pattern of flaw in bone matrices as Blender used to. And they likely come from 3ds max (I'm not 100% sure). So it is possible that the bone matrix calculation in 3ds max suffers from numerical instability too. And with this patch equipped, we'll get the leading accuracy in bone matrix calculation (since I can't see what algorithm can outperform mine regarding precision). Enjoy it. :)

tippisum commented

2014-03-30 20:04:29 +02:00

(Migrated from localhost:3001)

@karja

In a briefly view, I think the opportunity is not high.
#39326 comes from the float point deviation which we don't have much to do with.

Yes, it is possible to get higher precision calculating ordinary rotation matrix, but it will need a lot of work.
Seriously, if you guys surely want some accuracy boost without doing much coding or mathematical analysis, using double instead of float whenever possible is the quickest way.

@karja In a briefly view, I think the opportunity is not high. #39326 comes from the float point deviation which we don't have much to do with. Yes, it is possible to get higher precision calculating ordinary rotation matrix, but it will need a lot of work. Seriously, if you guys surely want some accuracy boost without doing much coding or mathematical analysis, using double instead of float whenever possible is the quickest way.

Bastien Montagne commented

2014-03-30 20:38:12 +02:00

@tippisum, thanks for clarification, makes sense indeed. :)

Bastien Montagne commented

2014-04-06 19:25:55 +02:00

Changed status from 'Open' to: 'Resolved'

Bastien Montagne closed this issue

2014-04-06 19:25:55 +02:00

Bastien Montagne commented

2014-04-06 19:25:55 +02:00

Applied to master as 07f8c5c3b6.

Many thanks again for the work & patch!

Applied to master as 07f8c5c3b6. Many thanks again for the work & patch!

Sign in to join this conversation.

No Label

Download

What's New

Blender Studio

Manual

Developers Blog

Documentation

Benchmark

Blender Conference

Development Fund

One-time Donations

New algorithm on computing bone matrix with improved accuracy #39470