Quaternions

Extentions implemented by Peter Meulbroek, mailto://meulbroek@wag.caltech.edu.  All extentions are covered by copywrite, but are released under the GNU public license (*copyleft *).   For a copy of the license agreement, check out http://www.fsf.org/copyleft/gpl.html

We will extend the package "Algebra`Quaternions`"

Introduction

This package extends the Quaternion package of Jason Kastner and Matthew Markert to add the functionality to use quaternions as surogates for rotation.  Quaternions can represent a rotation by defining the plane of the rotation and the angle of the rotation.  The plane is defined as normal to the quaternion vector, and the angle is defined as alpha, where the pure component of the quaternion is Cos[alpha/2].  The extensions add a rotation function, and the included documentation describes how to use quaternions to replace e.g., rotation matricies.

Several websites have been very helpful with understanding the graphical implications of quaternions .   They are:
The Matrix and Quaternion FAQ  mailto://"hexapod@netcom.com"
http: // www . martinb . com/maths/
                  quaternions/quaternions . htm--mailto://martinb@martinb . com
http: // www . gamasutra . com/
                        features/programming/19980703/quaternions_  
                01.  htm  --Nick Bobick,atnb@netcom . ca
http: // www . cs . berkeley . 
                      edu/~laura/cs184/quat/
              quaternion . html--Laura Downs (laura@cs . berkeley . edu)

Though the mathematics of the package is nice, it is not easy to use in the context of graphics.  I needed to play with quaternions as a replacement for rotation matrices (bad, bad , bad), and thus,  modified this package.  In this context, given a point {a,b,c} a quaternion {x,y,z,w} represents a rotation of α within the plane defined the point {a,b,c} normal to the vector {x,y,z}, where α is defined by w = Cos[(α/2)].  Most people solving a rotation problem use principle rotation matricies (rotations around each major axes) or euler angles (rotations around z-, x-, z-axes).  The problem with these methods are twofold:

1:  the system is overdetermined.  For example, the rotation matricies have too many elements (some of which are identical).  This means computational and storage expense, as well as a more insideous problem: drift.  Due to roundoff errors in computation, the values that are supposed to be identical can drift apart.  When this occurs, the matrix changes from a simple rotation to a shearing rotation, resulting in distortion.

2:  Both systems suffer from "gimble lock", where a rotation maps one principle axis onto another, resulting in a loss of a degree of freedom.

However, quaternions solve both these problems.  Any rotation of a point to another point occurs within a plane around a single axis, so no gimble lock.  All members of the quaternion are significant, so no drift.  Finally, as a bonus, it becomes easy to inturpolate between quaternions,hence we can do smooth animations using them.

Setup

Map[Needs,   {"Algebra`Quaternions`"}];

Note that we change the format of quaternions, so they print nicely.

Quaternion[1, 1, 1, 1]

Translations

To use quaternions in graphics, we need to be able to create them from a variety of sources.

To/From Vectors

Test the new translation rules

ToQuaternion[{1, 2, 3}]

ToVector[%]

First, from/to rotation matricies.  Note that these functions are inverse to each other modulo a constant

Matrix to Quaternions

A rotation may be converted back to a quaternion through the use of  the following algorithm:

Produce a Quaternion for a known rotation.  A rotation around the z-axis has the following form:

(({{Cos[θ], Sin[θ], 0, 0}, {-Sin[θ], Cos[θ], 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}))

We will test a rotation of π/2.  Note that this would tend to cause "gimble lock".

MatrixToQuaternion[(({{Cos[θ], Sin[θ], 0, 0}, {-Sin[θ], Cos[θ], 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}})) /. θ  ->  π/2]

This is correct. Cos[(α/2)] = ((Sqrt[2])/2) => α == (π/2);  Around the axis (0,0,z).  Note that the quaternion is normalized

Quaternions to Matrix

We now try to reverse this.

QuaternionToMatrix[Quaternion[1, 0, 0, 1]]  //  MatrixForm

MatrixToQuaternion[%]

Hm.  The problem is that the original Quaternion isn't normalized.  We get back the right FORM of the quaternion, but no longer symmetric.

Normalising Quaternions

Test the normalization function.  We can also assure ourselves that the To and From functions work

Normalize[Quaternion[1, 0, 0, 1]]

QuaternionToMatrix[%] //  MatrixForm

MatrixToQuaternion[%]

Bingo!

Quaternion Inverse

Quaternion[1, 1, 0, 1]  //  Conjugate

%  **  Quaternion[1, 1, 0, 1]

The conjugate is actually the inverse, modulo a constant.  we define a  new implementation of Inverse that takes this into account

Inverse[Quaternion[1, 1, 0, 1]]

%  **  Quaternion[1, 1, 0, 1]

We get back a unit Quaternion.

UnitQuaternionQ[%]

Rotations

The point of the extentions are to use them to predict and implement rotations.  Lets experiment on a geometric figure.

Map[Needs, { "Graphics`Polyhedra`"}];

we define a set of axes to make the rotations clearer.  The axes are: red=x, green=y, blue=z.

MyAxes  =  
      Show[Graphics3D[{Thickness[0.015], {Hue[0], 
              Line[{{-2, 0, 0},  {2, 0, 0}}]},  {Hue[0.3], 
              Line[{{0,  -2, 0},  {0,  2, 0}}]}, {Hue[0.6], 
              Line[{{0, 0, -2},  {0, 0, 2}}]}}]]

Define a function to plot any points generated by the functions that follow.  To help visualize the points, we'll connect them with polygons.  To make clear what's going on, we print the list of points that are graphed.

GraphPts[p : {__List},  opts___ ]  :=  
      Module[{pt,  pr},  
Print["The points are: ", 
          p]; 
pr  =  
          Permutations[
            p]; 
Join[
MapThread[
            Graphics3D[{Hue[#2],  PointSize[0.03], Point[#1]},  
                opts] &,  {p,  (Range[Length[p]] - 1) /Length[p]}], 
          Map[Graphics3D[Polygon[Rest[#1]],  opts] &,  
            pr ]
]
]

Rotations

Helper functions

Quaternions are defined (for some, odd reason) using 1/2 angles.  Thus, to easily implement them, and to avoid inverse trig functions, we'll use the 1/2 angle formulas In the functions below, cα  is the cosine of the angle alpha.  the functions return cos[α/2] and sin[α/2], resp.

Coshalf[α_] := (Sqrt[((1 + α)/2)])

Sinhalf[α_] := (Sqrt[((1 - α)/2)])

MyNorm[v : {x_, y_, z_}] := 
      Module[{n  =  Plus @@ ((v^2))}, v/(Sqrt[n]) /;  n > 0]

CreateQuaternion

We are given two points, and want to create the quaternion that defines the least rotation of one point to get to the other. To create the quaternion, we must define the rotation plane (using a normal vector), and the angle of rotatation.  If we assume that the points are defined by unit vectors, then the axes of rotation is the cross product, and the angle of rotation is the dot product (otherwise, it gets tricky).  Note that we define rotations around the origin.

A special note:  in the case the defining points form a line, the cross product fails, since an infinite number of planes are defined by a line running through a point.  In that case, we arbitrarily choose to rotate through an easy plane to visualize, the yz plane, defined by the x-axis

CreateQuaternion[start_List, 
        point_List] := 
Module[{x1 = MyNorm[start], 
          x2 = MyNorm[point], v, vp, cα}, 
v = 
          x1\[Cross]x2; 
        cα = 
          x1 . x2; 
vp = (Sqrt[Plus]) @@ (v^2); 
 (* if the points are linear,  
          define the rotation around the x -  
            axis *) 
v = 
          If[vp > 0, (v/vp),  {1, 0, 0}]; 
Quaternion @@ 
          Join[{Coshalf[cα]}, 
            v Sinhalf[cα]]
]

pts  =  Vertices[Tetrahedron];

Note that the default tetrahedron is vertical.  This leads to all sorts of issues, since spherical coordinates have a pole (no pun intended) at ϕ = 0.

Show[GraphPts[pts],  MyAxes]

Example 1:  rotation around the z axis

Create a quaternion that rotates {1,0,0} to the vector defined by {1,1,0}.  Note that we implicitly normalize the points, so that the rotation is a simple one.

CreateQuaternion[{1, 0, 0},  {1, 1, 0}]  // Simplify

Rotate the tetrahedra using this  quaternion

Map[Rotate[%,  #] &, pts]

Get back a set of vectors to plot

Map[ToVector,  %] //  Simplify

Plot both the original and the rotated version

Show[GraphPts[pts],  
      GraphPts[%,  ColorFunction  ->  GrayLevel], MyAxes]

Example 2:  rotation around an arbitrary axis

Create a quaternion that rotates {1,0,0} to the vector defined by {1,1,1}.  Note that we implicitly normalize the points, so that the rotation is a simple one.

CreateQuaternion[{1, 0, 0},  {1, 1, 1}]  // Simplify

Rotate the tetrahedron using this  quaternion

Map[Rotate[%,  #] &, pts]

Get back a set of vectors to plot

Map[ToVector,  %] //  Simplify

Plot both the original and the rotated version

Show[GraphPts[pts],  
      GraphPts[%,  ColorFunction  ->  GrayLevel], MyAxes]