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}]]]