Get rotation matrix from quaternion by quaternion.as_rotation_matrix()

Rotation matrix derived from quaternion by quaternion.as_rotation_matrix() The ratation matrix is not the same as neither function quat2dcm in MATLAB nor the regular formula of quaternion and rotation matrix written in the textbook. Actually, the two result matrix from quaternion.as_rotation_matrix() and quat2dcm in MATLAB are transposed each other.

Possible reason Looking into more details, I think the matrix should be

return np.array([
                [1 - 2*(q.y**2 + q.z**2),   2*(q.x*q.y + q.z*q.w),      2*(q.x*q.z - q.y*q.w)],
                [2*(q.x*q.y - q.z*q.w),     1 - 2*(q.x**2 + q.z**2),    2*(q.y*q.z + q.x*q.w)],
                [2*(q.x*q.z + q.y*q.w),     2*(q.y*q.z - q.x*q.w),      1 - 2*(q.x**2 + q.y**2)]
            ])

instead of the original code (from row 7 below)

def as_rotation_matrix(q):
if q.shape == () and not isinstance(q, np.ndarray):  # This is just a single quaternion
        n = q.norm()
        if n == 0.0:
            raise ZeroDivisionError("Input to `as_rotation_matrix({0})` has zero norm".format(q))
        elif abs(n-1.0) < _eps:  # Input q is basically normalized
            return np.array([
                [1 - 2*(q.y**2 + q.z**2),   2*(q.x*q.y - q.z*q.w),      2*(q.x*q.z + q.y*q.w)],
                [2*(q.x*q.y + q.z*q.w),     1 - 2*(q.x**2 + q.z**2),    2*(q.y*q.z - q.x*q.w)],
                [2*(q.x*q.z - q.y*q.w),     2*(q.y*q.z + q.x*q.w),      1 - 2*(q.x**2 + q.y**2)]
            ])
        else:  # Input q is not normalized
            return np.array([
                [1 - 2*(q.y**2 + q.z**2)/n,   2*(q.x*q.y - q.z*q.w)/n,      2*(q.x*q.z + q.y*q.w)/n],
                [2*(q.x*q.y + q.z*q.w)/n,     1 - 2*(q.x**2 + q.z**2)/n,    2*(q.y*q.z - q.x*q.w)/n],
                [2*(q.x*q.z - q.y*q.w)/n,     2*(q.y*q.z + q.x*q.w)/n,      1 - 2*(q.x**2 + q.y**2)/n]
            ])
    else:  # This is an array of quaternions
            ...

PS Maybe there are more in-depth reasons for the original code. But if I was wrong, please let me know. Thank you

Pretty much every implementation of anything having to do with rotations makes different choices for basically anything you can think of. Obviously, Euler angles are the worst in this way, but there are numerous choices to be made even with just quaternions and rotation matrices. For example, even the simple formula i*j*k = -1 is not universally accepted; some people use i*j*k = 1. Then there's the question of how to rotate a vector v with a quaternion Q: should the rotated vector be Q * v * Q.conjugate() or Q.conjugate() * v * Q? I've even seen disagreement about how a rotation matrix should be applied: does the matrix M rotate the vector v to M * v or v * M or even v * M.inverse()? Finally, there's the question of whether the transformation is "active" or "passive" — whether we're rotating a vector relative to a fixed basis, or rotating the basis in which a fixed vector is expressed.

There's no "right" answer to any of these questions; they're all just conventions we choose and stick with. In this package, I've chosen the first option in each of these questions. Your suggested version is just the transpose of my version, which should give you a hint that Matlab and I are just using "flipped" conventions relative to each other. Also note that the transpose of a rotation matrix is equal to its inverse (because it's an orthogonal matrix) — which is your hint that one of these conventions is to blame.

Looking at the documentation for quat2dcm, I see this:

The direction cosine matrix performs the coordinate transformation of a vector in inertial axes to a vector in body axes.

That's essentially saying that the output performs a passive transformation, which is not what this code does.

A little more googling leads me to another Matlab function: rotmat, which actually allows you to explicitly select the rotationType. Towards the bottom of that page they even show you the two matrices that could result from the two options from rotationType. The second of these — which they call the "point" rotation — agrees with my code, and my conventions.

If you want this code to produce the same results as quat2dcm, just remember that we are talking about opposite rotations, and use quaternion.as_rotation_matrix(Q.conjugate()).

But it's also a fair point that my conventions are not documented well enough. I'll add something to the docstrings.

moble / quaternion

Get rotation matrix from quaternion by quaternion.as_rotation_matrix() #145