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
AngleAxis from QuaternionBase , better precision
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Geometry
3.2
All All
: Normal Accuracy Problem
Assigned To: Nobody
:
Depends on:
Blocks: 1244
  Show dependency treegraph
 
Reported: 2016-09-26 09:21 UTC by Michael Norel
Modified: 2016-09-30 08:00 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."

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