Problem with precision in Euler Angles.

[DJuego]

I am working with Microsoft Visual Studio 2019 16.6.1 in Windows 10 x64 environment.

Well. There seems to be some problem in the numerical accuracy of Euler's new angle conversion feature.

For instance; When I convert, for example, roll=90º ,pitch=90º, yaw=0º to rotor, I get in rotor: roll=0º, pitch=89.972º, yaw=0º.I find a noticeable difference between 90º and 89,972º (0.028º) that can build up too quickly.

#include <klein/public/klein/klein.hpp>

#include <cmath>
#include <iostream>

const float M_PI = 3.141592653589793238462643383279502884; /* pi */

int main()
{
	float roll_1 = 90.0;
	float pitch_1 = 90.0;
	float yaw_1 = 0.0;

	std::cout << "INPUT: roll: " << roll_1 << " pitch: " << pitch_1 << " yaw: " << yaw_1 << std::endl;

	kln::euler_angles ea_1;
	ea_1.roll = roll_1 * (M_PI / 180.0f);
	ea_1.yaw = yaw_1 * (M_PI / 180.0f);
	ea_1.pitch = pitch_1 * (M_PI / 180.0f);

	kln::rotor r;
	r = kln::rotor(ea_1);
	kln::euler_angles ea_2 = r.as_euler_angles();
	
	float roll_2 = ea_2.roll * (180.0f / M_PI);
	float yaw_2 = ea_2.yaw * (180.0f / M_PI);
	float pitch_2 = ea_2.pitch * (180.0f / M_PI);

	std::cout << "OUTPUT: roll: " << roll_2 << " pitch: " << pitch_2 << " yaw: " << yaw_2 << std::endl;
	return 0;
}

I get:

INPUT: roll: 90 pitch: 90 yaw: 0
OUTPUT: roll: 0 pitch: 89.972 yaw: 0

DJuego

[Orbots]

Try defining a more exact constant for 180.0/pi rather than using M_PI. You will improve by about 4e-6 each time you convert to radians, and then subsequently propagate/compound that error through motor arithmetic.

julia> Float32(180.0/π) - 180.0f0/Float32(π)
3.8146973f-6

[jeremyong]

Can you try angles away from the 90 degree poles out of curiosity? Euler angles have an endemic precision issue around there so I'm wondering if there's some case I need to handle with branching.

[DJuego]

Thank you for your answers!

@Orbots It doesn't seem to be a precision problem with PI. I'm using 36 decimals. Now on a float, but earlier on a double. Plus I'm not doing any iterative operations that would accumulate error (only two radians<->degrees operations). .And yet the error is significant. What values would you use?

@jeremyong. Yes, indeed!.
Any samples;

INPUT: roll: 90 pitch: 90 yaw: 0
OUTPUT: roll: 0 pitch: 89.972 yaw: 0

INPUT: roll: 89 pitch: 89 yaw: 0
OUTPUT: roll: 89 pitch: 89.0002 yaw: -9.78401e-05

INPUT: roll: 88 pitch: 88 yaw: 0
OUTPUT: roll: 88 pitch: 88.0001 yaw: -4.89277e-05

INPUT: roll: 87 pitch: 87 yaw: 0
OUTPUT: roll: 87.0001 pitch: 87 yaw: 3.26266e-05

INPUT: roll: 86 pitch: 86 yaw: 0
OUTPUT: roll: 86.0001 pitch: 86.0001 yaw: 0

INPUT: roll: 85 pitch: 85 yaw: 0
OUTPUT: roll: 84.9999 pitch: 84.9999 yaw: 0

INPUT: roll: 84 pitch: 84 yaw: 0
OUTPUT: roll: 84 pitch: 84 yaw: 0

---

INPUT: roll: 180 pitch: 180 yaw: 0
OUTPUT: roll: 5.00896e-06 pitch: -5.00896e-06 yaw: 180

INPUT: roll: 179 pitch: 179 yaw: 0
OUTPUT: roll: -1 pitch: 1 yaw: 180

---

INPUT: roll: 270 pitch: 270 yaw: 0
OUTPUT: roll: 0 pitch: -89.972 yaw: 0

---

INPUT: roll: 360 pitch: 360 yaw: 0
OUTPUT: roll: 1.00179e-05 pitch: 1.00179e-05 yaw: 0

DJuego

[Orbots]

Your value of PI only is only precise out to 7 or 8 decimal places.

julia> Float32(3.141592653589793238462643383279502884) - 3.141592653589793238462643383279502884
8.742278012618954e-8

I'm assuming the error is accumulating in

        r = kln::rotor(ea_1);
	kln::euler_angles ea_2 = r.as_euler_angles();

Cheers.

[palazzol]

Hello,

I decided to step through this.

It seems like the issue is in creating the rotor. The calculation before the normalize() call sets the internal quaternion coefficients to:
(0.499999970, 0.499999970, 0.499999970,-0.49999970), but then the normalize() call is unable to adjust the result to (0.5,0.5,0.5,-0.5)

I'm not sure if this is a limitation of float precision, or the normalization technique.

-Frank

[DJuego]

Hello again! Long time no see! Any news regarding this problem? 😔

DJuego

[dobacetr]

Hello,

I decided to step through this.

It seems like the issue is in creating the rotor. The calculation before the normalize() call sets the internal quaternion coefficients to: (0.499999970, 0.499999970, 0.499999970,-0.49999970), but then the normalize() call is unable to adjust the result to (0.5,0.5,0.5,-0.5)

I'm not sure if this is a limitation of float precision, or the normalization technique.

-Frank

You are correct. The problem occurs when computing 1/sqrt(0.999999881) = 1.0000000595, which can not be represented in single point precision.
Consider using either higher precision or seperating reciprocal and sqrt.