In previous lessons, we introduced rotation matrices \( \mathbf{R} \in \mathrm{SO}(3) \), Euler angles, and axis–angle representation. Each has advantages and limitations: Euler angles are compact but suffer from singularities (gimbal lock), while rotation matrices are globally well-defined but redundant (nine parameters with orthonormality constraints).

Unit quaternions provide a smooth, globally nonsingular representation of orientation with four parameters and a single unit-length constraint. As a Lie group isomorphic to the unit 3-sphere \( \mathcal{S}^3 \subset \mathbb{R}^4 \), they are particularly convenient for:

• Interpolating orientations along shortest geodesics on \( \mathcal{S}^3 \),
• Composing rotations via quaternion multiplication,
• Avoiding Euler-angle singularities in control and trajectory generation.

Conceptually, we may view unit quaternions as “rotation points” on \( \mathcal{S}^3 \), with geodesic paths corresponding to constant-speed rotational motion in 3D space. The following diagram summarizes the role of quaternions in a typical robot orientation pipeline.

flowchart TD
  A["Orientation keyframes (R0, R1, ...)"] --> B["Convert each to unit quaternions q0, q1, ..."]
  B --> C["Choose interpolation rule on S3 (Slerp / Nlerp)"]
  C --> D["Evaluate interpolated quaternion q(t) for desired time"]
  D --> E["Convert q(t) to rotation matrix or Euler angles"]
  E --> F["Use as end-effector / joint orientation command"]
        

In this lesson we focus on mathematical structure: quaternion algebra, their interpretation as rotations, and how interpolation on \( \mathcal{S}^3 \) leads to algorithms like Slerp. Detailed conversion formulas between representations are deferred to the next lesson.

A quaternion is a 4-tuple consisting of a scalar part and a 3D vector part:

\[ q = \begin{bmatrix} q_0 \\ \mathbf{q}_v \end{bmatrix} = \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ q_3 \end{bmatrix}, \quad q_0 \in \mathbb{R},\; \mathbf{q}_v = (q_1,q_2,q_3)^\top \in \mathbb{R}^3. \]

The set of all quaternions forms a 4D real vector space with addition and scalar multiplication defined component-wise. The key operation is the Hamilton product (or quaternion multiplication). For two quaternions \( p = (p_0,\mathbf{p}_v) \) and \( q = (q_0,\mathbf{q}_v) \), their product \( p \otimes q \) is

\[ p \otimes q = \begin{bmatrix} p_0 q_0 - \mathbf{p}_v^\top \mathbf{q}_v \\ p_0 \mathbf{q}_v + q_0 \mathbf{p}_v + \mathbf{p}_v \times \mathbf{q}_v \end{bmatrix}. \]

The quaternion conjugate and norm are defined by

\[ \bar{q} = \begin{bmatrix} q_0 \\ -\mathbf{q}_v \end{bmatrix}, \qquad \lVert q \rVert^2 = q_0^2 + \lVert \mathbf{q}_v \rVert^2. \]

For nonzero \( q \), the inverse is \( q^{-1} = \bar{q} / \lVert q \rVert^2 \). A unit quaternion satisfies \( \lVert q \rVert = 1 \) and has inverse \( q^{-1} = \bar{q} \). The set of all unit quaternions is a 3D manifold:

\[ \mathcal{S}^3 = \left\{ q \in \mathbb{R}^4 \;\middle|\; \lVert q \rVert = 1 \right\}. \]

Under the Hamilton product, the set of unit quaternions is a Lie group. This group is a double cover of \( \mathrm{SO}(3) \): each rotation matrix corresponds to exactly two antipodal unit quaternions \( q \) and \( -q \).

Recall that any proper rotation in 3D can be parameterized by an axis–angle pair \( (\mathbf{u},\theta) \), where \( \mathbf{u} \in \mathbb{R}^3 \) is a unit vector and \( \theta \in \mathbb{R} \) is the rotation angle. The corresponding unit quaternion is

\[ q(\mathbf{u},\theta) = \begin{bmatrix} \cos(\theta / 2) \\ \mathbf{u} \,\sin(\theta / 2) \end{bmatrix} = \begin{bmatrix} \cos(\theta / 2) \\ u_x \sin(\theta / 2) \\ u_y \sin(\theta / 2) \\ u_z \sin(\theta / 2) \end{bmatrix}. \]

Proposition. If \( \lVert \mathbf{u} \rVert = 1 \), then \( q(\mathbf{u},\theta) \) is a unit quaternion.

Proof. We compute

\[ \lVert q(\mathbf{u},\theta) \rVert^2 = \cos^2(\theta/2) + \lVert \mathbf{u} \rVert^2 \sin^2(\theta/2) = \cos^2(\theta/2) + \sin^2(\theta/2) = 1. \]

Thus, axis–angle and unit quaternions are equivalent ways of describing elements of \( \mathrm{SO}(3) \), up to the sign ambiguity \( q \sim -q \).

To rotate a vector \( \mathbf{x} \in \mathbb{R}^3 \), we embed it into a pure quaternion

\[ x^q = \begin{bmatrix} 0 \\ \mathbf{x} \end{bmatrix}, \]

and define the rotated vector via conjugation:

\[ x^{q'} = q \otimes x^q \otimes q^{-1}, \qquad \mathbf{x}' = \operatorname{vec}(x^{q'}), \]

where \( \operatorname{vec}(\cdot) \) extracts the 3D vector part of the quaternion. A classical identity, obtained by expanding the Hamilton products, is

\[ \mathbf{x}' = \mathbf{x} + 2 q_0 (\mathbf{q}_v \times \mathbf{x}) + 2 \mathbf{q}_v \times (\mathbf{q}_v \times \mathbf{x}), \]

which coincides with the action of the corresponding rotation matrix \( \mathbf{R}(q) \) on \( \mathbf{x} \).

Double cover. Let \( q \) be a unit quaternion. Then

\[ (-q) \otimes x^q \otimes (-q)^{-1} = (-q) \otimes x^q \otimes (-\bar{q}) = q \otimes x^q \otimes \bar{q}, \]

so \( q \) and \( -q \) generate the same rotation on \( \mathbf{x} \). Hence each \( \mathbf{R} \in \mathrm{SO}(3) \) corresponds to two antipodal points on \( \mathcal{S}^3 \).

Endowing \( \mathcal{S}^3 \) with the metric inherited from \( \mathbb{R}^4 \), the geodesic (shortest path) between two unit quaternions \( q_0,q_1 \in \mathcal{S}^3 \) lies on the great circle that they span. The angle between them is

\[ \theta = \arccos\left( q_0^\top q_1 \right), \]

where \( q_0^\top q_1 \) is the Euclidean inner product in \( \mathbb{R}^4 \). Because \( q \sim -q \) represent the same rotation, we always choose the representative so that \( q_0^\top q_1 \ge 0 \), ensuring that \( \theta \in [0,\pi] \) corresponds to the shortest arc.

For robotics, this angle is directly related to the physical relative rotation between orientations. In fact, the relative rotation \( q_{\mathrm{rel}} = q_0^{-1} \otimes q_1 \) has an axis–angle representation \( (\mathbf{u},\theta) \), and the physical rotation angle in 3D space is \( \theta \).

Suppose we have two unit quaternions \( q_0,q_1 \in \mathcal{S}^3 \) representing the start and end orientations of a robot end-effector. We want a smooth interpolation \( q(t) \) such that:

• \( q(0) = q_0 \) and \( q(1) = q_1 \),
• \( q(t) \in \mathcal{S}^3 \) for all \( t \in [0,1] \),
• The path is a constant-speed geodesic on \( \mathcal{S}^3 \).

The spherical linear interpolation (Slerp) satisfies these properties. Let \( d = q_0^\top q_1 \) and, if needed, replace \( q_1 \) with \( -q_1 \) so that \( d \ge 0 \). Define \( \theta = \arccos(d) \). Then:

\[ \operatorname{slerp}(q_0,q_1,t) = \frac{\sin((1-t)\theta)}{\sin \theta}\, q_0 + \frac{\sin(t\theta)}{\sin \theta}\, q_1, \quad 0 \le t \le 1. \]

Boundary conditions. For \( t=0 \), the first coefficient becomes 1 and the second 0, so \( q(0)=q_0 \). For \( t=1 \), the roles are reversed and \( q(1)=q_1 \). For intermediate \( t \), the coefficients trace a great-circle arc on \( \mathcal{S}^3 \).

Lie-group view. We can also write Slerp using the exponential and logarithm maps of the unit-quaternion group, analogous to exponential coordinates of \( \mathrm{SO}(3) \):

\[ q_{\mathrm{rel}} = q_0^{-1} \otimes q_1, \qquad q(t) = q_0 \otimes \exp\!\left( t \, \log(q_{\mathrm{rel}}) \right). \]

The logarithm maps \( q_{\mathrm{rel}} \) to a tangent vector at the identity (an element of the Lie algebra), and the scaled exponential maps it back to the group along a geodesic.

Numerically, when \( \theta \) is extremely small, we can approximate Slerp via normalized linear interpolation (Nlerp), discussed in the next section.

flowchart TD
  I["Input: unit quaternions q0, q1, parameter t in [0,1]"] --> N["Normalize q0, q1"]
  N --> S["Compute dot = q0^T q1; if dot < 0, set q1 = -q1"]
  S --> T["theta = arccos(dot)"]
  T --> C["If theta is small: \nuse normalized linear interpolation"]
  T --> F["Else: compute Slerp coefficients \nand q(t) on S3"]
  C --> O["Return interpolated quaternion q(t)"]
  F --> O
        

Slerp is mathematically exact but involves trigonometric functions, which can be relatively expensive in embedded controllers. A simpler alternative is normalized linear interpolation (Nlerp):

\[ \operatorname{nlerp}(q_0,q_1,t) = \frac{(1-t)\,q_0 + t\,q_1}{\lVert (1-t)\,q_0 + t\,q_1 \rVert}, \quad 0 \le t \le 1. \]

Nlerp preserves end points and keeps quaternions on \( \mathcal{S}^3 \) after normalization. However:

• The trajectory is not exactly a geodesic on \( \mathcal{S}^3 \),
• Angular velocity magnitude is not constant along the path.

In many robotic applications (e.g., low-speed manipulator motion), Nlerp can be an acceptable approximation to Slerp. For high-performance tasks (e.g., fast camera motions, dynamic orientation control) Slerp is generally preferred.

In Python, we can either implement quaternion operations from scratch using numpy, or rely on libraries such as scipy.spatial.transform.Rotation or robotics toolboxes. Below is a minimal, self-contained implementation suitable for FK/trajectory labs in robotics.

import numpy as np

def quat_normalize(q):
    q = np.asarray(q, dtype=float)
    return q / np.linalg.norm(q)

def quat_conjugate(q):
    q = np.asarray(q, dtype=float)
    return np.array([q[0], -q[1], -q[2], -q[3]])

def quat_multiply(q1, q2):
    """Hamilton product q = q1 ⊗ q2."""
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    w = w1*w2 - x1*x2 - y1*y2 - z1*z2
    x = w1*x2 + x1*w2 + y1*z2 - z1*y2
    y = w1*y2 - x1*z2 + y1*w2 + z1*x2
    z = w1*z2 + x1*y2 - y1*x2 + z1*w2
    return np.array([w, x, y, z])

def quat_from_axis_angle(axis, theta):
    """Axis is a 3D vector, theta in radians."""
    axis = np.asarray(axis, dtype=float)
    axis = axis / np.linalg.norm(axis)
    half = 0.5 * theta
    s = np.sin(half)
    return np.array([np.cos(half), axis[0]*s, axis[1]*s, axis[2]*s])

def quat_rotate(q, v):
    """Rotate 3D vector v by unit quaternion q."""
    q = quat_normalize(q)
    v = np.asarray(v, dtype=float)
    vq = np.concatenate([[0.0], v])
    return quat_multiply(quat_multiply(q, vq), quat_conjugate(q))[1:]

def quat_slerp(q0, q1, t):
    """Spherical linear interpolation between q0 and q1."""
    q0 = quat_normalize(q0)
    q1 = quat_normalize(q1)
    dot = float(np.dot(q0, q1))

    # Use the shortest arc on S3
    if dot < 0.0:
        q1 = -q1
        dot = -dot

    # Clamp to handle numerical drift
    dot = np.clip(dot, -1.0, 1.0)
    theta = np.arccos(dot)

    # If the angle is tiny, fall back to normalized lerp
    if theta < 1e-6:
        q = (1.0 - t) * q0 + t * q1
        return quat_normalize(q)

    sin_theta = np.sin(theta)
    w0 = np.sin((1.0 - t) * theta) / sin_theta
    w1 = np.sin(t * theta) / sin_theta
    q = w0 * q0 + w1 * q1
    return quat_normalize(q)

# Example: interpolate between two axis-angle orientations
if __name__ == "__main__":
    q0 = quat_from_axis_angle([0, 0, 1], 0.0)
    q1 = quat_from_axis_angle([0, 0, 1], np.deg2rad(90.0))
    t = 0.5
    qm = quat_slerp(q0, q1, t)
    print("Mid quaternion:", qm)
    v = np.array([1.0, 0.0, 0.0])
    print("Rotated vector at t=0.5:", quat_rotate(qm, v))
      

In practice, you can wrap this module inside a higher-level kinematics library to manage orientation interpolation for trajectories in task space.

In C++, the Eigen library provides both a generic Quaternion class and conversion utilities to and from rotation matrices. Below is a lightweight example illustrating Slerp and vector rotation using Eigen::Quaterniond.

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

using Quaternion = Eigen::Quaterniond;
using Vec3 = Eigen::Vector3d;

Quaternion quatFromAxisAngle(const Vec3 &axis, double theta) {
    Vec3 u = axis.normalized();
    double half = 0.5 * theta;
    double s = std::sin(half);
    return Quaternion(std::cos(half), u.x() * s, u.y() * s, u.z() * s);
}

Quaternion quatSlerp(const Quaternion &q0_in,
                     const Quaternion &q1_in,
                     double t) {
    Quaternion q0 = q0_in.normalized();
    Quaternion q1 = q1_in.normalized();

    double dot = q0.w()*q1.w() + q0.vec().dot(q1.vec());

    // Ensure shortest path
    if (dot < 0.0) {
        q1.coeffs() = -q1.coeffs();
        dot = -dot;
    }

    dot = std::max(-1.0, std::min(1.0, dot));
    double theta = std::acos(dot);

    if (theta < 1e-6) {
        Quaternion q = Quaternion(
            (1.0 - t) * q0.w() + t * q1.w(),
            (1.0 - t) * q0.x() + t * q1.x(),
            (1.0 - t) * q0.y() + t * q1.y(),
            (1.0 - t) * q0.z() + t * q1.z()
        );
        return q.normalized();
    }

    double sin_theta = std::sin(theta);
    double w0 = std::sin((1.0 - t) * theta) / sin_theta;
    double w1 = std::sin(t * theta) / sin_theta;

    return Quaternion(
        w0 * q0.w() + w1 * q1.w(),
        w0 * q0.x() + w1 * q1.x(),
        w0 * q0.y() + w1 * q1.y(),
        w0 * q0.z() + w1 * q1.z()
    ).normalized();
}

int main() {
    Quaternion q0 = quatFromAxisAngle(Vec3(0, 0, 1), 0.0);
    Quaternion q1 = quatFromAxisAngle(Vec3(0, 0, 1), M_PI / 2.0);
    double t = 0.5;
    Quaternion qm = quatSlerp(q0, q1, t);

    Vec3 v(1.0, 0.0, 0.0);
    Vec3 v_rot = qm * v;  // Eigen applies rotation via q * v * q^{-1}

    std::cout << "qm = " << qm.coeffs().transpose() << std::endl;
    std::cout << "v_rot = " << v_rot.transpose() << std::endl;
    return 0;
}
      

This pattern is typical in robotics codebases: unit quaternions are used internally while interfaces may expose rotation matrices or Euler angles for compatibility with other subsystems.

Java has fewer standardized robotics libraries than C++ or MATLAB, but quaternions are easy to implement. The following example shows a simple Quaternion class with Slerp. It can be integrated into a kinematics module for Java-based robot simulators.

public final class Quaternion {
    public final double w, x, y, z;

    public Quaternion(double w, double x, double y, double z) {
        this.w = w; this.x = x; this.y = y; this.z = z;
    }

    public static Quaternion fromAxisAngle(double[] axis, double theta) {
        double norm = Math.sqrt(axis[0]*axis[0] +
                                axis[1]*axis[1] +
                                axis[2]*axis[2]);
        double ux = axis[0] / norm;
        double uy = axis[1] / norm;
        double uz = axis[2] / norm;
        double half = 0.5 * theta;
        double s = Math.sin(half);
        return new Quaternion(Math.cos(half), ux*s, uy*s, uz*s);
    }

    public Quaternion normalized() {
        double n = Math.sqrt(w*w + x*x + y*y + z*z);
        return new Quaternion(w / n, x / n, y / n, z / n);
    }

    public Quaternion conjugate() {
        return new Quaternion(w, -x, -y, -z);
    }

    public Quaternion multiply(Quaternion q) {
        double nw = w*q.w - x*q.x - y*q.y - z*q.z;
        double nx = w*q.x + x*q.w + y*q.z - z*q.y;
        double ny = w*q.y - x*q.z + y*q.w + z*q.x;
        double nz = w*q.z + x*q.y - y*q.x + z*q.w;
        return new Quaternion(nw, nx, ny, nz);
    }

    public double dot(Quaternion q) {
        return w*q.w + x*q.x + y*q.y + z*q.z;
    }

    public static Quaternion slerp(Quaternion q0_in,
                                   Quaternion q1_in,
                                   double t) {
        Quaternion q0 = q0_in.normalized();
        Quaternion q1 = q1_in.normalized();
        double dot = q0.dot(q1);

        // Ensure shortest arc
        if (dot < 0.0) {
            q1 = new Quaternion(-q1.w, -q1.x, -q1.y, -q1.z);
            dot = -dot;
        }

        dot = Math.max(-1.0, Math.min(1.0, dot));
        double theta = Math.acos(dot);

        if (theta < 1e-6) {
            // Nlerp fallback
            double w = (1.0 - t)*q0.w + t*q1.w;
            double x = (1.0 - t)*q0.x + t*q1.x;
            double y = (1.0 - t)*q0.y + t*q1.y;
            double z = (1.0 - t)*q0.z + t*q1.z;
            return new Quaternion(w, x, y, z).normalized();
        }

        double sinTheta = Math.sin(theta);
        double w0 = Math.sin((1.0 - t)*theta) / sinTheta;
        double w1 = Math.sin(t*theta) / sinTheta;

        Quaternion q = new Quaternion(
            w0*q0.w + w1*q1.w,
            w0*q0.x + w1*q1.x,
            w0*q0.y + w1*q1.y,
            w0*q0.z + w1*q1.z
        );
        return q.normalized();
    }

    public double[] rotateVector(double[] v) {
        // Convert v into a pure quaternion and apply q v q^{-1}
        Quaternion vq = new Quaternion(0.0, v[0], v[1], v[2]);
        Quaternion q = this.normalized();
        Quaternion res = q.multiply(vq).multiply(q.conjugate());
        return new double[]{res.x, res.y, res.z};
    }
}
      

This implementation can be extended with methods for converting to rotation matrices or Euler angles, which are covered in the next lesson.

MATLAB’s Robotics System Toolbox includes quaternion primitives (e.g., quaternion, rotmat, eul). Here we first implement basic Slerp manually, then comment on Simulink integration.

function q = quat_from_axis_angle(axis, theta)
    axis = axis(:) / norm(axis);
    half = 0.5 * theta;
    s = sin(half);
    q = [cos(half); axis * s];  % [w; x; y; z]
end

function q = quat_normalize(q)
    q = q(:);
    q = q / norm(q);
end

function q = quat_multiply(q1, q2)
    w1 = q1(1); x1 = q1(2); y1 = q1(3); z1 = q1(4);
    w2 = q2(1); x2 = q2(2); y2 = q2(3); z2 = q2(4);
    w = w1*w2 - x1*x2 - y1*y2 - z1*z2;
    x = w1*x2 + x1*w2 + y1*z2 - z1*y2;
    y = w1*y2 - x1*z2 + y1*w2 + z1*x2;
    z = w1*z2 + x1*y2 - y1*x2 + z1*w2;
    q = [w; x; y; z];
end

function v_rot = quat_rotate(q, v)
    q = quat_normalize(q);
    vq = [0; v(:)];
    q_conj = [q(1); -q(2:4)];
    res = quat_multiply(quat_multiply(q, vq), q_conj);
    v_rot = res(2:4);
end

function q = quat_slerp(q0, q1, t)
    q0 = quat_normalize(q0);
    q1 = quat_normalize(q1);
    dot = dot(q0, q1);

    if dot < 0
        q1 = -q1;
        dot = -dot;
    end

    dot = max(-1.0, min(1.0, dot));
    theta = acos(dot);

    if theta < 1e-6
        q = quat_normalize((1 - t)*q0 + t*q1);
        return;
    end

    sin_theta = sin(theta);
    w0 = sin((1 - t)*theta) / sin_theta;
    w1 = sin(t*theta) / sin_theta;
    q = quat_normalize(w0*q0 + w1*q1);
end

% Example script
q0 = quat_from_axis_angle([0;0;1], 0);
q1 = quat_from_axis_angle([0;0;1], pi/2);
t = 0.5;
qm = quat_slerp(q0, q1, t)
v = [1; 0; 0];
v_rot = quat_rotate(qm, v)
      

In Simulink, you can wrap quat_slerp inside a MATLAB Function block and connect it to blocks representing keyframe quaternions and interpolation parameters. Alternatively, the Robotics System Toolbox provides ready-made blocks for converting between quaternions and rotation matrices, which can be composed with custom Slerp logic.

Mathematica has built-in support for quaternions via the kernel package Quaternions` in some versions, but a direct symbolic implementation is straightforward and useful when deriving analytic expressions for robotic kinematics.

quatNormalize[{w_, x_, y_, z_}] := Module[{n},
  n = Sqrt[w^2 + x^2 + y^2 + z^2];
  {w, x, y, z} / n
];

quatFromAxisAngle[{ux_, uy_, uz_}, theta_] := Module[
  {u, half, s},
  u = {ux, uy, uz}/Sqrt[ux^2 + uy^2 + uz^2];
  half = theta/2;
  s = Sin[half];
  {Cos[half], u[[1]] s, u[[2]] s, u[[3]] s}
];

quatMultiply[{w1_, x1_, y1_, z1_}, {w2_, x2_, y2_, z2_}] :=
 Module[{w, x, y, z},
  w = w1 w2 - x1 x2 - y1 y2 - z1 z2;
  x = w1 x2 + x1 w2 + y1 z2 - z1 y2;
  y = w1 y2 - x1 z2 + y1 w2 + z1 x2;
  z = w1 z2 + x1 y2 - y1 x2 + z1 w2;
  {w, x, y, z}
 ];

quatConjugate[{w_, x_, y_, z_}] := {w, -x, -y, -z};

quatRotate[q_, {vx_, vy_, vz_}] := Module[
  {qn, vq, res},
  qn = quatNormalize[q];
  vq = {0, vx, vy, vz};
  res = quatMultiply[quatMultiply[qn, vq], quatConjugate[qn]];
  Rest[res]
];

quatSlerp[q0in_, q1in_, t_] := Module[
  {q0, q1, dot, theta, sinTheta, w0, w1, q},
  q0 = quatNormalize[q0in];
  q1 = quatNormalize[q1in];
  dot = q0.q1;

  If[dot < 0,
    q1 = -q1;
    dot = -dot;
  ];

  dot = Max[-1.0, Min[1.0, dot]];
  theta = ArcCos[dot];

  If[theta < 10^-6,
    q = (1 - t) q0 + t q1;
    Return[quatNormalize[q]];
  ];

  sinTheta = Sin[theta];
  w0 = Sin[(1 - t) theta]/sinTheta;
  w1 = Sin[t theta]/sinTheta;
  q = w0 q0 + w1 q1;
  quatNormalize[q]
];

(* Example usage *)
q0 = quatFromAxisAngle[{0, 0, 1}, 0];
q1 = quatFromAxisAngle[{0, 0, 1}, Pi/2];
qm = quatSlerp[q0, q1, 1/2];
vRot = quatRotate[qm, {1, 0, 0}]
      

Mathematica’s symbolic capabilities make it ideal for verifying properties of quaternion-based interpolation, such as constant angular velocity and equivalence with matrix exponential trajectories.

In this lesson we introduced unit quaternions as a compact, nonsingular representation of 3D orientation. We developed the Hamilton product, conjugation, and norm, and showed how axis–angle rotations map to unit quaternions on \( \mathcal{S}^3 \). We then characterized quaternion distance and geodesics and constructed spherical linear interpolation (Slerp), which yields constant-speed rotational motion along shortest arcs. We compared Slerp with normalized linear interpolation and implemented quaternion operations and interpolation in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

In the next lesson, we will derive explicit conversion formulas between quaternions, rotation matrices, Euler angles, and axis–angle, and discuss numerical issues when switching representations inside a robot kinematics pipeline.