Rotation interpolation: quaternions vs. Euler angles

Generic definitions

Turn off annoying warnings messages

Off[General::spell];
Off[General::spell1];

Load in necessary packages

<< Algebra`Quaternions`;
<< Graphics`Shapes`;
<< Geometry`Rotations`;
<< Graphics`Polyhedra`;

Quaternion conversion definitions

d[a_] := Sqrt[a . a];

\(\(quatSolve[r_] := 
	
        Module[{r11, r12, r13, \ r21, r22, r23, \ r31, r32, r33, \ q, \ t1, \ 
            t2, \ w}, 
		{r11, r12, r13, \ r21, r22, r23, \ r31, r32, 
              r33}\  = \ r\  // \ Flatten; 
		q\  = \ 
            Table[0. , \ {4}]; 
		t1\  = \ 
            Sqrt[\((\(\(({r11, r22, r33} . #)\) &\)\  /@ \ {{1, 1, 
                          1}, {1, \(-1\), \(-1\)}, {\(-1\), 
                          1, \(-1\)}, {\(-1\), \(-1\), 1}})\)/4 + 
                1/4]; 
		w\  = \ 
            First[Flatten[Position[t1, \ Max[t1]]]]; 
		q[\([w]\)]\  = \ 
            t1[\([w]\)]; 
		t2\  = \ 
            1/4 {r32\  - \ r23, \ r13\  - \ r31, \ r21\  - \ r12, \ 
                r12\  + \ r21, \ r13\  + \ r31, \ 
                r23\  + \ r32}; 
		Which[
			w\  == \ 
              1, 
				q[\([2]\)]\  = \ 
              t2[\([1]\)]/q[\([1]\)]; 
				q[\([3]\)]\  = \ 
              t2[\([2]\)]/q[\([1]\)]; 
				q[\([4]\)]\  = \ 
              t2[\([3]\)]/q[\([1]\)], 
			w\  == \ 
              2, 
				q[\([1]\)]\  = \ 
              t2[\([1]\)]/q[\([2]\)]; 
				q[\([3]\)]\  = \ 
              t2[\([4]\)]/q[\([2]\)]; 
				q[\([4]\)]\  = \ 
              t2[\([5]\)]/q[\([2]\)], 
			w\  == \ 
              3, 
				q[\([1]\)]\  = \ 
              t2[\([2]\)]/q[\([3]\)]; 
				q[\([2]\)]\  = \ 
              t2[\([4]\)]/q[\([3]\)]; 
				q[\([4]\)]\  = \ 
              t2[\([6]\)]/q[\([3]\)], 
			w\  == \ 
              4, 
				q[\([1]\)]\  = \ 
              t2[\([3]\)]/q[\([4]\)]; 
				q[\([2]\)]\  = \ 
              t2[\([5]\)]/q[\([4]\)]; 
				q[\([3]\)]\  = \ 
              t2[\([6]\)]/q[\([4]\)]]; 
		t1\ 
	];\)\)

\(\(quat2rot[q_] := 
	
        Module[{x, y, z, w}, 
		{w, \ x, y, z}\  = \ 
            List\  @@ \ q; 
		{1\  - \ 
                2 \((y^2\  + \ z^2)\), 
		2 \((x\ y\  - \ 
                    w\ z)\), 
		2 \((x\ z\  + \ 
                    w\ y)\), 
		2 \((x\ y\  + \ w\ z)\), 
		1\  - \ 
                2 \((x^2\  + \ z^2)\), 
		2 \((y\ z\  - \ 
                    w\ x)\), 
		2 \((x\ z\  - \ 
                    w\ y)\), 
		2 \((y\ z\  + \ w\ x)\), 
		1\  - \ 
                2 \((x^2\  + \ y^2)\)}\  // \ 
            Partition[#, \ 3] &];\)\)

\(\(eulerAngles[r_, \ 
          eps_:  10^\(-6\)] := 
		Module[{r11, r12, r13, \ r21, r22, 
            r23, \ r31, r32, r33, \ a, \ b, \ 
            g}, 
			{r11, r12, r13, \ r21, r22, r23, \ r31, r32, 
              r33}\  = \ 
            r\  // \ Flatten; 
			Which[
			Abs[\((r33\  - \ 
                    1)\)]\  < \ eps, 
				b\  = \ 0. ; 
				a\  = \ 
              0. ; 
				g\  = \ 
              ArcTan[r22, \ r12], 
			Abs[\((r33\  + \ 1)\)]\  < \ 
              eps, 
				b\  = \ Pi\  // \ N; 
				g\  = \ 
              0. ; 
				a\  = \ 
              ArcTan[r11, r21], 
		\ True, 
			\ b\  = \ 
              ArcTan[r33, \ Sqrt[r31^2\  + \ r32^2]]; 
			\ a\  = \ 
              ArcTan[\(-r32\)/Sin[b], \ r13/Sin[b]]; 
			\ g\  = \ 
              ArcTan[r23/Sin[b], \ r13/Sin[b]]]; 
		List[a, b, 
            g]];\)\)

Rotation generation definitions

colors  =  {RGBColor[1, 0, 0],  RGBColor[0, 1, 0],  
          RGBColor[0, 0, 1],  RGBColor[1, 0, 1]};

s1  =  {Polyhedron[
            Cube],  {{1, 1, 1, 1, 1}, {1, 1, 1, 1, 2}, {1, 1, 2, 1, 
              4}}};

s2  =   {Stellate[Polyhedron[Cube],  
            6.0],  {{1, 1, 1, 1, 1, 1}, {1, 1, 2, 1, 1, 1}, {1, 1, 3, 1, 1, 
              1}}};

\(\(rotationBoth[lst_, \ g_, \ parts_] := 
	
        Module[{p, \ plrng, \ pts, \ lst2, 
            estdist}, 
		p\  = \({Min[#], \ Max[#]} &\)\  /@ \ 
              Transpose[
                Partition[
                  Flatten[\((\(FullOptions[#, \ PlotRange] &\)\  /@ \ 
                        lst)\)], 6]]; 
		plrng\  = \ {{p[\([1, 1]\)], 
                p[\([2, 2]\)]}, {p[\([3, 1]\)], \ 
                p[\([4, 2]\)]}, {p[\([5, 1]\)], \ 
                p[\([6, 2]\)]}}; 
		pts\  = \ 
            Table[\(Graphics3D[{colors[\([Mod[i - 1, \ Length[colors]] + 
                            1]\)], \ PointSize[ .015], \ 
                      Point[#]}] &\)\  /@ 
								\((\(\(FullForm[#]\)[\([Sequence\  @@ \ parts[\([i]\)]]\)] &\)\  /@ \ 
                    lst)\), \ {i, \ 1, \ Length[parts]}]; 
		lst2\  = \ 
            Join[{lst}, pts]\  // \ Transpose; 
		estdist\  = \ 
            Plus\  @@ \ Flatten[
                Table[d[\((pts[\([j, i, 1, 3, 1]\)] - 
                        pts[\([j, i - 1, 1, 3, 1]\)])\)], {j, 1, 3}, {i, 2, 
                    Length[pts[\([j]\)]]}]]; 
		Print["Estimated distance traveled of selected points = ", \ estdist]; 
		{lst2, \ 
            plrng}];\)\)

\(\(rotEuler[start_, \ stop_, \ n_, \ {g_, \ parts_}] := 
	
        Module[{lst, \ i}, 
		lst\  = 
            Table[RotateShape[g, \ 
                Sequence\  @@ \ \((\((start\  + \ \((stop - start)\)*
                            i)\)\  // \ N)\)], \ {i, 0, 1, \ 
                Evaluate[1/n]}]; 
		rotationBoth[lst, \ g, \ 
            parts]];\)\)

\(\(rotQuat[start_, \ stop_, \ n_, \ {g_, \ parts_}] := 
	
        Module[{r1, r2, q1, q2, i, \ rlst, \ lst}, 
		r1\  = \ 
            RotationMatrix3D[Sequence\  @@ \ start]\ \  // \ 
              N; 
		r2\  = \ 
            RotationMatrix3D[Sequence\  @@ \ stop] // \ N; 
		q1\  = \ 
            Quaternion\  @@ \ quatSolve[r1]; 
		q2\  = \ 
            Quaternion\  @@ \ quatSolve[r2]; 
		rlst\  = \ 
            Table[quat2rot[\((q2\  ** \ Conjugate[q1])\)^\((i\  // \ 
                          N)\)\  ** \ q1], \ {i, \ 0, \ 1, 
                  Evaluate[1/n]}]\  // \ 
              Chop; 
		lst\  = \(RotateShape[g, 
                  Sequence\  @@ \ \((eulerAngles[#]\  // \ 
                        Chop)\)] &\)\  /@ \ rlst; 
		rotationBoth[lst, \ 
            g, \ parts]];\)\)

\(\(rotLerp[start_, \ stop_, \ n_, \ {g_, \ parts_}] := 
	
        Module[{r1, r2, q1, q2, i, \ rlst, \ lst}, 
		r1\  = \ 
            RotationMatrix3D[Sequence\  @@ \ start]\ \  // \ 
              N; 
		r2\  = \ 
            RotationMatrix3D[Sequence\  @@ \ stop] // \ N; 
		q1\  = \ 
            Quaternion\  @@ \ quatSolve[r1]; 
		q2\  = \ 
            Quaternion\  @@ \ quatSolve[r2]; 
		rlst\  = \ 
            Table[quat2rot[\((q2*i\  + \ q1*\((1 - i)\))\)/
                      Abs[\((q2*i\  + \ q1*\((1 - i)\))\)]\  // \ N], \ {i, \ 
                  0, \ 1, Evaluate[1/n]}]\  // \ 
              Chop; 
		lst\  = \(RotateShape[g, 
                  Sequence\  @@ \ \((eulerAngles[#]\  // \ 
                        Chop)\)] &\)\  /@ \ rlst; 
		rotationBoth[lst, \ 
            g, \ parts]];\)\)

\(\(vis[{lst_, \ 
            plrng_}] := 
		\ \(Show[#, \ PlotRange -> plrng, \ 
              Boxed -> False] &\)\  /@ \ lst;\)\)

\(\(vis2[{lst_, \ plrng_}] := 
		Show[
          Drop[lst\  // \ Transpose, \ 1]\  // \ Transpose, 
          PlotRange\  -> \ plrng, \ Boxed -> True];\)\)

Complete example

Rotation specifications (Z-X-Z Euler angles)

start  =  {Pi/200,  Pi/200,  Pi/200};
stop  =  {Pi,  Pi,  Pi};
n  =  30;

Rotation – Euler angle interpolation

tstEuler  =  rotEuler[start,  stop,  n,  s2];

Rotation – Quaternion interpolation (Lerp)

tstLerp  =  rotLerp[start,  stop,  n,  s2];

Rotation – Quaternion interpolation (Slerp)

tstSlerp  =  rotQuat[start,  stop,  n,  s2];

Visualization

Euler angle interpolation  (To animate, select second bracket from left, go to Cell menu and choose "Animate Graphics.")

vis[tstEuler];

Lerp interpolation (To animate, select second bracket from left, go to Cell menu and choose "Animate Graphics.")

vis[tstLerp];

Slerp interpolation (To animate, select second bracket from left, go to Cell menu and choose "Animate Graphics.")

vis[tstSlerp];

End point trajectories (Euler angles)

vis2[tstEuler];

End point trajectories (Lerp interpolation)

vis2[tstLerp];

End point trajectories (Slerp interpolation)

vis2[tstSlerp];