#ifndef INCLUDED_IMATHEULER_H
#define INCLUDED_IMATHEULER_H
#include "ImathMath.h"
#include "ImathVec.h"
#include "ImathQuat.h"
#include "ImathMatrix.h"
#include "ImathLimits.h"
#include <iostream>
namespace Imath {
#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
#pragma warning(disable:4244)
#endif
template <class T>
class Euler : public Vec3<T>
{
public:
using Vec3<T>::x;
using Vec3<T>::y;
using Vec3<T>::z;
enum Order
{
XYZ = 0x0101,
XZY = 0x0001,
YZX = 0x1101,
YXZ = 0x1001,
ZXY = 0x2101,
ZYX = 0x2001,
XZX = 0x0011,
XYX = 0x0111,
YXY = 0x1011,
YZY = 0x1111,
ZYZ = 0x2011,
ZXZ = 0x2111,
XYZr = 0x2000,
XZYr = 0x2100,
YZXr = 0x1000,
YXZr = 0x1100,
ZXYr = 0x0000,
ZYXr = 0x0100,
XZXr = 0x2110,
XYXr = 0x2010,
YXYr = 0x1110,
YZYr = 0x1010,
ZYZr = 0x0110,
ZXZr = 0x0010,
Legal = XYZ | XZY | YZX | YXZ | ZXY | ZYX |
XZX | XYX | YXY | YZY | ZYZ | ZXZ |
XYZr| XZYr| YZXr| YXZr| ZXYr| ZYXr|
XZXr| XYXr| YXYr| YZYr| ZYZr| ZXZr,
Min = 0x0000,
Max = 0x2111,
Default = XYZ
};
enum Axis { X = 0, Y = 1, Z = 2 };
enum InputLayout { XYZLayout, IJKLayout };
Euler();
Euler(const Euler&);
Euler(Order p);
Euler(const Vec3<T> &v, Order o = Default, InputLayout l = IJKLayout);
Euler(T i, T j, T k, Order o = Default, InputLayout l = IJKLayout);
Euler(const Euler<T> &euler, Order newp);
Euler(const Matrix33<T> &, Order o = Default);
Euler(const Matrix44<T> &, Order o = Default);
const Euler<T>& operator= (const Euler<T>&);
const Euler<T>& operator= (const Vec3<T>&);
static bool legal(Order);
void setXYZVector(const Vec3<T> &);
Order order() const;
void setOrder(Order);
void set(Axis initial,
bool relative,
bool parityEven,
bool firstRepeats);
void extract(const Matrix33<T>&);
void extract(const Matrix44<T>&);
void extract(const Quat<T>&);
Matrix33<T> toMatrix33() const;
Matrix44<T> toMatrix44() const;
Quat<T> toQuat() const;
Vec3<T> toXYZVector() const;
void angleOrder(int &i, int &j, int &k) const;
void angleMapping(int &i, int &j, int &k) const;
static float angleMod (T angle);
static void simpleXYZRotation (Vec3<T> &xyzRot,
const Vec3<T> &targetXyzRot);
static void nearestRotation (Vec3<T> &xyzRot,
const Vec3<T> &targetXyzRot,
Order order = XYZ);
void makeNear (const Euler<T> &target);
bool frameStatic() const { return _frameStatic; }
bool initialRepeated() const { return _initialRepeated; }
bool parityEven() const { return _parityEven; }
Axis initialAxis() const { return _initialAxis; }
protected:
bool _frameStatic : 1;
bool _initialRepeated : 1;
bool _parityEven : 1;
#if defined _WIN32 || defined _WIN64
Axis _initialAxis ;
#else
Axis _initialAxis : 2;
#endif
};
typedef Euler<float> Eulerf;
typedef Euler<double> Eulerd;
template<class T>
inline void
Euler<T>::angleOrder(int &i, int &j, int &k) const
{
i = _initialAxis;
j = _parityEven ? (i+1)%3 : (i > 0 ? i-1 : 2);
k = _parityEven ? (i > 0 ? i-1 : 2) : (i+1)%3;
}
template<class T>
inline void
Euler<T>::angleMapping(int &i, int &j, int &k) const
{
int m[3];
m[_initialAxis] = 0;
m[(_initialAxis+1) % 3] = _parityEven ? 1 : 2;
m[(_initialAxis+2) % 3] = _parityEven ? 2 : 1;
i = m[0];
j = m[1];
k = m[2];
}
template<class T>
inline void
Euler<T>::setXYZVector(const Vec3<T> &v)
{
int i,j,k;
angleMapping(i,j,k);
(*this)[i] = v.x;
(*this)[j] = v.y;
(*this)[k] = v.z;
}
template<class T>
inline Vec3<T>
Euler<T>::toXYZVector() const
{
int i,j,k;
angleMapping(i,j,k);
return Vec3<T>((*this)[i],(*this)[j],(*this)[k]);
}
template<class T>
Euler<T>::Euler() :
Vec3<T>(0,0,0),
_frameStatic(true),
_initialRepeated(false),
_parityEven(true),
_initialAxis(X)
{}
template<class T>
Euler<T>::Euler(typename Euler<T>::Order p) :
Vec3<T>(0,0,0),
_frameStatic(true),
_initialRepeated(false),
_parityEven(true),
_initialAxis(X)
{
setOrder(p);
}
template<class T>
inline Euler<T>::Euler( const Vec3<T> &v,
typename Euler<T>::Order p,
typename Euler<T>::InputLayout l )
{
setOrder(p);
if ( l == XYZLayout ) setXYZVector(v);
else { x = v.x; y = v.y; z = v.z; }
}
template<class T>
inline Euler<T>::Euler(const Euler<T> &euler)
{
operator=(euler);
}
template<class T>
inline Euler<T>::Euler(const Euler<T> &euler,Order p)
{
setOrder(p);
Matrix33<T> M = euler.toMatrix33();
extract(M);
}
template<class T>
inline Euler<T>::Euler( T xi, T yi, T zi,
typename Euler<T>::Order p,
typename Euler<T>::InputLayout l)
{
setOrder(p);
if ( l == XYZLayout ) setXYZVector(Vec3<T>(xi,yi,zi));
else { x = xi; y = yi; z = zi; }
}
template<class T>
inline Euler<T>::Euler( const Matrix33<T> &M, typename Euler::Order p )
{
setOrder(p);
extract(M);
}
template<class T>
inline Euler<T>::Euler( const Matrix44<T> &M, typename Euler::Order p )
{
setOrder(p);
extract(M);
}
template<class T>
inline void Euler<T>::extract(const Quat<T> &q)
{
extract(q.toMatrix33());
}
template<class T>
void Euler<T>::extract(const Matrix33<T> &M)
{
int i,j,k;
angleOrder(i,j,k);
if (_initialRepeated)
{
x = Math<T>::atan2 (M[j][i], M[k][i]);
Vec3<T> r (0, 0, 0);
r[i] = (_parityEven? -x: x);
Matrix44<T> N;
N.rotate (r);
N = N * Matrix44<T> (M[0][0], M[0][1], M[0][2], 0,
M[1][0], M[1][1], M[1][2], 0,
M[2][0], M[2][1], M[2][2], 0,
0, 0, 0, 1);
T sy = Math<T>::sqrt (N[j][i]*N[j][i] + N[k][i]*N[k][i]);
y = Math<T>::atan2 (sy, N[i][i]);
z = Math<T>::atan2 (N[j][k], N[j][j]);
}
else
{
x = Math<T>::atan2 (M[j][k], M[k][k]);
Vec3<T> r (0, 0, 0);
r[i] = (_parityEven? -x: x);
Matrix44<T> N;
N.rotate (r);
N = N * Matrix44<T> (M[0][0], M[0][1], M[0][2], 0,
M[1][0], M[1][1], M[1][2], 0,
M[2][0], M[2][1], M[2][2], 0,
0, 0, 0, 1);
T cy = Math<T>::sqrt (N[i][i]*N[i][i] + N[i][j]*N[i][j]);
y = Math<T>::atan2 (-N[i][k], cy);
z = Math<T>::atan2 (-N[j][i], N[j][j]);
}
if (!_parityEven)
*this *= -1;
if (!_frameStatic)
{
T t = x;
x = z;
z = t;
}
}
template<class T>
void Euler<T>::extract(const Matrix44<T> &M)
{
int i,j,k;
angleOrder(i,j,k);
if (_initialRepeated)
{
x = Math<T>::atan2 (M[j][i], M[k][i]);
Vec3<T> r (0, 0, 0);
r[i] = (_parityEven? -x: x);
Matrix44<T> N;
N.rotate (r);
N = N * M;
T sy = Math<T>::sqrt (N[j][i]*N[j][i] + N[k][i]*N[k][i]);
y = Math<T>::atan2 (sy, N[i][i]);
z = Math<T>::atan2 (N[j][k], N[j][j]);
}
else
{
x = Math<T>::atan2 (M[j][k], M[k][k]);
Vec3<T> r (0, 0, 0);
r[i] = (_parityEven? -x: x);
Matrix44<T> N;
N.rotate (r);
N = N * M;
T cy = Math<T>::sqrt (N[i][i]*N[i][i] + N[i][j]*N[i][j]);
y = Math<T>::atan2 (-N[i][k], cy);
z = Math<T>::atan2 (-N[j][i], N[j][j]);
}
if (!_parityEven)
*this *= -1;
if (!_frameStatic)
{
T t = x;
x = z;
z = t;
}
}
template<class T>
Matrix33<T> Euler<T>::toMatrix33() const
{
int i,j,k;
angleOrder(i,j,k);
Vec3<T> angles;
if ( _frameStatic ) angles = (*this);
else angles = Vec3<T>(z,y,x);
if ( !_parityEven ) angles *= -1.0;
T ci = Math<T>::cos(angles.x);
T cj = Math<T>::cos(angles.y);
T ch = Math<T>::cos(angles.z);
T si = Math<T>::sin(angles.x);
T sj = Math<T>::sin(angles.y);
T sh = Math<T>::sin(angles.z);
T cc = ci*ch;
T cs = ci*sh;
T sc = si*ch;
T ss = si*sh;
Matrix33<T> M;
if ( _initialRepeated )
{
M[i][i] = cj; M[j][i] = sj*si; M[k][i] = sj*ci;
M[i][j] = sj*sh; M[j][j] = -cj*ss+cc; M[k][j] = -cj*cs-sc;
M[i][k] = -sj*ch; M[j][k] = cj*sc+cs; M[k][k] = cj*cc-ss;
}
else
{
M[i][i] = cj*ch; M[j][i] = sj*sc-cs; M[k][i] = sj*cc+ss;
M[i][j] = cj*sh; M[j][j] = sj*ss+cc; M[k][j] = sj*cs-sc;
M[i][k] = -sj; M[j][k] = cj*si; M[k][k] = cj*ci;
}
return M;
}
template<class T>
Matrix44<T> Euler<T>::toMatrix44() const
{
int i,j,k;
angleOrder(i,j,k);
Vec3<T> angles;
if ( _frameStatic ) angles = (*this);
else angles = Vec3<T>(z,y,x);
if ( !_parityEven ) angles *= -1.0;
T ci = Math<T>::cos(angles.x);
T cj = Math<T>::cos(angles.y);
T ch = Math<T>::cos(angles.z);
T si = Math<T>::sin(angles.x);
T sj = Math<T>::sin(angles.y);
T sh = Math<T>::sin(angles.z);
T cc = ci*ch;
T cs = ci*sh;
T sc = si*ch;
T ss = si*sh;
Matrix44<T> M;
if ( _initialRepeated )
{
M[i][i] = cj; M[j][i] = sj*si; M[k][i] = sj*ci;
M[i][j] = sj*sh; M[j][j] = -cj*ss+cc; M[k][j] = -cj*cs-sc;
M[i][k] = -sj*ch; M[j][k] = cj*sc+cs; M[k][k] = cj*cc-ss;
}
else
{
M[i][i] = cj*ch; M[j][i] = sj*sc-cs; M[k][i] = sj*cc+ss;
M[i][j] = cj*sh; M[j][j] = sj*ss+cc; M[k][j] = sj*cs-sc;
M[i][k] = -sj; M[j][k] = cj*si; M[k][k] = cj*ci;
}
return M;
}
template<class T>
Quat<T> Euler<T>::toQuat() const
{
Vec3<T> angles;
int i,j,k;
angleOrder(i,j,k);
if ( _frameStatic ) angles = (*this);
else angles = Vec3<T>(z,y,x);
if ( !_parityEven ) angles.y = -angles.y;
T ti = angles.x*0.5;
T tj = angles.y*0.5;
T th = angles.z*0.5;
T ci = Math<T>::cos(ti);
T cj = Math<T>::cos(tj);
T ch = Math<T>::cos(th);
T si = Math<T>::sin(ti);
T sj = Math<T>::sin(tj);
T sh = Math<T>::sin(th);
T cc = ci*ch;
T cs = ci*sh;
T sc = si*ch;
T ss = si*sh;
T parity = _parityEven ? 1.0 : -1.0;
Quat<T> q;
Vec3<T> a;
if ( _initialRepeated )
{
a[i] = cj*(cs + sc);
a[j] = sj*(cc + ss) * parity,
a[k] = sj*(cs - sc);
q.r = cj*(cc - ss);
}
else
{
a[i] = cj*sc - sj*cs,
a[j] = (cj*ss + sj*cc) * parity,
a[k] = cj*cs - sj*sc;
q.r = cj*cc + sj*ss;
}
q.v = a;
return q;
}
template<class T>
inline bool
Euler<T>::legal(typename Euler<T>::Order order)
{
return (order & ~Legal) ? false : true;
}
template<class T>
typename Euler<T>::Order
Euler<T>::order() const
{
int foo = (_initialAxis == Z ? 0x2000 : (_initialAxis == Y ? 0x1000 : 0));
if (_parityEven) foo |= 0x0100;
if (_initialRepeated) foo |= 0x0010;
if (_frameStatic) foo++;
return (Order)foo;
}
template<class T>
inline void Euler<T>::setOrder(typename Euler<T>::Order p)
{
set( p & 0x2000 ? Z : (p & 0x1000 ? Y : X),
!(p & 0x1),
!!(p & 0x100),
!!(p & 0x10));
}
template<class T>
void Euler<T>::set(typename Euler<T>::Axis axis,
bool relative,
bool parityEven,
bool firstRepeats)
{
_initialAxis = axis;
_frameStatic = !relative;
_parityEven = parityEven;
_initialRepeated = firstRepeats;
}
template<class T>
const Euler<T>& Euler<T>::operator= (const Euler<T> &euler)
{
x = euler.x;
y = euler.y;
z = euler.z;
_initialAxis = euler._initialAxis;
_frameStatic = euler._frameStatic;
_parityEven = euler._parityEven;
_initialRepeated = euler._initialRepeated;
return *this;
}
template<class T>
const Euler<T>& Euler<T>::operator= (const Vec3<T> &v)
{
x = v.x;
y = v.y;
z = v.z;
return *this;
}
template<class T>
std::ostream& operator << (std::ostream &o, const Euler<T> &euler)
{
char a[3] = { 'X', 'Y', 'Z' };
const char* r = euler.frameStatic() ? "" : "r";
int i,j,k;
euler.angleOrder(i,j,k);
if ( euler.initialRepeated() ) k = i;
return o << "("
<< euler.x << " "
<< euler.y << " "
<< euler.z << " "
<< a[i] << a[j] << a[k] << r << ")";
}
template <class T>
float
Euler<T>::angleMod (T angle)
{
angle = fmod(T (angle), T (2 * M_PI));
if (angle < -M_PI) angle += 2 * M_PI;
if (angle > +M_PI) angle -= 2 * M_PI;
return angle;
}
template <class T>
void
Euler<T>::simpleXYZRotation (Vec3<T> &xyzRot, const Vec3<T> &targetXyzRot)
{
Vec3<T> d = xyzRot - targetXyzRot;
xyzRot[0] = targetXyzRot[0] + angleMod(d[0]);
xyzRot[1] = targetXyzRot[1] + angleMod(d[1]);
xyzRot[2] = targetXyzRot[2] + angleMod(d[2]);
}
template <class T>
void
Euler<T>::nearestRotation (Vec3<T> &xyzRot, const Vec3<T> &targetXyzRot,
Order order)
{
int i,j,k;
Euler<T> e (0,0,0, order);
e.angleOrder(i,j,k);
simpleXYZRotation(xyzRot, targetXyzRot);
Vec3<T> otherXyzRot;
otherXyzRot[i] = M_PI+xyzRot[i];
otherXyzRot[j] = M_PI-xyzRot[j];
otherXyzRot[k] = M_PI+xyzRot[k];
simpleXYZRotation(otherXyzRot, targetXyzRot);
Vec3<T> d = xyzRot - targetXyzRot;
Vec3<T> od = otherXyzRot - targetXyzRot;
T dMag = d.dot(d);
T odMag = od.dot(od);
if (odMag < dMag)
{
xyzRot = otherXyzRot;
}
}
template <class T>
void
Euler<T>::makeNear (const Euler<T> &target)
{
Vec3<T> xyzRot = toXYZVector();
Vec3<T> targetXyz;
if (order() != target.order())
{
Euler<T> targetSameOrder = Euler<T>(target, order());
targetXyz = targetSameOrder.toXYZVector();
}
else
{
targetXyz = target.toXYZVector();
}
nearestRotation(xyzRot, targetXyz, order());
setXYZVector(xyzRot);
}
#if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
#pragma warning(default:4244)
#endif
}
#endif