Decomposition to rotation around arbitrary axis

Question

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).

1. 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.
2. 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?
3. 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?

Thanks for helping

Answer 1

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.

Answer 2

First, a lot of background to get everybody up to speed.

It is typical to use 4-vectors (implicit, because the last component is always 1, and not usually stored) to describe 3D vectors, $$\vec{p} = \left [ \begin{matrix} x \\ y \\ z \\ 1 \end{matrix} \right ]$$ and corresponding 4×4 transformation matrices: $$\mathbf{M} = \left [ \begin{matrix} X_x & Y_x & Z_x & T_x \\ X_y & Y_y & Z_y & T_y \\ X_z & Y_z & Z_z & T_z \\ 0 & 0 & 0 & 1 \end{matrix} \right ]$$ To transform a vector, we pre-multiply on the left: $$\vec{p}' = \mathbf{M}\vec{p}$$ I named the entries in the transformation matrix this way, because the coordinate axis unit vectors are $$\hat{e}_x = \left [ \begin{matrix} 1 \\ 0 \\ 0 \\ 1 \end{matrix} \right ], \; \hat{e}_y = \left [ \begin{matrix} 0 \\ 1 \\ 0 \\ 1 \end{matrix} \right ], \; \hat{e}_z = \left [ \begin{matrix} 0 \\ 0 \\ 1 \\ 1 \end{matrix} \right ] $$ and transforming them yields $$\mathbf{M}\hat{e}_x = \left [ \begin{matrix} X_x+T_x \\ X_y+T_y \\ X_z+T_z \\ 1 \end{matrix} \right ], \; \mathbf{M}\hat{e}_y = \left [ \begin{matrix} Y_x+T_x \\ Y_y+T_y \\ Y_z+T_z \\ 1 \end{matrix} \right ], \; \mathbf{M}\hat{e}_z = \left [ \begin{matrix} Z_x+T_x \\ Z_y+T_y \\ Z_z+T_z \\ 1 \end{matrix} \right ]$$ and In other words, the three first columns (with fourth component changed to $1$) define the $x$, $y$, and $z$ axes in the rotated coordinate system, with the fourth column being the translation vector.

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).

Wikipedia article on Rotation matrix contains an example how to construct the rotation part from an unit axis vector $\hat{u}$ (unit meaning its length is $1 = u_x^2+u_y^2+u_z^2$) and rotation angle $\theta$. I've used shorthand $s = \sin\theta$, $c = \cos\theta$, and $k = 1 - \cos\theta$: $$\mathbf{M} = \left [ \begin{matrix} u_x u_x k + c \; \; & u_x u_y k - u_z s & u_x u_z k + u_y s & 0 \\ u_y u_x k + u_z s & u_y u_y k + c \; \; & u_y u_z k - u_x s & 0 \\ u_z u_x k - u_y s & u_z u_y k + u_x s & u_z u_z k + c \; \; & 0 \\ 0 & 0 & 0 & 1 \end{matrix}\right ]$$

---

1. 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.

2. 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.

The trick is that a versor can also be described as $$Q = \begin{cases} w = \cos\left(\frac{\theta}{2}\right) \\ x = u_x \sin\left(\frac{\theta}{2}\right) \\ y = u_y \sin\left(\frac{\theta}{2}\right) \\ z = u_z \sin\left(\frac{\theta}{2}\right)\end{cases}$$

Note that if there is no rotation, you should get $w=1$ and $x=y=z=0$. If so, there is no rotation axis, and no reason to proceed. (If you do proceed, you'll end up dividing by zero.)

In other words, $$\hat{u} = \left ( \frac{x^2}{x^2 + y^2 + z^2}, \; \frac{y^2}{x^2 + y^2 + z^2}, \; \frac{z^2}{x^2 + y^2 + z^2} \right )$$ and $$\theta = 2 \operatorname{acos}(w)$$

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.