The quaternion
a_0 + i a_1 + j b_0 + k b_1
(underscores indicate subscripts) can be represented by the matrix
(a,b) ,: (-+b),+a NB. b is b_0 + i b_1, monadic + is conjugate
Then the determinant gives the square of the magnitude, matrix
multiplication gives the product, matrix inverse the reciprocal, matrix
addition the sum.
Kip Murray
On Sat, 19 Apr 2008, Chris Burke wrote:
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
Kip Murray
[EMAIL PROTECTED]
http://www.math.uh.edu/~km
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
