How can you extract orientation from a transformation matrix?

Question

I have a 4x4 transformation matrix M, and I want to find out the shape of a sphere when transformed by M. (Sphere is at the origin and has radius 1.)

I know I can find the center by just multiplying M by (0,0,0,1).

However, the radius becomes a problem as M can squash and rotate the sphere. How can I find out the new radius(es) of the resulting ellipsoid? Is there any way to find out the orientation?

More specifically, I need to know the size of the bounding sphere which would enclose the transformed sphere. In other words, what is the maximum of |M*V - M*(0,0,0,1)|, where V is a unit vector (a point on the original sphere).

Answer 1

Mathematically, the quantity you're asking about is called the operator norm. Unfortunately, there's no simple formula for it. If it's a fully general affine transformation - for instance, if it could have an arbitrary combination of rotations and nonuniform scales, in any order - then I'm afraid there's nothing for it but to use singular value decomposition. If you apply SVD to your matrix then the largest singular value will be the maximum radius of the resulting ellipsoid. The other singular values will also be its other two radii, and the SVD procedure can also extract the orientation of the axes for you.

Implementing SVD is not for the faint of heart, as it involves finding eigenvalues. If all you want is the singular values themselves, they are the square roots of the eigenvalues of M^T * M. So if you have a 3x3 eigenvalue solver handy, or you don't mind writing one, you can use that. If you want to extract the orientations of the axes as well, then it gets more involved as you have to find eigenvectors too. On that Wikipedia article there is a list of links to libraries for doing SVD, one of which you may be able to use in your project.

If the form of your matrix is restricted in such a way that nonuniform scale happens at most once and is the first transform applied, i.e. is rightmost when you're using column vectors, then you can simplify this to just look at the lengths of the transformed axis vectors. In that case alone - i.e. a single nonuniform scale followed by any sequence of rotations, reflections, and uniform scales - looking at just the axis vectors will give you the right answer.

Answer 2

Maybe extract scale factors from the matrix and then use max value of its components. Using SRT (Scale-Rotation-Translation) matrix you can do this like that:

glm::mat4 m = ...;
// Extract col vectors of the matrix
glm::vec3 col1(m[0][0], m[0][1], m[0][2]);
glm::vec3 col2(m[1][0], m[1][1], m[1][2]);
glm::vec3 col3(m[2][0], m[2][1], m[2][2]);
//Extract the scaling factors
glm::vec3 scaling;
scaling.x = glm::length(col1);
scaling.y = glm::length(col2);
scaling.z = glm::length(col3);

float scaleFactor = MAX(scaling.x, MAX(scaling.y, scaling.z));

(based on http://wklej.org/id/950061/ - the name is decomposeTRS and not decomposeSRT because I use names besed on order which matrices are multiplied in OpenGL).

Now you can multiply original sphere radius by scaleFactor and you have your bounding sphere.