Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I've wrote a function that will extract euler angle from matrix according to Matlab's function in AeroSpace toolbox's dcm2angle package file. But not quite right.

 threeaxisrot[r11_, r12_, r21_, r31_, r32_, r11a_, r12a_] :=  Module[{r1, r2, r3},  (*        % find angles for rotations about X, Y, and Z axes*)          r1 = ArcTan[ r11/ r12];          r2 = ArcSin[ r21 ];          r3 = ArcTan[ r31/ r32];   If[Abs[r21] >= 1, Print["HA-ha"]; r1 = ArcTan[r11a/r12a];    r2 = ArcSin[r21]; r3 = 0];  {r1, r2, r3}  ]       

Matrix2EulerAngleZXY[matrix0_] :=  Module[{dcm = matrix0},   threeaxisrot[-dcm[[2, 1]], dcm[[2, 2]], dcm[[2, 3]], -dcm[[1, 3]],    dcm[[3, 3]], dcm[[1, 2]], dcm[[1, 1]]]]

Matlab package directory in my computer is

"F:\\Program Files\\MATLAB\\toolbox\\aero\\aero\\dcm2angle.m"

matrixList =  Table[RotationMatrix[i, {1, 0, 0}], {i, 0, 2 \[Pi], \[Pi]/6}]//N;

Matrix2EulerAngleTrackExpression // Definition

(*
    Matrix2EulerAngleTrackExpression[matrix0_,type_:ZXY]:=Module[{matrix=matrix0[[1;;3,1;;3]]}
,ttt=ToExpression[StringCases[StringReplace[ppp=MATLink`MEvaluate[[yaw, pitch, roll] = dcm
2angle(<>matrix2Matlab[N[Chop[matrix,1/10^5]]]<>,'<>type<>')],x:(NumberString~~e-~~NumberS
tring):>0],NumberString]]]
*)

Result1


This is the Matlab's result(good one).

Matrix2EulerAngleTrackExpression /@ matrixList

(*
    {{0,0,0},{0,-0.5236,0},{0,-1.0472,0},{0,-1.5708,0},{-3.1416,-1.0472,-3.1416},{-3.1416,-0.5
236,-3.1416},{-3.1416,0,-3.1416},{-3.1416,0.5236,-3.1416},{-3.1416,1.0472,-3.1416},{0,1.57
08,0},{0,1.0472,0},{0,0.5236,0},{0,0,0}}
*)

Matrix2EulerAngleZXY /@ matrixList

(*
    {{0.,0.,0.},{0.,-0.523599,0.},{0.,-1.0472,0.},{0.,-1.5708,0},{0.,-1.0472,0.},{0.,-0.523599
,0.},{0.,0.,0.},{0.,0.523599,0.},{0.,1.0472,0.},{0.,1.5708,0},{0.,1.0472,0.},{0.,0.523599,
0.},{0.,0.,0.}}
*)

Where goes wrong?

Do you have a robust euler angle function/package?

share|improve this question
    
I use the MATLink package to calculate the angles before. –  HyperGroups Oct 30 '13 at 15:30
    
Apart from the other answers in this: mathematica.stackexchange.com/q/29924/131, David Park´s Presentations package seems to deal with Euler angles: mathematica.stackexchange.com/a/30163/131 –  Yves Klett Oct 30 '13 at 15:52

2 Answers 2

More in-depth explanation of the math can be found at: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToEuler/

float yaw,pitch,roll=0.0f,epsilon=0.0001f;
if (m[1][0] > (1.0f-epsilon)) { // singularity at north pole
    yaw = atan2(m[0][2],m[2][2]);
    pitch = M_PI / 2.0f;
} else if (m[1][0] < -(1.0f-epsilon)) { // singularity at south pole
    yaw = atan2(m[0][2],m[2][2]);
    pitch = M_PI / -2.0f;
} else {
    yaw = atan2(-m[2][0],m[0][0]);
    pitch = asin(m[1][0]);
    roll = atan2(-m[1][2],m[1][1]);
}
return Vector3(yaw,pitch,roll);

You can of course switch handedness by transposing all the m[x][y] accesses.

share|improve this answer
    
Not exactly a Mathematica implementation... –  Yves Klett Oct 30 '13 at 17:36
    
He cross posted this to GameDev asking for a C/C++ implementation specifically. I had answered there but the moderators asked me to move it here. –  MickLH Oct 30 '13 at 17:42
1  
Talk about convolution :-) –  Yves Klett Oct 30 '13 at 18:02
    
Yes, I asked here because I'd like a mathematica version. Can you make it the same result as my sample? –  HyperGroups Oct 31 '13 at 1:30
3  
You have to do that much on your own, come on. I've even tested it against the transpose of the matlab spec. –  MickLH Oct 31 '13 at 2:11

Let's call the two answers that you call the good one and the bad one q1 and q2. The difference is:

q1-q2
{{0., 0., 0.}, {0., 1.*10^-6, 0.}, {0., 0., 0.}, {0., 0., 0}, {3.1416, 0., 3.1416}, {3.1416, 1.*10^-6, 3.1416}, {3.1416, 0.,   3.1416}, {3.1416, -1.*10^-6, 3.1416}, {3.1416, 0., 3.1416}, {0., 0., 0}, {0., 0., 0.}, {0., -1.*10^-6, 0.}, {0., 0., 0.}}

Looking through this, all the terms are either roundoff error (10^-6) or are Pi. So all you need do to make them match is to replace the near-Pi values with 0 (or vice versa).

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.