Re: Coordinate Transformations
- To: mathgroup at smc.vnet.net
- Subject: [mg19910] Re: Coordinate Transformations
- From: Dave Grosvenor <dag at hplb.hpl.hp.com>
- Date: Tue, 21 Sep 1999 02:22:46 -0400
- Organization: Hewlett Packard Laboratories, Bristol, UK
- References: <7rsitq$3r3@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Organization: Hewlett Packard Laboratories, Bristol, UK
Jason Rupert wrote:
>
> Please Reply to: rupertj at email.uah.edu
Hi Jason
This is my solution for the problem I think you have.
I will post some mathematica code to help you but you still need to do
some work.
Problem
--------
Calculate the rotation transformation between two right-handed
coordinate systems.
The appraoch is in two stages
i) Pick on one pair of corresponding axes, determine a transform which
will align the z-axes and the x-y planes. This is done by rotating about
the cross product of the corresponding z-axes, by the angle between the
two vectors defining the z-axes. Call this the intermediate coordinate
system where the z-axis are aligned but the x-y axes may be an arbitrary
rotation different.
ii) Now align the x-y planes of the two cordinate systems. This is done
by rotating about the z-axis (in the target and intermediate coordinate
system), by the angle between the intermediate x -axis and the target
x-axis. This is where the assumption that both coordinate systems are
right handed gets used.
QED
There is loads of scope for going wrong (signs of angles,cross products
handedness of coordinate system,etc..) but it should be obvious when you
get it wrong. So you have plenty to do!
So making a general routine is a pain.
Rotation about an axis and point
----------------------------------
Approach is to define a 3D coordinate transformation by specifying axis
of rotation (a point {0,0,0} and a vector), plus an angle of rotation.
This is standard (see Graphic Gems I "Rotation Tools", Michael E.
Pique.p466.)
see yuz
Dave Grosvenor
Mathematica code
=================
The two stages for the approach are something like below:-
i) Mat =
rotatePointVector[{0,0,0},angleBetweenVectors[{ax,ay,az},{bx,by,bz}],
Cross[{ax,ay,az},{bx,by,bz}]]
where a and b are the vector of the corresponding axes
ii) let {icx,icy,icz} = deHomogenise[Mat . homogenise[{cx,cy,cz]]
where homogenise[{x_,y_,z_}]:= {x,y,z,1}
deHomogenise[{x_,y_,z_,w_}]:= {x/w,y/w,z/w}
so now we rotate about
Mat1 =
rotatePointVector[{0,0,0},angleBetweenVectors[{icx,icy,icz},{dx,dy,dz}],
Cross[{icx,icy,icz},{dx,dy,dz}]]
{dx,dy,dz} is the corresponding axis in the target
Angle between Vector
---------------------
These expressions are derived by dot and cross product for vectors.
magnitude[v_]:= Sqrt[v.v]
sineBetweenVectors[u_,v_]:=
magnitude[Cross[u,v]]/(magnitude[u] magnitude[v])
cosineBetweenVectors[u_,v_]:= Dot[u,v]/(magnitude[u] magnitude[v])
angleBetweenVectors[u_,v_]:=With[{acos=cosineBetweenVectors[u,v],
asin=sineBetweenVectors[u,v]},
ArcTan[acos,asin]]
Rotation 4x4 matrix routine
-----------------------------
rotate routine:- from graphics gems, this routine is more general than
required (it allows a translate to an arbitrary point)--but is slow as
it does not pre-multiply the translation.
Rotate3x3[{x_,y_,z_},\[Theta]_]:=
With[{c=Cos[\[Theta]],s=Sin[\[Theta]],t=1-Cos[\[Theta]]},
{{t x^2 + c,t x y + s z, t x z - s y},
{t x y - s z, t y^2 + c,t y z + s x},
{t x z + s y, t y z - s x, t z^2 + c}}]
Now a function to uplift the rotate matrix to a 4x4 transformation
uplift[{{a00_,a01_,a02_},{a10_,a11_,a12_},{a20_,a21_,a22_}}]:=
{{a00,a01,a02,0},{a10,a11,a12,0},{a20,a21,a22,0},{0,0,0,1}}
In[29]:=
MatrixForm[uplift[Rotate3x3[{x,y,z},theta]]]
Out[29]//MatrixForm=
\!\(\*
TagBox[
RowBox[{"(", GridBox[{
{\(x\^2\ \((1 - Cos[theta])\) + Cos[theta]\),
\(x\ y\ \((1 - Cos[theta])\) + z\ Sin[theta]\),
\(x\ z\ \((1 - Cos[theta])\) - y\ Sin[theta]\), "0"},
{\(x\ y\ \((1 - Cos[theta])\) - z\ Sin[theta]\),
\(y\^2\ \((1 - Cos[theta])\) + Cos[theta]\),
\(y\ z\ \((1 - Cos[theta])\) + x\ Sin[theta]\), "0"},
{\(x\ z\ \((1 - Cos[theta])\) + y\ Sin[theta]\),
\(y\ z\ \((1 - Cos[theta])\) - x\ Sin[theta]\),
\(z\^2\ \((1 - Cos[theta])\) + Cos[theta]\), "0"},
{"0", "0", "0", "1"}
}], ")"}],
(MatrixForm[ #]&)]\)
Now we define a function to define a translate transformation as a 4x4
matrix
In[231]:=
translate[{tx_,ty_,tz_}] := {{1,0,0,tx},{0,1,0,ty},{0,0,1,tz},{0,0,0,1}}
In[232]:=
translate[{ax,ay,az}]
Out[232]=
{{1,0,0,ax},{0,1,0,ay},{0,0,1,az},{0,0,0,1}}
Conceptually we apply a transformation to move the origin to the desired
point and then perform the rotate and then we change the coordinate
system
back to undo the translate. We calculate and simplify this 4x4
transformation
In[32]:=
FullSimplify[translate[{tx,ty,tz}] . uplift[Rotate3x3[{x,y,z},theta]] .
translate[{-tx,-ty,-tz}]]
Out[32]=
\!\({{x\^2 + Cos[theta] - x\^2\ Cos[theta],
x\ y - x\ y\ Cos[theta] + z\ Sin[theta],
x\ z - x\ z\ Cos[theta] - y\ Sin[theta],
\((tx\ \((\(-1\) + x\^2)\) + x\ \((ty\ y + tz\ z)\))\)\
\((\(-1\) + Cos[theta])\) + \((tz\ y - ty\ z)\)\ Sin[theta]},
{
x\ y - x\ y\ Cos[theta] - z\ Sin[theta],
y\^2 + Cos[theta] - y\^2\ Cos[theta],
y\ z - y\ z\ Cos[theta] + x\ Sin[theta],
\((tx\ x\ y + ty\ \((\(-1\) + y\^2)\) + tz\ y\ z)\)\
\((\(-1\) + Cos[theta])\) + \((\(-tz\)\ x + tx\ z)\)\
Sin[theta]}, {
x\ z - x\ z\ Cos[theta] + y\ Sin[theta],
y\ z - y\ z\ Cos[theta] - x\ Sin[theta],
z\^2 + Cos[theta] - z\^2\ Cos[theta],
\((\((tx\ x + ty\ y)\)\ z + tz\ \((\(-1\) + z\^2)\))\)\
\((\(-1\) + Cos[theta])\) + \((ty\ x - tx\ y)\)\ Sin[theta]},
{0,
0, 0, 1}}\)
In[33]:=
MatrixForm[%%]
Out[33]//MatrixForm=
\!\(\*
TagBox[
RowBox[{"(", GridBox[{
{"1", "0", "0", "ax"},
{"0", "1", "0", "ay"},
{"0", "0", "1", "az"},
{"0", "0", "0", "1"}
}], ")"}],
(MatrixForm[ #]&)]\)
Then (lazilly) just define the desired routine as the un-simplified
matrix multiplications -- the assumption is that the vector defining the
axis of rotation is a unit vector. .
Clear[rotatePointVector]
rotatePointVector[{tx_,ty_,tz_},theta_,{x0_,y0_,z0_}]:=
With[{unit = unitVector[{x0,y0,z0}]},
With[{x=unit[[1]],y=unit[[2]],z=unit[[3]]},
translate[{tx,ty,tz}] . uplift[Rotate3x3[{x,y,z},theta]] .
translate[{-tx,-ty,-tz}]]]