2

I have the following code:

#define GLM_ENABLE_EXPERIMENTAL

#include <iostream>
#include <glm/glm.hpp>
#include <glm/ext.hpp>

// combines 2 XYZ euler angles given in degrees
glm::vec3   EulerCombine(const glm::vec3& first, const glm::vec3& second)
{
        glm::vec3 output;
        glm::mat4 m1 = glm::eulerAngleXYZ(glm::radians(first.x), glm::radians(first.y), glm::radians(first.z));
        glm::mat4 m2 = glm::eulerAngleXYZ(glm::radians(second.x), glm::radians(second.y), glm::radians(second.z));
        glm::extractEulerAngleXYZ(m2 * m1, output.x, output.y,  output.z);
        return glm::degrees(output);
}

// applies the XYZ euler rotation on p
glm::vec3   EulerRotate(const glm::vec3& p, const glm::vec3& euler)
{
        glm::vec3 output;
        output = glm::rotateX(p, glm::radians(euler.x));
        output = glm::rotateY(output, glm::radians(euler.y));
        output = glm::rotateZ(output, glm::radians(euler.z));
        return output;
}

int     main(void)
{
        glm::vec3 euler1(30, 20, 90); // euler angles in degrees
        glm::vec3 euler2(20, 30, 10);
        glm::vec3 euler3 = EulerCombine(euler1, euler2);
        glm::vec3 p(-10, 7, 23);

        glm::vec3 result1 = EulerRotate(EulerRotate(p, euler1), euler2);
        glm::vec3 result2 = EulerRotate(p, euler3);

        std::cout << result1.x << " " << result1.y << " " << result1.z << std::endl;
        std::cout << result2.x << " " << result2.y << " " << result2.z << std::endl;
}

but it prints:

17.9056 -6.99702 17.5622
17.2369 4.15094 19.0699

which is an issue because both result1 and result2 (the vectors being printed) should be the same.

I feel like my EulerCombine function has to be correct because its almost exclusively glm, but then that would mean EulerRotate was incorrect - but I thought an XYZ Euler angle meant you rotate it on X axis followed by Y axis followed by Z axis.

Whats going on?

Nicol Bolas
  • 378,677
  • 53
  • 635
  • 829
Theo Walton
  • 973
  • 6
  • 17
  • Correct order of rotations - everytimes easy to stumble... IMHO, you have to apply in `EulerRotate()` `glm::rotateX(glm::rotateY(glm::rotateZ(p, glm::radians(euler.z)), glm::radians(euler.y)), glm::radians(euler.x))`. If I only could find some code in our S/W where I did it right... (...I would write an answer.) – Scheff's Cat Jul 22 '18 at 09:06
  • _I thought an XYZ Euler angle meant you rotate it on X axis followed by Y axis followed by Z axis._ Turn it around: rotation of p about Z axis, rotation of p' about Y axis, rotation of p'' about X axis. Rotations are not commutative (in opposition to translations) - hence the order is so important (and makes rotations somehow complicated). – Scheff's Cat Jul 22 '18 at 09:27

1 Answers1

4

As OP already suspected: the implementation of EulerRotate() is wrong. It should be instead:

glm::vec3 EulerRotate(const glm::vec3& p, const glm::vec3& euler)
{
  return glm::rotateX(glm::rotateY(glm::rotateZ(p, glm::radians(euler.z)),
    glm::radians(euler.y)), glm::radians(euler.x));
}

which can also be written (more readable) as:

glm::vec3 EulerRotate(const glm::vec3& p, const glm::vec3& euler)
{
  const glm::vec3 p1 = glm::rotateZ(p, glm::radians(euler.z));
  const glm::vec3 p2 = glm::rotateY(p1, glm::radians(euler.y));
  const glm::vec3 p3 = glm::rotateX(p2, glm::radians(euler.x));
  return p3;
}

The transformation matrix corresponding to a Euler angles can be split like this:

MrXYZ = MrX · MrY · MrZ

and hence:

P' = MrXYZ · P = MrX · MrY · MrZ · P

This is

P' = P rotated about Z, rotated about Y, rotated about X.

As simple as it sounds – to do this correctly in practical work drives me crazy from time to time. Rotations in 3D space are not commutative → the order is very important.

I prepared a sample to demonstrate this. As I had no glm at hand (and was not willing to install it on my side), I used own code (I once prepared for another sample) and tried to resemble OPs code as close as possible:

#include <iostream>
#include "linMath.h"

double radians(double);
Vec3 degrees(Vec3 angles);
Mat4x4 eulerAngleXYZ(double rX, double rY, double rZ);
void extractEulerAngles(const Mat4x4 &mat, double &rX, double &rY, double &rZ);
Vec3 rotateX(const Vec3 &p, double angle);
Vec3 rotateY(const Vec3 &p, double angle);
Vec3 rotateZ(const Vec3 &p, double angle);

// combines 2 XYZ euler angles given in degrees
Vec3 eulerCombine(const Vec3 &first, const Vec3 &second)
{
  const Mat4x4 mat1 = eulerAngleXYZ(radians(first.x), radians(first.y), radians(first.z));
  const Mat4x4 mat2 = eulerAngleXYZ(radians(second.x), radians(second.y), radians(second.z));
  Vec3 output;
  extractEulerAngles(mat2 * mat1, output.x, output.y, output.z);
  return degrees(output);
}

// applies the XYZ euler rotation on p
Vec3 eulerRotate(const Vec3 &p, const Vec3 &euler)
{
#ifndef FIX // Theo:
  Vec3 output;
  output = rotateX(p, radians(euler.x));
  output = rotateY(output, radians(euler.y));
  output = rotateZ(output, radians(euler.z));
  return output;
#else // Dirk:
  return rotateX(rotateY(rotateZ(p, radians(euler.z)), radians(euler.y)), radians(euler.x));
#endif // FIX
}

int main()
{
  Vec3 euler1(30, 20, 90); // euler angles in degrees
  Vec3 euler2(20, 30, 10);
  Vec3 euler3 = eulerCombine(euler1, euler2);
  Vec3 p(-10, 7, 23);
  
  Vec3 result1 = eulerRotate(eulerRotate(p, euler1), euler2);
  Vec3 result2 = eulerRotate(p, euler3);

  std::cout << result1.x << " " << result1.y << " " << result1.z << std::endl;
  std::cout << result2.x << " " << result2.y << " " << result2.z << std::endl;
  // done
  return 0;
}

double radians(double angle) { return degToRad(angle); }

Vec3 degrees(Vec3 angles)
{
  return Vec3(radToDeg(angles.x), radToDeg(angles.y), radToDeg(angles.z));
}

Mat4x4 eulerAngleXYZ(double rX, double rY, double rZ)
{
  return Mat4x4(InitRotX, rX) * Mat4x4(InitRotY, rY) * Mat4x4(InitRotZ, rZ);
}

void extractEulerAngles(const Mat4x4 &mat, double &rX, double &rY, double &rZ)
{
  decompose(mat, RotX, RotY, RotZ, rX, rY, rZ);
}

Vec3 rotateX(const Vec3 &p, double angle)
{
  const Vec4 p_ = Mat4x4(InitRotX, angle) * Vec4(p, 1.0);
  return Vec3(p_.x, p_.y, p_.z);
}

Vec3 rotateY(const Vec3 &p, double angle)
{
  const Vec4 p_ = Mat4x4(InitRotY, angle) * Vec4(p, 1.0);
  return Vec3(p_.x, p_.y, p_.z);
}

Vec3 rotateZ(const Vec3 &p, double angle)
{
  const Vec4 p_ = Mat4x4(InitRotZ, angle) * Vec4(p, 1.0);
  return Vec3(p_.x, p_.y, p_.z);
}

First I tried the code with the eulerRotate() function which resembles OPs transformation order and got the following output:

17.9056 -6.99702 17.5622
17.2369 4.15094 19.0698

It seems that I reproduced the issue of OP correctly. Now, I defined FIX to use the correct transformation order and got the following output:

12.1019 -22.7589 3.68476
12.1019 -22.7589 3.68476

Live Demo on Wandbox

Now, both computations provide the same result as expected.

I thought an XYZ Euler angle meant you rotate it on X axis followed by Y axis followed by Z axis.

Not quite. As shown above, it has to be interpreted from right to left.

Scheff's Cat
  • 16,517
  • 5
  • 25
  • 45