diff options
author | Gael Guennebaud <g.gael@free.fr> | 2013-06-09 23:14:45 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2013-06-09 23:14:45 +0200 |
commit | c98fd7a6cae853f1ca8570994ae9ba3c13e9c4bd (patch) | |
tree | 6280d9b1bb53c14e203f066be5c56854bf82c025 | |
parent | e04b59929e45564dbbc80f4e1482892458983ac0 (diff) |
Fix bug #609: avoid if statement and improve consistency of eulerAngles method
-rw-r--r-- | Eigen/src/Geometry/EulerAngles.h | 58 | ||||
-rw-r--r-- | test/geo_eulerangles.cpp | 58 |
2 files changed, 85 insertions, 31 deletions
diff --git a/Eigen/src/Geometry/EulerAngles.h b/Eigen/src/Geometry/EulerAngles.h index 216307706..3f6ecc6d9 100644 --- a/Eigen/src/Geometry/EulerAngles.h +++ b/Eigen/src/Geometry/EulerAngles.h @@ -27,12 +27,18 @@ namespace Eigen { * * AngleAxisf(ea[1], Vector3f::UnitX()) * * AngleAxisf(ea[2], Vector3f::UnitZ()); \endcode * This corresponds to the right-multiply conventions (with right hand side frames). + * + * The returned angles are in the ranges [0:pi]x[0:pi]x[-pi:pi]. + * + * \sa class AngleAxis */ template<typename Derived> inline Matrix<typename MatrixBase<Derived>::Scalar,3,1> MatrixBase<Derived>::eulerAngles(Index a0, Index a1, Index a2) const { using std::atan2; + using std::sin; + using std::cos; /* Implemented from Graphics Gems IV */ EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Derived,3,3) @@ -44,39 +50,53 @@ MatrixBase<Derived>::eulerAngles(Index a0, Index a1, Index a2) const const Index i = a0; const Index j = (a0 + 1 + odd)%3; const Index k = (a0 + 2 - odd)%3; - + if (a0==a2) { - Scalar s = Vector2(coeff(j,i) , coeff(k,i)).norm(); - res[1] = atan2(s, coeff(i,i)); - if (s > epsilon) + res[0] = atan2(coeff(j,i), coeff(k,i)); + if((odd && res[0]<0) || ((!odd) && res[0]>0)) { - res[0] = atan2(coeff(j,i), coeff(k,i)); - res[2] = atan2(coeff(i,j),-coeff(i,k)); + res[0] = (res[0] > 0) ? res[0] - M_PI : res[0] + M_PI; + Scalar s2 = Vector2(coeff(j,i), coeff(k,i)).norm(); + res[1] = -atan2(s2, coeff(i,i)); } else { - res[0] = Scalar(0); - res[2] = (coeff(i,i)>0?1:-1)*atan2(-coeff(k,j), coeff(j,j)); + Scalar s2 = Vector2(coeff(j,i), coeff(k,i)).norm(); + res[1] = atan2(s2, coeff(i,i)); } - } + + // With a=(0,1,0), we have i=0; j=1; k=2, and after computing the first two angles, + // we can compute their respective rotation, and apply its inverse to M. Since the result must + // be a rotation around x, we have: + // + // c2 s1.s2 c1.s2 1 0 0 + // 0 c1 -s1 * M = 0 c3 s3 + // -s2 s1.c2 c1.c2 0 -s3 c3 + // + // Thus: m11.c1 - m21.s1 = c3 & m12.c1 - m22.s1 = s3 + + Scalar s1 = sin(res[0]); + Scalar c1 = cos(res[0]); + res[2] = atan2(c1*coeff(j,k)-s1*coeff(k,k), c1*coeff(j,j) - s1 * coeff(k,j)); + } else { - Scalar c = Vector2(coeff(i,i) , coeff(i,j)).norm(); - res[1] = atan2(-coeff(i,k), c); - if (c > epsilon) - { - res[0] = atan2(coeff(j,k), coeff(k,k)); - res[2] = atan2(coeff(i,j), coeff(i,i)); + res[0] = atan2(coeff(j,k), coeff(k,k)); + Scalar c2 = Vector2(coeff(i,i), coeff(i,j)).norm(); + if((odd && res[0]<0) || ((!odd) && res[0]>0)) { + res[0] = (res[0] > 0) ? res[0] - M_PI : res[0] + M_PI; + res[1] = atan2(-coeff(i,k), -c2); } else - { - res[0] = Scalar(0); - res[2] = (coeff(i,k)>0?1:-1)*atan2(-coeff(k,j), coeff(j,j)); - } + res[1] = atan2(-coeff(i,k), c2); + Scalar s1 = sin(res[0]); + Scalar c1 = cos(res[0]); + res[2] = atan2(s1*coeff(k,i)-c1*coeff(j,i), c1*coeff(j,j) - s1 * coeff(k,j)); } if (!odd) res = -res; + return res; } diff --git a/test/geo_eulerangles.cpp b/test/geo_eulerangles.cpp index 9bf149d2a..26456beee 100644 --- a/test/geo_eulerangles.cpp +++ b/test/geo_eulerangles.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -12,22 +12,18 @@ #include <Eigen/LU> #include <Eigen/SVD> -template<typename Scalar> void eulerangles(void) +template<typename Scalar> void check_all_var(const Matrix<Scalar,3,1>& ea) { typedef Matrix<Scalar,3,3> Matrix3; typedef Matrix<Scalar,3,1> Vector3; - typedef Quaternion<Scalar> Quaternionx; typedef AngleAxis<Scalar> AngleAxisx; - - Scalar a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); - Quaternionx q1; - q1 = AngleAxisx(a, Vector3::Random().normalized()); - Matrix3 m; - m = q1; - + #define VERIFY_EULER(I,J,K, X,Y,Z) { \ - Vector3 ea = m.eulerAngles(I,J,K); \ - VERIFY_IS_APPROX(m, Matrix3(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z()))); \ + Matrix3 m(AngleAxisx(ea[0], Vector3::Unit##X()) * AngleAxisx(ea[1], Vector3::Unit##Y()) * AngleAxisx(ea[2], Vector3::Unit##Z())); \ + Vector3 eabis = m.eulerAngles(I,J,K); \ + Matrix3 mbis(AngleAxisx(eabis[0], Vector3::Unit##X()) * AngleAxisx(eabis[1], Vector3::Unit##Y()) * AngleAxisx(eabis[2], Vector3::Unit##Z())); \ + VERIFY_IS_APPROX(m, mbis); \ + if(I!=K || ea[1]!=0) VERIFY_IS_APPROX(ea, eabis); \ } VERIFY_EULER(0,1,2, X,Y,Z); VERIFY_EULER(0,1,0, X,Y,X); @@ -45,6 +41,44 @@ template<typename Scalar> void eulerangles(void) VERIFY_EULER(2,1,2, Z,Y,Z); } +template<typename Scalar> void eulerangles(void) +{ + typedef Matrix<Scalar,3,3> Matrix3; + typedef Matrix<Scalar,3,1> Vector3; + typedef Array<Scalar,3,1> Array3; + typedef Quaternion<Scalar> Quaternionx; + typedef AngleAxis<Scalar> AngleAxisx; + + Scalar a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI)); + Quaternionx q1; + q1 = AngleAxisx(a, Vector3::Random().normalized()); + Matrix3 m; + m = q1; + + Vector3 ea = m.eulerAngles(0,1,2); + check_all_var(ea); + ea = m.eulerAngles(0,1,0); + check_all_var(ea); + + ea = (Array3::Random() + Array3(1,1,0))*M_PI*Array3(0.5,0.5,1); + check_all_var(ea); + + ea[2] = ea[0] = internal::random<Scalar>(0,M_PI); + check_all_var(ea); + + ea[0] = ea[1] = internal::random<Scalar>(0,M_PI); + check_all_var(ea); + + ea[1] = 0; + check_all_var(ea); + + ea.head(2).setZero(); + check_all_var(ea); + + ea.setZero(); + check_all_var(ea); +} + void test_geo_eulerangles() { for(int i = 0; i < g_repeat; i++) { |