Services & ResourcesWolfram Forums
 MathGroup Archive
 1998 January February March April May June July August September October November December

Re: Rotate3D bug solution

• To: mathgroup@smc.vnet.net
• Subject: [mg10382] Re: Rotate3D bug solution
• From: sidles@u.washington.edu (John Sidles)
• Date: Mon, 12 Jan 1998 04:10:12 -0500
• Organization: University of Washington, Seattle
• References: <199712250838.DAA15159@smc.vnet.net.> <34A499DD.5E16@physics.uwa.edu.au> <684v6o$q8r@smc.vnet.net> <68csvb$aa4@smc.vnet.net>

In article <68csvb$aa4@smc.vnet.net>, Selwyn Hollis <shollis@peachnet.campus.mci.net> wrote: >Mark Evans wrote: > >> Paul Abbott wrote: >> > >> > The use of Eulerian angles for specifying rotations in 3D IS >> standard >> > (at least in maths and physics -- especially in quantum mechanics, >> > crystallography, and angular momemntum theory). Note that the >> eulerian >> > angle parametrization avoids the singularities that arise in other >> > parametrizations. >> > >> >> In the same sense, you could say that sea shells are standard legal >> tender if you live in a certain part of the world. >> >> Paul is right that there is nothing technically wrong with this kind >> of >> rotation. My point was that Mathematica packages should be written >> for >> a wider audience. It seems intuitive that the most common >> understanding of a rotation matrix is one that rotates sequentially >> about each of the three coordinate axes. The fact that Mathematica >> does not offer this rotation by default is a slip-up in my mind. > >Right on, Mark! > >I've never heard of Eulerian angles before encountering them in >Mathematica. Maybe they're the usual tricks-of-the-trade to a few >quantum physicists, but you'll be hard-pressed to find a reference to >them in any but the most esoteric mathematics literature. > Well, the latitude and longitude coordinates on the globe are (essentially) Eulerian coordinates -- so they're not *too* esoteric. Here's a very important and useful, yet simple, theorem which anyone working with rotations needs to know -- it explains why seemingly inequivalent conventions are actually precisely equivalent. Let$R(v)$be a function which computes the three-by-three matrix associated with a three-vector$v$, where the direction of$v$gives the axis of rotation, and the magnitude of$v$gives the angle of rotation. (Exercise: program this useful utility function in Mathematica, together with its inverse function v(R) -- answer given at end!). Let$U$be an arbitrary rotation matrix, with$U^{t}$the matrix transpose of$U$. Since$U$is a rotation,$U^{t}$is the matrix inverse of$U$, i.e.,$U^{t} U = I$. Then for any$U$and$v$, here's the key theorem! U R(v) U^{t} = R(U v) To see how this clarifies the literature on rotations, suppose Textbook A defines Euler matrices in terms of rotation angles${\theta, \phi,
\psi}$about fixed unit axes${\hat{n}_1,\hat{n}_2,\hat{n}_3}\$ as
follows

R(\theta, \phi, \psi)
\edef R(\theta \hat{n}_1) R(\phi \hat{n}_2) R(\psi \hat{n}_3)
\edef R_1 R_2 R_3

Now we use the above theorem to rewrite this in terms of *moving* axes.

R(\theta,\phi,\psi) = R(\psi R1 R2 \hat{n_3})
R(\phi R1 \hat{n_2}) R(\theta \hat{n_1})

Cool! Its the same three angles, but now applied in the *opposite*
order, and about moving instead of fixed axes.  Yet the final  matrix
is the same.  And it is perfectly reasonable for Textbook B to adopt
this moving-axis convention to define Euler angles.

Given all this ambiguity, I no longer use Euler angles when doing
calculations involving rotations.  It is just too easy to  confuse the
various signs and conventions!

A much safer strategy, which I recommend, is to simply code the
function R(v) and its inverse v(R) as utility routines in whatever
language you prefer.  You will find that these two functions,  plus
ordinary matrix multiplication, suffice for *any* calculation
involving rotations.  No more Euler angle torture!  No more sine and
cosine functions with obscure arguments!

Here's some Mathematica code for R(v):

rotationMatrix[x_] := Block[
{angle,xHat,pPar,pPerp,pOrthog},

angle = Sqrt[x.x]//N;
If[angle<10^-6,Return[id]];

xHat = x/angle;

pPar = Outer[Times,xHat,xHat];

pPerp =DiagonalMatrix[{1.0,1.0,1.0}] - pPar;

pOrthog = {
{0.0,xHat[[3]],-xHat[[2]]},
{-xHat[[3]],0.0,xHat[[1]]},
{xHat[[2]],-xHat[[1]],0.0}
};

pPar + Cos[angle]* pPerp + Sin[angle]*pOrthog
]

It is left as an exercise to (a) figure out how the above works, and (b)
code the inverse function!  Coding the inverse function  is quite a
nontrivial exercise.  Hint: (a) determine the cosine() from the trace
of R, then determine the sine() and the the axis of rotation from the
antisymmetric part of R.

This works for all rotations *except* for angles near Pi, for which the
axis of rotation should be set to the (unique) eigenvector of the
symmetric part of R that has unit  eigenvalue (this is because sin(Pi)
= 0).

A final refinement: the above algorithm assumes that R is orthogonal.
But what if R is just *close* to orthogonal, but has some accumulated
numerical imprecision?  Even worse, you can bet that some future user
of your v(R) routine will hand it an R matrix that is grossly
non-orthonormal!  The right thing to do, therefore, is  to condition R
before calculating v, by (a) calculating the singular  value
decomposition for R, then (b) adjusting all the singular values  to
have unit magnitude.  This will yield a cleaned-up exactly orthogonal
R which (formally) is the orthogonal matrix that is closest to the
input matrix in the least-mean squares sense.  Better issue an error
message if the singular values are far from unity -- because this
indicates abuse of the inverse routine!

The Mathematica routines SingularValue[] and Eigensystem[] are well
suited to the above tasks -- it's tougher in "C".

I wrote this up at length because I have a student who is generating
animations in Mathematica and POV-Ray -- he might as well learn to do
things the easy way!

Happy rotating ... JAS



• Prev by Date: Re: Math Problem in Mathematica
• Next by Date: Re: writing out a general n-th order iterative eqn.
• Prev by thread: Re: Rotate3D bug solution
• Next by thread: RE: Re: Rotate3D bug solution