This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
Bug 1312 - AngleAxis from QuaternionBase , better precision
Summary: AngleAxis from QuaternionBase , better precision
Status: RESOLVED FIXED
Alias: None
Product: Eigen
Classification: Unclassified
Component: Geometry (show other bugs)
Version: 3.2
Hardware: All All
: Normal Accuracy Problem
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks: 1244
  Show dependency treegraph
 
Reported: 2016-09-26 09:21 UTC by Michael Norel
Modified: 2019-12-04 16:21 UTC (History)
4 users (show)



Attachments

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.

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