Rosario 3D

My developer blog

Category Archives: Maths

Base quaternion math

Difficult: basic, let’s keep it simple.

\hat{q}=(q_v, q_w) = i q_x + j q_y + k q_z + q_w

A quaternion is a tuple of four number that can be used to represent a rotation. They are composed of number q_w for the real part and a imaginary vector q_v.

The alternatives

You can say, “I already represent rotation with Euler angles and matrix, why I need other maths to do things that I already know?“.

Well, Euler angles are very simple, they only use three numbers but they can give you headaches in some cases, also while summing two rotation is quite simple, rotate a point with Euler angle require expensive trigonometric function, the really bad part about Euler angles is that different combination of angles can give you the same final result and the final result depends on the order in witch you combine the rotation.

Matrix are very generic, a 3×3 matrix can represent rotation, scale and shearing, combining two transformation is more expensive, but rotating a point quite easy. The bad part about rotation matrix is that given a matrix is not easy understand the orientation of the object. Also they use nine numbers, sixteen if you use homogeneous matrices, this can be a problem if you have few memory (for example when you want to pass a lot of parameter to a shader).

Quaternion operation

Quaternion can save you a lot of time, the basic idea is to store a quaternion to represent the orientation of your object and then use them to compute the matrix to send to openGL or if you are using shader you can use quaternion directly.

It can be strange but rotation in 4D are linear, so you can threat quaternion as a simple polynom:

Addition:

\hat{q}+\hat{r}= (q_v+r_v, q_w+r_w) = (i(q_x+r_x)+j(q_y+r_y)+k(q_z+r_z)+q_w+r_w)

Multiplication (to solve the equation remember that i^2=j^2=k^2=-1 and ij=-ji=k, ki=-ik=j, ij=-ji=k ):

\hat{q}\hat{r}=(iq_x+jq_y+kq_z+q_w)(ir_x+jr_y+kr_z+r_w)\newline=i(q_z r_z -q_z r_y + r_w q_z + q_w r_x)\newline+j(q_z r_x -q_x r_z + r_w q_y + q_w r_y)\newline+k(q_x r_y -q_y r_x + r_w q_z + q_w r_z)\newline+ q_w r_w -q_x r_x - q_y r_y - q_z r_z\newline=(q_v\times r_v + r_w q_v + q_w r_v, q_w r_w - q_v \cdot r_v)

While adding to quaternion is not meaningful for rotation, the product of two quaternion represent the composed rotation.

We can define the conjugate like this:

\hat{q}^\star = (q_v, q_w)^\star = (-q_v, q_w)

Conjugate quaternion represent opposite rotation, rotate a point with a quaternion and then rotate again with his conjugate put back in the place.

The norm is defined like the length of the vector.

n(\hat{q}) = \sqrt{\hat{q}\hat{q}^\star} = \sqrt{\hat{q}^\star\hat{q}} = \sqrt{q_v \cdot q_v +q_w^2} = \sqrt{q_x^2 + q_y^2 + q_z^2 + q_q^2}

Now, we would like to define the inverse operator, so we can derived a rotation that multiplied by another give 1 as result. Let’s try. Remember the norm definition? n(\hat{q}) = \sqrt{\hat{q}\hat{q}^\star} I try to square elevate both sides and I get n(\hat{q})^2 = \hat{q}\hat{q}^\star so

\dfrac{\hat{q}\hat{q}^\star }{ n(\hat{q})^2 } \rightarrow \hat{q}^-1 = \dfrac{\hat{q}^\star}{n(\hat{q})^2}

With quaternion compute the inverse is easy as compute a square radix and make a division. Even better is that we will always work with normalized quaternion, so the length will always be one. In this happy case the conjugate is equivalent to the inverse.

Nothing to scare till now, I present only definition, they are quite useless if we didn’t find any practical application, going out with some friend speaking how beautiful quaternion are is pretty lame.

Well as I wrote before, quaternion represent rotation so they can be used to rotate vertexes. If you have a point p and a normalized quaternion \hat{q} you can rotate the point with:

\hat{q} p \hat{q}^-1\newline =(q_v, q_w) p (-q_v, q_w)\newline =(q_v\times p + q_w p, -q_v \cdot p) (-q_v, q_w)\newline

Using quaternions multiplications the real part became:

-q_w (p \cdot q_v) -(q_v \times p + q_w \cdot p) \cdot (-q_v)=-q_w (p \cdot q_v) +q_w (p \cdot q_v) = 0

q_v\times p\cdot q_v is zero cause q_v\times p is perpendicular to q_v

While the imaginary part became:

(q_v\times p + q_w\cdot p)\times (-q_v) + (q_w\cdot p + q_v\times p)\cdot q_v + (q_v\cdot p) q_v\newline = q_v\times (q_w P + q_v\times p)+q_w^2 p +q_w (q_v\times p) + (q_v\cdot p) q_v\newline = q_w (q_v\times p) + q_v\times (q_v\times p) + q_w^2 p + q_w (q_v\times p) + (q_v\cdot p) q_v

Lagrange’s formula say a \times (b \times c) = b (a \cdot c) - c (a \cdot b) so we can simplify:

2 q_w (q_v\times p) + q_v (q_v\cdot p) - p (q_v\cdot q_v) + q_w^2 p +  (q_v\cdot p) q_v\newline = 2 q_w (q_v\times p) + 2 q_v (q_v\cdot p) + p  (q_w^2 - q_v\cdot q_v)

There are the coordinate of the rotated point.

If you want to rotate your object by an \theta angle around an axis v you can create create your quaternion like this:

\hat{q}= (\sin\left(\dfrac{\theta}{2}\right)v, \cos\left(\dfrac{\theta}{2}\right))

So the axis is related to the imaginary part and the angle represent somewhat the real part, and you can combine multiple quaternion to set complex pose. Now that you have quaternion and we can compose them, you can find useful to convert them to a matrix. Given a quaternion \hat{q} = (w x y z) this is the matrix that represent the same rotation:

\left(\begin{array}{ccc}1-2 (y^2+z^2 )&2 (xy-wz )&2 (wy+zx )\\2 (xy+sz )&1-2 (x^2+z^2 )&2 (yz-wx )\\2 (xz-qy )&2 (qx+yz )&1-2 (x^2+y^2 )\end{array}\right)

Quaternion are very fascinating, but they look quite complex, we really need them? Well the best part of quaternion is that you can interpolate them without to much problem and you get a really smooth orientation change. There are different technique to interpolate quaternion with different result, I’ll show them on another post. Stay tuned.

Vector class

Difficulty: medium (a lot of template)

I want to write an article to speak  about quaternion and how beautiful they are to implement orientation and  keyframe interpolation. But before speaking about quaternion is better to start from the base and create a class for vector.

Vector class

When I speak about vector I refer to maths vector, not an array.  We will develop general version of a vector class that can be used to record our geometry data or other multidimensional data.
This implementation is inspired by the nVidia vector class that you can find with the nVidia SDK, is one of the best Vector class implementation I found. Let’s try to make it better.

Requirement:

  1. Speed: Vector class should be very fast, in a 3d program we will use a lot of vector maths and the speed of this class is crucial for the whole program.
  2. Small size: I also want my class to be compact, I will use the vector class to store the vertex attribute, I don’t want any other data in the middle.
  3. Readability: I want to access to the y coordinate with the .y notation

For the first two reason it’s not wise to use virtual function in our vector class. Virtual function are (slightly) slower, but the most important is that using virtual function will insert an hidden pointer that will increase the size considerably.

Let’s start with the first part of the class, first of all I’ll not create a class but a struct, cause I want all my data to be public, yep I can make the data private and then use getX, getY, setX, setY… but why? This will make only code less readable. A vector is a very simple data that don’t have anything to hide.

template<typename T>
struct Vec4Container
{
   static unsigned const len = 4;
   union{
      struct { T x, y, z, w; };
      struct { T r, g, b, a; };
      struct { T s, t, u, v; };
      T data_[len];
   };
};

I let you immagine how the Vec3Container and Vec2Container is be made.
I have used a template so my class is type independent, now I can create a vector of float, int and so on. I have used union, so I can access to my value with the .x .y .z notation or if my vec represent a color I can use .r .g .b .a or stuv for texture coordinate. OpenGL use s, t, r, q but we already used r in rgb.

Here the real vector class, the one that will perform the computation

template<class container, class T>
struct Vec : public container
{
public:
  static size_t size() { return container::len; }
  Vec() { }
  explicit Vec(const T *v) { for(size_t i = 0; i < size(); i++) container::data_[i]=v[i]; }
  explicit Vec(const T &v) { for(size_t i = 0; i < size(); i++) container::data_[i]=v; }
  Vec(const Vec &vv) { for(size_t i = 0; i < size(); i++) container::data_[i]=vv.data_[i]; }
...

This struct derive from container, so we can access the same data. I also added some constructor,  the first one is the default constructor, I don’t want to initialize my vector if I don’t explicitly needed, then I have a constructor from an array, a constructor from a value. Notice that I have added the explicit keyword, cause I don’t want that a single number can be interpreted as a Vec, this can lead to strange behaviors and confusion. I have added a size operator that return the length of the vector, useful for cycles. Note that this struct don’t add any overhead. The only problem is that cause we are using template C++ don’t know where data_ came from, and we must specify it every time.

We need one more constructor, is very useful to specify a vector like this b(4, 3) or w(1, 3, 5), we need some other classes

template<class T>
struct Vec3 : public Vec<Vec3Container<T> >
{
   Vec3() {};
   explicit Vec3(const T &v) : Vec<Vec3Container<T>, T>(v){ }
   explicit Vec3(const T *v) : Vec<Vec3Container<T>, T>(v){ }
   Vec3(const T& v0, const T& v1, const T& v2){
      Vec<Vec3Container<T>, T>::x = v0;
      Vec<Vec3Container<T>, T>::y = v1;
      Vec<Vec3Container<T>, T>::z = v2;
   }
};

this class will declare only the constructors. The Vec2 and Vec4 are very similar.

Now we can define some basic operation for Vec:

...
   T& operator[] (int i){
      return container::data_[i];
   }
   const T& operator[] (int i) const {
      return container::data_[i];
   }
   Vec& operator +=(const Vec &rop){
      for(size_t i = 0; i < size(); i++)
         container::data_[i]+= rop.data_[i];
      return *this;
   }
   Vec& operator -=(const Vec &rop){
      for(size_t i = 0; i < size(); i++)
         container::data_[i]-= rop.data_[i];
      return *this;
   }
...

There is no problem with sum or subtraction, but what we have to do with multiplication?
We can define different type of multiplication for vector, the dot product, the cross product, pairwise multiplication or scalar multiplication, every operation deserve a * operator. For clarity sake I’ll use a named method for dot and cross product and I’ll reserve the * operator for pairwise and scalar product.

...
   Vec& operator *=(const T &rop){
      for(size_t i = 0; i < size(); i++)
         container::data_[i]*= rop;
      return *this;
   }
   Vec& operator *=(const Vec &rop){
      for(size_t i = 0; i < size(); i++)
         container::data_[i]*= rop.data_[i];
      return *this;
   }
...

Dot product is a binary operation, so it’s better to make it extern to the class, same thing for cross product.

template<class container, typename T>
T dot(const Vec<container, T>& rOp, const Vec<container, T>& lOp)
{
  T ret(0);
  for(size_t i = 0; i < container::len; i++)
     ret += rOp[i] * lOp[i];
  return ret;
}

We can specialize the cross funtion to make is more efficient

template<typename T>
Vec3<T> cross(const Vec3<T>& rOp, const Vec3<T&>& lOp)
{
    Vec3 res;
    res.x = rOp.y * lOp.z - rOp.z * lOp.y;
    res.y = rOp.x * lOp.z - rOp.z * lOp.x;
    res.z = rOp.x * lOp.y - rOp.y * lOp.x;
    return res;
}

Other function like normalize, length and binary operator are trivial, you can implement them as exercise. 😛
Final touch:

typedef Vec3 Vec3f;
typedef Vec3 Vec3i;
typedef Vec3 Vec3ui;
typedef Vec2 Vec2f;
typedef Vec2 Vec2i;
...

Edit: I notice that a lot of compiler don’t like template of template, so I modify the code a little bit.