1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. Eigen itself is part of the KDE project. 3// 4// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> 5// 6// This Source Code Form is subject to the terms of the Mozilla 7// Public License v. 2.0. If a copy of the MPL was not distributed 8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10// no include guard, we'll include this twice from All.h from Eigen2Support, and it's internal anyway 11 12namespace Eigen { 13 14template<typename Other, 15 int OtherRows=Other::RowsAtCompileTime, 16 int OtherCols=Other::ColsAtCompileTime> 17struct ei_quaternion_assign_impl; 18 19/** \geometry_module \ingroup Geometry_Module 20 * 21 * \class Quaternion 22 * 23 * \brief The quaternion class used to represent 3D orientations and rotations 24 * 25 * \param _Scalar the scalar type, i.e., the type of the coefficients 26 * 27 * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of 28 * orientations and rotations of objects in three dimensions. Compared to other representations 29 * like Euler angles or 3x3 matrices, quatertions offer the following advantages: 30 * \li \b compact storage (4 scalars) 31 * \li \b efficient to compose (28 flops), 32 * \li \b stable spherical interpolation 33 * 34 * The following two typedefs are provided for convenience: 35 * \li \c Quaternionf for \c float 36 * \li \c Quaterniond for \c double 37 * 38 * \sa class AngleAxis, class Transform 39 */ 40 41template<typename _Scalar> struct ei_traits<Quaternion<_Scalar> > 42{ 43 typedef _Scalar Scalar; 44}; 45 46template<typename _Scalar> 47class Quaternion : public RotationBase<Quaternion<_Scalar>,3> 48{ 49 typedef RotationBase<Quaternion<_Scalar>,3> Base; 50 51public: 52 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,4) 53 54 using Base::operator*; 55 56 /** the scalar type of the coefficients */ 57 typedef _Scalar Scalar; 58 59 /** the type of the Coefficients 4-vector */ 60 typedef Matrix<Scalar, 4, 1> Coefficients; 61 /** the type of a 3D vector */ 62 typedef Matrix<Scalar,3,1> Vector3; 63 /** the equivalent rotation matrix type */ 64 typedef Matrix<Scalar,3,3> Matrix3; 65 /** the equivalent angle-axis type */ 66 typedef AngleAxis<Scalar> AngleAxisType; 67 68 /** \returns the \c x coefficient */ 69 inline Scalar x() const { return m_coeffs.coeff(0); } 70 /** \returns the \c y coefficient */ 71 inline Scalar y() const { return m_coeffs.coeff(1); } 72 /** \returns the \c z coefficient */ 73 inline Scalar z() const { return m_coeffs.coeff(2); } 74 /** \returns the \c w coefficient */ 75 inline Scalar w() const { return m_coeffs.coeff(3); } 76 77 /** \returns a reference to the \c x coefficient */ 78 inline Scalar& x() { return m_coeffs.coeffRef(0); } 79 /** \returns a reference to the \c y coefficient */ 80 inline Scalar& y() { return m_coeffs.coeffRef(1); } 81 /** \returns a reference to the \c z coefficient */ 82 inline Scalar& z() { return m_coeffs.coeffRef(2); } 83 /** \returns a reference to the \c w coefficient */ 84 inline Scalar& w() { return m_coeffs.coeffRef(3); } 85 86 /** \returns a read-only vector expression of the imaginary part (x,y,z) */ 87 inline const Block<const Coefficients,3,1> vec() const { return m_coeffs.template start<3>(); } 88 89 /** \returns a vector expression of the imaginary part (x,y,z) */ 90 inline Block<Coefficients,3,1> vec() { return m_coeffs.template start<3>(); } 91 92 /** \returns a read-only vector expression of the coefficients (x,y,z,w) */ 93 inline const Coefficients& coeffs() const { return m_coeffs; } 94 95 /** \returns a vector expression of the coefficients (x,y,z,w) */ 96 inline Coefficients& coeffs() { return m_coeffs; } 97 98 /** Default constructor leaving the quaternion uninitialized. */ 99 inline Quaternion() {} 100 101 /** Constructs and initializes the quaternion \f$ w+xi+yj+zk \f$ from 102 * its four coefficients \a w, \a x, \a y and \a z. 103 * 104 * \warning Note the order of the arguments: the real \a w coefficient first, 105 * while internally the coefficients are stored in the following order: 106 * [\c x, \c y, \c z, \c w] 107 */ 108 inline Quaternion(Scalar w, Scalar x, Scalar y, Scalar z) 109 { m_coeffs << x, y, z, w; } 110 111 /** Copy constructor */ 112 inline Quaternion(const Quaternion& other) { m_coeffs = other.m_coeffs; } 113 114 /** Constructs and initializes a quaternion from the angle-axis \a aa */ 115 explicit inline Quaternion(const AngleAxisType& aa) { *this = aa; } 116 117 /** Constructs and initializes a quaternion from either: 118 * - a rotation matrix expression, 119 * - a 4D vector expression representing quaternion coefficients. 120 * \sa operator=(MatrixBase<Derived>) 121 */ 122 template<typename Derived> 123 explicit inline Quaternion(const MatrixBase<Derived>& other) { *this = other; } 124 125 Quaternion& operator=(const Quaternion& other); 126 Quaternion& operator=(const AngleAxisType& aa); 127 template<typename Derived> 128 Quaternion& operator=(const MatrixBase<Derived>& m); 129 130 /** \returns a quaternion representing an identity rotation 131 * \sa MatrixBase::Identity() 132 */ 133 static inline Quaternion Identity() { return Quaternion(1, 0, 0, 0); } 134 135 /** \sa Quaternion::Identity(), MatrixBase::setIdentity() 136 */ 137 inline Quaternion& setIdentity() { m_coeffs << 0, 0, 0, 1; return *this; } 138 139 /** \returns the squared norm of the quaternion's coefficients 140 * \sa Quaternion::norm(), MatrixBase::squaredNorm() 141 */ 142 inline Scalar squaredNorm() const { return m_coeffs.squaredNorm(); } 143 144 /** \returns the norm of the quaternion's coefficients 145 * \sa Quaternion::squaredNorm(), MatrixBase::norm() 146 */ 147 inline Scalar norm() const { return m_coeffs.norm(); } 148 149 /** Normalizes the quaternion \c *this 150 * \sa normalized(), MatrixBase::normalize() */ 151 inline void normalize() { m_coeffs.normalize(); } 152 /** \returns a normalized version of \c *this 153 * \sa normalize(), MatrixBase::normalized() */ 154 inline Quaternion normalized() const { return Quaternion(m_coeffs.normalized()); } 155 156 /** \returns the dot product of \c *this and \a other 157 * Geometrically speaking, the dot product of two unit quaternions 158 * corresponds to the cosine of half the angle between the two rotations. 159 * \sa angularDistance() 160 */ 161 inline Scalar eigen2_dot(const Quaternion& other) const { return m_coeffs.eigen2_dot(other.m_coeffs); } 162 163 inline Scalar angularDistance(const Quaternion& other) const; 164 165 Matrix3 toRotationMatrix(void) const; 166 167 template<typename Derived1, typename Derived2> 168 Quaternion& setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b); 169 170 inline Quaternion operator* (const Quaternion& q) const; 171 inline Quaternion& operator*= (const Quaternion& q); 172 173 Quaternion inverse(void) const; 174 Quaternion conjugate(void) const; 175 176 Quaternion slerp(Scalar t, const Quaternion& other) const; 177 178 template<typename Derived> 179 Vector3 operator* (const MatrixBase<Derived>& vec) const; 180 181 /** \returns \c *this with scalar type casted to \a NewScalarType 182 * 183 * Note that if \a NewScalarType is equal to the current scalar type of \c *this 184 * then this function smartly returns a const reference to \c *this. 185 */ 186 template<typename NewScalarType> 187 inline typename internal::cast_return_type<Quaternion,Quaternion<NewScalarType> >::type cast() const 188 { return typename internal::cast_return_type<Quaternion,Quaternion<NewScalarType> >::type(*this); } 189 190 /** Copy constructor with scalar type conversion */ 191 template<typename OtherScalarType> 192 inline explicit Quaternion(const Quaternion<OtherScalarType>& other) 193 { m_coeffs = other.coeffs().template cast<Scalar>(); } 194 195 /** \returns \c true if \c *this is approximately equal to \a other, within the precision 196 * determined by \a prec. 197 * 198 * \sa MatrixBase::isApprox() */ 199 bool isApprox(const Quaternion& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const 200 { return m_coeffs.isApprox(other.m_coeffs, prec); } 201 202protected: 203 Coefficients m_coeffs; 204}; 205 206/** \ingroup Geometry_Module 207 * single precision quaternion type */ 208typedef Quaternion<float> Quaternionf; 209/** \ingroup Geometry_Module 210 * double precision quaternion type */ 211typedef Quaternion<double> Quaterniond; 212 213// Generic Quaternion * Quaternion product 214template<typename Scalar> inline Quaternion<Scalar> 215ei_quaternion_product(const Quaternion<Scalar>& a, const Quaternion<Scalar>& b) 216{ 217 return Quaternion<Scalar> 218 ( 219 a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z(), 220 a.w() * b.x() + a.x() * b.w() + a.y() * b.z() - a.z() * b.y(), 221 a.w() * b.y() + a.y() * b.w() + a.z() * b.x() - a.x() * b.z(), 222 a.w() * b.z() + a.z() * b.w() + a.x() * b.y() - a.y() * b.x() 223 ); 224} 225 226/** \returns the concatenation of two rotations as a quaternion-quaternion product */ 227template <typename Scalar> 228inline Quaternion<Scalar> Quaternion<Scalar>::operator* (const Quaternion& other) const 229{ 230 return ei_quaternion_product(*this,other); 231} 232 233/** \sa operator*(Quaternion) */ 234template <typename Scalar> 235inline Quaternion<Scalar>& Quaternion<Scalar>::operator*= (const Quaternion& other) 236{ 237 return (*this = *this * other); 238} 239 240/** Rotation of a vector by a quaternion. 241 * \remarks If the quaternion is used to rotate several points (>1) 242 * then it is much more efficient to first convert it to a 3x3 Matrix. 243 * Comparison of the operation cost for n transformations: 244 * - Quaternion: 30n 245 * - Via a Matrix3: 24 + 15n 246 */ 247template <typename Scalar> 248template<typename Derived> 249inline typename Quaternion<Scalar>::Vector3 250Quaternion<Scalar>::operator* (const MatrixBase<Derived>& v) const 251{ 252 // Note that this algorithm comes from the optimization by hand 253 // of the conversion to a Matrix followed by a Matrix/Vector product. 254 // It appears to be much faster than the common algorithm found 255 // in the litterature (30 versus 39 flops). It also requires two 256 // Vector3 as temporaries. 257 Vector3 uv; 258 uv = 2 * this->vec().cross(v); 259 return v + this->w() * uv + this->vec().cross(uv); 260} 261 262template<typename Scalar> 263inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const Quaternion& other) 264{ 265 m_coeffs = other.m_coeffs; 266 return *this; 267} 268 269/** Set \c *this from an angle-axis \a aa and returns a reference to \c *this 270 */ 271template<typename Scalar> 272inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const AngleAxisType& aa) 273{ 274 Scalar ha = Scalar(0.5)*aa.angle(); // Scalar(0.5) to suppress precision loss warnings 275 this->w() = ei_cos(ha); 276 this->vec() = ei_sin(ha) * aa.axis(); 277 return *this; 278} 279 280/** Set \c *this from the expression \a xpr: 281 * - if \a xpr is a 4x1 vector, then \a xpr is assumed to be a quaternion 282 * - if \a xpr is a 3x3 matrix, then \a xpr is assumed to be rotation matrix 283 * and \a xpr is converted to a quaternion 284 */ 285template<typename Scalar> 286template<typename Derived> 287inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const MatrixBase<Derived>& xpr) 288{ 289 ei_quaternion_assign_impl<Derived>::run(*this, xpr.derived()); 290 return *this; 291} 292 293/** Convert the quaternion to a 3x3 rotation matrix */ 294template<typename Scalar> 295inline typename Quaternion<Scalar>::Matrix3 296Quaternion<Scalar>::toRotationMatrix(void) const 297{ 298 // NOTE if inlined, then gcc 4.2 and 4.4 get rid of the temporary (not gcc 4.3 !!) 299 // if not inlined then the cost of the return by value is huge ~ +35%, 300 // however, not inlining this function is an order of magnitude slower, so 301 // it has to be inlined, and so the return by value is not an issue 302 Matrix3 res; 303 304 const Scalar tx = Scalar(2)*this->x(); 305 const Scalar ty = Scalar(2)*this->y(); 306 const Scalar tz = Scalar(2)*this->z(); 307 const Scalar twx = tx*this->w(); 308 const Scalar twy = ty*this->w(); 309 const Scalar twz = tz*this->w(); 310 const Scalar txx = tx*this->x(); 311 const Scalar txy = ty*this->x(); 312 const Scalar txz = tz*this->x(); 313 const Scalar tyy = ty*this->y(); 314 const Scalar tyz = tz*this->y(); 315 const Scalar tzz = tz*this->z(); 316 317 res.coeffRef(0,0) = Scalar(1)-(tyy+tzz); 318 res.coeffRef(0,1) = txy-twz; 319 res.coeffRef(0,2) = txz+twy; 320 res.coeffRef(1,0) = txy+twz; 321 res.coeffRef(1,1) = Scalar(1)-(txx+tzz); 322 res.coeffRef(1,2) = tyz-twx; 323 res.coeffRef(2,0) = txz-twy; 324 res.coeffRef(2,1) = tyz+twx; 325 res.coeffRef(2,2) = Scalar(1)-(txx+tyy); 326 327 return res; 328} 329 330/** Sets *this to be a quaternion representing a rotation sending the vector \a a to the vector \a b. 331 * 332 * \returns a reference to *this. 333 * 334 * Note that the two input vectors do \b not have to be normalized. 335 */ 336template<typename Scalar> 337template<typename Derived1, typename Derived2> 338inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b) 339{ 340 Vector3 v0 = a.normalized(); 341 Vector3 v1 = b.normalized(); 342 Scalar c = v0.eigen2_dot(v1); 343 344 // if dot == 1, vectors are the same 345 if (ei_isApprox(c,Scalar(1))) 346 { 347 // set to identity 348 this->w() = 1; this->vec().setZero(); 349 return *this; 350 } 351 // if dot == -1, vectors are opposites 352 if (ei_isApprox(c,Scalar(-1))) 353 { 354 this->vec() = v0.unitOrthogonal(); 355 this->w() = 0; 356 return *this; 357 } 358 359 Vector3 axis = v0.cross(v1); 360 Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2)); 361 Scalar invs = Scalar(1)/s; 362 this->vec() = axis * invs; 363 this->w() = s * Scalar(0.5); 364 365 return *this; 366} 367 368/** \returns the multiplicative inverse of \c *this 369 * Note that in most cases, i.e., if you simply want the opposite rotation, 370 * and/or the quaternion is normalized, then it is enough to use the conjugate. 371 * 372 * \sa Quaternion::conjugate() 373 */ 374template <typename Scalar> 375inline Quaternion<Scalar> Quaternion<Scalar>::inverse() const 376{ 377 // FIXME should this function be called multiplicativeInverse and conjugate() be called inverse() or opposite() ?? 378 Scalar n2 = this->squaredNorm(); 379 if (n2 > 0) 380 return Quaternion(conjugate().coeffs() / n2); 381 else 382 { 383 // return an invalid result to flag the error 384 return Quaternion(Coefficients::Zero()); 385 } 386} 387 388/** \returns the conjugate of the \c *this which is equal to the multiplicative inverse 389 * if the quaternion is normalized. 390 * The conjugate of a quaternion represents the opposite rotation. 391 * 392 * \sa Quaternion::inverse() 393 */ 394template <typename Scalar> 395inline Quaternion<Scalar> Quaternion<Scalar>::conjugate() const 396{ 397 return Quaternion(this->w(),-this->x(),-this->y(),-this->z()); 398} 399 400/** \returns the angle (in radian) between two rotations 401 * \sa eigen2_dot() 402 */ 403template <typename Scalar> 404inline Scalar Quaternion<Scalar>::angularDistance(const Quaternion& other) const 405{ 406 double d = ei_abs(this->eigen2_dot(other)); 407 if (d>=1.0) 408 return 0; 409 return Scalar(2) * std::acos(d); 410} 411 412/** \returns the spherical linear interpolation between the two quaternions 413 * \c *this and \a other at the parameter \a t 414 */ 415template <typename Scalar> 416Quaternion<Scalar> Quaternion<Scalar>::slerp(Scalar t, const Quaternion& other) const 417{ 418 static const Scalar one = Scalar(1) - machine_epsilon<Scalar>(); 419 Scalar d = this->eigen2_dot(other); 420 Scalar absD = ei_abs(d); 421 422 Scalar scale0; 423 Scalar scale1; 424 425 if (absD>=one) 426 { 427 scale0 = Scalar(1) - t; 428 scale1 = t; 429 } 430 else 431 { 432 // theta is the angle between the 2 quaternions 433 Scalar theta = std::acos(absD); 434 Scalar sinTheta = ei_sin(theta); 435 436 scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta; 437 scale1 = ei_sin( ( t * theta) ) / sinTheta; 438 if (d<0) 439 scale1 = -scale1; 440 } 441 442 return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs()); 443} 444 445// set from a rotation matrix 446template<typename Other> 447struct ei_quaternion_assign_impl<Other,3,3> 448{ 449 typedef typename Other::Scalar Scalar; 450 static inline void run(Quaternion<Scalar>& q, const Other& mat) 451 { 452 // This algorithm comes from "Quaternion Calculus and Fast Animation", 453 // Ken Shoemake, 1987 SIGGRAPH course notes 454 Scalar t = mat.trace(); 455 if (t > 0) 456 { 457 t = ei_sqrt(t + Scalar(1.0)); 458 q.w() = Scalar(0.5)*t; 459 t = Scalar(0.5)/t; 460 q.x() = (mat.coeff(2,1) - mat.coeff(1,2)) * t; 461 q.y() = (mat.coeff(0,2) - mat.coeff(2,0)) * t; 462 q.z() = (mat.coeff(1,0) - mat.coeff(0,1)) * t; 463 } 464 else 465 { 466 int i = 0; 467 if (mat.coeff(1,1) > mat.coeff(0,0)) 468 i = 1; 469 if (mat.coeff(2,2) > mat.coeff(i,i)) 470 i = 2; 471 int j = (i+1)%3; 472 int k = (j+1)%3; 473 474 t = ei_sqrt(mat.coeff(i,i)-mat.coeff(j,j)-mat.coeff(k,k) + Scalar(1.0)); 475 q.coeffs().coeffRef(i) = Scalar(0.5) * t; 476 t = Scalar(0.5)/t; 477 q.w() = (mat.coeff(k,j)-mat.coeff(j,k))*t; 478 q.coeffs().coeffRef(j) = (mat.coeff(j,i)+mat.coeff(i,j))*t; 479 q.coeffs().coeffRef(k) = (mat.coeff(k,i)+mat.coeff(i,k))*t; 480 } 481 } 482}; 483 484// set from a vector of coefficients assumed to be a quaternion 485template<typename Other> 486struct ei_quaternion_assign_impl<Other,4,1> 487{ 488 typedef typename Other::Scalar Scalar; 489 static inline void run(Quaternion<Scalar>& q, const Other& vec) 490 { 491 q.coeffs() = vec; 492 } 493}; 494 495} // end namespace Eigen 496