New user self-registration is currently disabled. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Bug 1312 - AngleAxis from QuaternionBase , better precision
 Summary: AngleAxis from QuaternionBase , better precision
 Status: RESOLVED FIXED Eigen Unclassified Geometry 3.2 All All Normal Accuracy Problem Nobody 1244 Show dependency tree / graph

 Reported: 2016-09-26 09:21 UTC by Michael Norel 2016-09-30 08:00 UTC (History) 4 users (show) chtz gael.guennebaud hauke.heibel jacob.benoit.1

Attachments

 Michael Norel 2016-09-26 09:21:27 UTC ```1. m_angle = Scalar(2)*acos((min)((max)(Scalar(-1),q.w()),Scalar(1))); Angle can be computed with better precision using atan2(sqrt(n2), q.w()), without check for "1" overflow. 2. Also, for better precision , it is good to switch sign of input quaternion to have positive q.w(), it keeps rotation angle in [0, PI]. 3. The check "if (n2 < NumTraits::dummy_precision()*NumTraits::dummy_precision())" can be changed to strong comparison with "zero", and provide same relative precision for small rotations without side effects.``` Michael Norel 2016-09-28 13:49:52 UTC `Sorry , this is dublicate, i posted this year ago. But still , not fixed in 3.xx.` Gael Guennebaud 2016-09-29 11:46:54 UTC ```hm, the current implementation is: using std::atan2; Scalar n = q.vec().norm(); if(n::epsilon()) n = q.vec().stableNorm(); if (n > Scalar(0)) { m_angle = Scalar(2)*atan2(n, q.w()); m_axis = q.vec() / n; } else { m_angle = Scalar(0); m_axis << Scalar(1), Scalar(0), Scalar(0); } return *this; So 1. and 3. were already "fixed". Regarding 2, I would prefer to keep it like this so that angle-axis -> quat -> angle-axis returns exactly the same axis.``` Michael Norel 2016-09-29 13:29:06 UTC ```Yes, my mistake , i checked "Eigen 3.2.9" instead of current "Eigen 3.3-rc1". Thank for correction, ill try to avoid such mistakes. About 2) "that angle-axis -> quat -> angle-axis returns exactly the same axis." If this is not important issue, a can demonstrate the loose of accuracy. double small_angle = 1.0e-16; // any value smaller epsilon double a1 = atan2(small_angle, -1.0);// == PI - small_angle; double a2 = atan2(small_angle, 1.0);// == small_angle; double a1Check = acos(-1.0) - a1;// a1Check is zero here :( For all angles less then epsilon * PI, rotation angle will be lost. So, current conversion: - loose absolute precision for small angles. - breaks relative precision for small angles. fix can be simple ... Scalar w = q.w(); if(w < 0) { w=-w; n=-n; } if (n != Scalar(0)) ... but introduce additional branch. Here, i prefer precision on full range of angles over performance. Anyway, it is not a big issue.``` Gael Guennebaud 2016-09-29 20:37:30 UTC `Sorry, you're right, and in my angleaxis-quat->angleaxis scenario, w will be positive anyway (unless the angle is outside [-2pi,2pi])` Gael Guennebaud 2016-09-29 20:44:20 UTC ```(In reply to Gael Guennebaud from comment #4) > (unless the angle is outside [-2pi,2pi]) I meant [-pi,pi].``` Gael Guennebaud 2016-09-29 21:24:22 UTC `Fixed: https://bitbucket.org/eigen/eigen/commits/5998f71df011` Michael Norel 2016-09-30 08:00:22 UTC ```As i see from patch, code looks template template AngleAxis& AngleAxis::operator=(const QuaternionBase& q) { using std::atan2; using std::abs; Scalar n = q.vec().norm(); if(n::epsilon()) n = q.vec().stableNorm(); if (n != Scalar(0)) { m_angle = Scalar(2)*atan2(n, std::abs(q.w())); if(q.w() < 0) n = -n; m_axis = q.vec() / n; } else { m_angle = Scalar(0); m_axis << Scalar(1), Scalar(0), Scalar(0); } return *this; } 1. "n" is positive here , and (n > Scalar(0)) still can be used 2. m_angle = Scalar(2)*atan2(n, std::abs(q.w())); returns value in [0,pi], so this comment should be changed "angle is in the [-pi,pi] range."```

 Note You need to log in before you can comment on or make changes to this bug.