3D rotations with 3D Gibbs vectors in place of 4D quaternions

[Deleted User]

#12726

August 3, 2017

Recently I discovered that Gibbs vectors multiply and rotate exactly like quaternions (because they are quaternions), but using only 3 values instead of 4.

Here's C code to try them yourself, it can be compiled with no compile parameters in both Clang and Visual Studio:

#include <windows.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define deg(x) ((x) * 57.295779513)
#define rad(x) ((x) *  0.017453292)

typedef struct {
  float x, y, z;
} v3;

typedef struct {
  float x, y, z, w;
} v4;

static inline float dot3(v3 v) {
  return v.x * v.x +
         v.y * v.y +
         v.z * v.z;
}

static inline float dot4(v4 v) {
  return v.x * v.x +
         v.y * v.y +
         v.z * v.z +
         v.w * v.w;
}

static inline float length3(v3 v) {
  return sqrt(dot3(v));
}

static inline float length4(v4 v) {
  return sqrt(dot4(v));
}

static inline v3 normalize3(v3 v) {
  float l = length3(v);
  if (l != 0) {
    v.x /= l;
    v.y /= l;
    v.z /= l;
  }
  return v;
}

static inline v4 normalize4(v4 v) {
  float l = length4(v);
  if (l != 0) {
    v.x /= l;
    v.y /= l;
    v.z /= l;
    v.w /= l;
  }
  return v;
}

static inline v4 qmul(v4 a, v4 b) {
  v4 o = {
    a.x * b.w + b.x * a.w + (a.y * b.z - b.y * a.z),
    a.y * b.w + b.y * a.w + (a.z * b.x - b.z * a.x),
    a.z * b.w + b.z * a.w + (a.x * b.y - b.x * a.y),
    a.w * b.w - (a.x * b.x + a.y * b.y + a.z * b.z)
  };
  return o;
}

static inline v3 gmul(v3 A, v3 B) {
  v4 a = (v4){A.x, A.y, A.z, 1};
  v4 b = (v4){B.x, B.y, B.z, 1};
  v4 o = {
    a.x * b.w + b.x * a.w + (a.y * b.z - b.y * a.z),
    a.y * b.w + b.y * a.w + (a.z * b.x - b.z * a.x),
    a.z * b.w + b.z * a.w + (a.x * b.y - b.x * a.y),
    a.w * b.w - (a.x * b.x + a.y * b.y + a.z * b.z)
  };
  return (v3){
    o.x / o.w,
    o.y / o.w,
    o.z / o.w
  };
}

static inline v3 qrot(v3 v, v4 q) {
  v4 o = qmul(qmul(q, (v4){v.x, v.y, v.z, 0}), (v4){-q.x, -q.y, -q.z, q.w});
  return (v3){o.x, o.y, o.z};
}

static inline v3 grot(v3 v, v3 g) {
  v3 o = gmul(gmul(g, v), (v3){-g.x, -g.y, -g.z});
  return o;
}

static inline v4 qaxisrad(v3 axis, float rad) {
  float s = (float)sin(rad / 2.0);
  float c = (float)cos(rad / 2.0);
  v4 o;
  o.x = s * axis.x;
  o.y = s * axis.y;
  o.z = s * axis.z;
  o.w = c;
  return o;
}

static inline v3 gaxisrad(v3 axis, float rad) {
  float t = (float)tan(rad / 2.0);
  v3 o;
  o.x = t * axis.x;
  o.y = t * axis.y;
  o.z = t * axis.z;
  return o;
}

static inline v3 qaxis(v4 q) {
  return normalize3((v3){q.x, q.y, q.z});
}

static inline v3 gaxis(v3 g) {
  return normalize3(g);
}

static inline float qrad(v4 q) {
  return atan2(length3((v3){q.x, q.y, q.z}), q.w) * 2.0;
}

static inline float grad(v3 g) {
  return atan(length3(g)) * 2.0;
}

static inline v3 qtog(v4 q) {
  return (v3){
    q.x / q.w,
    q.y / q.w,
    q.z / q.w
  };
}

static inline v4 gtoq(v3 g) {
  return qaxisrad(gaxis(g), grad(g));
}

static inline float lerp(float a, float b, float t) {
  return (1 - t) * a + t * b;
}

static inline v4 qlerp(v4 q1, v4 q2, float t) {
  v4 o;
  o.x = lerp(q1.x, q2.x, t);
  o.y = lerp(q1.y, q2.y, t);
  o.z = lerp(q1.z, q2.z, t);
  o.w = lerp(q1.w, q2.w, t);
  return normalize4(o);
}

static inline v3 glerp(v3 g1, v3 g2, float t) {
  v3 o;
  o.x = lerp(g1.x, g2.x, t);
  o.y = lerp(g1.y, g2.y, t);
  o.z = lerp(g1.z, g2.z, t);
  return o;
}

static inline v3 gderp(v3 g1, v3 g2, float t) {
  return qtog(qlerp(gtoq(g1), gtoq(g2), t));
}

int main() {
  srand(time(NULL));

  v3 axis_1 = normalize3((v3){1, 1, 0});
  v3 axis_2 = normalize3((v3){0, 1, 1});
  float rad_1 = 10.0;
  float rad_2 = 10.0;

  {
    printf("Comparison:\n");
    printf("\n");

    v4 q1 = qaxisrad(axis_1, rad(rad_1));
    v4 q2 = qaxisrad(axis_2, rad(rad_2));
    v4 q = qmul(q1, q2);

    v3 g1 = gaxisrad(axis_1, rad(rad_1));
    v3 g2 = gaxisrad(axis_2, rad(rad_2));
    v3 g = gmul(g1, g2);

    printf("q: %.4f %.4f %.4f %.4f\n", q.x, q.y, q.z, q.w);
    printf("g: %.4f %.4f %.4f\n", g.x, g.y, g.z);
    printf("\n");

    for (int i = 0; i < 10; i += 1) {
      v3 point = {
        (((float)rand() / RAND_MAX) * 2.0 - 1.0) * 100.0,
        (((float)rand() / RAND_MAX) * 2.0 - 1.0) * 100.0,
        (((float)rand() / RAND_MAX) * 2.0 - 1.0) * 100.0
      };
      v3 Q = qrot(point, q);
      v3 G = grot(point, g);
      printf("Q: %.4f %.4f %.4f\n", Q.x, Q.y, Q.z);
      printf("G: %.4f %.4f %.4f,\t", G.x, G.y, G.z);
      printf("error: %.4f %.4f %.4f\n", fabs(Q.x - G.x),
                                        fabs(Q.y - G.y),
                                        fabs(Q.z - G.z));
      printf("\n");
    }

    v3 q_axis = qaxis(q);
    v3 g_axis = gaxis(g);
    float q_angle = deg(qrad(q));
    float g_angle = deg(grad(g));

    printf("q_axis:  %.4f %.4f %.4f\n", q_axis.x, q_axis.y, q_axis.z);
    printf("g_axis:  %.4f %.4f %.4f\n", g_axis.x, g_axis.y, g_axis.z);
    printf("q_angle: %.4f\n", q_angle);
    printf("g_angle: %.4f\n", g_angle);
    printf("\n");
  }

  {
    printf("Interpolation:\n");
    printf("\n");

    float deg_diff = 90;

    v3 g1 = gaxisrad((v3){0, 1, 0}, rad(0));
    v3 g2 = gaxisrad((v3){0, 1, 0}, rad(deg_diff));

    v4 q1 = qaxisrad(gaxis(g1), grad(g1));
    v4 q2 = qaxisrad(gaxis(g2), grad(g2));

    printf("q1: %.4f %.4f %.4f %.4f\n", q1.x, q1.y, q1.z, q1.w);
    printf("q2: %.4f %.4f %.4f %.4f\n", q2.x, q2.y, q2.z, q2.w);
    printf("\n");

    printf("g1: %.4f %.4f %.4f\n", g1.x, g1.y, g1.z);
    printf("g2: %.4f %.4f %.4f\n", g2.x, g2.y, g2.z);
    printf("\n");

    v3 q1_axis = qaxis(q1);
    v3 q2_axis = qaxis(q2);
    float q1_angle = deg(qrad(q1));
    float q2_angle = deg(qrad(q2));

    v3 g1_axis = gaxis(g1);
    v3 g2_axis = gaxis(g2);
    float g1_angle = deg(grad(g1));
    float g2_angle = deg(grad(g2));

    printf("q1_axis:  %.4f %.4f %.4f\n", q1_axis.x, q1_axis.y, q1_axis.z);
    printf("q2_axis:  %.4f %.4f %.4f\n", q2_axis.x, q2_axis.y, q2_axis.z);
    printf("q1_angle: %.4f\n", q1_angle);
    printf("q2_angle: %.4f\n", q2_angle);
    printf("\n");

    printf("g1_axis:  %.4f %.4f %.4f\n", g1_axis.x, g1_axis.y, g1_axis.z);
    printf("g2_axis:  %.4f %.4f %.4f\n", g2_axis.x, g2_axis.y, g2_axis.z);
    printf("g1_angle: %.4f\n", g1_angle);
    printf("g2_angle: %.4f\n", g2_angle);
    printf("\n");

    v3 P = normalize3((v3){
      ((float)rand() / RAND_MAX) * 2.0 - 1.0,
      ((float)rand() / RAND_MAX) * 2.0 - 1.0,
      ((float)rand() / RAND_MAX) * 2.0 - 1.0
    });

    for (int i = 0, size = 10; i <= size; i += 1) {

      float t = i / (float)size;

      v4 q = qlerp(q1, q2, t);
      v3 g = glerp(g1, g2, t);

      v3 q_axis = qaxis(q);
      v3 g_axis = gaxis(g);

      float q_angle = deg(qrad(q));
      float g_angle = deg(grad(g));

      printf("[%.4f] q_axis: %.4f %.4f %.4f, ", t, q_axis.x, q_axis.y, q_axis.z);
      printf("angle: %.4f\n", q_angle);
      printf("[%.4f] g_axis: %.4f %.4f %.4f, ", t, g_axis.x, g_axis.y, g_axis.z);
      printf("angle: %.4f,\t", g_angle);
      printf("error: %.4f %.4f %.4f, %.4f\n", fabs(q_axis.x - g_axis.x),
                                              fabs(q_axis.y - g_axis.y),
                                              fabs(q_axis.z - g_axis.z),
                                              fabs(q_angle  - g_angle));
      v3 Q = qrot(P, q);
      v3 G = grot(P, g);
      printf("[%.4f] Q: %.4f %.4f %.4f\n", t, Q.x, Q.y, Q.z);
      printf("[%.4f] G: %.4f %.4f %.4f, ", t, G.x, G.y, G.z);
      printf("error: %.4f %.4f %.4f\n", fabs(Q.x - G.x),
                                        fabs(Q.y - G.y),
                                        fabs(Q.z - G.z));
      printf("\n");
    }
  }

  {
    printf("Compression:\n");
    printf("\n");

    for (int i = 0; i < 10; i += 1) {
      v3 axis = normalize3((v3){
        ((float)rand() / RAND_MAX) * 2.0 - 1.0,
        ((float)rand() / RAND_MAX) * 2.0 - 1.0,
        ((float)rand() / RAND_MAX) * 2.0 - 1.0
      });
      float rad = rad(((float)rand() / RAND_MAX) * 360.0);
      v4 q = qaxisrad(axis, rad);
      v3 g_comp = qtog(q);
      v4 q_uncomp = gtoq(g_comp);

      v3 q_axis = qaxis(q);
      float q_angle = deg(qrad(q));

      v3 q_uncomp_axis = qaxis(q_uncomp);
      float q_uncomp_angle = deg(qrad(q_uncomp));

      printf("q_axis:         %.4f %.4f %.4f\n", q_axis.x, q_axis.y, q_axis.z);
      printf("q_angle:        %.4f\n", q_angle);
      printf("q_uncomp_axis:  %.4f %.4f %.4f\n", q_uncomp_axis.x, q_uncomp_axis.y, q_uncomp_axis.z);
      printf("q_uncomp_angle: %.4f\n", q_uncomp_angle);
      printf("\n");
    }
  }

  {
    printf("Direction:\n");
    printf("\n");

    v4 q1 = qaxisrad((v3){-0.5774, -0.5774, -0.5774}, rad(1));
    printf("q1: %.4f %.4f %.4f %.4f\n", q1.x, q1.y, q1.z, q1.w);
    v3 q1_axis = qaxis(q1);
    float q1_angle = deg(qrad(q1));
    printf("q1_axis:  %.4f %.4f %.4f\n", q1_axis.x, q1_axis.y, q1_axis.z);
    printf("q1_angle: %.4f\n", q1_angle);
    printf("\n");

    v4 q2 = qaxisrad((v3){0.5774, 0.5774, 0.5774}, rad(359));
    printf("q2: %.4f %.4f %.4f %.4f\n", q2.x, q2.y, q2.z, q2.w);
    v3 q2_axis = qaxis(q2);
    float q2_angle = deg(qrad(q2));
    printf("q2_axis:  %.4f %.4f %.4f\n", q2_axis.x, q2_axis.y, q2_axis.z);
    printf("q2_angle: %.4f\n", q2_angle);
    printf("\n");

    v3 g1 = gaxisrad((v3){-0.5774, -0.5774, -0.5774}, rad(1));
    printf("g1: %.4f %.4f %.4f\n", g1.x, g1.y, g1.z);
    v3 g1_axis = gaxis(g1);
    float g1_angle = deg(grad(g1));
    printf("g1_axis:  %.4f %.4f %.4f\n", g1_axis.x, g1_axis.y, g1_axis.z);
    printf("g1_angle: %.4f\n", g1_angle);
    printf("\n");

    v3 g2 = gaxisrad((v3){0.5774, 0.5774, 0.5774}, rad(359));
    printf("g2: %.4f %.4f %.4f\n", g2.x, g2.y, g2.z);
    v3 g2_axis = gaxis(g2);
    float g2_angle = deg(grad(g2));
    printf("g2_axis:  %.4f %.4f %.4f\n", g2_axis.x, g2_axis.y, g2_axis.z);
    printf("g2_angle: %.4f\n", g2_angle);
    printf("\n");

    v3 P = {
      (((float)rand() / RAND_MAX) * 2.0 - 1.0) * 100.0,
      (((float)rand() / RAND_MAX) * 2.0 - 1.0) * 100.0,
      (((float)rand() / RAND_MAX) * 2.0 - 1.0) * 100.0
    };
    v3 Q = qrot(P, q2);
    v3 G = grot(P, g2);
    printf("P: %.4f %.4f %.4f\n", P.x, P.y, P.z);
    printf("Q: %.4f %.4f %.4f\n", Q.x, Q.y, Q.z);
    printf("G: %.4f %.4f %.4f\n", G.x, G.y, G.z);
    printf("\n");
  }
}

Basically, Gibbs vectors encode quaternion's scalar part to the vector part with divisions on all vector elements. That's it.

You can also try them in GpuLib's Instancing and MRT example, `build_vs2015.bat` and `build_vs2017.bat` are available.

Edited by Procedural on August 12, 2017, 6:04pm

Mārtiņš Možeiko

#12727

August 3, 2017

Hmm... This looks very interesting!

If I understand correctly this cannot represent 180° rotations because of division by 0, right? Probably this is not issue for some cases, but with this limitation it cannot completely replace matrix.

In your example this works because tan value is approximated to 13245402.0 which is "very large" number (but not infinity). I'm wondering if for some cases tan will return ±∞ and all calculations further will be incorrect... Maybe it will never happen because no floating point value will get close enough to ±π/2.

Edited by Mārtiņš Možeiko on August 3, 2017, 7:01am

[Deleted User]

#12728

August 3, 2017

mmozeiko
but with this limitation it cannot completely replace matrix.

In the real world, in the actual reality of real time graphics, when we multiply at least 2 orientations, quaternions will always be cheaper than matrices, so even if Gibbs vectors fail somehow you can always fallback to them.

mmozeiko
I'm wondering if for some cases tan will return ±∞ and all calculations further will be incorrect...

Right, in that case I would write a custom tan that won't return infinity.

mmozeiko
Maybe it will never happen because no floating point value will get close enough to ±π/2.

This is correct too, in GpuLib example I never saw any issues like that.
What are the chances of a camera rotation or animation hitting one of the axes exactly at floating point 0b0000_0000_0000_0000_0000_0000_0000_0000 or 0b0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000? Close to 0. :)

ratchetfreak

#12729

August 3, 2017

there is also the chance that a.w * b.w - (a.x * b.x + a.y * b.y + a.z * b.z) = 1-dot(A, B) will end up 0.

This is a very likely spot for catastrophic cancellation which will ruin your accuracy.

Mārtiņš Možeiko

#12730

August 3, 2017

Procedural
What are the chances of a camera rotation or animation hitting one of the axes exactly at floating point 0b0000_0000_0000_0000_0000_0000_0000_0000 or 0b0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000? Close to 0. :)

The problem is not about hitting value exactly 0 or 180 or any other. The problem with floats is that if you get close enough to some discontinuity which can happen with much higher probability, the calculations will become very unstable. The results produced will become very large or ∞ or NaN. And later calculations will propagate this error further. As a result you'll get weird glitches in rendering or sudden jumps in rotation which could be super annoying to debug or eliminate.

Edited by Mārtiņš Možeiko on August 3, 2017, 9:12am

[Deleted User]

#12731

August 3, 2017

mmozeiko
if you get close enough

Isn't any divide in any math code ever suffers from the same problem?

I will not guard my code against divide by 0 anywhere other than in vector length procedure... Will you?

Ginger Bill

#12732

August 3, 2017

With unit quaternions, there is no worry with dividing by zero. However, with Rodrigues/Gibbs Parameters (RPs) as Mārtiņš has pointed out, 180° rotations cannot be represented, even in floating pointer.

Modified Rodrigues parameters (MRPs) can be used to solve this issue but it only pushes the problem to 360°.

These objects are an exponential mapping of the 3D rotations (exponential map of so(3), the Lie algebra of SO(3)). To represent a 3D rotation, you need at a minimum 4 degrees of freedom (pieces of information): the axis of rotation (x, y, z) and the angle of rotation (θ). Quaternions encode this with no information loss as they have 4 values.

RPs rely on the fact that the axis of rotation is unit length thus you can "lose" a degree of freedom. Even with this fact, you lose other niceties as you have found.

Quaternions seem to be the best solution to the problem for representing 3D rotations on a computer.

[Deleted User]

#12733

August 3, 2017

gingerBill
To represent a 3D rotation, you need at a minimum 4 degrees of freedom (pieces of information): the axis of rotation (x, y, z) and the angle of rotation (θ). Quaternions encode this with no information loss as they have 4 values.

The code is updated for extracting an axis and angle from a Gibbs vector.

Mārtiņš Možeiko

#12734

August 3, 2017

Procedural
mmozeiko
if you get close enough

Isn't any divide in any math code ever suffers from the same problem?

I will not guard my code against divide by 0 anywhere other than in vector length procedure... Will you?

That's exactly my point - you don't guard against those because there is no solution. Nor in vector normalization, nor for this rotation.

For vectors normalization this is not a problem, because you cannot normalize vector without length. You design your algorithms / functionality in a way that doesn't produce these invalid vectors. And if they occur, that's a bug or impossible situation.

For example, let's say you need normalized direction from some object to light source to use in lighting equation. In real life this will always be non-0 vector when you subtract object position from light position. Sure in case the positions are the same (or almost the same) the vector will be non-normalizeable (the result will be huge or NaN). But that's ok, in such case it is ok to have a glitch. Because in real-world this never happens - light cannot be at same position as object. And if it is, then nobody knows what to expect result to be. So the software can do whatever it wants, no need to "guard" against divide by 0.

But in Gibbs rotations the algorithm is designed to not work for 180° or in case ratchetfreak mentions. But these cases are perfectly normal, the expected result is strictly defined, but they will misbehave.

Nevertheless, I like the idea. I'm sure this representation can work just fine for some use cases. Just not the 100% replacement for rotation matrix/quaternion.

Edited by Mārtiņš Možeiko on August 3, 2017, 4:57pm

[Deleted User]

#12745

August 4, 2017

mmozeiko
But in Gibbs rotations the algorithm is designed to not work for 180° or in case ratchetfreak mentions.

ratchetfreak's problem is not a problem for unit quaternions (and Gibbs vectors are of unit length) as Ginger Bill mentioned, you can check it yourself by printing dot product in `gmul` procedure.

As of 180°, the only solution to this problem is:

Deal With It™ :)

I mean, the context is everything. Mouse / gamepad camera rotations? 180° rotations are not possible here. Character animation? If you can think of a character that can rotate its limbs for more than 180° in one operation? I personally can't. I don't know where else you can use >180° rotations, leave a comment below which real world example actually require this. In the context of gamedev I can't think of any!

Ginger Bill

#12746

August 4, 2017

It actually is a problem. RPs can be zero unlike quaternions.

Other than curiosity, I see no benefit to using RPs over quaternions for practical use. RPs cause more problems than they solve. Other than memory footprint, they have no advantage over quaternions.

[Deleted User]

#12747

August 4, 2017

gingerBill
It actually is a problem. RPs can be zero unlike quaternions.

It is not, until proven otherwise. I can't make dot product return 1, if you know how let us know.

gingerBill
Other than curiosity, I see no benefit to using RPs over quaternions for practical use. RPs cause more problems than they solve.

By more problems you mean one? Saving space is not practical enough for you?

gingerBill
Other than memory footprint, they have no advantage over quaternions.

This is one MASSIVE advantage if we talk space for hundreds of thousands of vectors, and one 180° in one op problem no one even proved to be real so far, because in the real world, in animation, people always lerp.

Ginger Bill

#12748

August 4, 2017

First off, how do you encode an zero rotation, half a turn or a full turn around a particular axis? At these particular rotations (depending on the particular encoding), it removes all the information about its axis. With quaternions, all angles of rotation are unique and still keep all the information unlike RPs where the same rotation can be represented in multiple ways.

You lose information when you reduce the degrees of freedom.

If saving space is that much of a concern, multiplying the axis by its angle will achieve the same thing whilst also removing most of the angle rotation problems but then the calculations become a little more costly. Another way to save space is to use a smaller type like 32-bit float or a 16-bit float. But personally, I don't think that saving one extra component outweighs the costs you get from the extra instructions needed to calculate the products.

Show me the evidence that the memory saving cost is worth it over the performance one in the real world and not the theoretical one.

Edited by Ginger Bill on August 4, 2017, 10:56am

ratchetfreak

#12749

August 4, 2017

Procedural

ratchetfreak's problem is not a problem for unit quaternions (and Gibbs vectors are of unit length) as Ginger Bill mentioned, you can check it yourself by printing dot product in `gmul` procedure.

eh no. The entire point of Gibbs is that the length encodes the rotation angle, so they are not of unit length.

Procedural

As of 180°, the only solution to this problem is:

Deal With It™ :)

I mean, the context is everything. Mouse / gamepad camera rotations? 180° rotations are not possible here. Character animation? If you can think of a character that can rotate its limbs for more than 180° in one operation? I personally can't. I don't know where else you can use >180° rotations, leave a comment below which real world example actually require this. In the context of gamedev I can't think of any!

Rotations are composed from a base orientation. So any time the object is 180° in the other direction the Gibbs will break down.

For example take a car on a flat surface turning in a circle.

At some point it needs to be rotated 180° compared to the starting position. At that point your car will glitch out and depending on how you created the turning transform (iteratively or closed function) will add inaccuracies for the future. So it may recover or it will remain with a glitched transform until it gets reset.

You can't lerp or slerp Gibbs and get the same result as the quaternions unless you normalize them to quaternions. So if you want space savings then just drop the w component and use sqrt(1-dot(q3, q3)) to get it back when needed.

Edited by ratchetfreak on August 4, 2017, 11:02am

Ginger Bill

#12750

August 4, 2017

Quaternions are used for storing an orientation, not just rotations. RPs cannot do this. I'll repeat this again:

You lose information when you reduce the degrees of freedom.

Exponential maps in general have this problem.

Matrix: Easily composed; cannot interpolate
Euler Angle: Not easily composed; you can sometimes interpolate them; can joint limit; has Gimbal lock
Angle Axis: Not easily composed; can interpolate; can joint limit; avoids gimbal lock
Quaternion: Easily composed; can interpolate; can joint limit (but difficult); avoids gimbal lock
Exponential Map: Easily composed; cannot interpolate; cannot joint limit easily; can have gimbal lock

Edited by Ginger Bill on August 4, 2017, 11:12am