Relative motion between Quaternions

I'm working on Quaternions. I wanna find out the relative motion between quaternions in MATLAB.

For example, I've

QuatA = [1283842729  -831786718  1290396116  1271562329] and  
QuatB = [1313190417  1274713520  -845615766  1306114950]

How can I find how much QuatB is rotated with respect to QuatA?

Any help is appreciated.

Thanks.

[James Tursa]

What do these quaternions represent?

[BAN]

Here the quaternions represent the motion.

[James Tursa]

Angular motion? You mention rotation, which normally means I would have expected unit quaternions, but your quaternions are not unit. Hence I am trying to understand what exactly they represent.

[BAN]

Yeah, I'm looking for angular motion. I'm sorry for not adding the proper tags. I'll edit.

[BAN]

We can easily get the unit quaternions by normalizing it. I don't have any problem converting it to unit quaternions. But I don't understand how to get the relative motion.

Answer 1

[James Tursa]

In general, for unit quaternions you would multiply the conjugate of one times the other, and then extract the angle and rotation axis from that. Be sure to do a quaternion multiply, not a regular multiply. Also be sure the quaternion multiply routine you are using assumes the scalar is in the same spot (1st or 4th) as the quaternions you have. You would have to divide your quaternions above by their respective norms to apply this technique, but not knowing what these quaternions represent I can't advise if this is the correct procedure or not.

[BAN]

I'm sorry. Can you please explain it.

[James Tursa]

E.g., assuming quaternion scalar is in the 1st spot to match the Aerospace Toolbox:

QuatA = [1283842729  -831786718  1290396116  1271562329];
QuatB = [1313190417  1274713520  -845615766  1306114950];
QuatAu = QuatA / norm(QuatA);
QuatBu = QuatB / norm(QuatB);
QuatDiff = quatmultiply(quatconj(QuatAu),QuatBu);

Then QuatDiff will be equivalent to [cos(ang/2) sin(ang/2)*eigx]

Where ang is the angle you are looking for, and eigx is a 1x3 unit vector about which the rotation is taking place. So this can recover the angle:

ang = 2 * atan2(norm(QuatDiff(2:4)),QuatDiff(1));

Again, this all assumes that normalizing the quaternions has meaning in your setup, and that the scalar is in the 1st spot.

[BAN]

Everything looks good James. Sorry to bother you so much. I'm new to this.

I've a small question. You've mentioned that, "QuatDiff will be equivalent to [cos(ang/2) sin(ang/2)*eigx]"

Here QuatDiff is a 1x4 vector. How could eigx be a 1x3 vector.

[James Tursa]

Sorry I was not more clear. QuatDiff(1) is equal to cos(ang/2), and QuatDiff(2:4) is equal to sin(ang/2)*eigx. So, only the last three elements of the 1x4 vector are used for the sin(ang/2)*eigx part. I could have written it like this (with a comma):

[ cos(ang/2) , sin(ang/2)*eigx ]

[BAN]

It worked out. Everything looks great.

Thank you for your support.

[manikya valli]

I am having similar problem. I want to find difference between reference quaternion and measured quaternion. I want to verify by converting into euler angles. But the difference quaternion when converted back to euler angles is not expected euler angles. Here is the code: q1=angle2quat(deg2rad(30),deg2rad(40),deg2rad(50));

q2=angle2quat(deg2rad(10),deg2rad(10),deg2rad(10));

q3=quatmultiply(quatconj(q2),q1);

[r,p,y]=quat2angle(q3);

rad2deg(r);rad2deg(p);rad2deg(y); result I am getting is

r,p,y=11.6921,33.0794,34.6602

but I have to get 30-10,40-10,50-10=20,30,40 respectively. Why I am not getting this result?

[Aitor Burdaspar]

Hello Manikya,

I have exactly the same problem as yours. Did you finally solve it?

Best regards,

[James Tursa]

There is a fundamental misunderstanding of how to work with quaternions and Euler Angle sequences. In the first place, Euler Angle sequences do not behave linearly. That is, if you have a (yaw1,pitch1,roll1) sequence followed by a (yaw2,pitch2,roll2) sequence, you should not be expecting the result to be the same as a (yaw1+yaw2,pitch1+pitch2,roll1+roll2) sequence. Sequential rotations simply do not behave that way. So the subtractions you are doing (30-10, 40-10, 50-10) to get your "expected" result is simply not true with regards to sequential rotations.

In particular, I will use a couple of example coordinate frames to illustrate the issue.

Suppose q1 and q2 are both ECI_to_BODY quaternions. I.e.,

q1 = ECI_to_BODY1 (BODY1 being the BODY frame orientation for q1)

q2 = ECI_to_BODY2 (BODY2 being the BODY frame orientation for q2)

Then

q3 = conj(ECI_to_BODY2) * ECI_to_BODY1

= BODY2_to_ECI * ECI_to_BODY1

= BODY2_to_BODY1

So q3 represents a rotation between the two end BODY frames for the quaternions. Note that q3 is not another ECI_to_BODY quaternion. Since it is not another ECI_to_BODY quaternion, you should not expect the angles associated with it to be the same as what you get from a straight subtraction of your q1 and q2 angles.

Side note: The default order for the angles is y,p,r and not r,p,y.

Answer 2

[Mo]

How is this method used for a matrix? It seems "norm(A)" for a quaternion matrix is different. I appreciate any help.

[Luca Pozzi]

If A is an nx4 matrix (a 'column' of quaternions) the following code should work:

for i=1:size(A,1) % for each row of A
    Anorm(i,:)=norm(A(i,:)); % dividing a row by its norm
end

norm(A) returns the norm of the whole matrix A, the for cycle here above allows you to normalize every single row of the matrix (i.e. every quaternion) instead. Hope this helps!