58 #ifndef ACG_QUATERNION_HH
59 #define ACG_QUATERNION_HH
65 #include "Matrix4x4T.hh"
80 template <
class Scalar>
85 #define W Base::data()[0]
86 #define X Base::data()[1]
87 #define Y Base::data()[2]
88 #define Z Base::data()[3]
99 QuaternionT(Scalar _w=1.0, Scalar _x=0.0, Scalar _y=0.0, Scalar _z=0.0)
100 : Vec4(_w, _x, _y, _z) {}
104 : Vec4(0, _p[0], _p[1], _p[2]) {}
108 : Vec4(_p[0], _p[1], _p[2], _p[3]) {}
114 Scalar theta = 0.5 * _angle;
115 Scalar sin_theta = sin(theta);
117 X = sin_theta * _axis[0];
118 Y = sin_theta * _axis[1];
119 Z = sin_theta * _axis[2];
123 template <
class MatrixT>
134 {
return Quaternion(W, -X, -Y, -Z); }
144 {
return Quaternion(W*_q.W - X*_q.X - Y*_q.Y - Z*_q.Z,
145 W*_q.X + X*_q.W + Y*_q.Z - Z*_q.Y,
146 W*_q.Y - X*_q.Z + Y*_q.W + Z*_q.X,
147 W*_q.Z + X*_q.Y - Y*_q.X + Z*_q.W); }
152 {
return *
this = *
this * _q; }
156 template <
class Vec3T>
159 Quaternion q = *
this * Quaternion(0,_v[0],_v[1],_v[2]) *
conjugate();
160 return Vec3T(q[1], q[2], q[3]);
167 if (fabs(W) > 0.999999)
174 _angle = 2.0 * acos(W);
175 _axis = Vec3(X, Y, Z).normalize();
185 ww(W*W), xx(X*X), yy(Y*Y), zz(Z*Z), wx(W*X),
186 wy(W*Y), wz(W*Z), xy(X*Y), xz(X*Z), yz(Y*Z);
190 m(0,0) = ww + xx - yy - zz;
191 m(1,0) = 2.0*(xy + wz);
192 m(2,0) = 2.0*(xz - wy);
194 m(0,1) = 2.0*(xy - wz);
195 m(1,1) = ww - xx + yy - zz;
196 m(2,1) = 2.0*(yz + wx);
198 m(0,2) = 2.0*(xz + wy);
199 m(1,2) = 2.0*(yz - wx);
200 m(2,2) = ww - xx - yy + zz;
202 m(0,3) = m(1,3) = m(2,3) = m(3,0) = m(3,1) = m(3,2) = 0.0;
238 m(0,0) = W; m(0,1) = -X; m(0,2) = -Y; m(0,3) = -Z;
239 m(1,0) = X; m(1,1) = W; m(1,2) = Z; m(1,3) = -Y;
240 m(2,0) = Y; m(2,1) = -Z; m(2,2) = W; m(2,3) = X;
241 m(3,0) = Z; m(3,1) = Y; m(3,2) = -X; m(3,3) = W;
250 m(0,0) = W; m(0,1) = -X; m(0,2) = -Y; m(0,3) = -Z;
251 m(1,0) = X; m(1,1) = W; m(1,2) = -Z; m(1,3) = Y;
252 m(2,0) = Y; m(2,1) = Z; m(2,2) = W; m(2,3) = -X;
253 m(3,0) = Z; m(3,1) = -Y; m(3,2) = X; m(3,3) = W;
259 template<
class MatrixT>
262 Scalar trace = _rot(0,0) + _rot(1,1) + _rot(2,2);
265 Scalar s = 0.5 / sqrt(trace + 1.0);
267 X = ( _rot(2,1) - _rot(1,2) ) * s;
268 Y = ( _rot(0,2) - _rot(2,0) ) * s;
269 Z = ( _rot(1,0) - _rot(0,1) ) * s;
273 if ( _rot(0,0) > _rot(1,1) && _rot(0,0) > _rot(2,2) )
275 Scalar s = 2.0 * sqrt( 1.0 + _rot(0,0) - _rot(1,1) - _rot(2,2));
276 W = (_rot(2,1) - _rot(1,2) ) / s;
278 Y = (_rot(0,1) + _rot(1,0) ) / s;
279 Z = (_rot(0,2) + _rot(2,0) ) / s;
282 if (_rot(1,1) > _rot(2,2))
284 Scalar s = 2.0 * sqrt( 1.0 + _rot(1,1) - _rot(0,0) - _rot(2,2));
285 W = (_rot(0,2) - _rot(2,0) ) / s;
286 X = (_rot(0,1) + _rot(1,0) ) / s;
288 Z = (_rot(1,2) + _rot(2,1) ) / s;
292 Scalar s = 2.0 * sqrt( 1.0 + _rot(2,2) - _rot(0,0) - _rot(1,1) );
293 W = (_rot(1,0) - _rot(0,1) ) / s;
294 X = (_rot(0,2) + _rot(2,0) ) / s;
295 Y = (_rot(1,2) + _rot(2,1) ) / s;
306 Scalar theta( n.norm());
307 Scalar sin_theta = sin(theta);
308 Scalar cos_theta = cos(theta);
311 n *= sin_theta/theta;
315 return Quaternion( cos_theta, n[0], n[1], n[2]);
324 if( w > 1.0) w = 1.0;
325 else if( w < -1.0) w = -1.0;
327 Scalar theta_half = acos(w);
330 Scalar n_norm( n.norm());
333 n *= theta_half/n_norm;
337 return Quaternion( 0, n[0], n[1], n[2]);
349 std::cerr <<
"quaternion : " << (*this) << std::endl;
350 std::cerr <<
"length : " << this->norm() << std::endl;
351 std::cerr <<
"axis, angle: " << axis <<
", " << angle*180.0/M_PI <<
"\n";
352 std::cerr <<
"rot matrix :\n";
353 std::cerr << m << std::endl;
364 typedef QuaternionT<float> Quaternionf;
365 typedef QuaternionT<double> Quaterniond;
371 #endif // ACG_QUATERNION_HH defined
Quaternion logarithm() const
quaternion logarithm (for unit quaternions)
Namespace providing different geometric functions concerning angles.
void identity()
identity rotation
Quaternion operator*(const Quaternion &_q) const
quaternion * quaternion
Matrix rotation_matrix() const
cast to rotation matrix
Quaternion invert() const
invert quaternion
Matrix left_mult_matrix() const
get matrix for mult from left (q*p = Qp)
QuaternionT(Vec3 _axis, Scalar _angle)
construct from rotation axis and angle (in radians)
QuaternionT(const Vec3 &_p)
construct from 3D point (pure imaginary quaternion)
Vec3T rotate(const Vec3T &_v) const
rotate vector
Quaternion conjugate() const
conjugate quaternion
QuaternionT(const Vec4 &_p)
construct from 4D vector
void axis_angle(Vec3 &_axis, Scalar &_angle) const
get rotation axis and angle (only valid for unit quaternions!)
Quaternion & operator*=(const Quaternion &_q)
quaternion *= quaternion
QuaternionT(Scalar _w=1.0, Scalar _x=0.0, Scalar _y=0.0, Scalar _z=0.0)
construct from 4 Scalars (default)
void init_from_matrix(const MatrixT &_rot)
get quaternion from rotation matrix
QuaternionT(const MatrixT &_rot)
construct from rotation matrix (only valid for rotation matrices!)
Matrix right_mult_matrix() const
get matrix for mult from right (p*q = Qp)
Quaternion exponential() const
quaternion exponential (for unit quaternions)