Quaternion Initialization and Basic Arithmetic

In computer graphics and robotics, quaternions are the preferred method for handling 3D rotations due to their robustness against gimbal lock and computational efficiency. The Eigen library provides a comprehensive Quaternion class template to faiclitate these operations.

To begin working with quaternions, include the appropriate module and declare instances using the template class Quaternion<typename Scalar>. Common typedefs include Quaternionf for float and Quaterniond for double precision.

#include <iostream>
#include <Eigen/Dense>
#include <cmath>

int main() {
    // Initialize quaternions (w, x, y, z)
    Eigen::Quaternionf qA(1.0f, 0.0f, 0.0f, 0.0f);
    Eigen::Quaternionf qB(0.5f, 0.5f, 0.5f, 0.5f);

    std::cout << "Quaternion A: " << qA.coeffs().transpose() << std::endl;
    std::cout << "Quaternion B: " << qB.coeffs().transpose() << std::endl;

    return 0;
}

Component-wise Addition and Scalar Scaling

Unlike multiplication, quaternion addition is rarely used in rotation contexts but is mathematical valid as a vector addition in R4. Scalar scaling multiplies every component of the quaternion by a scalar value.

struct QuaternionUtils {
    static Eigen::Quaternionf add(const Eigen::Quaternionf& a, const Eigen::Quaternionf& b) {
        // Add components directly
        return Eigen::Quaternionf(
            a.w() + b.w(), 
            a.x() + b.x(), 
            a.y() + b.y(), 
            a.z() + b.z()
        );
    }

    static Eigen::Quaternionf scale(const Eigen::Quaternionf& q, float scalar) {
        return Eigen::Quaternionf(
            q.w() * scalar, 
            q.x() * scalar, 
            q.y() * scalar, 
            q.z() * scalar
        );
    }
};

Quaternion Multiplication

Quaternion multiplication (the Hamilton product) is the core operation for combining rotations. While Eigen overloads the * operator, understanding the underlying math is useful for debugging and optimization.

Method 1: Native Eigen Operator

The most efficient method is using the library's built-in operator, which utilizes SIMD optimizations where available.

// Native multiplication
Eigen::Quaternionf productNative = qA * qB;

Method 2: Explicit Component Formula

To manually implement the multiplication, we calculate the scalar and vector parts. Given $q_1 = (w_1, \vec{v}_1)$ and $q_2 = (w_2, \vec{v}_2)$, the product $q_{res}$ is defined as:

• $w_{res} = w_1 w_2 - \vec{v}_1 \cdot \vec{v}_2$
• $\vec{v}_{res} = w_1 \vec{v}_2 + w_2 \vec{v}_1 + \vec{v}_1 \times \vec{v}_2$

We can implement this explicitly using component accessors:

Eigen::Quaternionf manualMultiplication(const Eigen::Quaternionf& a, const Eigen::Quaternionf& b) {
    float newW = a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z();
    float newX = a.w() * b.x() + a.x() * b.w() + a.y() * b.z() - a.z() * b.y();
    float newY = a.w() * b.y() - a.x() * b.z() + a.y() * b.w() + a.z() * b.x();
    float newZ = a.w() * b.z() + a.x() * b.y() - a.y() * b.x() + a.z() * b.w();

    return Eigen::Quaternionf(newW, newX, newY, newZ);
}

Applying Quaternion Rotation to 3D Points

To rotate a 3D point $p$ using a quaternion $q$, we represent the point as a "pure" quaternion $p_q = (0, \vec{p})$. The rotated point $p'$ is calculated using the sandwich product:

$p' = q \cdot p_q \cdot q^{-1}$

For unit quaternions (which represent rotation), the inverse $q^{-1}$ is equivalent to the conjugate $q^*$.

The following example demonstrates rotating a vector by 90 degrees around the Z-axis. We will construct the rotation quaternion from an axis and an angle, apply the transformation, and verify the result.

#include <iostream>
#include <Eigen/Dense>
#include <cmath>

int main() {
    // 1. Define the point to rotate (e.g., along the X-axis)
    Eigen::Vector3f point(1.0f, 0.0f, 0.0f);
    std::cout << "Original Point: " << point.transpose() << std::endl;

    // 2. Create a pure quaternion from the point (w=0)
    Eigen::Quaternionf pQuat;
    pQuat.w() = 0.0f;
    pQuat.vec() = point;

    // 3. Define the rotation: 90 degrees (PI/2) around Z-axis (0, 0, 1)
    float angle = M_PI / 2.0f;
    float halfAngle = angle / 2.0f;
    float s = std::sin(halfAngle);
    float c = std::cos(halfAngle);

    Eigen::Vector3f axis(0.0f, 0.0f, 1.0f);
    Eigen::Quaternionf rotation(c, axis.x() * s, axis.y() * s, axis.z() * s);
    
    // Normalize to ensure it is a valid rotation quaternion
    rotation.normalize();

    std::cout << "Rotation Quaternion: " << rotation.coeffs().transpose() << std::endl;

    // 4. Perform the rotation: p' = q * p * q.inverse()
    Eigen::Quaternionf rotatedQuat = rotation * pQuat * rotation.inverse();

    // 5. Extract the rotated vector
    Eigen::Vector3f rotatedPoint = rotatedQuat.vec();
    std::cout << "Rotated Point:   " << rotatedPoint.transpose() << std::endl;

    // Logic Check: Rotating (1,0,0) by 90 deg around Z should yield (0,1,0)
    // Due to floating point precision, values might be very close to 0.0 or 1.0
    
    return 0;
}