/Users/barrand/private/dev/softinex/old/inexlib-1.2/inlib/inlib/qrot

Go to the documentation of this file.

// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file inlib.license for terms.

#ifndef inlib_qrot
#define inlib_qrot

// rotation done with quaternion.

#include "a4"
#include "a3"

namespace inlib {

template <class T>
class qrot {
public:
  qrot()
  :m_quat(0,0,0,1)   //zero rotation around the positive Z axis.
  {}
  qrot(const a3::vec<T>& a_axis,T a_radians){
    if(!set_value(a_axis,a_radians)) {} //FIXME : throw
  }
  //qrot(const a4::sqm<T>& a_matrix){
  //  set_value(a_matrix);
  //}
  virtual ~qrot(){}
public:
  qrot(const qrot& a_from)
  :m_quat(a_from.m_quat)
  {}
  qrot& operator=(const qrot& a_from){
    m_quat = a_from.m_quat;
    return *this;
  }
protected:
  qrot(T a_q0,T a_q1,T a_q2,T a_q3)
  :m_quat(a_q0,a_q1,a_q2,a_q3)
  {
    if(!m_quat.normalize()) {} //FIXME throw
  }

public:
  qrot& operator*=(const qrot& a_q) {
    //Multiplies the quaternions.
    //Note that order is important when combining quaternions with the
    //multiplication operator.
    // Formula from <http://www.lboro.ac.uk/departments/ma/gallery/quat/>

    T tx = m_quat.v0();
    T ty = m_quat.v1();
    T tz = m_quat.v2();
    T tw = m_quat.v3();

    T qx = a_q.m_quat.v0();
    T qy = a_q.m_quat.v1();
    T qz = a_q.m_quat.v2();
    T qw = a_q.m_quat.v3();

    m_quat.set_value(qw*tx + qx*tw + qy*tz - qz*ty,
                     qw*ty - qx*tz + qy*tw + qz*tx,
                     qw*tz + qx*ty - qy*tx + qz*tw,
                     qw*tw - qx*tx - qy*ty - qz*tz);
    m_quat.normalize();
    return *this;
  }

  bool operator==(const qrot& a_r) const {
    return m_quat.equal(a_r.m_quat);
  }
  bool operator!=(const qrot& a_r) const {
    return !operator==(a_r);
  }
  qrot operator*(const qrot& a_r) const {
    qrot tmp(*this);
    tmp *= a_r;
    return tmp;
  }

  bool invert(){
    T length = m_quat.length();
    if(length==T()) return false;

    // Optimize by doing 1 div and 4 muls instead of 4 divs.
    T inv = one() / length;

    m_quat.set_value(-m_quat.v0() * inv,
                     -m_quat.v1() * inv,
                     -m_quat.v2() * inv,
                      m_quat.v3() * inv);

    return true;
  }

  bool inverse(qrot& a_r) const {
    //Non-destructively inverses the rotation and returns the result.
    T length = m_quat.length();
    if(length==T()) return false;

    // Optimize by doing 1 div and 4 muls instead of 4 divs.
    T inv = one() / length;
  
    a_r.m_quat.set_value(-m_quat.v0() * inv,
                         -m_quat.v1() * inv,
                         -m_quat.v2() * inv,
                          m_quat.v3() * inv);
  
    return true;
  }


  bool set_value(const a3::vec<T>& a_axis,T a_radians) {
    // Reset rotation with the given axis-of-rotation and rotation angle.
    // Make sure axis is not the null vector when calling this method.
    // From <http://www.automation.hut.fi/~jaro/thesis/hyper/node9.html>.
    if(a_axis.length()==T()) return false;
    m_quat.v3(::cos(a_radians/2));
    T sineval = ::sin(a_radians/2);
    a3::vec<T> a = a_axis;
    a.normalize();
    m_quat.v0(a.v0() * sineval);
    m_quat.v1(a.v1() * sineval);
    m_quat.v2(a.v2() * sineval);
    return true;
  }

  bool value(a3::vec<T>& a_axis,T& a_radians) const {
    //WARNING : can fail.
    if( (m_quat.v3()<minus_one()) || (m_quat.v3()> one()) ){ 
      a_axis.set_value(0,0,1);
      a_radians = 0;
      return false;
    }

    a_radians = ::acos(m_quat.v3()) * 2;
    T sineval = ::sin(a_radians/2);

    if(sineval==T()) { //???
      a_axis.set_value(0,0,1);
      a_radians = 0;
      return false;
    }
    a_axis.set_value(m_quat.v0()/sineval,
                     m_quat.v1()/sineval,
                     m_quat.v2()/sineval);
    return true;
  }

/*
  void set_value(const a4::sqm<T>& a_m) {
    //Set the rotation from the components of the given matrix.
    T scalerow = a_m.value(0,0) + a_m.value(1,1) + a_m.value(2,2);
    if (scalerow > T()) {
      T s = fsqrt(scalerow + a_m.value(3,3));
      m_quat.v3(s * 0.5f);
      s = 0.5f / s;

      m_quat.v0((a_m.value(1,2) - a_m.value(2,1)) * s);
      m_quat.v1((a_m.value(2,0) - a_m.value(0,2)) * s);
      m_quat.v2((a_m.value(0,1) - a_m.value(1,0)) * s);
    } else {
      unsigned int i = 0;
      if (a_m.value(1,1) > a_m.value(0,0)) i = 1;
      if (a_m.value(2,2) > a_m.value(i,i)) i = 2;

      unsigned int j = (i+1)%3;
      unsigned int k = (j+1)%3;

      T s = fsqrt((a_m.value(i,i) - (a_m.value(j,j) + a_m.value(k,k))) + a_m.value(3,3));
  
      m_quat.set_value(i,s * 0.5f);
      s = 0.5f / s;
  
      m_quat.v3((a_m.value(j,k) - a_m.value(k,j)) * s);
      m_quat.set_value(j,(a_m.value(i,j) + a_m.value(j,i)) * s);
      m_quat.set_value(k,(a_m.value(i,k) + a_m.value(k,i)) * s);
    }

    if (a_m.value(3,3)!=1.0f) {
      m_quat.multiply(1.0f/fsqrt(a_m.value(3,3)));
    }
  }
*/

  void value(a4::sqm<T>& matrix) const {
    //Return this rotation in the form of a matrix.
    const T x = m_quat.v0();
    const T y = m_quat.v1();
    const T z = m_quat.v2();
    const T w = m_quat.v3();
    // z = w + x * i + y * j + z * k
  
    // first row :
    matrix.set_value(0,0,w*w + x*x - y*y - z*z);
    matrix.set_value(0,1,2*x*y - 2*w*z);
    matrix.set_value(0,2,2*x*z + 2*w*y);
    matrix.set_value(0,3,0);
  
    // second row :
    matrix.set_value(1,0,2*x*y + 2*w*z);
    matrix.set_value(1,1,w*w - x*x + y*y - z*z);
    matrix.set_value(1,2,2*y*z - 2*w*x);
    matrix.set_value(1,3,0);
  
    // third row :
    matrix.set_value(2,0,2*x*z - 2*w*y);
    matrix.set_value(2,1,2*y*z + 2*w*x);
    matrix.set_value(2,2,w*w - x*x - y*y + z*z);
    matrix.set_value(2,3,0);
  
    // fourth row :
    matrix.set_value(3,0,0);
    matrix.set_value(3,1,0);
    matrix.set_value(3,2,0);
    matrix.set_value(3,3,w*w + x*x + y*y + z*z);
  }

  void mult_vec(const a3::vec<T>& a_in,a3::vec<T>& a_out) const {
    const T x = m_quat.v0();
    const T y = m_quat.v1();
    const T z = m_quat.v2();
    const T w = m_quat.v3();
  
    // first row :
    T v0 =     (w*w + x*x - y*y - z*z) * a_in.v0()
              +        (2*x*y - 2*w*z) * a_in.v1()
              +        (2*x*z + 2*w*y) * a_in.v2();
  
    T v1 =             (2*x*y + 2*w*z) * a_in.v0()
              +(w*w - x*x + y*y - z*z) * a_in.v1()
              +        (2*y*z - 2*w*x) * a_in.v2();
  
    T v2 =             (2*x*z - 2*w*y) * a_in.v0()
              +        (2*y*z + 2*w*x) * a_in.v1()
              +(w*w - x*x - y*y + z*z) * a_in.v2();
  
    a_out.set_value(v0,v1,v2);
  }

public: //for io::streamer
  const a4::vec<T>& quat() const {return m_quat;}
  a4::vec<T>& quat() {return m_quat;}

protected:
  static T one() {return T(1);}
  static T minus_one() {return T(-1);}

protected:
  a4::vec<T> m_quat;

public:
  //NOTE : don't handle a static object because of mem balance.
  //static const qrot<double>& identity() {
  //  static const qrot<double> s_v(0,0,0,1);
  //  return s_v;
  //}

};

}

#endif