This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

Bug 1312

Summary: AngleAxis from QuaternionBase , better precision
Product: Eigen Reporter: Michael Norel <minorlogic>
Component: GeometryAssignee: Nobody <eigen.nobody>
Status: RESOLVED FIXED    
Severity: Accuracy Problem CC: chtz, gael.guennebaud, hauke.heibel, jacob.benoit.1
Priority: Normal    
Version: 3.2   
Hardware: All   
OS: All   
Whiteboard:
Bug Depends on:    
Bug Blocks: 1244    

Description 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<Scalar>::dummy_precision()*NumTraits<Scalar>::dummy_precision())"
can be changed to strong comparison with "zero", and provide same relative precision for small rotations without side effects.
Comment 1 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.
Comment 2 Gael Guennebaud 2016-09-29 11:46:54 UTC
hm, the current implementation is:

  using std::atan2;
  Scalar n = q.vec().norm();
  if(n<NumTraits<Scalar>::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.
Comment 3 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.
Comment 4 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])
Comment 5 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].
Comment 6 Gael Guennebaud 2016-09-29 21:24:22 UTC
Fixed: https://bitbucket.org/eigen/eigen/commits/5998f71df011
Comment 7 Michael Norel 2016-09-30 08:00:22 UTC
As i see from patch, code looks 

template<typename Scalar>
template<typename QuatDerived>
AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived>& q)
{
  using std::atan2;
  using std::abs;
  Scalar n = q.vec().norm();
  if(n<NumTraits<Scalar>::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."
Comment 8 Nobody 2019-12-04 16:21:26 UTC
-- GitLab Migration Automatic Message --

This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.

You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/1312.