Euler: Fix gimbal lock precision

[TheophileMot]

When setting Euler angles from a rotation matrix, imprecise calculations cause observable glitches near gimbal lock.

This happens because the cutoff threshold 0.99999 is not close enough to 1. Note that arcsin(0.99999) ≃ 89.74°, so there is a range of half a degree (i.e., a quarter of a degree on either side of 90°) that is problematic.

---

Example demonstrating the problem: https://jsfiddle.net/theophilemot/yhxte5vd/12/

Every frame, with probability 0.5, the object's Euler angles are converted to a quaternion and back; if the calculation is done properly, this should have no visual consequences. Near gimbal lock, the object will jitter.

[WestLangley]

When setting Euler angles from a rotation matrix...

Set object.quaternion from the rotation matrix, instead.

This happens because the cutoff threshold 0.99999 is not close enough to 1.

You have changed the threshold to 1.0000. I'm not sure of the consequences of doing that.

Let's first see if setting the quaternion solves your use case before pursuing this further.

[Mugen87]

This happens because the cutoff threshold 0.99999 is not close enough to 1.

The book 3D Math Primer for Graphics and Game Development (2nd Ed) suggest the following C++ code for converting a matrix to euler angles (chapter 8.7.2).

// Assume the matrix is stored in these variables:
float m11,m12,m13; 
float m21,m22,m23; 
float m31,m32,m33;

// We will compute the Euler angle values in radians and store them here:
float h,p,b;

// Extract pitch from m32, being careful for domain errors with
// asin (). We could have values slightly out of range due to
// floating point arithmetic . 

float sp = −m32;

if (sp <= −1.0f) {
	p = −1.570796f; // −pi/2
} else if (sp >= 1.0f) {
	p = 1.570796f; // pi/2
}else{
	p = asin(sp);
}

// Check for the Gimbal lock case, giving a slight tolerance 
// for numerical imprecision
if (fabs(sp) > 0.9999f) {
	// We are looking straight up or down.
	// Slam bank to zero and just set heading b = 0.0f;
	h = atan2(−m13, m11);
}else{
	// Compute heading from m13 and m33
	h = atan2(m31, m33);
	// Compute bank from m21 and m22
	b = atan2(m12, m22);
}

They assume head/pitch/bank order which is essentially Y/X/Z. Also you have to transpose the matrix so it matches to three.js.

In any event, they use a similar tolerance value 0.9999 compared to three.js. Probably for good reasons. I think we should not deviate from this reference if the consequences of a value of 1 are not clear.

[TheophileMot]

@WestLangley @mrdoob
Thank you for your suggestion. In many cases, I agree that it would be good for people to work around the issue. The fact remains, however, that the calculations are broken near gimbal lock.

I've put together another demo, a unit test of sorts, to show which angles are incorrectly calculated. The same demo shows that the proposed changes are accurate. I would be curious to know your thoughts on this.

Demo: https://jsfiddle.net/theophilemot/7xfrz9je/

[bhouston]

@Mugen87 wrote:

In any event, they use a similar tolerance value 0.9999 compared to three.js. Probably for good reasons. I think we should not deviate from this reference if the consequences of a value of 1 are not clear.

You are aware that the code you are referencing is using float32 where we operate in doubles. Thus at least one should increase it from 0.9999f to something appropriate for doubles like 0.99999999.

I asked @TheophileMot to contribute this to Three.JS as we have made the same change in our private Three.JS to fix a real bug. I do not mind if it isn't merged in but we are doing this out of good intentions.

[Mugen87]

0.99999999

That seems a better value than 1. I'm fine with this change. And your explanation makes of course sense. Thanks!

[TheophileMot]

@Mugen87
Another factor to consider is that the expansion of arcsin(x) near x=1 is

arcsin(1 - x) = pi/2 - sqrt(2x) + ...

The presence of the square root means that for each digit of precision we want in the output, we need twice as much precision in the threshold. In other words, changing the threshold to 0.9999999 would cut the error by ten. (This would be acceptable, in my opinion, as an alternative to 1.)

[Mugen87]

Let's try this. Feel free to reopen the PR 😊

[WestLangley]

@bhouston @TheophileMot Increasing the constant to a value other than 1 is fine with me. But as I said, it is better to set the object quaternion in your app instead.

[bhouston]

@WestLangley You are completely right about the quaternions. The issue is Clara.io, like Max and Maya only saves Euler angles in its user accessible transforms.

[TheophileMot]

@WestLangley Sure, it sounds like we have a compromise on 0.9999999. That way the error will only occur for angles within ±0.026° of the pole instead of ±0.26°, which should eliminate virtually all noticeable glitches. I've updated the PR.

[mrdoob]

Thanks!