In 3d, I have a $4\times4$ matrix $M$, which has only a rotation part and a translation part. In other words, I can compute $X'=RX+T$ ( with $R$ a $3\times3$ rotation matrix, $T$ a vector for the translation, $X$ a point and $X'$ its image).

Can I write the $4\times4$ matrix $M$ as a rotation around an arbitrary axis (not necessarily through origin)? I would write it as $X'=R(X-X_a)+X_a$ with $X_a$ a point on the axis, and $R$ the same rotation matrix.

If yes, I would like to get one point $X_a$, i.e. a point on the axis. My first thought was to solve the system $RX_a+T=X_a$. It's a singular system because we look for an axis. So I set $z=0$ and solve for $x$ and $y$ to get one point. Is there a mistake somewhere?

I thought about another way to get this axis, but I don't understand why it does not work. I take two random points, $X_1$ and $X_2$ and compute their image $X_1'$ and $X_2'$ by $M$. I assume the axis is on the plane perpendicular to $X_1X_1'$ through the middle of $X_1X_1'$. Indeed, $X_1$ would be at the same distance to the axis than $X_1'$. I assume this also for $X_2X_2'$. So I get the axis with the intersection between both planes. What is wrong in this approach?

$\begingroup$I would like to clarify: I know that we can build a matrix around an arbitrary line and how to do it. My question is: I have an input matrix M (with only rotation and translation, no scale, no shear), is there a line u+tv and an angle theta, such that the matrix of the rotation around this line is M? And how to find the line and theta?$\endgroup$
– AlexJul 26 '16 at 9:49

2 Answers
2

A rotation about an arbitrary line $\mathbf u+t\mathbf v$ in $\mathbb R^3$ can be decomposed into a translation to the line, a rotation about an axis through the origin, and a translation back. In homogeneous coordinates, this can be represented by the $4\times4$ matrix $$M=\left[\begin{array}{c|c}I&\mathbf u\\\hline 0&1\end{array}\right]\left[\begin{array}{c|c}R_{\mathbf v,\theta}&0\\\hline 0&1\end{array}\right]\left[\begin{array}{c|c}I&-\mathbf u\\\hline 0&1\end{array}\right]=\left[\begin{array}{c|c}R_{\mathbf v,\theta}&(I-R_{\mathbf v,\theta})\mathbf u\\\hline 0&1\end{array}\right].$$ It’s clear from inspection that this is equivalent to a rotation about a parallel axis that goes through the origin followed by a translation. You have a matrix of the form $$\left[\begin{array}{c|c}R_{\mathbf v,\theta}&\mathbf w\\\hline 0&1\end{array}\right]\tag{*}$$ so recovering the axis of rotation is simply a matter of solving $(I-R_{\mathbf v,\theta})\mathbf u=\mathbf w$ for $\mathbf u$. Unless $\theta$ is an integer multiple of $2\pi$, the matrix $I-R_{\mathbf v,\theta}$ has rank 2, so the solution set—if it exists—will be a line. Note that this equation’s having a solution is a necessary condition for a matrix of the form (*) to represent a rotation. Also, $(I-R_{\mathbf v,\theta})\mathbf u$ is the difference between $\mathbf u$ and its image under $R_{\mathbf v,\theta}$ and so is orthogonal to $\mathbf v$. This gives us another necessary condition: $\mathbf v\cdot\mathbf w=0$.

You can also work directly with the $4\times4$ matrix $M$: $1$ is an eigenvalue and $M$ is the identity on its eigenspace. For most angles $\ker(I-M)$ will be two-dimensional, and so correspond to a line in $\mathbb R^3$. If you use the usual method of row-reduction to find a basis for the kernel, one of the resulting vectors will have $0$ for its last coordinate—it’s a point at infinity. You can take the first three coordinates of this vector as the direction vector $\mathbf v$, while the other vector in the kernel basis can be taken as the fixed vector $\mathbf u$.

The rotation angle can be recovered from either the full matrix $M$ via $\operatorname{tr}M=2+2\cos\theta$, or from the submatrix $R$ via $\operatorname{tr}R=1+2\cos\theta$. These identities come from the facts that the trace of a matrix is equal to the sum of its eigenvalues (with the appropriate number of repetitions), and that the eigenvalues of these matrices are $1$, $e^{i\theta}$ and $e^{-i\theta}$.

Addition: Your idea of using two points and their images to define a pair of planes that intersect along the axis of rotation is sound. It seems like a lot more work than recovering it directly from the matrix to me, though.

As an example, let $\mathbf u=(1,0,-3)^T$, $\mathbf v=(1,-1,-1)^T$ and $\theta = \pi/6$. The rotation submatrix is easily obtained via Rodrigues’ formula and other methods. After multiplying everything out, the resulting matrix is $$M=\begin{bmatrix}
{1+\sqrt3\over3}&{-1+\sqrt3\over3}&-\frac13&{-1-\sqrt3\over3} \\
-\frac13&{1+\sqrt3\over3}&{1-\sqrt3\over3}&{4-3\sqrt3\over3} \\
{-1+\sqrt3\over3}&\frac13&{1+\sqrt3\over3}&{-5+2\sqrt3\over3} \\
0&0&0&1\end{bmatrix}.$$ $\operatorname{tr}M=2+\sqrt3$, so $\cos\theta=\frac{\sqrt3}2$ and so $\theta=\frac\pi6$. Row-reducing $I-M$ produces $$\begin{bmatrix}1&0&1&2\\0&1&-1&-3\\0&0&0&0\\0&0&0&0\end{bmatrix}$$ from which we can read that the rotation axis is $t(-1,1,1)^T+(-2,3,0)^T=(-(t+2),t+3,t)^T$. Setting $t=-3$ gives us the original vector $\mathbf u$. Similarly, row-reducing $\left[\begin{array}{c|c}I-R&\mathbf w\end{array}\right]$ produces $$\left[\begin{array}{ccc|c}1&0&1&-2\\0&1&-1&3\\0&0&0&0\end{array}\right]$$ from which we can see that $\mathbf u=(-(z+2),z+3,z)^T$ as before.

$\begingroup$Thanks a lot. It helped me to understand some things. I played a bit with Scilab, and I'm able to reproduce your example. However, I tried with randomly generated matrices (still a rotation and a translation). I highlight the fact that I only have one post-rotation (@Nominal Animal's definition). My input matrix is not build as a composition of a translation, then rotation, then translation back. With this kind of matrix, is your method still valid? Is there a rotation around a line matching any post-rotation matrices?$\endgroup$
– AlexJul 26 '16 at 9:32

$\begingroup$@Alex The first part of my answer shows that for a matrix in the form (*) to represent a rotation about some axis, the equation $(I-R)\mathbf u=\mathbf w$ must have a solution. For randomly-generated matrices, I would expect that most would fail this test. Of course, you also have to check that $R$ is actually a rotation, since the first equation above holds for any submatrix $R$.$\endgroup$
– amdJul 26 '16 at 16:15

$\begingroup$Here is a Scilab script for your solution: link. I made two tests: one with a matrix built in the form (*), and the other with a randomly-generated rotation and a random translation. The first test always works, but the second tests always fails. It provide the rotation axis $\mathbf{u}$, but not the vector $\mathbf{v}$ for the line $\mathbf{u}+\mathbf{v}t$.$\endgroup$
– AlexJul 27 '16 at 6:54

$\begingroup$@Alex It’s not at all surprising that a randomly-generated matrix doesn’t represent a rotation about some axis. Not every matrix of the form (*) does. $(I-R_{\mathbf v})\mathbf u$ is orthogonal to $\mathbf v$ (think about it). So $\mathbf v\cdot\mathbf w=0$ is another necessary condition for this matrix to represent a rotation. It’s pretty unlikely that a random vector is going to be orthogonal to the axis of the (also-random) rotation submatrix.$\endgroup$
– amdJul 27 '16 at 7:35

$\begingroup$Thanks to your last comment, I was able to understand why it's not possible. Even if the condition $\mathbf{v}\cdot\mathbf{w}$ is not sufficient, it proves that the general case is not possible. I tried the script with a randomly-generated translation orthogonal to $\mathbf{v}$ and it works well. Thanks for your help$\endgroup$
– AlexJul 27 '16 at 8:53

Note that if the upper left 3×3 part contains a pure rotation (no shearing, stretching, or skewing), it is orthogonal, and therefore its inverse is its transpose.

The translation is done after the rotation; I call this post-rotation translation. If you have a pre-rotation translation -- say, you want to move the origin, i.e. the point we rotate around, so you pre-translate by the negated vector to the centerpoint --, you simply apply the rotation matrix to the pre-rotation translation vector $\vec{t}$, i.e.
$$\vec{T} = \left [ \begin{matrix} T_x \\ T_y \\ T_z \\ 1 \end{matrix} \right ] = \left [ \begin{matrix} X_x t_x + Y_x t_y + Z_x t_z \\ X_y t_x + Y_y t_y + Z_y t_z \\ X_z t_x + Y_z t_y + Z_z t_z \\ 1 \end{matrix} \right ]$$
If $\vec{t}$ is the center of the rotation, and you want to apply a pre-translation by $-\vec{t}$ and a post-translation by $\vec{t}$, to keep the center of rotation at $\vec{t}$, the combined post-translation is
$$\vec{T} = \left [ \begin{matrix} T_x \\ T_y \\ T_z \\ 1 \end{matrix} \right ] = \left [ \begin{matrix} t_x - X_x t_x - Y_x t_y - Z_x t_z \\ t_y - X_y t_x - Y_y t_y - Z_y t_z \\ X_z t_x - Y_z t_y - Z_z t_z \\ 1 \end{matrix} \right ]$$
When constructing transformations, you work out the rotation part first. You can pick any point as the center of rotation, or even pick one point in the coordinate system before the rotation and another in after the rotation. After you have the rotation worked out, you can work out the post-rotation translation (as a sum of post-rotation translation(s) and rotated pre-rotation translations).

Can I write M as a rotation around an arbitrary axis, not necessarily through the origin?

Yes. First, find the direction of the axis, and save it as unit vector ${u}$. With this, you can construct the rotation part of $\mathbf{M}$.

After you have the rotation part, you can trivially add the translations to the post-translation part (rightmost column) of $\mathbf{M}$) using the method I showed above: substract the rotated center vector from the center vector.

How to find one point $\vec{a}$ on the axis?

If $\mathbf{M}$ describes a rotation (around any point or axis), then indeed
$$\vec{c} = \mathbf{M} \vec{c}$$
is true for all points $\vec{c}$ on the axis.

This is a system of three linear equations in three variables:
$$\begin{cases}
c_x = X_x c_x + Y_x c_y + Z_x c_z + T_x \\
c_y = X_y c_x + Y_y c_y + Z_y c_z + T_y \\
c_z = X_z c_x + Y_z c_y + Z_z c_z + T_z
\end{cases}$$
If there is a solution (there is no solution if there is no rotation but there is some translation), then there are an infinite number of solutions, because the axis is infinite. (If $\mathbf{M}$ is an identity matrix, then there is no rotation nor translation, and all points $\vec{c}$ fulfill the three equations.)

If you choose $c_z = 0$, then you'll find a solution only if the axis intersects the $z=0$ plane.

3.

I don't really follow, but I think I get your idea.

Let's say you have two random points, $\vec{p}$ and $\vec{q}$. They are transformed to $\vec{p}' = \mathbf{M}\vec{p}$ and $\vec{q}' = \mathbf{M}\vec{q}$.

Line through $\vec{p}$ and $\vec{p}'$ is parallel to a plane that is perpendicular to the axis of rotation. Similarly, line through $\vec{q}$ and $\vec{q}'$ is parallel to a plane that is perpendicular to the axis of rotation. (Consider $\vec{p}$ and $\vec{q}$ as points on two separate wheels, maybe of different radii, rotating on the same axis. The two wheels' planes are parallel, and perpendicular to the axis.)

If the two vectors $(\vec{p}'-\vec{p})$ and $(\vec{q}'-\vec{q})$ are not collinear, then you should find the direction of the axis using their cross product,
$$\vec{a} = (\vec{p}' - \vec{p}) \times (\vec{q}' - \vec{q})$$
Note that $\vec{a}$ is not an unit vector, and that it is only parallel to the rotation axis. You'll usually want to normalize it:
$$\hat{u} = \frac{\vec{a}}{\lVert\vec{a}\rVert}$$
where
$$\lVert\vec{a}\rVert = \sqrt{a_x^2 + a_y^2 + a_z^2}$$
as usual.

However, there is a much better method instead, one that reliably gives you the direction of the axis $\hat{u}$ (i.e. an unit vector parallel to the axis) and the angle $\theta$: convert the rotation part of the matrix to a versor (unit quaternion).

If you use the method described at the Wikipedia Rotation matrix page, you can obtain $w$ (the scalar part) and $x, y, z$ (the vector part) for the quaternion corresponding to the 3×3 rotation matrix. Note that $w^2+x^2+y^2+z^2 = 1$ for versors/rotation quaternions.

When you do have $\hat{u}$, set the coordinate with the highest coefficient in magnitude in $\hat{u}$ to zero, and solve the other two. This finds the point where the axis intersects the coordinate plane ($x=0$, $y=0$, or $z=0$) with the largest angle to the axis.

$\begingroup$$(\vec{p}'-\vec{p})\times(\vec{q}'-\vec{q})$ produces a vector that is parallel to the rotation axis, but not the rotation axis itself that the OP wants.$\endgroup$
– amdJul 26 '16 at 6:46

$\begingroup$@amd: Right. I used "the direction of the axis" because of that, but it wasn't enough. I added an explicit note about it.$\endgroup$
– Nominal AnimalJul 26 '16 at 12:35