Jump to content

  • Log In with Google      Sign In   
  • Create Account

Converting Euler angles to DCM and back again


Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

  • You cannot reply to this topic
11 replies to this topic

#1 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 21 January 2010 - 11:27 AM

Hi all, First off, I apologise for the length of this query. I've double checked everything, and tried to be as clear as possible. It may even turn out to be an error with Wikipedia's article, so some other good can come out of this if we find the problem. I'm trying to combine two Euler angle sets (named say, "E1" and "E2", comprising x,y,z angles for each of course) by converting them to DCM (Direct cosine matrix), multiplying the matrices and then converting them back to Euler angles again (say, "E3"). The calcs are going screwy somewhere, but I'm not sure where. I'm only using one source for now - Wikipedia, to prevent possible convention mixup of Euler angle types: http://en.wikipedia.org/wiki/Rotation_representation_(mathematics) Before I show my working out, maybe there's a more direct way of combining Euler angles without converting to DCM etc. Unfortunately, it's not as easy as merely adding the angles by calculating:
E3.x = E1.angx + E2.angx
E3.y = E1.angy + E2.angy
E3.z = E1.angz + E2.angz
...because that will only work if we rotate on the x axis, then the y next, and then the z axis last. We may want to rotate the y axis first etc. Hence why we 'need to' (afaik) convert to DCM and back. Below is my entire working out. At the end of the post is the entire code, which is fully compilable. Sometimes I will display simple pseudo-code for speed, (and partially because it'll be a lot quicker me to do, and also shorten what I write). The variables E1 and E2 are the Euler angles (each branch off into three components: angx, angy, and angz angles) Okay so first we convert from Euler angle set 1 (E1) to DCM. (for reference see http://en.wikipedia.org/wiki/Rotation_representation_(mathematics)#Euler_angles_.E2.86.92_DCM ) sx = sin(E1.angx); sy = sin(E1.angy); sz = sin(E1.angz); cx = cos(E1.angx); cy = cos(E1.angy); cz = cos(E1.angz); The three matrices for x,y,z are as follows:
[1      0       0 ] [cy    0    -sy] [cz    sz    0]
[0      cx      sx] [0     1    0  ] [-sz   cz    0]
[0      -sx     cx] [sy    0    cy ] [0     0     1]
We now multiply these matrices, and end up with the matrix ABC (convert E1 to DCM matrix):
[ (cy*cz)                 (cy*sz)                 (-sy)   ]	
[ (sx*sy*cz + cx*-sz)     (sx*sy*sz + cx*cz)      (sx*cy) ]	
[ (cx*sy*cz+(-sx*-sz))    (cx*sy*sz + -sx*cz)     (cx*cy) ]
I'll assign this matrix to temporary variables:
a1 = (cy*cz);                   a2=(cy*sz);                     a3=(-sy);
b1 = (sx*sy*cz + cx*-sz);       b2=(sx*sy*sz + cx*cz);          b3=(sx*cy);
c1 = (cx*sy*cz+(-sx*-sz));      c2=(cx*sy*sz + -sx*cz);         c3=(cx*cy);
Okay, now we do the exactly the same with the other set of Euler angles (E2), and end up with: sx = sin(E2.angx); sy = sin(E2.angy); sz = sin(E2.angz); cx = cos(E2.angx); cy = cos(E2.angy); cz = cos(E2.angz); ....creating this matrix MNO (convert E2 to DCM matrix):
m1 = (cy*cz);                   m2=(cy*sz);                     m3=(-sy);
n1 = (sx*sy*cz + cx*-sz);       n2=(sx*sy*sz + cx*cz);          n3=(sx*cy);
o1 = (cx*sy*cz+(-sx*-sz));      o2=(cx*sy*sz + -sx*cz);         o3=(cx*cy);
Next, we multiply the ABC matrix with the MNO matrix, ending up with:
[ (a1*m1 + a2*n1 + a3*o1)   (a1*m2 + a2*n2 + a3*o2)    (a1*m3 + a2*n3 + a3*o3) ]
[ (b1*m1 + b2*n1 + b3*o1)   (b1*m2 + b2*n2 + b3*o2)    (b1*m3 + b2*n3 + b3*o3) ]
[ (c1*m1 + c2*n1 + c3*o1)   (c1*m2 + c2*n2 + c3*o2)    (c1*m3 + c2*n3 + c3*o3) ]
And finally we convert some of these entries back to a single set of Euler angles using the formula given at Wikipedia here: http://en.wikipedia.org/wiki/Rotation_representation_(mathematics)#DCM_.E2.86.92_Euler_angles ):
E3.angx =   acos ( (c1*m3 + c2*n3 + c3*o3)) ;
E3.angy = - atan2( (a1*m3 + a2*n3 + a3*o3)  ,  (b1*m3 + b2*n3 + b3*o3));
E3.angz =   atan2( (c1*m1 + c2*n1 + c3*o1)  ,  (c1*m2 + c2*n2 + c3*o2) );
The above three assigns are the only ones I'm not sure about. But even assuming I mixed up say, E3.angx with E3.angy above, the results as shown below would still show any combintation to still be garbage. Okay finished. Unfortunately, as said above, when I enter some numbers for the initial two Euler angle sets, I end up with garbage. For example, if I try some very simple ones: INPUT: E1.angx=0, E1.angy=0.05, E1.angz=0 E2.angx=0, E2.angy=0, E2.angz=0 ... we would expect E3 to be the same as E1 (since the angles are all adding zero in E2). However, we get: OUTPUT: E3.angx=0.05, E3.angy=0.5pi, E3.angz=0.5pi Let's try different ones (E2 angles will stay at zero as before): INPUT: E1.angx=0.05, E1.angy=0, E1.angz=0 OUTPUT: E3.angx=0.05, E3.angy=0, E3.angz=1pi INPUT: E1.angx=0, E1.angy=0, E1.angz=0.05 OUTPUT: E3.angx=0, E3.angy=0, E3.angz=0 Not just slightly wrong, those answers aren't even consistently wrong. What's going on? Here's the full compilable code for testing purposes:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/timeb.h>
#include <time.h>
#include <math.h>

void combineEuler(double E1angx, double E1angy, double E1angz,
		 double E2angx, double E2angy, double E2angz )
{

    double sx = sin(E1angx);    double sy = sin(E1angy);    double sz = sin(E1angz);
    double cx = cos(E1angx);    double cy = cos(E1angy);    double cz = cos(E1angz);

    // ABC matrix (convert E1 to DCM):
    double a1 = (cy*cz);                    double a2=(cy*sz);              double a3=(-sy)     ;
    double b1 = (sx*sy*cz + cx*-sz);        double b2=(sx*sy*sz + cx*cz);   double b3=(sx*cy)   ;
    double c1 = (cx*sy*cz+(-sx*-sz));       double c2=(cx*sy*sz + -sx*cz);  double c3=(cx*cy)   ;

    // Okay, now we do the exactly the same with the other set of Euler angles (E2), and end up with:
    sx = sin(E2angx);   sy = sin(E2angy);   sz = sin(E2angz);
    cx = cos(E2angx);   cy = cos(E2angy);   cz = cos(E2angz);

    //matrix MNO (convert E2 to DCM):
    double m1 = (cy*cz);                double m2=(cy*sz);                  double m3=(-sy) ;
    double n1 = (sx*sy*cz + cx*-sz);    double n2=(sx*sy*sz + cx*cz);       double n3=(sx*cy)   ;
    double o1 = (cx*sy*cz+(-sx*-sz));   double o2=(cx*sy*sz + -sx*cz);      double o3=(cx*cy)   ;

    /*Next, we multiply the ABC matrix with the MNO matrix, ending up with:
    [ (a1*m1 + a2*n1 + a3*o1)   (a1*m2 + a2*n2 + a3*o2) (a1*m3 + a2*n3 + a3*o3)  ]
    [ (b1*m1 + b2*n1 + b3*o1)   (b1*m2 + b2*n2 + b3*o2) (b1*m3 + b2*n3 + b3*o3) ]
    [ (c1*m1 + c2*n1 + c3*o1)   (c1*m2 + c2*n2 + c3*o2) (c1*m3 + c2*n3 + c3*o3)  ]*/

    double E3angx =    acos( (c1*m3 + c2*n3 + c3*o3)) ;
    double E3angy = - atan2( (a1*m3 + a2*n3 + a3*o3)  ,  (b1*m3 + b2*n3 + b3*o3));
    double E3angz =   atan2( (c1*m1 + c2*n1 + c3*o1)  ,  (c1*m2 + c2*n2 + c3*o2));

    printf("E3xang: %f ",E3angx);
    printf(" E3yang: %f ",E3angy);
    printf(" E3zang: %f\n", E3angz);
    
}

int main(int argc, char *argv[])
{
    //    Output should be same as first 3 parameters of input.
    combineEuler(0,0.05,0  ,0,0,0);        // OUTPUT (WRONG): E3.angx=0.05, E3.angy=0.5pi, E3.angz=0.5pi
    combineEuler(0.05,0,0  ,0,0,0);        // OUTPUT (WRONG): E3.angx=0.05, E3.angy=0,     E3.angz=1pi
    combineEuler(0,0,0.05  ,0,0,0);        // OUTPUT (WRONG): E3.angx=0,    E3.angy=0,     E3.angz=0

}
[Edited by - TwinbeeUK on January 21, 2010 6:37:25 PM]

Sponsor:

#2 scgames   Members   -  Reputation: 1977

Like
0Likes
Like

Posted 21 January 2010 - 01:18 PM

Converting to matrices and back seems as reasonable a way as any of 'combining' two Euler-angle rotations. May I ask why you need to do this though?

I think you touched on this in your post, but could you perhaps restate in summary what leads you to believe that the results are incorrect? Keep in mind that a) the output angles won't necessarily 'look like' any sort of simple sum of the input angles, and b) multiple Euler-angle sets can represent the same orientation, so you may not always get the results you expect (in fact, you may get different results depending on what conversion code you use).

#3 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 21 January 2010 - 10:31 PM

Quote:
Converting to matrices and back seems as reasonable a way as any of 'combining' two Euler-angle rotations. May I ask why you need to do this though?

If there's an easier, quicker, terser, less error-prone way, I'd love to hear. It's hard to believe that one has to go through all those calculations to do what I need to do - there must be a shortcut, with a few short sin/cos commands.

It's for a raytracer I'm making. The angles are initially Euler (as at least the user will want to enter data like this), and the conversion to DCM was a last resort.

Quote:
I think you touched on this in your post, but could you perhaps restate in summary what leads you to believe that the results are incorrect?

Sure. In a nutshell, I'm 'adding' angles of 0,0,0 to a positive number say 0.05,0,0. As you can see I'm testing the simplest case by simply adding zeroed angles to the original. The result should stay the same, i.e. 0.05,0,0. It's not even close as you can see above.

Quote:
Keep in mind that a) the output angles won't necessarily 'look like' any sort of simple sum of the input angles

Yes, agreed. I would imagine the numbers can often change from very simple to transcendental even.

Quote:
and b) multiple Euler-angle sets can represent the same orientation, so you may not always get the results you expect

Yes, agreed.

[Edited by - TwinbeeUK on January 22, 2010 5:31:14 AM]

#4 scgames   Members   -  Reputation: 1977

Like
0Likes
Like

Posted 21 January 2010 - 11:26 PM

It appears that you're initially building the orientation matrices using XYZ (or maybe ZYX) order, but are then extracting the angles using ZXZ order. As such, it stands to reason that the results would not be as expected. In any case, I don't think the article you linked to provides the particular conversion algorithm that you're looking for.

One thing you might do is check out geometrictools.com. The math code there includes functions for converting between matrices and Euler angles using a variety of axis orders.

#5 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 22 January 2010 - 04:09 AM

Wow you're right - thank you!! I think I saw the mention of 'ZXZ' and thought they must meant some other order of XYZ, but now I think about, it makes sense that you can still represent all angles by not even touching the Y axis.

Here is the aforementioned document that lists all conversions (all 12 of them):
http://www.geometrictools.com/Documentation/EulerAngles.pdf

I'm now having a different problem now (rotation A followed by B works, but not the other way around), but I'll work on it a lot before I write again.

In the mean time, is it possible to combine Euler set 1 and 2, without touching DCM?

#6 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 24 January 2010 - 05:00 AM

Okay, I couldn't figure it out, despite studying quite a few websites and videos all over the net. Obviously, there's issues with ambiguity (multiple Euler systems), and also it's hard to tell often if demonstrations are using local or world rotation systems (the former are where the axes are fixed to the object and rotate with it......, and the latter is where the axes are fixed in space independent of the current object's rotation).

As before, I've gone the extra mile to make it clear what I need to do. Below is a bunch of diagrams showing how my raytracer renders a colour cube at various angles. Everything now is based from SECTION 2.1 - "Factor as RxRyRz" from:
http://www.geometrictools.com/Documentation/EulerAngles.pdf

Unfortunately, diagram E in the picture below is wrong. It should look like F but doesn't. (the yellow numbers use degrees from 0 to 360 for clarity).

Diagram D is basically rotB followed by rotC and it's fine.
Diagram E is basically rotC followed by rotB and it's incorrect.




Here is the function I'm using, nice and compact, but there's presumably a mistake somewhere (I double checked from the geometrictool.com website though):


// Reference from: http://www.geometrictools.com/Documentation/EulerAngles.pdf
// SEE SECTION 2.1 - "Factor as RxRyRz"
//
euler camadd(euler E1, euler E2) {
euler E3;

double sx,sy,sz,cx,cy,cz;

sx = sin(E1.angx); sy = sin(E1.angy); sz = sin(E1.angz);
cx = cos(E1.angx); cy = cos(E1.angy); cz = cos(E1.angz);
double a1 = cy*cz; double a2 = -cy*sz; double a3 = sy;
double b1 = cz*sx*sy + cx*sz; double b2 = cx*cz-sx*sy*sz; double b3 = -cy*sx;
double c1 = -cx*cz*sy+ sx*sz; double c2 = cz*sx + cx*sy*sz; double c3 = cx*cy;


sx = sin(E2.angx); sy = sin(E2.angy); sz = sin(E2.angz);
cx = cos(E2.angx); cy = cos(E2.angy); cz = cos(E2.angz);
double m1 = cy*cz; double m2 = -cy*sz; double m3 = sy;
double n1 = cz*sx*sy + cx*sz; double n2 = cx*cz-sx*sy*sz; double n3 = -cy*sx;
double o1 = -cx*cz*sy+ sx*sz; double o2 = cz*sx + cx*sy*sz; double o3 = cx*cy;


double r00=(a1*m1 + a2*n1 + a3*o1); double r01=(a1*m2 + a2*n2 + a3*o2); double r02=(a1*m3 + a2*n3 + a3*o3);
double r10=(b1*m1 + b2*n1 + b3*o1); double r11=(b1*m2 + b2*n2 + b3*o2); double r12=(b1*m3 + b2*n3 + b3*o3);
double r20=(c1*m1 + c2*n1 + c3*o1); double r21=(c1*m2 + c2*n2 + c3*o2); double r22=(c1*m3 + c2*n3 + c3*o3);


if (r02 < +1)
{
if (r02 > -1) {
E3.angy = asin(r02);
E3.angx = atan2(-r12,r22);
E3.angz = atan2(-r01,r00);
}
else // r02 = -1
{
// Not a unique solution: thetaZ - thetaX = atan2(r10,r11)
E3.angy = -0.5*PI;
E3.angx = -atan2(r10,r11);
E3.angz = 0;
}
}
else // r02 = +1
{
// Not a unique solution: thetaZ + thetaX = atan2(r10,r11)
E3.angy = 0.5*PI;
E3.angx = atan2(r10,r11);
E3.angz = 0;
}

return E3;
}


[Edited by - TwinbeeUK on January 24, 2010 5:00:23 PM]

#7 scgames   Members   -  Reputation: 1977

Like
0Likes
Like

Posted 24 January 2010 - 10:34 AM

I'd really have to dig into that code in order to find the error (if there is one). I do have a couple of other suggestions though. (By the way, are you coding in straight C? Or is this C++?)

My first question is, do you have a math library or math code base of some sort that you're using? I would assume you do (since you've written a raytracer), but I also see that you're multiplying matrices 'by hand' in your Euler-angle functions, which is something I'd expect your math library would handle.

In any case, the the first thing I would do here is to break things down into smaller steps that can be tested and debugged separately:

1. Instead of multiplying matrices by hand at random places in your code, you should have a matrix multiplication function that does this for you. Once you've written this function and tested it to make sure it's correct, you can then more or less eliminate that as a source of error in your code.

2. Next, Euler-to-matrix and matrix-to-Euler conversions should each have their own function. Again, these functions should be tested in isolation to make sure they work.

3. With these functions in place, the problem you're trying to solve should basically reduce to the following pseudocode:
matrix m1 = matrix_from_euler(euler1);
matrix m2 = matrix_from_euler(euler2);
matrix m3 = m1 * m2; // Or maybe m2 * m1
return euler_from_matrix(m3);


#8 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 24 January 2010 - 11:54 AM

Hi, I usually code in C for speed, though will use C++ stuff sometimes.

I don't use any math libraries (apart from my own, and the basic trig from math.h), so the raytracer is entirely hand-cooked by myself. I'm tending to think that implementing full global illumination and translucent surfaces will be easier than this ;-)

You're right about breaking it up into separate functions (exactly as you've described). I would usually do (and probably intend to do) that, but this is the first time that I've started using matrices in my code. I'm 99% sure they're right though - I have checked them quite thoroughly. The code is basically just taken from SECTION 2.1 - "Factor as RxRyRz" from the document at geometrictools.com

For now, can you definitely say that the diagrams I have printed are what we would expect to see, and that we should be seeing F instead of E?

[Edited by - TwinbeeUK on January 24, 2010 6:54:10 PM]

#9 scgames   Members   -  Reputation: 1977

Like
0Likes
Like

Posted 24 January 2010 - 12:12 PM

Quote:
I don't use any math libraries (apart from my own, and the basic sin/cos from math.h), so the raytracer is entirely hand-cooked by myself.
Just to be clear, I wasn't suggesting that you use somebody else's code necessarily. I'm just saying that if you haven't already done so, you should factor out common operations such as matrix multiplication in one way or another (whether that means using an existing library, or just adding said functions to your own math module).
Quote:
For now, can you definitely say that the diagrams I have printed are what we would expect to see, and that we should be seeing F instead of E?
Yes, it does look that way - given the values in your example, F would be the expected result.

#10 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 24 January 2010 - 12:44 PM

SOLVED! It occurred to me that diagram E is a simple Z rotation of F (check in a image editor to show this).

When I use:
-0.30409,-0.1666, 0.1959

instead of:
-0.30409,-0.1666, -0.1959

...we have success!

A couple of other unlucky factors contributed to me not spotting this sooner. Thanks everyone for your help! If there's a super easy way of adjusting the code to using the local object's axes instead of the global axes I'm currently using, I'd be interested to hear. But I'm happy for now anyway :)

#11 Sneftel   Senior Moderators   -  Reputation: 1781

Like
0Likes
Like

Posted 24 January 2010 - 12:52 PM

XYZ-ordered local rotations are identical to ZYX-ordered global rotations. Similar for other axis orderings.

#12 TwinbeeUK   Members   -  Reputation: 100

Like
0Likes
Like

Posted 24 January 2010 - 01:00 PM

Wow, thanks yet again - I thought it may be something like that. Funnily enough D would swap with F in the diagrams to show this also.




Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.



PARTNERS