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.
)
=================================================================
-Patrick
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
