"Visualizing Quaternions"
Andrew J. Hanson
Indiana University
2004-2005

Master Utility File for  quatbook.tex  Mathematica files

Utilities for plotting

<< Graphics`ParametricPlot3D`

MakePolygons[vl_List] := 
      Module[{l = vl, l1 = Map[RotateLeft, vl], meshxx}, 
        meshxx = {l, l1, RotateLeft[l1], 
            RotateLeft[l]}; 
meshxx = 
          Map[Drop[#, -1] &, meshxx, {1}]; 
meshxx = 
          Map[Drop[#, -1] &, meshxx, {2}]; 
Polygon /@ 
          Transpose[Map[Flatten[#, 1] &, meshxx]]]

MakeWireFrame[vl_List] := 
      Module[{l = vl, l1 = Map[RotateLeft, vl], meshxx, polys}, 
        meshxx = {l, l1, RotateLeft[l1], 
            RotateLeft[l]}; 
meshxx = 
          Map[Drop[#, -1] &, meshxx, {1}]; 
meshxx = 
          Map[Drop[#, -1] &, meshxx, {2}]; 
polys  = 
          Polygon /@ 
            Transpose[
              Map[Flatten[#, 1] &, 
                meshxx]]; 
polys /. (Polygon[
                pts_] :> 
              Line[Append[pts, First[pts]]])]

wireframeParPlot3D[
        eqn_, {u_, u0_, u1_, du_}, {v_, v0_, v1_, 
          dv_}]  := 
Module[{plot, 
          graphics}, 
plot  =  
          ParametricPlot3D[eqn, {u, u0, u1, du}, {v, v0, v1, dv}, 
            DisplayFunction -> Identity]; 
graphics = 
          plot /. (Polygon[pts_] :> 
                Line[Append[pts, First[pts]]]); 
Show[
          graphics, DisplayFunction -> $DisplayFunction]]

meshviewDOT[file_, ptlist_, ptsize_]  := 
Export[
        file, 
Join[{{"SDOT"}, {ToString[
                ptsize]}, {ToString[Length[ptlist]]}}, 
          ptlist], 
"TSV"]

meshviewDOTC[file_, ptlist_, color4_, 
        ptsize_]  := 
Export[
        file, 
Join[{{"SDOTC"}, 
            color4, {ToString[ptsize]}, {ToString[Length[ptlist]]}}, 
          ptlist], 
"TSV"]

\(sphere[nth_:  16, nph_:  16, r_:  1.0]\  := \ 
      Module[{}, 
		\ \ Table[{r\ Cos[th]\ Sin[ph], r\ Sin[th]\ Sin[ph], \ 
              r\ Cos[ph]} // N, 
			\ \ \ {th, 0, 2\ Pi, 
            2\ \ Pi/nth}, \ {ph, \(-\ Pi\), \ Pi, 2\ \ Pi/nph}]]\)

wiresphere[ nth_:  16, nph_:  16, r_:  1.0]  :=  
      Graphics3D[
        MakePolygons[
            sphere[nth, nph, r]] /. (Polygon[pts_]  :> 
              Line[Append[pts, First[pts]]])]

\(drawframe[frame3dmat_, 
        xform_]\  := 
Module[{lines}, 
lines\  = \ 
          Map[Line[{{0, 0, 0}, #}] &, 
            Transpose[\((xform\  . \ frame3dmat)\)]]; 
\ 
        Show[Graphics3D[AbsoluteThickness[3]], 
          Graphics3D[{RGBColor[1, 0, 0], 
              lines[\([1]\)]}], 
Graphics3D[{RGBColor[0, 1, 
                0], lines[\([2]\)]}], 
Graphics3D[{RGBColor[
                0, 0, 1], lines[\([3]\)]}], 
Axes -> 
            True, \ AxesLabel -> {"X axis", \ "Y axis", \ "Z axis"}, 
PlotRange -> {{\(-1\), 1}, {\(-1\), 
                1}, {\(-1\), 1}}]\ 
]\)

Quaternion Library:
   q[[1]]        -->      q0 = Cos[th/2]
   q[[2,3,4]]] -->  qvec = Sin[th/2] nhat

qprod[q_List,p_List] :=
{q[[1]]*p[[1]] - q[[2]]*p[[2]] - q[[3]]*p[[3]] - q[[4]]*p[[4]],
 q[[1]]*p[[2]] + q[[2]]*p[[1]] + q[[3]]*p[[4]] - q[[4]]*p[[3]],
 q[[1]]*p[[3]] + q[[3]]*p[[1]] + q[[4]]*p[[2]] - q[[2]]*p[[4]],
 q[[1]]*p[[4]] + q[[4]]*p[[1]] + q[[2]]*p[[3]] - q[[3]]*p[[2]]} /;
   Length[q] == Length[p] == 4

qinverse[q_List] := {q[[1]], - q[[2]], -q[[3]], - q[[4]]}

normalize[
        lst_List] := Module[{norm  =  lst  .  lst, 
          eps = 10^(-14)}, 
If[
          NumberQ[norm], 
If[norm < eps,  lst,  
            N[lst/Sqrt[norm]] ], 
lst/
            Sqrt[norm]]]

proj4 [pt4D_,angle_:0] := Chop[{pt4D[[1]],pt4D[[2]], 
                   Cos[angle]*pt4D[[3]] + Sin[angle]*pt4D[[4]]}]//N

ROTATIONS and quaternions:
q*v*qbar  with v= (0,vec) a "pure" quaternion
gives EACH COLUMN of a rotation matrix:

QuatToRot[q_List] :=
Module[{q0= q[[1]], q1=q[[2]], q2 = q[[3]], q3 = q[[4]]},
   Module[{d23 = 2 q2 q3, d1 = 2 q0 q1,
          d31 = 2 q3 q1, d2 = 2 q0 q2,
          d12 = 2 q1 q2, d3 = 2 q0 q3,
          q0sq = q0^2, q1sq = q1^2, q2sq = q2^2, q3sq = q3^2},
          {{q0sq + q1sq - q2sq - q3sq, d12 - d3, d31 + d2},
           {d12 + d3, q0sq - q1sq + q2sq - q3sq, d23 - d1},
           {d31 - d2, d23 + d1, q0sq - q1sq - q2sq + q3sq}} ]]

Warning: these are the reverse of conventional argument order:
change order once all is consistent to (angle, nhat).

makeQRot[n_List:{0,0,1},angle_:0] :=
     Module[{c = Cos[angle/2], s = Sin[angle/2]},
     {c, n[[1]]*s, n[[2]]*s, n[[3]]*s}//N]

makeRotation[n_List:{0,0,1},angle_:0] :=
Module[{c = Cos[angle]//N, s = Sin[angle]//N, cm},
cm = 1 - c;
{{c + cm*n[[1]]^2, cm*n[[2]]*n[[1]] - s*n[[3]], cm*n[[3]]*n[[1]] + s*n[[2]]},
{cm*n[[1]]*n[[2]] + s*n[[3]], c + cm*n[[2]]^2,  cm*n[[3]]*n[[2]] - s*n[[1]]},
{ cm*n[[1]]*n[[3]] - s*n[[2]], 
              cm*n[[2]]*n[[3]] + s*n[[1]], c + cm*n[[3]]^2}}//N]

RotToQuat[mat_List] :=
   Module[{q0,q1,q2,q3,trace,s,t1,t2,t3},
         trace = Sum[mat[[i,i]],{i,1,3}];
         If[trace > 0,
             s = Sqrt[trace + 1.0];
               q0 = s/2;
               s = 1/(2 s);
               q1 = (mat[[3,2]] - mat[[2,3]])*s;
               q2 = (mat[[1,3]] - mat[[3,1]])*s;
               q3 = (mat[[2,1]] - mat[[1,2]])*s,
            If[(mat[[1,1]] >= mat[[2,2]]) &&
                (mat[[1,1]] >= mat[[3,3]]), (* i=0,  j = 1, k = 2 *)
                  s = Sqrt[mat[[1,1]] - mat[[2,2]] - mat[[3,3]] + 1.0];
                  q1 = s/2; s = 1/(2 s);
                  q0 = (mat[[3,2]] - mat[[2,3]])*s;
                  q2= (mat[[2,1]] + mat[[1,2]])*s;
                  q3 = (mat[[1,3]] + mat[[3,1]])*s,
               If[(mat[[1,1]] < mat[[2,2]]) &&
                (mat[[1,1]] >= mat[[3,3]]), (* i=1,  j = 2, k = 0 *)
                  s = Sqrt[mat[[2,2]] - mat[[3,3]] - mat[[1,1]] + 1.0];
                  q2 = s/2; s = 1/(2 s);
                  q0 = (mat[[1,3]] - mat[[3,1]])*s;
                  q3= (mat[[3,2]] + mat[[2,3]])*s;
                  q1 = (mat[[2,1]] + mat[[1,2]])*s,
             (* Else: (mat[[1,1]] < mat[[2,2]]) && (mat[[1,1]] < mat[[3,3]])  *)
             (* i=2,  j = 0, k = 1 *)
            s = Sqrt[mat[[3,3]] - mat[[1,1]] - mat[[2,2]] + 1.0];
                  q3 = s/2; s = 1/(2 s);
                  q0 = (mat[[2,1]] - mat[[1,2]])*s;
                  q1= (mat[[1,3]] + mat[[3,1]])*s;
                  q2 = (mat[[3,2]] + mat[[2,3]])*s]]];
                  normalize[N[{q0,q1,q2,q3}]]]

RotToQuatSym[mat_List] :=
   Module[{q0,q1,q2,q3,trace,s,t1,t2,t3},
         trace = Sum[mat[[i,i]],{i,1,3}];
             s = Sqrt[trace + 1.0];
               q0 = s/2;
               s = 1/(2 s);
               q1 = (mat[[3,2]] - mat[[2,3]])*s;
               q2 = (mat[[1,3]] - mat[[3,1]])*s;
               q3 = (mat[[2,1]] - mat[[1,2]])*s;
                {q0,q1,q2,q3}]

Distances in quaternion arcs:

QLength[curveOnSphere_List] :=
  Sum[ArcCos[curveOnSphere[[i]] . curveOnSphere[[i+1]]],
        {i,1,Length[curveOnSphere]-1}]

Arrays of Quaternion frames need these:

ForceCloseQFrames[qfrms_List] :=
  Module[{lastfrm = qfrms[[1]], thisfrm},
    Table[thisfrm = qfrms[[i]];
          lastfrm =
             If[thisfrm . lastfrm >= 0, thisfrm, - thisfrm],
             {i,1,Length[qfrms]}]]

ForceClose2DQFrames[qrows_List] :=
  Module[{lastrow = qrows[[1]], thisrow},
    Table[thisrow = ForceCloseQFrames[qrows[[i]]];
          lastrow = 
             If[thisrow[[1]] . lastrow[[1]] >= 0, thisrow, - thisrow],
              {i,1,Length[qrows]}]]

ForceClose3DQFrames[qplanes_List] := 
      Module[{lastplane = qplanes[[1]], thisplane}, 
        Table[thisplane = 
            ForceClose2DQFrames[
              qplanes[[i]]]; 
lastplane = 
            If[thisplane[[1]] . lastplane[[1]] ≥ 0, 
              thisplane, -thisplane], {i, 1, 
            Length[qplanes]}]]

4D rotation double quaternions

QQTo4DRot[p_List,q_List] :=
  Module[{q0= q[[1]], q1=q[[2]], q2 = q[[3]], q3 = q[[4]],
         p0= p[[1]], p1=p[[2]], p2 = p[[3]], p3 = p[[4]]},
{{p0*q0 + p1*q1 + p2*q2 + p3*q3, p1*q0 - p0*q1 - p3*q2 + p2*q3, 
            p2*q0 + p3*q1 - p0*q2 - p1*q3, p3*q0 - p2*q1 + p1*q2 - p0*q3}, 
 {(-(p1*q0) + p0*q1 - p3*q2 + p2*q3), (p0*q0 + p1*q1 - p2*q2 - p3*q3), 
     (-(p3*q0) + p2*q1 + p1*q2 - p0*q3), (p2*q0 + p3*q1 + p0*q2 + p1*q3)}, 
 {(-(p2*q0) + p3*q1 + p0*q2 - p1*q3), (p3*q0 + p2*q1 + p1*q2 + p0*q3),
     (p0*q0 - p1*q1 + p2*q2 - p3*q3), (-(p1*q0) - p0*q1 + p3*q2 + p2*q3)}, 
 {(-(p3*q0) - p2*q1 + p1*q2 + p0*q3), (-(p2*q0) + p3*q1 - p0*q2 + p1*q3), 
     (p1*q0 + p0*q1 + p3*q2 + p2*q3), (p0*q0 - p1*q1 - p2*q2 + p3*q3)}}]

Spline Library is here

Lerp[p0_List, p1_List, t_]  :=  (1 - t)*p0  +  t*p1

Slerp[p0_List, p1_List, t_]  :=  
      Module[{costh  =  p0  .  p1 // Chop,  th,  
          sinth}, 
  If[costh  >  0.0, 
          costh  =  Chop[costh  - 1. ]  + 1. ]; 
  If[
          costh  <  0.0, 
          costh  =  Chop[costh  + 1. ]  - 1. ]; 
 
        th  =  N[ArcCos[costh]]; 
 
        sinth  =  N[Sin[th]]; 
If[
          sinth  ==  0,  
 (1 - t)*p0  +  
            t*p1, 
 (Sin[th*(1 - t)]/sinth)*
              p0  + (Sin[th*t]/sinth)*p1 ]]

Spline and Spherical Splines

\(CRspline[p0_, p1_, p2_, p3_, t_]\  := 
	
      Lerp[Lerp[Lerp[p0, p1, t + 1], 
          Lerp[p1, p2, t], \((t + 1)\)/
            2. ], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Lerp[Lerp[p1, p2, t], 
          Lerp[p2, p3, t - 1], \ \ \ \ \ \ \ \ \ \ \ t/
            2. ], \ \ \ \ \ \ \ \ \ t]\)

\(BZspline[p0_, p1_, p2_, p3_, t_]\  := 
	
      Lerp[Lerp[Lerp[p0, p1, t], Lerp[p1, p2, t], 
          t], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Lerp[Lerp[p1, p2, t], 
          Lerp[p2, p3, t], t], t]\)

\(UBspline[p0_, p1_, p2_, p3_, t_]\  := 
	
      Lerp[Lerp[Lerp[p0, p1, \((t + 2)\)/3. ], 
          Lerp[p1, p2, \((t + 1)\)/3. ], \((t + 1)\)/
            2. ], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Lerp[
          Lerp[p1, p2, \((t + 1)\)/3. ], 
          Lerp[p2, p3, 
            t/3. ], \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t/2. ], 
        t]\)

\(SCRspline[p0_, p1_, p2_, p3_, t_]\  := 
	
      Slerp[Slerp[Slerp[p0, p1, t + 1], 
          Slerp[p1, p2, t], \((t + 1)\)/
            2. ], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Slerp[Slerp[p1, p2, t], 
          Slerp[p2, p3, t - 1], \ \ \ \ \ \ \ \ \ \ \ t/
            2. ], \ \ \ \ \ \ \ \ \ t]\)

\(SBZspline[p0_, p1_, p2_, p3_, t_]\  := 
	
      Slerp[Slerp[Slerp[p0, p1, t], Slerp[p1, p2, t], 
          t], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Slerp[Slerp[p1, p2, t], 
          Slerp[p2, p3, t], t], t]\)

\(SUBspline[p0_, p1_, p2_, p3_, t_]\  := 
	
      Slerp[Slerp[Slerp[p0, p1, \((t + 2)\)/3. ], 
          Slerp[p1, p2, \((t + 1)\)/3. ], \((t + 1)\)/
            2. ], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ Slerp[
          Slerp[p1, p2, \((t + 1)\)/3. ], 
          Slerp[p2, p3, 
            t/3. ], \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ t/2. ], 
        t]\)

squad[x0_, x1_, x2_, x3_, 
        t_]  := (1 - 2  t (1  -  t))*((1 - t) x0  +  
              t x3)  +  
        2  t*(1 - t) ((1 - t) x1  +  t x2)

Usage:

wiresphere[4, 4]

\(\(\(Module[{sspl, sp0\  = normalize[\ {\(-1\), 0, 0}], 
        sp1\  = \ 
          normalize[\ {\(- .75\), 1, \(- .25\)}], 
\ 
        sp2\  = \ normalize[{ .75, 1,  .5}], sp3 = normalize[{1, 0, 0}]}, 
      sspl\  = \ 
        ParametricPlot3D[
          Evaluate[SCRspline[sp0, sp1, sp2, sp3, t]], {t, 0, 1}, 

          PlotRange -> {{\(-1.1\), 1.1}, {\(-1.1\), 1.1}, {\(-1.1\), 
                1.1}}, 
	DisplayFunction\  -> \ Identity]; 

      Show[Graphics3D[
          RGBColor[0,  .7, 
            1]], 
wiresphere[], 
Graphics3D[RGBColor[0, 0, 0]], 
        Graphics3D[Thickness[ .015]], sspl, 
	
        Graphics3D[{RGBColor[0, 0, 0], PointSize[0.04], 
            Map[Point[#1] &, {sp0, sp1, sp2, sp3}]}], 

        DisplayFunction\  -> $DisplayFunction\ ]\ ]\)\(\ \)\)\)

TUBING Library is here: patterned after Gray's code

\(\(\(\(kappa[alpha_]\)[t_] := 
      Simplify[Factor[
                Cross[D[alpha[tt], tt], D[alpha[tt], {tt, 2}]] . 
                  Cross[D[alpha[tt], tt], D[alpha[tt], {tt, 2}]]]]^\((1/2)\)/
          Simplify[Factor[D[alpha[tt], tt] . D[alpha[tt], tt]]]^\((3/2)\) /. 
        tt -> t\)\(
\)
    \)\)

\(\(\(\(tau[alpha_]\)[t_] := 
      Simplify[Det[{D[alpha[tt], tt], D[alpha[tt], {tt, 2}], 
                D[alpha[tt], {tt, 3}]}]]/
          Simplify[
            Factor[Cross[D[alpha[tt], tt], D[alpha[tt], {tt, 2}]] . 
                Cross[D[alpha[tt], tt], D[alpha[tt], {tt, 2}]]]] /. 
        tt -> t\)\(
\)\(
\)
    \)\)

tangent[alpha_][t_] := 
      Module[{deriv = D[alpha[tt], tt]}, 
        deriv/Sqrt[deriv . deriv] /. tt -> t]

binormal[alpha_][t_] := 
Module[{binorm= Cross[
  D[alpha[tt],tt],D[alpha[tt],{tt,2}]]},
  binorm/Sqrt[binorm.binorm]/.tt->t]

normal[alpha_][t_] := Cross[binormal[alpha][t],
               tangent[alpha][t]]

tubecurve[gamma_][r_][t_,theta_] := Simplify[gamma[t] +
    r*(Cos[theta]*normal[gamma][t] +
       Sin[theta]*binormal[gamma][t])]

Core functions for tubed curve object definitions are here:

torusknot[a_, b_, c_][p_, q_][
        t_] := {(a + b*Cos[q*t])*Cos[p*t], (a + b*Cos[q*t])*Sin[p*t], 
        c*Sin[q*t]}

\(\(\(Trefoil[t_] := \(\(torusknot[8, 3, 5]\)[2, 3]\)[t]\)\(
\)
    \)\)

\(\(\(\(helix[a_, b_]\)[t_] := {a*Cos[t], a*Sin[t], b*t}\)\(
\)
    \)\)

Usage:

ParametricPlot3D[helix[1.0,  .1*Pi][t] // Evaluate, {t, 0, 4  Pi}, 
      AspectRatio -> 1]

ParametricPlot3D[tubecurve[helix[2, 0.5]][ .5][t, theta] // 
        Evaluate, {theta, 0, 2*Pi}, {t, 0, 4*Pi}, 
      PlotPoints -> {10, 40}, Axes -> None]