17

I'm looking for a way to convert direction vector (X,Y,Z) into Euler angles (heading, pitch, bank). I know that direction vector by itself is not enough to get the bank angle, so there's also another so-called Up vector.

Having direction vector (X,Y,Z) and up vector (X,Y,Z) how do I convert that into Euler angles?

Kromster
  • 6,665
  • 7
  • 55
  • 98
  • 7
    This question appears to be off-topic because it is about mathematics it should be moved to http://math.stackexchange.com/. – Bas Bossink Feb 07 '14 at 09:07
  • 2
    @BasBossink I agree, but unfortunately the question is too old to be moved now. – jrh Jun 28 '18 at 19:47

2 Answers2

26

Let's see if I understand correctly. This is about the orientation of a rigid body in three dimensional space, like an air plane during flight. The nose of that airplane points towards the direction vector

D=(XD,YD,ZD) .

Towards the roof is the up vector

U=(XU,YU,ZU) .

Then heading H would be the direction vector D projected onto the earth surface: 

H=(XD,YD,0) ,

with an associated angle 

angle_H=atan2(YD,XD) .

Pitch P would be the up/down angle of the nose with respect to the horizon, if the direction vector D is normalized you get it from

ZD=sin(angle_P)

resulting in 

angle_P=asin(ZD) .

Finally, for the bank angle we consider the direction of the wings, assuming the wings are perpendicular to the body. If the plane flies straight towards D, the wings point perpendicular to D and parallel to the earth surface:

W0 = ( -YD, XD, 0 )

This would be a bank angle of 0. The expected Up Vector would be perpendicular to W0 and perpendicular to D

U0 = W0 × D

with × denoting the cross product. U equals U0 if the bank angle is zero, otherwise the angle between U and U0 is the bank angle angle_B, which can be calculated from

cos(angle_B) = Dot(U0,U) / abs(U0) / abs(U)
sin(angle_B) = Dot(W0,U) / abs(W0) / abs(U) .

Here 'abs' calculates the length of the vector. From that you get the bank angle as

angle_B = atan2( Dot(W0,U) / abs(W0), Dot(U0,U) / abs(U0) ) .

The normalization factors cancel each other if U and D are normalized.

pentadecagon
  • 4,317
  • 1
  • 16
  • 21
  • 1
    I know this is an old post but it's useful for me too. I'd like to know what is the meaning of `abs(W0) * abs(U0)`. Does the `abs()` function return the lenght of the vector? And then, shouldn't the `atan2` function take two arguments? x and y? In the example is passed olny one paramater and I can't understand why. – zeb Mar 18 '15 at 23:05
  • 5
    atan2 assumes two arguments, I suppose that `angle_B = atan2( Dot(W0,U) / Dot(U0,U) / abs(W0) * abs(U0) ) .` should be rather `angle_B = atan2( Dot(W0,U) , Dot(U0,U) / abs(W0) * abs(U0) ) .` – Samuel Jun 12 '15 at 13:34
  • I voted this answer down since the meaning of abs() is not clear and as Samuel said atan2 expects two parameters, not just one. – padmalcom Oct 21 '17 at 21:16
  • Great algorithm. It solved my problem too :) That problem is part of a thesis however. Therefore do you happen to have a "source" where this is explained in a similar way, so I can reference it? – DragonGamer Feb 04 '19 at 11:08
  • I think `angle_B` is wrong for the edgecase where `D=(0,0,1)` or `D=(0,0,-1)`. `W0` and `U0` both end up as `(0,0,0)`. I'm not currently sure how to address that. – Robadob Mar 13 '21 at 15:27
  • Implementing this in GLSL I've also found that `angle_H` comes out wrong for the same edgecase. However, that's due to GLSL's `atan(0,0)` having undefined behaviour, so findings may vary. – Robadob Mar 13 '21 at 18:16
5

we need three vectors: X1, Y1, Z1 of local coordinate system (LCS) expressed in terms of world coordinate system (WCS). The code below presents how to calculate three Euler angles based on these 3 vectors.

#include <math.h> 
#include <float.h> 

#define PI 3.141592653589793 
/**
 * @param X1x
 * @param X1y
 * @param X1z X1 vector coordinates
 * @param Y1x
 * @param Y1y
 * @param Y1z Y1 vector coordinates
 * @param Z1x
 * @param Z1y
 * @param Z1z Z1 vector coordinates
 * @param pre precession rotation
 * @param nut nutation rotation
 * @param rot intrinsic rotation
 */
void lcs2Euler(
        double X1x, double X1y, double X1z,
        double Y1x, double Y1y, double Y1z,
        double Z1x, double Z1y, double Z1z,
        double *pre, double *nut, double *rot) {
    double Z1xy = sqrt(Z1x * Z1x + Z1y * Z1y);
    if (Z1xy > DBL_EPSILON) {
        *pre = atan2(Y1x * Z1y - Y1y*Z1x, X1x * Z1y - X1y * Z1x);
        *nut = atan2(Z1xy, Z1z);
        *rot = -atan2(-Z1x, Z1y);
    }
    else {
        *pre = 0.;
        *nut = (Z1z > 0.) ? 0. : PI;
        *rot = -atan2(X1y, X1x);
    }
}
4pie0
  • 27,469
  • 7
  • 70
  • 110
  • 1
    When you say 3 vectors (named X Y Z for some reason?) you mean the origin (0,0,0) the direction and the up vector? Which is which in your answer? – Kromster Feb 07 '14 at 09:12
  • X,Y,Z are the LCS coordinate system axises direction vectors in [GCS]. Z is your direction and the two can be obtain by vector multiplication of Z and your up vector – Spektre Feb 07 '14 at 09:15
  • I don't understand that. Can we please switch to simple GCS for clarity sake? – Kromster Feb 07 '14 at 09:20
  • 1
    LCS is coordinate system of your 'camera' or plane or whatever viewing/moving in its Z direction. (LCS Z axis) = (Z direction vector) = (your direction vector). Orthogonal systems has 3 axises not one so you need 2 perpendicular vectors to it. If your up vector is perpendicular to direction then (Y-axis vector) = (your up vector) and then the X axis is easy X = Z * Y or X = Y * Z ... – Spektre Feb 07 '14 at 09:26
  • this code snippet is taken from here: http://geom3d.com/data/documents/Calculation=20of=20Euler=20angles.pdf – christian Jan 19 '20 at 15:36