\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{epsfig}

\textheight= 8in \textwidth=5.5in
\parindent =0pt
\hsize = 5in
%\voffset = 1 in
%\hoffset= 1 in
\parindent = 0pt
\parskip = 0.2cm
\tolerance = 10000
\input{../../newcomm}

\begin{document}
\bibliographystyle{plain}

\title{Quaternions}
    \author {Gerald Farin}
\maketitle
%**************************
\section{Motivation}
%**************************
For animation, objects have to be moved through space. If the
object's orientation changes during such a motion, a rotation is
invariably involved. A typical problem is that the object is only
defined at several key positions, and the intermediate ones have
to be computed. In other words, we have to be able to interpolate
between rotations. This is not trivial: if $R_1$ and $R_2$ are two
rotation matrices, then the blend
$$
    Q=(1-t)R_1 + tR_2
$$
will, in general, not be a rotation matrix since we cannot
guarantee that $Q^{\rm T}Q=I$ as is required for rotation
matrices. A different entity, called a quaternion, may also be
used to describe rotations, but in addition is capable of handling
the interpolation problem.
%***************************
\section{Introduction}
%***************************
Let $R$ be rotation matrix, such that
$$
    R^TR=I.
$$
First we show that $R$ has an eigenvalue 1, i.e., there is a unit
vector $\uu$, corresponding to this eigenvalue,  such that
$$
    A\uu=\uu.
$$
Assume $\lambda$ is the eigenvector corresponding to $\uu$. Then
$$
    (\lambda \uu)^{\rm T}(\lambda \uu)=\lambda^2 \uu^{\rm T}\uu =
    \uu^{\rm T}R^{\rm T}R\uu = \uu^{\rm T}\uu = 1,
$$
hence $\lambda^2 = \pm 1$. To show that in fact $+1$ is an
eigenvector requires some more work.

Thus in order to find the axis of a general rotation matrix, all
we have to do is to find its eigenvector corresponding to the
eigenvalue $+1$. The angle $\Phi$ by which $\vv$ is rotated around
$\uu$ to yield $\vv'$ can be computed from
$$
     {\rm trace}(R)=1+2\cos\Phi.
$$
\hrule {\bf Example}~~ Let a rotation matrix be given by
$$
    R=\begin{bmatrix}
        \tau^2 & \tau^2-\tau & \tau^2 + \tau \\
        \tau^2 + \tau & \tau^2 & \tau^2 - \tau \\
        \tau^2 - \tau & \tau^2 + \tau & \tau^2
        \end{bmatrix}
$$
with $\tau = 1/\sqrt{3}$. The matrix has an eigenvector
$\uu=[\tau,\tau,\tau]^{\rm T}$ corresponding to the eigenvalue 1:
$$
\begin{bmatrix}
        \tau^2 & \tau^2-\tau & \tau^2 + \tau \\
        \tau^2 + \tau & \tau^2 & \tau^2 - \tau \\
        \tau^2 - \tau & \tau^2 + \tau & \tau^2
        \end{bmatrix}
        \begin{bmatrix} \tau \\ \tau  \\ \tau \end{bmatrix}
            =
        \begin{bmatrix}
            \tau^3 + \tau^3 - \tau^2 + \tau^3 + \tau^2\\
            \tau^3 + \tau^2 + \tau^3 + \tau^3 -\tau^2 \\
            \tau^3 -\tau^2 + \tau^3 + \tau^2 + \tau^3
        \end{bmatrix}
            =
        \begin{bmatrix} \tau \\ \tau  \\ \tau \end{bmatrix}
$$
We have ${\rm trace}(R) = 3 \tau^2 = 1$. Thus $\cos\Phi = 0$ and
thus $\Phi = 90$.
 \hrule

If $\vv' = R \vv$, then $\vv'$ is in the plane which contains
$\vv$ and has $\uu$ as its normal. To see this, note that
$(\uu^{\rm T}\vv)\uu$ is the perpendicular projection of $\vv$
onto $\uu$, marked by a solid circle in Figure \ref{fig:rotv}. If
we can show that the projection of $\vv'$ onto $\uu$ is the same
vector, our claim is proved. That projection is given by
$(\uu^{\rm T}\vv')\uu$. But (noting that $R\uu=\uu$):
$$
    (\uu^{\rm T}\vv')\uu = (\uu^{\rm T}R^{\rm T}R\vv)\uu =(\uu^{\rm
    T}\vv)\uu,
$$
thus proving our claim.

 \bef
    \epsfig{file=rotate.ps, height = 2 in}
    \caption{The vector $\vv$ is rotated into $\vv'$ around axis $\uu$.}
   \label{fig:rotv}
\eef
 Conversely, we may ask how to compute the rotation matrix $R$
if we are given the axis vector $\uu$ and the rotation angle
$\Phi$. This matrix turns out to be (check any linear algebra
text)
 \beq
    \begin{bmatrix}
    u_1^2 + c(1-u_1^2) & (1-c)u_1u_2-su_3   &  (1-c)u_1u_3 + su_2\\
    (1-c)u_1u_2+su_3   & u_2^2+c(1-u_2^2)   &  (1-c)u_2u_3-su_1  \\
    (1-c)u_1u_3-su_2   & (1-c)u_2u_3 + su_1 &  u_3^2+c(1-u_3^2)
    \end{bmatrix},
    \label{eq:rotmat}
\eeq
 where $c=\cos\Phi, s=\sin\Phi$. For a sanity check, verify
that ${\rm trace}(R)=1+2c$!

The essential information defining a rotation is given by $\uu$
and $\Phi$. This information is collected in the {\em quaternion}
\beq
    {\hat \uu} =
    \begin{bmatrix}
        s\uu \\
            c
    \end{bmatrix}
        \label{eq:quat} \eeq

Clearly we can rebuild the matrix $R$ from this information. But
instead of working with matrices, some applications related to
animation are formulated much more elegantly using quaternions.

%*******************************
\section{Definitions}\label{sec:intro}
%*******************************
Since quaternions were originally invented as a generalization of
complex numbers, we now turn to a formulation which is more
general than the above. Now, a quaternion is defined as
$$
    \hqq=\begin{bmatrix}
    \xx\\
    w
    \end{bmatrix}
$$
where
$$
    \xx = ix + jy + kz
$$
Note that on one hand, we write $\xx$ as a vector, on the other
hand as a generalized complex number. It's not a totally bad idea
to think of $i,j,k$ as the unit direction vectors of $\eee$, but
they are better thought of as generalized complex numbers:
$$
    i^2=j^2=k^2=-1; ~~ ij=k, jk=i, ki=j.
$$
We also have $ij=-ji$ etc.

The product $\hqq_1\hqq_2$ of two quaternions may again expressed
as a quaternion. It is found by multiplying
$$
    (ix_1 + jy_1 + kz_1 + w_1)(ix_2 + jy_2 + kz_2 + w_2)
$$
which gives
$$
    \hqq_1\hqq_2=
    \begin{bmatrix}
        w_1 x_2 +x_1w_2 + y_1 z_2 - z_1 y_2\\
        w_1 y_2 -x_1 z_2 + y_1 w_2 + z_1 x_2\\
        w_1 z_2 + x_1 y_2 -y_1 x_2 + z_1 w_2\\
        w_1w_2 - x_1x_2 -y_1 y_2 - z_1z_2
    \end{bmatrix}.
$$
It will be helpful to rewrite this in matrix form:
$$
   \hqq_1\hqq_2=
    \begin{bmatrix}
        w_1 & -z_1 & y_1 & x_1\\
        z_1 & w_1 & -x_1 & y_1\\
        -y_1& x_1 & w_1 & z_1 \\
        -x_1 & -y_1 & -z_1 & w_1
    \end{bmatrix}
    \begin{bmatrix} x_2 \\ y_2 \\ z_2 \\ w_2 \end{bmatrix},
$$
which we abbreviate as
\beq
    \hqq_1\hqq_2= M_1\hqq_2.
    \label{eq:mult1}
\eeq
 There is a second possibility:
$$
   \hqq_1\hqq_2=
    \begin{bmatrix}
        w_2 & z_2 & -y_2 & x_2\\
        -z_2 & w_2 & x_2 & y_2 \\
        y_2 & -x_2 & w_2 & z_2 \\
        -x_2 & -y_2 & -z_2 & w_2
    \end{bmatrix}
    \begin{bmatrix} x_1 \\ y_1 \\ z_1 \\ w_1 \end{bmatrix}
$$
which we abbreviate to
\beq
    \hqq_1\hqq_2= M_2\hqq_1.
    \label{eq:mult2}
\eeq

This multiplication is clearly not commutative!

Yet a third possibility for the computation of $\hqq_1\hqq_2$ is
given by
$$
    \hqq_1\hqq_2=\begin{bmatrix}
                    \xx_1 \wedge \xx_2 + w_2\xx_1 + w_1\xx_2 \\
                        w_1w_2 - \xx_1\xx_2
                \end{bmatrix}
$$

 Other operations
on quaternions are most easily understood by knowing they are
equivalent to matrices. For example, addition is defined by
$$
    \hqq + \hpp =
    \begin{bmatrix}
        \qq_v + \pp_v\\
        q_w + p_w
    \end{bmatrix}.
$$
Another property is linearity:
$$
    \hqq(\alpha \hpp_1 + \beta \pp_2) = \alpha \hqq\hpp_1 + \beta
    \hqq\hpp_2.
$$


 Since the $i,j,k$ parts are imaginary numbers, we may
define the complex conjugate $\hqq^*$ of a quaternion $\hqq$.
First, recall the standard complex conjugate of a complex number:
$$
    (a+ib)^*=a-ib.
$$

For quaternions, the same pattern holds:
$$
    \hqq^*=\begin{bmatrix}
        -\xx\\
        w
        \end{bmatrix}.
$$
The {\em norm} (length) of a complex number is given by
$|a+ib|^2=(a+ib)(a+ib)^*=(a+ib)(a-ib)=a^2-aib+aib-i^2b=a^2+b^2$.
The same holds for quaternions:
$$
    \| \hqq \|^2 =\hqq \hqq^*.
$$
Note that $\hqq\hqq^*$ yields a real number because of the
properties of the imaginary units in quaternions.

 A {\em unit quaternion} is one with $\| \hqq \|^2 =1$ and
ties in with the quaternions introduced in Section
\ref{sec:intro}. It is thus of the form of (\ref{eq:unitquat}):
\beq
    \hqq=\begin{bmatrix}
        \sin \Phi \uu\\
        \cos\Phi
        \end{bmatrix}
        \label{eq:unitquat}
\eeq
 We see that $\hqq\hqq^* = s^2(u_1^2 + u_2^2 + u_3^2) + c^2 =
1$ since we assumed that $\uu$ is a unit vector. In what follows,
we will only deal with unit quaternions.

Finally, the identity quaternion:
$$
    {\hat \ii}=\begin{bmatrix}
            {\bf 0}\\
                1
            \end{bmatrix}.
$$
It is often convenient to simply treat this quaternion as the
scalar 1.

We may also define the inverse $\hqq^{-1}$ of $\hqq$, defined by
$$
    \hqq^{-1}\hqq = \hqq \hqq^{-1} = 1.
$$
It is given by
$$
    \hqq^{-1}=\frac{\hqq^*}{\| \hqq \|^2},
$$
and is verified by
$$
    \frac{\hqq^{-1}\hqq^*}{\| \hqq \|^2}=1.
$$

%*****************************************************************
\section{Quaternions and Matrices}
%*****************************************************************

Rotating a vector $\vv$ around an axis $\uu$ by an angle $2\Phi$
may be expressed by $\vv'=R\vv$, where $R$ is given by
(\ref{eq:rotmat}). In terms of quaternions, the corresponding
equation is
$$
    \vv' = \hqq \vv \hqq^{-1} = \hqq \vv \hqq^*
$$
This makes sense since we may write $\vv$ in quaternion form with
a zero real part. This operation corresponds to a rotation by
$2\Phi$ if we set $w=\cos \Phi$.

Invoking (\ref{eq:mult1}) and (\ref{eq:mult2}), we obtain
$$
    \hqq \vv \hqq^* = \hqq(M_2^* \vv) = M_1M_2^*\vv.
$$
Thus $M_1M_2^*$ is the desired rotation matrix $R$. It is
explicitly given by
$$ R =
    \begin{bmatrix}
        w^2 +x^2 -y^2-z^2 & 2xy -2wz & 2xz+2wy & 0\\
        2xy + 2wz& w^2-x^2+y^2-z^2 & 2yz-2wx & 0\\
        2xz-2wy& 2yz+2wx& w^2 - x^2 - y^2 +z^2 & 0\\
        & 0 & 0 & x^2 + y^2 + z^2 +w^2
    \end{bmatrix}
$$

\hrule {\bf Example} A rotation around the $z-$ axis by 90 degrees
is given by the unit quaternion
$$
    \hqq = [0, 0, \frac{1}{2}\sqrt{2}, \frac{1}{2}\sqrt{2}]^{\rm
    T}.
$$
If we use the above to calculate the corresponding rotation
matrix, we find
$$
    R=\begin{bmatrix}
    0 & -1 & 0 & 0\\
    1 & 0  & 0 & 0\\
    0 & 0 &  1 & 0\\
    0 & 0 & 0 & 1
    \end{bmatrix}.
$$
\hrule

We may concatenate quaternions as follows. Transform a vector
$\vv$ using quaternion $\hqq_1$, then transform the result using a
quaternion $\hqq_2$. We have
$$
    \vv' = \hqq_2(\hqq_1 \vv \hqq_1^*)\hqq_2^*
$$
This may easily be shown to yield
$$
    \vv' = (\hqq_2\hqq_1) \vv (\hqq_2\hqq_1)^*.
$$

Another application is this: suppose you have two unit vectors
$\vv_1$ and $\vv_2$, and you are interested in the rotation which
takes $\vv_1$ to $\vv_2$. That rotation will have the axis
$$
    \uu=\frac{\vv_1 \wedge \vv_2}{\|\vv_1 \wedge \vv_2 \|};
$$
 its angle $\Phi$ will be determined by $\cos 2\Phi = \vv_1\vv_2$.
 Note also that $\sin 2\Phi =\|\vv_1 \wedge \vv_2 \|$.
The quaternion $\hqq$ which realizes this rotation is then given
by
\begin{equation}
    \hqq=\begin{bmatrix}
            \sin \Phi \uu  \\
            \cos \Phi
        \end{bmatrix}
        =
        \frac{1}{\sqrt{2}}
        \begin{bmatrix}
            \sqrt{ 1- \vv_1\vv_2} \uu \\
            \sqrt{1 + \vv_1\vv_2}
        \end{bmatrix}
\end{equation}
where we had to use these trig identities:
$$
    \cos \Phi = \frac{1}{\sqrt {2}}\sqrt{1+\cos 2\Phi},
    \sin \Phi = \frac{1}{\sqrt {2}}\sqrt{1-\sin 2\Phi}.
$$
\hrule {\bf Example:} Let
$$
    \vv_1=\begin{bmatrix} 1 \\ 0 \\ 0  \end{bmatrix},
    \vv_2=\begin{bmatrix} 0 \\ 1 \\ 0  \end{bmatrix},
    \uu=\begin{bmatrix} 0 \\ 0 \\ 1  \end{bmatrix}.
$$
Thus
$$
    \hqq=
    \frac{1}{\sqrt{2}}
    \begin{bmatrix}
        0 \\ 0 \\ 1 \\ 1
    \end{bmatrix}
$$
which is the quaternion from the previous example in this section.
\hrule

%*****************************************************************
\section{Interpolation using quaternions}
%*****************************************************************

Before we discuss quaternions, let us first solve a simple 3D
problem. Given are two unit vectors $\vv_0$ and $\vv_1$, forming
an angle $\alpha$. We would like to rotate $\vv_0$ into $\vv_1$,
creating a sequence of unit vectors $\vv(t)$ which, for $t=0$,
starts out at $\vv_0$ and for $t=\alpha$ ends at $\vv_1$. We write
$\vv(t)$ as
\beq
    \vv(t)= f(t)\vv_0 + g(t)\vv_1.
    \label{eq:defv}
\eeq
The functions $f$ and $g$ do not sum to unity -- $\vv(t)$
would not be a unit vector then. We have
$$
    \| \vv_0 \wedge \vv_1 \|= \sin \alpha ~~{\rm and} ~~ \| \vv_0 \wedge
    \vv(t) \|= \sin t.
$$
Substituting (\ref{eq:defv}) into this yields
$$
    \vv(t) = \frac{\sin(\alpha - t)}{\sin(alpha)} + \frac{\sin t}{\sin
    \alpha}.
$$
For an illustration, see Fig. \ref{fig:slerp}.

\bef
    \epsfig{file=slerp.ps, height = 2 in}
    \caption{Rotation interpolation in 3D.}
   \label{fig:slerp}
\eef


 When we interpolate
between unit quaternions, the same principle applies. While
$\vv(t)$ above traced out a unit circle, quaternions, being 4D
vectors, trace out the 4D unit sphere.


Suppose $\hqq_0$ and $\hqq_1$ are two unit quaternions, such that
at time $ t=0$ we assume orientation $\hqq_0$ and at time $t=1$,
we assume $\hqq_1$. What is a meaningful orientation at an
arbitrary time $t$? The answer is given by {\em spherical linear
interpolation} or {\em slerp}:
$$
    \hqq(t)= \frac{\sin((1-t)\Phi)}{\sin \Phi} \hqq_0 + \frac{\sin(t\Phi)}{\sin
    \Phi}\hqq_1.
$$
where $\cos \Phi=\hqq_0 \cdot \hqq_1$ (dot product, not quaternion
multiplication).
\end{document}
