See also:
open '~system/packages/math/quatern.ijs'
J. Patrick Harrington wrote:
> On Fri, 18 Apr 2008, Dominic Genest wrote:
>> Hi,
>>
>> Are there phrases for quaternions ?
>> Mainly, I need to convert from quaternions to rotation matrices.
>>
>> Thanks
>>
>> Dom
>>
> I've used quaternions on occasion, and find two sets of routines
> that differ depending on the way you store an array of 4-element
> quaternions: (4 n) or (n 4). Here are both. I think they work but
> no warranty... In the following,
>
> mm=: +/ . * NB. matrix multiply
>
> NB. Quaternion maths. NB. Four-element quaternion is scalar followed by
> 3-vector
> NB. Should work on (n 4) arrays of quaternions.
>
> L=: +/"1 &. (*:"_) NB. Gives the lengths of vector arrays
> L2=: [: +/"1 *:"_ NB. squares of lengths
> dot=: [: +/"1 * NB. dot products of arrays of vectors
> pm1=: 1&|."1 NB. 1st cyclic permutation
> pm2=: 2&|."1 NB. 2nd cyclic permutation
> pm12=: 4 : '(pm1 x)*pm2 y'
> cross=: pm12 - pm12~ NB. cross products of vectors
> NB. cross=: pm12~ - pm12 NB. for left-handed coordinate system
>
> NB. Multiply quaternions: q1 QM q2 --> q3
> QM=: 4 : 0
> s1=. 0{"1 x
> v1=. }."1 x
> s2=. 0{"1 y
> v2=. }."1 y
> s=. (s1*s2) - v1 dot v2
> v=. (s1*v2)+(s2*v1)+(v1 cross v2)
> s ,&.|: v NB. catenate under transpose
> )
>
> NB. Quaternion for rotation of points by theta about unit vector u
> NB. (Rotates *points*; use -theta to rotate coordinate axes)
> NB. Usage: u Qrot theta --> q then q Pr P --> P'
> Qrot=: 4 : 0
> u=. x
> s=. cos -:y
> f=. sin -:y
> s ,&.|: (f*u) NB. catenate under transpose
> )
>
> NB. The conjuate of the quaternions: Qb q --> q_bar
> Qb=: 3 : '(1 _1 _1 _1)*"1 y'
>
> NB. The inverse of quaternions: Qi q --> q^(-1)
> Qi=: Qb % L2
>
> NB. Rotate point using quaternion q
> NB. Usage: q Pr (x y z) --> (x' y' z')
> Pr=: 4 : 0
> qb=. Qb q=. x
> }."1 q QM (0 ,&.|: y) QM qb
> )
>
> norm=: ] % L NB. normalize to unit length
>
> NB. Rotate points P (P[n,3]) about some vector u NB. by angle theta: (0
> 1 0; 0.25p1) Qrtn P --> P'
> NB. (theta in radians; u need not be a unit vector)
> Qrtn=: 4 : 0
> 'u theta'=. x
> ((norm u) Qrot theta) & Pr"1 y
> )
>
> NB. Usage: DC Trans P --> P', where DC = direction cosines of z-axis
> NB. P point in scattering frame, P' point in fixed system
> Trans=: 4 : 0
> NB. rotation axis: u=, (x y z) cross (0 0 1), normalized
> u=. norm (1 0 2){"1 (_1 1 0)*"1 x theta=. -acos 2{"1 x NB. rotate point
> back to fixed axes
> q=. u Qrot theta NB. the rotation quaternion
> q Pr y NB. rotate P in scattering frame to P' in fixed
> )
>
> NB. Form 3x3 rotation matrix which effects the same transformation as
> NB. the quaternion q. Thus, Quat_2_Mat q --> M ; M mm P --> P'
> Quat_2_Mat=: 3 : 0
> 'w x yy z'=. |: y NB. seperate components of q
> n=. #x2=. *:x NB. n is number of quaternions
> y2=. *:yy
> z2=. *:z
> a=. (1- +: y2+z2),(+:(x*yy)-(w*z)),(+:(x*z)+(w*yy))
> b=. (+:(x*yy)+(w*z)),(1- +: x2+z2),(+:(yy*z)-(w*x))
> c=. (+:(x*z)-(w*yy)),(+:(yy*z)+(w*x)),(1- +:x2+y2)
> if. n=1 do.
> (3,3)$ a,b,c
> else.
> |: (3,3,n)$ |. a,b,c
> end.
> )
>
> ================================================================
> NB. Quaternion maths. NB. Four-element quaternion is scalar followed by
> 3-vector
> NB. Should work on (4 n) arrays of quaternions.
>
> L=: +/ &. (*:"_) NB. Gives the lengths of vector arrays
> L2=: [: +/ *:"_ NB. squares of lengths
> dot=: [: +/ * NB. dot products of arrays of vectors: [4 n] dot [4 n]
> dot2=: [: +/ ]* [: ,[ NB. 1 vector dotted into array: [4 1] dot2 [4 n]
> NB. pm1=: 1&|. NB. 1st cyclic permutation
> NB. pm2=: 2&|. NB. 2nd cyclic permutation
> NB. pm12=: 4 : '(pm1 x)*pm2 y'
> NB. cross=: pm12 - pm12~ NB. cross products of vectors
> NB. cross=: pm12~ - pm12 NB. for left-handed coordinate system
> cross=: ((1: |.[)*(_1: |. ]))-((_1: |.[)*(1:|.]))
>
> NB. Multiply quaternions: q1 QM q2 --> q3
> QM=: 4 : 0
> s1=. 0{ x
> v1=. }. x
> s2=. 0{ y
> v2=. }. y
> s=. (s1*s2) - v1 dot v2
> s,(s1*"1 v2)+(s2*"1 v1)+(v1 cross v2)
> )
>
> NB. Quaternion for rotation of points by theta about unit vector u
> NB. (Rotates *points*; use -theta to rotate coordinate axes)
> NB. Usage: u Qrot theta --> q then q Pr P --> P'
> NB. u may be [3 n] array and theta [n] array.
> Qrot=: 4 : 0
> s=. cos -: y
> f=. sin -: y
> s, (f*"1 x) )
>
> NB. The conjuate of the quaternions: Qb q --> q_bar
> Qb=: 3 : '(1 _1 _1 _1)* y'
>
> NB. The inverse of quaternions: Qi q --> q^(-1)
> Qi=: Qb %"1 L2
>
> NB. Rotate point using quaternion q
> NB. Usage: q Pr P --> P' where q is [4 n] array of quaternions,
> NB. P is [3 n] array of point coordinates, and P' also [3 n]
> NB. Pr=: 4 : '}. x QM (0,y) QM Qb x'
> Pr=: [: }. [ QM (0"_ , ]) QM [: Qb [
>
> norm=: ] %"1 L NB. normalize to unit length
>
> NB. Rotate points P (P[n,3]) about some vector u NB. by angle theta: (0
> 1 0; 0.25p1) Qrtn P --> P'
> NB. (theta in radians; u need not be a unit vector)
> Qrtn=: 4 : 0
> 'u theta'=. x
> ((norm u) Qrot theta) & Pr y
> )
>
> NB. Usage: DC Trans P --> P', where DC = direction cosines of z-axis
> NB. P point in scattering frame, P' point in fixed system
> Trans=: 4 : 0
> NB. rotation axis: u=, (x y z) cross (0 0 1), normalized
> u=. norm (1 0 2){ (_1 1 0)* x theta=. -acos 2{ x NB. rotate point
> back to fixed axes
> q=. u Qrot theta NB. the rotation quaternion
> q Pr y NB. rotate P in scattering frame to P' in fixed
> )
>
> NB. Form 3x3 rotation matrix which effects the same transformation as
> NB. the quaternion q. Thus, Quat_2_Mat q --> M ; M mm P --> P'
> Quat_2_Mat=: 3 : 0
> 'w x yy z'=. |: y NB. seperate components of q
> n=. #x2=. *:x NB. n is number of quaternions
> y2=. *:yy
> z2=. *:z
> a=. (1- +: y2+z2), (+:(x*yy)-(w*z)), (+:(x*z)+(w*yy))
> b=. (+:(x*yy)+(w*z)), (1 - +: x2+z2), (+:(yy*z)-(w*x))
> c=. (+:(x*z)-(w*yy)), (+:(yy*z)+(w*x)), (1- +:x2+y2)
> if. n=1 do.
> (3,3)$ a,b,c
> else.
> |: (3,3,n)$ |. a,b,c
> end.
> )
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
