devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merged eigen/eigen into default

Browse Source
This commit is contained in:
Benoit Steiner
2016-11-26 11:28:25 -08:00
parent 7318daf887 67b2c41f30
commit 9f8fbd9434
17 changed files with 551 additions and 477 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
unsupported/Eigen/CXX11/src/Tensor/README.md
View File

@@ -1737,11 +1737,9 @@ TODO
## Representation of scalar values
Scalar values are often represented by tensors of size 1 and rank 1. It would be
more logical and user friendly to use tensors of rank 0 instead. For example
Tensor<T, N>::maximum() currently returns a Tensor<T, 1>. Similarly, the inner
product of 2 1d tensors (through contractions) returns a 1d tensor. In the
future these operations might be updated to return 0d tensors instead.
Scalar values are often represented by tensors of size 1 and rank 0.For example
Tensor<T, N>::maximum() currently returns a Tensor<T, 0>. Similarly, the inner
product of 2 1d tensors (through contractions) returns a 0d tensor.
## Limitations

2
unsupported/Eigen/CXX11/src/Tensor/TensorDimensions.h
View File

@@ -33,7 +33,7 @@ namespace Eigen {
namespace internal {
template<std::size_t n, typename Dimension> struct dget {
static const std::size_t value = get<n, Dimension>::value;
static const std::ptrdiff_t value = get<n, Dimension>::value;
};

4
unsupported/Eigen/CXX11/src/util/CXX11Meta.h
View File

@@ -123,6 +123,10 @@ template<typename a, typename... as> struct get<0, type_lis
template<typename T, int n, T a, T... as> struct get<n, numeric_list<T, a, as...>> : get<n-1, numeric_list<T, as...>> {};
template<typename T, T a, T... as> struct get<0, numeric_list<T, a, as...>> { constexpr static T value = a; };
template<std::size_t n, typename T, T a, T... as> constexpr inline const T array_get(const numeric_list<T, a, as...>& l) {
return get<(int)n, numeric_list<T, a, as...>>::value;
}
/* always get type, regardless of dummy; good for parameter pack expansion */
template<typename T, T dummy, typename t> struct id_numeric { typedef t type; };

257
unsupported/Eigen/src/EulerAngles/EulerAngles.h
View File

@@ -12,11 +12,6 @@
namespace Eigen
{
/*template<typename Other,
int OtherRows=Other::RowsAtCompileTime,
int OtherCols=Other::ColsAtCompileTime>
struct ei_eulerangles_assign_impl;*/
/** \class EulerAngles
*
* \ingroup EulerAngles_Module
@@ -36,7 +31,7 @@ namespace Eigen
* ### Rotation representation and conversions ###
*
* It has been proved(see Wikipedia link below) that every rotation can be represented
* by Euler angles, but there is no singular representation (e.g. unlike rotation matrices).
* by Euler angles, but there is no single representation (e.g. unlike rotation matrices).
* Therefore, you can convert from Eigen rotation and to them
* (including rotation matrices, which is not called "rotations" by Eigen design).
*
@@ -55,33 +50,27 @@ namespace Eigen
* Additionally, some axes related computation is done in compile time.
*
* #### Euler angles ranges in conversions ####
* Rotations representation as EulerAngles are not single (unlike matrices),
* and even have infinite EulerAngles representations.<BR>
* For example, add or subtract 2*PI from either angle of EulerAngles
* and you'll get the same rotation.
* This is the general reason for infinite representation,
* but it's not the only general reason for not having a single representation.
*
* When converting some rotation to Euler angles, there are some ways you can guarantee
* the Euler angles ranges.
* When converting rotation to EulerAngles, this class convert it to specific ranges
* When converting some rotation to EulerAngles, the rules for ranges are as follow:
* - If the rotation we converting from is an EulerAngles
* (even when it represented as RotationBase explicitly), angles ranges are __undefined__.
* - otherwise, alpha and gamma angles will be in the range [-PI, PI].<BR>
* As for Beta angle:
* - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2].
* - otherwise:
* - If the beta axis is positive, the beta angle will be in the range [0, PI]
* - If the beta axis is negative, the beta angle will be in the range [-PI, 0]
*
* #### implicit ranges ####
* When using implicit ranges, all angles are guarantee to be in the range [-PI, +PI],
* unless you convert from some other Euler angles.
* In this case, the range is __undefined__ (might be even less than -PI or greater than +2*PI).
* \sa EulerAngles(const MatrixBase<Derived>&)
* \sa EulerAngles(const RotationBase<Derived, 3>&)
*
* #### explicit ranges ####
* When using explicit ranges, all angles are guarantee to be in the range you choose.
* In the range Boolean parameter, you're been ask whether you prefer the positive range or not:
* - _true_ - force the range between [0, +2*PI]
* - _false_ - force the range between [-PI, +PI]
*
* ##### compile time ranges #####
* This is when you have compile time ranges and you prefer to
* use template parameter. (e.g. for performance)
* \sa FromRotation()
*
* ##### run-time time ranges #####
* Run-time ranges are also supported.
* \sa EulerAngles(const MatrixBase<Derived>&, bool, bool, bool)
* \sa EulerAngles(const RotationBase<Derived, 3>&, bool, bool, bool)
*
* ### Convenient user typedefs ###
*
* Convenient typedefs for EulerAngles exist for float and double scalar,
@@ -103,7 +92,7 @@ namespace Eigen
*
* More information about Euler angles: https://en.wikipedia.org/wiki/Euler_angles
*
* \tparam _Scalar the scalar type, i.e., the type of the angles.
* \tparam _Scalar the scalar type, i.e. the type of the angles.
*
* \tparam _System the EulerSystem to use, which represents the axes of rotation.
*/
@@ -111,8 +100,11 @@ namespace Eigen
class EulerAngles : public RotationBase<EulerAngles<_Scalar, _System>, 3>
{
public:
typedef RotationBase<EulerAngles<_Scalar, _System>, 3> Base;
/** the scalar type of the angles */
typedef _Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
/** the EulerSystem to use, which represents the axes of rotation. */
typedef _System System;
@@ -146,67 +138,56 @@ namespace Eigen
public:
/** Default constructor without initialization. */
EulerAngles() {}
/** Constructs and initialize Euler angles(\p alpha, \p beta, \p gamma). */
/** Constructs and initialize an EulerAngles (\p alpha, \p beta, \p gamma). */
EulerAngles(const Scalar& alpha, const Scalar& beta, const Scalar& gamma) :
m_angles(alpha, beta, gamma) {}
/** Constructs and initialize Euler angles from a 3x3 rotation matrix \p m.
*
* \note All angles will be in the range [-PI, PI].
*/
template<typename Derived>
EulerAngles(const MatrixBase<Derived>& m) { *this = m; }
// TODO: Test this constructor
/** Constructs and initialize an EulerAngles from the array data {alpha, beta, gamma} */
explicit EulerAngles(const Scalar* data) : m_angles(data) {}
/** Constructs and initialize Euler angles from a 3x3 rotation matrix \p m,
* with options to choose for each angle the requested range.
/** Constructs and initializes an EulerAngles from either:
* - a 3x3 rotation matrix expression(i.e. pure orthogonal matrix with determinant of +1),
* - a 3D vector expression representing Euler angles.
*
* If positive range is true, then the specified angle will be in the range [0, +2*PI].
* Otherwise, the specified angle will be in the range [-PI, +PI].
*
* \param m The 3x3 rotation matrix to convert
* \param positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \param positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \param positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
*/
* \note If \p other is a 3x3 rotation matrix, the angles range rules will be as follow:<BR>
* Alpha and gamma angles will be in the range [-PI, PI].<BR>
* As for Beta angle:
* - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2].
* - otherwise:
* - If the beta axis is positive, the beta angle will be in the range [0, PI]
* - If the beta axis is negative, the beta angle will be in the range [-PI, 0]
*/
template<typename Derived>
EulerAngles(
const MatrixBase<Derived>& m,
bool positiveRangeAlpha,
bool positiveRangeBeta,
bool positiveRangeGamma) {
System::CalcEulerAngles(*this, m, positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma);
}
explicit EulerAngles(const MatrixBase<Derived>& other) { *this = other; }
/** Constructs and initialize Euler angles from a rotation \p rot.
*
* \note All angles will be in the range [-PI, PI], unless \p rot is an EulerAngles.
* If rot is an EulerAngles, expected EulerAngles range is __undefined__.
* (Use other functions here for enforcing range if this effect is desired)
* \note If \p rot is an EulerAngles (even when it represented as RotationBase explicitly),
* angles ranges are __undefined__.
* Otherwise, alpha and gamma angles will be in the range [-PI, PI].<BR>
* As for Beta angle:
* - If the system is Tait-Bryan, the beta angle will be in the range [-PI/2, PI/2].
* - otherwise:
* - If the beta axis is positive, the beta angle will be in the range [0, PI]
* - If the beta axis is negative, the beta angle will be in the range [-PI, 0]
*/
template<typename Derived>
EulerAngles(const RotationBase<Derived, 3>& rot) { *this = rot; }
EulerAngles(const RotationBase<Derived, 3>& rot) { System::CalcEulerAngles(*this, rot.toRotationMatrix()); }
/** Constructs and initialize Euler angles from a rotation \p rot,
* with options to choose for each angle the requested range.
*
* If positive range is true, then the specified angle will be in the range [0, +2*PI].
* Otherwise, the specified angle will be in the range [-PI, +PI].
*
* \param rot The 3x3 rotation matrix to convert
* \param positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \param positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \param positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
*/
template<typename Derived>
EulerAngles(
const RotationBase<Derived, 3>& rot,
bool positiveRangeAlpha,
bool positiveRangeBeta,
bool positiveRangeGamma) {
System::CalcEulerAngles(*this, rot.toRotationMatrix(), positiveRangeAlpha, positiveRangeBeta, positiveRangeGamma);
}
/*EulerAngles(const QuaternionType& q)
{
// TODO: Implement it in a faster way for quaternions
// According to http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
// we can compute only the needed matrix cells and then convert to euler angles. (see ZYX example below)
// Currently we compute all matrix cells from quaternion.
// Special case only for ZYX
//Scalar y2 = q.y() * q.y();
//m_angles[0] = std::atan2(2*(q.w()*q.z() + q.x()*q.y()), (1 - 2*(y2 + q.z()*q.z())));
//m_angles[1] = std::asin( 2*(q.w()*q.y() - q.z()*q.x()));
//m_angles[2] = std::atan2(2*(q.w()*q.x() + q.y()*q.z()), (1 - 2*(q.x()*q.x() + y2)));
}*/
/** \returns The angle values stored in a vector (alpha, beta, gamma). */
const Vector3& angles() const { return m_angles; }
@@ -246,90 +227,48 @@ namespace Eigen
return inverse();
}
/** Constructs and initialize Euler angles from a 3x3 rotation matrix \p m,
* with options to choose for each angle the requested range (__only in compile time__).
/** Set \c *this from either:
* - a 3x3 rotation matrix expression(i.e. pure orthogonal matrix with determinant of +1),
* - a 3D vector expression representing Euler angles.
*
* If positive range is true, then the specified angle will be in the range [0, +2*PI].
* Otherwise, the specified angle will be in the range [-PI, +PI].
*
* \param m The 3x3 rotation matrix to convert
* \tparam positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \tparam positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \tparam positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
*/
template<
bool PositiveRangeAlpha,
bool PositiveRangeBeta,
bool PositiveRangeGamma,
typename Derived>
static EulerAngles FromRotation(const MatrixBase<Derived>& m)
{
EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Derived, 3, 3)
EulerAngles e;
System::template CalcEulerAngles<
PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma, _Scalar>(e, m);
return e;
}
/** Constructs and initialize Euler angles from a rotation \p rot,
* with options to choose for each angle the requested range (__only in compile time__).
*
* If positive range is true, then the specified angle will be in the range [0, +2*PI].
* Otherwise, the specified angle will be in the range [-PI, +PI].
*
* \param rot The 3x3 rotation matrix to convert
* \tparam positiveRangeAlpha If true, alpha will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \tparam positiveRangeBeta If true, beta will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* \tparam positiveRangeGamma If true, gamma will be in [0, 2*PI]. Otherwise, in [-PI, +PI].
* See EulerAngles(const MatrixBase<Derived, 3>&) for more information about
* angles ranges output.
*/
template<
bool PositiveRangeAlpha,
bool PositiveRangeBeta,
bool PositiveRangeGamma,
typename Derived>
static EulerAngles FromRotation(const RotationBase<Derived, 3>& rot)
template<class Derived>
EulerAngles& operator=(const MatrixBase<Derived>& other)
{
return FromRotation<PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma>(rot.toRotationMatrix());
}
/*EulerAngles& fromQuaternion(const QuaternionType& q)
{
// TODO: Implement it in a faster way for quaternions
// According to http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
// we can compute only the needed matrix cells and then convert to euler angles. (see ZYX example below)
// Currently we compute all matrix cells from quaternion.
// Special case only for ZYX
//Scalar y2 = q.y() * q.y();
//m_angles[0] = std::atan2(2*(q.w()*q.z() + q.x()*q.y()), (1 - 2*(y2 + q.z()*q.z())));
//m_angles[1] = std::asin( 2*(q.w()*q.y() - q.z()*q.x()));
//m_angles[2] = std::atan2(2*(q.w()*q.x() + q.y()*q.z()), (1 - 2*(q.x()*q.x() + y2)));
}*/
/** Set \c *this from a rotation matrix(i.e. pure orthogonal matrix with determinant of +1). */
template<typename Derived>
EulerAngles& operator=(const MatrixBase<Derived>& m) {
EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(Derived, 3, 3)
EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename Derived::Scalar>::value),
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
System::CalcEulerAngles(*this, m);
internal::eulerangles_assign_impl<System, Derived>::run(*this, other.derived());
return *this;
}
// TODO: Assign and construct from another EulerAngles (with different system)
/** Set \c *this from a rotation. */
/** Set \c *this from a rotation.
*
* See EulerAngles(const RotationBase<Derived, 3>&) for more information about
* angles ranges output.
*/
template<typename Derived>
EulerAngles& operator=(const RotationBase<Derived, 3>& rot) {
System::CalcEulerAngles(*this, rot.toRotationMatrix());
return *this;
}
// TODO: Support isApprox function
/** \returns \c true if \c *this is approximately equal to \a other, within the precision
* determined by \a prec.
*
* \sa MatrixBase::isApprox() */
bool isApprox(const EulerAngles& other,
const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const
{ return angles().isApprox(other.angles(), prec); }
/** \returns an equivalent 3x3 rotation matrix. */
Matrix3 toRotationMatrix() const
{
// TODO: Calc it faster
return static_cast<QuaternionType>(*this).toRotationMatrix();
}
@@ -347,6 +286,15 @@ namespace Eigen
s << eulerAngles.angles().transpose();
return s;
}
/** \returns \c *this with scalar type casted to \a NewScalarType */
template <typename NewScalarType>
EulerAngles<NewScalarType, System> cast() const
{
EulerAngles<NewScalarType, System> e;
e.angles() = angles().template cast<NewScalarType>();
return e;
}
};
#define EIGEN_EULER_ANGLES_SINGLE_TYPEDEF(AXES, SCALAR_TYPE, SCALAR_POSTFIX) \
@@ -379,8 +327,29 @@ EIGEN_EULER_ANGLES_TYPEDEFS(double, d)
{
typedef _Scalar Scalar;
};
// set from a rotation matrix
template<class System, class Other>
struct eulerangles_assign_impl<System,Other,3,3>
{
typedef typename Other::Scalar Scalar;
static void run(EulerAngles<Scalar, System>& e, const Other& m)
{
System::CalcEulerAngles(e, m);
}
};
// set from a vector of Euler angles
template<class System, class Other>
struct eulerangles_assign_impl<System,Other,4,1>
{
typedef typename Other::Scalar Scalar;
static void run(EulerAngles<Scalar, System>& e, const Other& vec)
{
e.angles() = vec;
}
};
}
}
#endif // EIGEN_EULERANGLESCLASS_H

184
unsupported/Eigen/src/EulerAngles/EulerSystem.h
View File

@@ -18,7 +18,7 @@ namespace Eigen
namespace internal
{
// TODO: Check if already exists on the rest API
// TODO: Add this trait to the Eigen internal API?
template <int Num, bool IsPositive = (Num > 0)>
struct Abs
{
@@ -36,6 +36,12 @@ namespace Eigen
{
enum { value = Axis != 0 && Abs<Axis>::value <= 3 };
};
template<typename System,
typename Other,
int OtherRows=Other::RowsAtCompileTime,
int OtherCols=Other::ColsAtCompileTime>
struct eulerangles_assign_impl;
}
#define EIGEN_EULER_ANGLES_CLASS_STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1]
@@ -69,7 +75,7 @@ namespace Eigen
*
* You can use this class to get two things:
* - Build an Euler system, and then pass it as a template parameter to EulerAngles.
* - Query some compile time data about an Euler system. (e.g. Whether it's tait bryan)
* - Query some compile time data about an Euler system. (e.g. Whether it's Tait-Bryan)
*
* Euler rotation is a set of three rotation on fixed axes. (see \ref EulerAngles)
* This meta-class store constantly those signed axes. (see \ref EulerAxis)
@@ -80,7 +86,7 @@ namespace Eigen
* signed axes{+X,+Y,+Z,-X,-Y,-Z} are supported:
* - all axes X, Y, Z in each valid order (see below what order is valid)
* - rotation over the axis is supported both over the positive and negative directions.
* - both tait bryan and proper/classic Euler angles (i.e. the opposite).
* - both Tait-Bryan and proper/classic Euler angles (i.e. the opposite).
*
* Since EulerSystem support both positive and negative directions,
* you may call this rotation distinction in other names:
@@ -90,7 +96,7 @@ namespace Eigen
* Notice all axed combination are valid, and would trigger a static assertion.
* Same unsigned axes can't be neighbors, e.g. {X,X,Y} is invalid.
* This yield two and only two classes:
* - _tait bryan_ - all unsigned axes are distinct, e.g. {X,Y,Z}
* - _Tait-Bryan_ - all unsigned axes are distinct, e.g. {X,Y,Z}
* - _proper/classic Euler angles_ - The first and the third unsigned axes is equal,
* and the second is different, e.g. {X,Y,X}
*
@@ -112,9 +118,9 @@ namespace Eigen
*
* \tparam _AlphaAxis the first fixed EulerAxis
*
* \tparam _AlphaAxis the second fixed EulerAxis
* \tparam _BetaAxis the second fixed EulerAxis
*
* \tparam _AlphaAxis the third fixed EulerAxis
* \tparam _GammaAxis the third fixed EulerAxis
*/
template <int _AlphaAxis, int _BetaAxis, int _GammaAxis>
class EulerSystem
@@ -138,14 +144,16 @@ namespace Eigen
BetaAxisAbs = internal::Abs<BetaAxis>::value, /*!< the second rotation axis unsigned */
GammaAxisAbs = internal::Abs<GammaAxis>::value, /*!< the third rotation axis unsigned */
IsAlphaOpposite = (AlphaAxis < 0) ? 1 : 0, /*!< weather alpha axis is negative */
IsBetaOpposite = (BetaAxis < 0) ? 1 : 0, /*!< weather beta axis is negative */
IsGammaOpposite = (GammaAxis < 0) ? 1 : 0, /*!< weather gamma axis is negative */
IsOdd = ((AlphaAxisAbs)%3 == (BetaAxisAbs - 1)%3) ? 0 : 1, /*!< weather the Euler system is odd */
IsEven = IsOdd ? 0 : 1, /*!< weather the Euler system is even */
IsAlphaOpposite = (AlphaAxis < 0) ? 1 : 0, /*!< whether alpha axis is negative */
IsBetaOpposite = (BetaAxis < 0) ? 1 : 0, /*!< whether beta axis is negative */
IsGammaOpposite = (GammaAxis < 0) ? 1 : 0, /*!< whether gamma axis is negative */
IsTaitBryan = ((unsigned)AlphaAxisAbs != (unsigned)GammaAxisAbs) ? 1 : 0 /*!< weather the Euler system is tait bryan */
// Parity is even if alpha axis X is followed by beta axis Y, or Y is followed
// by Z, or Z is followed by X; otherwise it is odd.
IsOdd = ((AlphaAxisAbs)%3 == (BetaAxisAbs - 1)%3) ? 0 : 1, /*!< whether the Euler system is odd */
IsEven = IsOdd ? 0 : 1, /*!< whether the Euler system is even */
IsTaitBryan = ((unsigned)AlphaAxisAbs != (unsigned)GammaAxisAbs) ? 1 : 0 /*!< whether the Euler system is Tait-Bryan */
};
private:
@@ -180,127 +188,99 @@ namespace Eigen
static void CalcEulerAngles_imp(Matrix<typename MatrixBase<Derived>::Scalar, 3, 1>& res, const MatrixBase<Derived>& mat, internal::true_type /*isTaitBryan*/)
{
using std::atan2;
using std::sin;
using std::cos;
using std::sqrt;
typedef typename Derived::Scalar Scalar;
typedef Matrix<Scalar,2,1> Vector2;
res[0] = atan2(mat(J,K), mat(K,K));
Scalar c2 = Vector2(mat(I,I), mat(I,J)).norm();
if((IsOdd && res[0]<Scalar(0)) || ((!IsOdd) && res[0]>Scalar(0))) {
if(res[0] > Scalar(0)) {
res[0] -= Scalar(EIGEN_PI);
}
else {
res[0] += Scalar(EIGEN_PI);
}
res[1] = atan2(-mat(I,K), -c2);
const Scalar plusMinus = IsEven? 1 : -1;
const Scalar minusPlus = IsOdd? 1 : -1;
const Scalar Rsum = sqrt((mat(I,I) * mat(I,I) + mat(I,J) * mat(I,J) + mat(J,K) * mat(J,K) + mat(K,K) * mat(K,K))/2);
res[1] = atan2(plusMinus * mat(I,K), Rsum);
// There is a singularity when cos(beta) == 0
if(Rsum > 4 * NumTraits<Scalar>::epsilon()) {// cos(beta) != 0
res[0] = atan2(minusPlus * mat(J, K), mat(K, K));
res[2] = atan2(minusPlus * mat(I, J), mat(I, I));
}
else if(plusMinus * mat(I, K) > 0) {// cos(beta) == 0 and sin(beta) == 1
Scalar spos = mat(J, I) + plusMinus * mat(K, J); // 2*sin(alpha + plusMinus * gamma
Scalar cpos = mat(J, J) + minusPlus * mat(K, I); // 2*cos(alpha + plusMinus * gamma)
Scalar alphaPlusMinusGamma = atan2(spos, cpos);
res[0] = alphaPlusMinusGamma;
res[2] = 0;
}
else {// cos(beta) == 0 and sin(beta) == -1
Scalar sneg = plusMinus * (mat(K, J) + minusPlus * mat(J, I)); // 2*sin(alpha + minusPlus*gamma)
Scalar cneg = mat(J, J) + plusMinus * mat(K, I); // 2*cos(alpha + minusPlus*gamma)
Scalar alphaMinusPlusBeta = atan2(sneg, cneg);
res[0] = alphaMinusPlusBeta;
res[2] = 0;
}
else
res[1] = atan2(-mat(I,K), c2);
Scalar s1 = sin(res[0]);
Scalar c1 = cos(res[0]);
res[2] = atan2(s1*mat(K,I)-c1*mat(J,I), c1*mat(J,J) - s1 * mat(K,J));
}
template <typename Derived>
static void CalcEulerAngles_imp(Matrix<typename MatrixBase<Derived>::Scalar,3,1>& res, const MatrixBase<Derived>& mat, internal::false_type /*isTaitBryan*/)
static void CalcEulerAngles_imp(Matrix<typename MatrixBase<Derived>::Scalar,3,1>& res,
const MatrixBase<Derived>& mat, internal::false_type /*isTaitBryan*/)
{
using std::atan2;
using std::sin;
using std::cos;
using std::sqrt;
typedef typename Derived::Scalar Scalar;
typedef Matrix<Scalar,2,1> Vector2;
res[0] = atan2(mat(J,I), mat(K,I));
if((IsOdd && res[0]<Scalar(0)) || ((!IsOdd) && res[0]>Scalar(0)))
{
if(res[0] > Scalar(0)) {
res[0] -= Scalar(EIGEN_PI);
}
else {
res[0] += Scalar(EIGEN_PI);
}
Scalar s2 = Vector2(mat(J,I), mat(K,I)).norm();
res[1] = -atan2(s2, mat(I,I));
}
else
{
Scalar s2 = Vector2(mat(J,I), mat(K,I)).norm();
res[1] = atan2(s2, mat(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
const Scalar plusMinus = IsEven? 1 : -1;
const Scalar minusPlus = IsOdd? 1 : -1;
Scalar s1 = sin(res[0]);
Scalar c1 = cos(res[0]);
res[2] = atan2(c1*mat(J,K)-s1*mat(K,K), c1*mat(J,J) - s1 * mat(K,J));
const Scalar Rsum = sqrt((mat(I, J) * mat(I, J) + mat(I, K) * mat(I, K) + mat(J, I) * mat(J, I) + mat(K, I) * mat(K, I)) / 2);
res[1] = atan2(Rsum, mat(I, I));
// There is a singularity when sin(beta) == 0
if(Rsum > 4 * NumTraits<Scalar>::epsilon()) {// sin(beta) != 0
res[0] = atan2(mat(J, I), minusPlus * mat(K, I));
res[2] = atan2(mat(I, J), plusMinus * mat(I, K));
}
else if(mat(I, I) > 0) {// sin(beta) == 0 and cos(beta) == 1
Scalar spos = plusMinus * mat(K, J) + minusPlus * mat(J, K); // 2*sin(alpha + gamma)
Scalar cpos = mat(J, J) + mat(K, K); // 2*cos(alpha + gamma)
res[0] = atan2(spos, cpos);
res[2] = 0;
}
else {// sin(beta) == 0 and cos(beta) == -1
Scalar sneg = plusMinus * mat(K, J) + plusMinus * mat(J, K); // 2*sin(alpha - gamma)
Scalar cneg = mat(J, J) - mat(K, K); // 2*cos(alpha - gamma)
res[0] = atan2(sneg, cneg);
res[2] = 0;
}
}
template<typename Scalar>
static void CalcEulerAngles(
EulerAngles<Scalar, EulerSystem>& res,
const typename EulerAngles<Scalar, EulerSystem>::Matrix3& mat)
{
CalcEulerAngles(res, mat, false, false, false);
}
template<
bool PositiveRangeAlpha,
bool PositiveRangeBeta,
bool PositiveRangeGamma,
typename Scalar>
static void CalcEulerAngles(
EulerAngles<Scalar, EulerSystem>& res,
const typename EulerAngles<Scalar, EulerSystem>::Matrix3& mat)
{
CalcEulerAngles(res, mat, PositiveRangeAlpha, PositiveRangeBeta, PositiveRangeGamma);
}
template<typename Scalar>
static void CalcEulerAngles(
EulerAngles<Scalar, EulerSystem>& res,
const typename EulerAngles<Scalar, EulerSystem>::Matrix3& mat,
bool PositiveRangeAlpha,
bool PositiveRangeBeta,
bool PositiveRangeGamma)
{
CalcEulerAngles_imp(
res.angles(), mat,
typename internal::conditional<IsTaitBryan, internal::true_type, internal::false_type>::type());
if (IsAlphaOpposite == IsOdd)
if (IsAlphaOpposite)
res.alpha() = -res.alpha();
if (IsBetaOpposite == IsOdd)
if (IsBetaOpposite)
res.beta() = -res.beta();
if (IsGammaOpposite == IsOdd)
if (IsGammaOpposite)
res.gamma() = -res.gamma();
// Saturate results to the requested range
if (PositiveRangeAlpha && (res.alpha() < 0))
res.alpha() += Scalar(2 * EIGEN_PI);
if (PositiveRangeBeta && (res.beta() < 0))
res.beta() += Scalar(2 * EIGEN_PI);
if (PositiveRangeGamma && (res.gamma() < 0))
res.gamma() += Scalar(2 * EIGEN_PI);
}
template <typename _Scalar, class _System>
friend class Eigen::EulerAngles;
template<typename System,
typename Other,
int OtherRows,
int OtherCols>
friend struct internal::eulerangles_assign_impl;
};
#define EIGEN_EULER_SYSTEM_TYPEDEF(A, B, C) \

50
unsupported/Eigen/src/Polynomials/Companion.h
View File

@@ -75,7 +75,7 @@ class companion
void setPolynomial( const VectorType& poly )
{
const Index deg = poly.size()-1;
m_monic = -1/poly[deg] * poly.head(deg);
m_monic = Scalar(-1)/poly[deg] * poly.head(deg);
//m_bl_diag.setIdentity( deg-1 );
m_bl_diag.setOnes(deg-1);
}
@@ -107,8 +107,8 @@ class companion
* colB and rowB are repectively the multipliers for
* the column and the row in order to balance them.
* */
bool balanced( Scalar colNorm, Scalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB );
bool balanced( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, RealScalar& colB, RealScalar& rowB );
/** Helper function for the balancing algorithm.
* \returns true if the row and the column, having colNorm and rowNorm
@@ -116,8 +116,8 @@ class companion
* colB and rowB are repectively the multipliers for
* the column and the row in order to balance them.
* */
bool balancedR( Scalar colNorm, Scalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB );
bool balancedR( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, RealScalar& colB, RealScalar& rowB );
public:
/**
@@ -139,10 +139,10 @@ class companion
template< typename _Scalar, int _Deg >
inline
bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB )
bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, RealScalar& colB, RealScalar& rowB )
{
if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
else
{
//To find the balancing coefficients, if the radix is 2,
@@ -150,29 +150,29 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
// \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
// then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
// and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
rowB = rowNorm / radix<Scalar>();
colB = Scalar(1);
const Scalar s = colNorm + rowNorm;
rowB = rowNorm / radix<RealScalar>();
colB = RealScalar(1);
const RealScalar s = colNorm + rowNorm;
while (colNorm < rowB)
{
colB *= radix<Scalar>();
colNorm *= radix2<Scalar>();
colB *= radix<RealScalar>();
colNorm *= radix2<RealScalar>();
}
rowB = rowNorm * radix<Scalar>();
rowB = rowNorm * radix<RealScalar>();
while (colNorm >= rowB)
{
colB /= radix<Scalar>();
colNorm /= radix2<Scalar>();
colB /= radix<RealScalar>();
colNorm /= radix2<RealScalar>();
}
//This line is used to avoid insubstantial balancing
if ((rowNorm + colNorm) < Scalar(0.95) * s * colB)
if ((rowNorm + colNorm) < RealScalar(0.95) * s * colB)
{
isBalanced = false;
rowB = Scalar(1) / colB;
rowB = RealScalar(1) / colB;
return false;
}
else{
@@ -182,21 +182,21 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
template< typename _Scalar, int _Deg >
inline
bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB )
bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, RealScalar& colB, RealScalar& rowB )
{
if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
else
{
/**
* Set the norm of the column and the row to the geometric mean
* of the row and column norm
*/
const _Scalar q = colNorm/rowNorm;
const RealScalar q = colNorm/rowNorm;
if( !isApprox( q, _Scalar(1) ) )
{
rowB = sqrt( colNorm/rowNorm );
colB = Scalar(1)/rowB;
colB = RealScalar(1)/rowB;
isBalanced = false;
return false;
@@ -219,8 +219,8 @@ void companion<_Scalar,_Deg>::balance()
while( !hasConverged )
{
hasConverged = true;
Scalar colNorm,rowNorm;
Scalar colB,rowB;
RealScalar colNorm,rowNorm;
RealScalar colB,rowB;
//First row, first column excluding the diagonal
//==============================================

18
unsupported/Eigen/src/Polynomials/PolynomialSolver.h
View File

@@ -99,7 +99,7 @@ class PolynomialSolverBase
*/
inline const RootType& greatestRoot() const
{
std::greater<Scalar> greater;
std::greater<RealScalar> greater;
return selectComplexRoot_withRespectToNorm( greater );
}
@@ -108,7 +108,7 @@ class PolynomialSolverBase
*/
inline const RootType& smallestRoot() const
{
std::less<Scalar> less;
std::less<RealScalar> less;
return selectComplexRoot_withRespectToNorm( less );
}
@@ -213,7 +213,7 @@ class PolynomialSolverBase
bool& hasArealRoot,
const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
{
std::greater<Scalar> greater;
std::greater<RealScalar> greater;
return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
}
@@ -236,7 +236,7 @@ class PolynomialSolverBase
bool& hasArealRoot,
const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
{
std::less<Scalar> less;
std::less<RealScalar> less;
return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
}
@@ -259,7 +259,7 @@ class PolynomialSolverBase
bool& hasArealRoot,
const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
{
std::greater<Scalar> greater;
std::greater<RealScalar> greater;
return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
}
@@ -282,7 +282,7 @@ class PolynomialSolverBase
bool& hasArealRoot,
const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
{
std::less<Scalar> less;
std::less<RealScalar> less;
return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
}
@@ -327,7 +327,7 @@ class PolynomialSolverBase
* However, almost always, correct accuracy is reached even in these cases for 64bit
* (double) floating types and small polynomial degree (<20).
*/
template< typename _Scalar, int _Deg >
template<typename _Scalar, int _Deg>
class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
{
public:
@@ -337,7 +337,9 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
typedef EigenSolver<CompanionMatrixType> EigenSolverType;
typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
ComplexEigenSolver<CompanionMatrixType>,
EigenSolver<CompanionMatrixType> >::type EigenSolverType;
public:
/** Computes the complex roots of a new polynomial. */