|
Sandro <sandro2109@aol.com> wrote in message <9fa07a19-86f8-4b95-b717-8e1294ff7349@i76g2000hsf.googlegroups.com>...
> Hey guys,
>
> I am given a quaternion, which rotates three orthogonal vectors (the
> coordinate system basically).
>
> And now I have to determine the minimal angle over which to rotate the
> new system to coincide with the old one
>
> Here is how far I have come:
>
> %convert to rotation matrix
> R = quat2dcm(q);
>
> %determine rotation axis
> if (abs(det(R - eye(3))) == 0)
> [V, D] = eig(R);
> if (norm(V(:, 1)) == 1)
> v = V(:, 1);
> elseif (norm(V(:, 2)) == 1)
> v = V(:, 2);
> elseif (norm(V(:, 3)) == 1);
> v = V(:, 3);
> end
> else
> v = null(R - eye(3));
> end
>
> %find another vector orthogonal to v
> if (v(3) == 0) && (v(2) == 0)
> w = cross(v, [v(2) v(1) v(3)])';
> else
> w = cross(v, [v(1) v(3) v(2)])';
> end
>
> ang = acos(dot(w, R*w) / (norm(w) * norm(R*w));
>
> When testing the results I find that sometimes I have to use the
> negative rotation axis to get the proper rotation. The angle seems to
> fit anyway, but could anyone explain why that is?
> And more importantly, when I have to go the eigenvector route, I
> pretty much always get complex eigenvectors, which don't seem to go
> together well with the rest of the plot. Either the rotation in the
> end doesn't fit or there might also be not eigenvector of the norm
> 1...
>
> How do I do it properly?
>
> Thx a lot
>
> Sandro
If you are starting from a quaternion then you can get the rotation angle and axis easily from the quaternion itself without converting to direction cosine matrix and using the eig function. The vector part of the quaternion contains the rotation axis times sin(angle / 2), and using an atan2 function on the magnitude of the vector part and the scalar part, which is cos(angle / 2), will yield half the rotation angle. To ensure you get the minimum single axis rotation angle, make sure the scalar part is non-negative first. e.g., assuming the scalar is in the first element (I presume your quat2dcm is from the Aerospace Toolbox and this is the convention for that function):
if q(1) < 0
p = -q;
else
p = q;
end
n = norm(p(2:4));
if n == 0
ang = 0;
v = [1 0 0];
else
ang = 2 * atan2(n,p(1));
v = p(2:4) / n;
end
James Tursa
|