Harold Serrano

View Original

Developing a Math Engine in C++: Implementing Quaternions

Harold Serrano

Quaternions are used in computer graphics when a 3D character rotation is involved. Quaternions allows a character to rotate about multiple axis simultaneously instead of sequentially, as matrix rotation allows.

Why use Quaternions to rotate a 3D character when matrices can do the same job? There are two reasons why Quaternions are preferred in computer graphics:

  1. Matrix rotations suffer from what is known as Gimbal Lock.
  2. Quaternions consume less memory and are faster to compute than matrices.

Gimbal lock is the loss of one degree of freedom in a three-dimensional, three-gimbal mechanism that occurs when the axes of two of the three gimbals are driven into a parallel configuration, "locking" the system into rotation in a degenerate two-dimensional space.

At first, I was thinking to simply develop our math engine with vectors and matrices. But decided to add a Quaternion class to our engine as well. Quaternions, although difficult to understand, are simple to use.

So, let's implement a Quaternion class in the math engine.

Quaternions

Quaternions are composed of a scalar and a vector. They are represented in various different ways such as:

where S is a scalar number and V is a vector representing an axis.

Let's implement a Quaternion class. Download the math engine and create a new C++ class file. Call it R4DQuaternion. Since we are creating a C++ class in an iOS environment, change the .hpp and .cpp file to .h and .mm, respectively.

In the header file, declare the constructors and data members of the quaternion.

namespace R4DEngine {

    class R4DQuaternion{

    public:
        //scalar
        float s;

        //vector
        R4DVector3n v;

        //Constructor
        R4DQuaternion();

        //Constructor
        R4DQuaternion(float uS,R4DVector3n& uV);

        //Destructor
        ~R4DQuaternion();

        //Copy Constructor
        R4DQuaternion(const R4DQuaternion & value);

        //Copy Constructor
        R4DQuaternion& operator=(const R4DQuaternion& value);
  };
}

In the .mm file implement the following construction methods:

R4DQuaternion::R4DQuaternion(float uS,R4DVector3n& uV):s(uS),v(uV){}

R4DQuaternion::R4DQuaternion(const R4DQuaternion & value){

    s=value.s;
    v=value.v;

};

R4DQuaternion& R4DQuaternion::operator=(const R4DQuaternion& value){

    s=value.s;
    v=value.v;

    return *this;
};

To create an instance of a quaternion in the engine we do the following:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create a quaternion
R4DEngine::R4DQuaternion q(90,axis);

Adding and Subtracting

Quaternions can be added and subtracted among themselves. In quaternion addition and subtraction, the corresponding scalar and vectors from each quaternion are added/subtracted, respectively. Mathematically, quaternion addition/subtraction is represented as follows:

For example, the addition of two quaternions is calculated as follows:

Addition and subtraction is implemented in C++ as follows:

void R4DQuaternion::operator+=(const R4DQuaternion& q){

    s+=q.s;
    v+=q.v;

}

R4DQuaternion R4DQuaternion::operator+(const R4DQuaternion& q)const{

    float scalar=s+q.s;
    R4DVector3n imaginary=v+q.v;

    return R4DQuaternion(scalar,imaginary);
}

//substract
void R4DQuaternion::operator-=(const R4DQuaternion& q){

    s-=q.s;
    v-=q.v;
}

R4DQuaternion R4DQuaternion::operator-(const R4DQuaternion& q)const{

    float scalar=s-q.s;
    R4DVector3n imaginary=v-q.v;

    return R4DQuaternion(scalar,imaginary);
}

If you have been reading my blog for a while, you would have noticed that I overload C++ arithmetic operators often. I do this, because it makes working with the math engine quite simple. Read more about it here.

In the engine we can add and subtract quaternions as follows:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create two quaternions
R4DEngine::R4DQuaternion q(90,axis);
R4DEngine::R4DQuaternion p(10,axis);

//Add quaternions
R4DEngine::R4DQuaternion nQ=q+p;

//Add quaternions
q+=p;

Product

Quaternions can be multiplied among themselves. Mathematically, quaternion multiplication is represented as follows:

Quaternion multiplication is implemented as follows:

//multiply
void R4DQuaternion::operator*=(const R4DQuaternion& q){

    (*this)=multiply(q);
}

R4DQuaternion R4DQuaternion::operator*(const R4DQuaternion& q)const{

    float scalar=s*q.s - v.dot(q.v);

    R4DVector3n imaginary=q.v*s + v*q.s + v.cross(q.v);

    return R4DQuaternion(scalar,imaginary);

}

R4DQuaternion R4DQuaternion::multiply(const R4DQuaternion& q)const{

    float scalar=s*q.s - v.dot(q.v);

    R4DVector3n imaginary=q.v*s + v*q.s + v.cross(q.v);

    return R4DQuaternion(scalar,imaginary);

}

In the engine, we can multiply two quaternions as follows:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create two quaternions
R4DEngine::R4DQuaternion q(90,axis);
R4DEngine::R4DQuaternion p(10,axis);

//Quaternion multiplication
R4DEngine::R4DQuaternion nQ=q*p;

//Quaternion multiplication
q-=p;

//Quaternion multiplication
R4DEngine::R4DQuaternion nP=q.multiply(p);

Scalar Multiplication

Just like vectors and matrices, a quaternion can be multiplied by a scalar. Mathematically, this is represented as follows:

where k is a scalar real number.

Quaternion scalar multiplication is done as follows in C++:

void R4DQuaternion::operator*=(const float value){

    s*=value;
    v*=value;
}


R4DQuaternion R4DQuaternion::operator*(const float value)const{

    float scalar=s*value;
    R4DVector3n imaginary=v*value;

    return R4DQuaternion(scalar,imaginary);

}

Multiplying a quaternion by a scalar in the engine is done as follows:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create a quaternion
R4DEngine::R4DQuaternion q(90,axis);

//scalar multiplication
R4DEngine::R4DQuaternion nQ=q*3;

//scalar multiplication
q*=3;

Norm

Mathematically, the norm of a quaternion is represented as follows:

For example, the norm of quaterion q is calcuated as follows:

This can be implemented as follows in C++:

//norm
float R4DQuaternion::norm(){

    float scalar=s*s;
    float imaginary=v*v;

    return sqrt(scalar+imaginary);
}

To get a quaternion's magnitude in the engine we do the following:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create a quaternion
R4DEngine::R4DQuaternion q(90,axis);

//Magnitude (Norm)
float magnitude=q.norm();

Unit Norm

The unit norm of a quaternion, also known as Normalized quaternion, is calculated as follows:

For example, normalizing the quaternion q is calculated as follows:

The normalization of a quaternion is implemented as follows:

//unit-norm
void R4DQuaternion::normalize(){

    if (norm()!=0) {

        float normValue=1/norm();

        s*=normValue;
        v*=normValue;
    }
}

To normalize a quaternion in the engine we call the following function:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create a quaternion
R4DEngine::R4DQuaternion q(90,axis);

//Normalize
q.normalize();

The type of unit-norm quaternion we will use for rotations are of the form:

And is implemented as:

//Unit-Norm Quaternion (Special Form)
void R4DQuaternion::convertToUnitNormQuaternion(){

    float angle=DegreesToRad(s);
    v.normalize();
    s=cosf(angle*0.5);
    v=v*sinf(angle*0.5);

}

I will talk more about it later in this article.

Conjugate

The conjugate of a quaternion is very important in computing the inverse of a quaternion. Mathematically, the conjugate of quaternion q is calculated as follows:

The conjugate is implemented as follows:

//conjugate
R4DQuaternion R4DQuaternion::conjugate(){

    float scalar=s;
    R4DVector3n imaginary=v*(-1);

    return R4DQuaternion(scalar,imaginary);
}

A quaternion's conjugate is calculated by calling the following function:

//Create an axis vector
R4DEngine::R4DVector3n axis(1,0,0);

//create a quaternion
R4DEngine::R4DQuaternion q(90,axis);

//Quaternion conjugate
R4DEngine::R4DQuaternion conjugate=q.conjugate();

Inverse

The ability to divide in quaternion arithmetic is very important. To be able to divide a quaternion, you need its inverse. By definition the inverse quaternion satisfies the following:

The inverse of a quaternion is defined as:

The inverse of a quaternion is done as follows:

//inverse
R4DQuaternion R4DQuaternion::inverse(){

    float absoluteValue=norm();
    absoluteValue*=absoluteValue;
    absoluteValue=1/absoluteValue;

    R4DQuaternion conjugateValue=conjugate();

    float scalar=conjugateValue.s*absoluteValue;
    R4DVector3n imaginary=conjugateValue.v*absoluteValue;

    return R4DQuaternion(scalar,imaginary);
}

To get the inverse of a quatertion in the engine, we do the following:

//create a quaternion
R4DEngine::R4DQuaternion q(90,axis);

//Quaternion inverse
R4DEngine::R4DQuaternion conjugate=q.inverse();

Rotating a Vector

There are two types of quaternions:

  • Real Quaternions
  • Pure Quaternions

A real quaternion has a Zero vector term. Think of them as real numbers:

A pure quaternion has a Zero scalar term. Think of them as vectors:

A vector rotation occurs when you multiply a pure quaternion p by a unit-norm quaternion q and its inverse:

Where p represents your vector in quaternion form. And q represents a unit-norm quaternion in the form of:

Rotating a vector using quaternions is implemented as follows:

R4DVector3n R4DVector3n::rotateVectorAboutAngleAndAxis(float uAngle, R4DVector3n& uAxis){

    //convert our vector to a pure quaternion
    R4DQuaternion p(0,(*this));

    //normalize the axis
    uAxis.normalize();

    //create the real quaternion
    R4DQuaternion q(uAngle,uAxis);

    //convert quaternion to unit norm quaternion
    q.convertToUnitNormQuaternion();

    //Get the inverse of the quaternion
    R4DQuaternion qInverse=q.inverse();

    //rotate the quaternion
    R4DQuaternion rotatedVector=q*p*qInverse;

    //return the vector part of the quaternion
    return rotatedVector.v;

}

Notice that we convert the vector into a pure quaternion and convert the real quaternion into the special unit-norm form.

Let's do a simple test. Let's rotate a vector about the x-axis by 90 degrees.

In main.mm test the following:

int main(int argc, const char * argv[]) {

    //Vector to rotate
    R4DEngine::R4DVector3n v(0,1,0);

    //axis of rotation: x-axis
    R4DEngine::R4DVector3n axis(1,0,0);

    //rotate vector v about axis by 90 degrees
    R4DEngine::R4DVector3n rotatedVector=v.rotateVectorAboutAngleAndAxis(90, axis);

    //print result
    rotatedVector.show();

    return 0;
}

Build the project using Command+B. The result should show the vector:

(0,1,0)

Being rotated to:

(0,0,1)

As shown here:

Source Code

The complete source code can be found here. Feel free to download it and play with it.