Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
1998
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

RE: Re: Rotate3D bug solution



I followed with great interest all the discussions about Eulerian
angles, I  have not got much to say about them, (the name is the same
in French  litterature) but I remember having spent some time on the
following code,  which uses a generalization of Eulerian angles in
unspecified dimension.  The code uses Eulerian angles, and mostly shows
the way in which these  angles are defined:

(*
this module generates a random point on a sphere of center "center" and 
radius "radius".
The algorithm uses generalized polar angles as random numbers and 
transforms the polar
coordinates in cartesian co ordinates. *)

randomPointOnASphere[center_,radius_]:=Module[
	{dim,rp,rc,x},
	dim=Length[center];
	x[i_Integer,pc_]:=pc[[1]] Switch[i,
		1,Product[Sin[pc[[j]] ],{j,2,dim}],
		_,Product[Sin[pc[[k]] ],{k,i+1,dim}] Cos[pc[[i]]]
		];
	rp=Flatten[
		{radius,Random[Real,{0,2
Pi//N}],Table[Random[Real,{0,Pi//N}],{dim-2}]}];
	rc=Table[x[i,rp],{i,dim}];
	Return[rc+center]
	]
(*  end *)


Eulerian angles are here called polar angles. I'm not a mathematician,
so  the terminology might be wrong. Anyway this code returns random
points on a  sphere, random in the sense of uniform density on the
surface of the  sphere.

Hope this helps,

-----------------------------------------------

Jean-Marie THOMAS
Conseil et Audit en Ingenierie de Calcul jmthomas@cybercable.tm.fr
+33 (0)3 88 32 93 64
www.cybercable.tm.fr/~jmthomas
=======================


-----Message d'origine-----
De:	John Sidles [SMTP:sidles@u.washington.edu] Date:	lundi 12 janvier
1998 10:10
A:	mathgroup@smc.vnet.net
Objet:	[mg10382] Re: Rotate3D bug solution

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: Why doesn't this work ?
  • Next by Date: Re: Numeric overflow
  • Prev by thread: Re: Rotate3D bug solution
  • Next by thread: Eigenvectors