Angle between 2 quaternions

How do we calculate angle between 2 quaternions? For example, what is the angle between x = ( 0.968, 0.008, -0.008, 0.252) and y = (0.382, 0.605, 0.413, 0.563)? How can I do it in matlab?

[Pawel Jastrzebski]

I think the identical problem has already been answered in a following thread:

https://uk.mathworks.com/matlabcentral/answers/150894-calculate-angles-between-two-intersecting-lines-using-the-slopes

[Patrícia Falcão]

This is about quaternions, not lines. It doesn't work, but thanks anyway!

[Patrícia Falcão]

I don't have matlab 2018

[Rik]

What is the angle between two points? You need to define that first, or explain the mathematical definition, as I couldn't find the definition with a quick search.

Answer 1

[James Tursa]

Assuming these represent attitude rotations from one coordinate frame to another, if you are simply asking what is the minimum rotation to take you from one quaternion to the other, you simply multiply one quaternion by the conjugate of the other and then pick off the rotation angle of the resulting quaternion.

But we really need to know what these quaternions represent, and what angle you are trying to recover, before we know what you want.

E.g., suppose x and y represent ECI->BODY rotation quaternions, and you want to know the minimum rotation angle that would take you from the x BODY position to the y BODY position. Then you could do this:

>> x = [ 0.968, 0.008, -0.008, 0.252]; x = x/norm(x); % ECI->BODY1
>> y = [ 0.382, 0.605,  0.413, 0.563]; y = y/norm(y); % ECI->BODY2
>> z = quatmultiply(quatconj(x),y) % BODY1->BODY2
z =
    0.5132    0.6911    0.2549    0.4405
>> a = 2*acosd(z(4)) % min angle rotation from BODY1 to BODY2
a =
  127.7227

But, again, these calculations are dependent on how I have the quaternions defined. Your specific case may be different.

[Silas Waxter]

The accepted answer here seems to work better.

[James Tursa]

Agreed for matching MATLAB. I should have noted this above. The Answer above assumes the scalar part is in the 4th element (does not match MATLAB toolboxes). The Answer in the other post assumes the scalar part is in the 1st element (which matches MATLAB toolboxes). They are both valid conventions and both are widely used in industry. Which is the correct method to use for a given application will depend on how the user has his/her quaternions defined.

[Brent Raiteri]

As both quatmultiply and quatconj require that the scalar is in the first column, I think that you need to use the following equation instead:

a = 2*acosd(z(1)) %correct

rater than:

a = 2*acosd(z(4)) %incorrect

I double-checked this answer using the code available here. Hopefully this is helpful to someone else. Below is the complete code I used to obtain answers from both sets of equations:

% Above
x = [ 0.968, 0.008, -0.008, 0.252]; 
y = [ 0.382, 0.605,  0.413, 0.563];
z = quatmultiply(quatconj(x),y); 
a = 2*atan2d(norm(z(2:4)),z(1))
%or
a2 = 2*acosd(z(1))

% Code from here, note that x and y become Q1 and Q2 respectively to be consistent with that code

Q1 = x;
Q2 = y;
%conj(Q1) did not change the signs of the vector of Q1
%convert Q1 and Q2 to quaternions
Q1 = quaternion(Q1);
Q2 = quaternion(Q2);
Q12 = conj(Q1) * Q2;
%convert back to 4 columns
Q12 = compact(Q12);
angle = 2*atan2d(norm(Q12(2:4)),Q12(1))
%or
angle2 = 2*acosd(Q12(1))
a =
       118.25
a2 =
       118.25
angle =
       118.25
angle2 =
       118.25

[James Tursa]

@Brent Raiteri Yes, correct. Looking back at this post I see I was using my own scalar-last functions for the reply. Which in retrospect is of course confusing given that the MATLAB functions are the same name. Good catch.

Answer 2

[Erik Blake]

Just as with vectors, the cosine of the rotation angle between two quaternions can be calculated as the dot product of the two quaternions divided by the 2-norm of the both quaternions. Normalization by the 2-norms is not required if the quaternions are unit quaternions (as is often the case when describing rotations).

As with vectors, the dot product is calculated by summing the products of the four elements of the quaternion.

Note that this calculation yields the full rotation angle, not the half-angle as when converting from quaternions to rotation vectors.

[James Tursa]

What you are describing is just the math for the scalar part of the quaternion multiply I described. This still only recovers the half angle of the rotation since it is the same calculation. I.e., the scalar part of the quatmultiply(quatconj(x),y) is this:

xs*ys - dot(-xv,yv) = xs*ys + dot(xv,yv)

And this is the same as the dot(x,y) operation you describe.

So if the quaternions represent two coordinate system transformations (my assumption), this result contains cos(half rotation angle), so acos( ) of this will recover the half rotation angle, not the full rotation angle. Or maybe I am misunderstanding you?

To be more robust for numerical issues, I could have done this:

z = quatmultiply(quatconj(x),y)

angle = 2 * atan2(norm(zv),zs)

where zv is the vector part of z and zs is the scalar part of z (either first or last element depending on convention).