For context, I am trying to apply the same rotations from one set of 3D points to another set of 3D points. Concretely, I am trying to model the movements of human joints from a video on a 3D human body model.
From what I gathered, I need to calculate the Euler angle between two 3D points around the third 3D point (they all represent joint locations in the video), so I can apply the same rotation to my 3D model.
So, my problem boils down to that I want to calculate the Euler angles between X, Y and Z axes (in radians) between these two 3D points:
(224, 129, 76)
(225, 129, 78)
rotating around this 3D point:
(181, 151, 67)
How do I do this? I am aiming to implement this in software.
Best Answer
Denote the two points as $$ P_1:(x_1,y_1,z_1), \quad P_2: (x_2,y_2,z_2) $$ and the center as $P_0: (x_0,y_0,z_0)$. First of all, note that the rotation matrix corresponding to the rotation about $P_0$ will take the vector $v_1 = P_1 - P_0$ to the vector $v_2 = P_2 - P_0$.
Note that a suitable rotation will only exist if $\|v_1\| = \|v_2\|$, i.e. both points are equidistant from the center of the rotation. If that is the case, then it suffices to find rotation matrix $R$ satisfying $Rv_1 = v_2$. Once this is done, the desired affine rotation is the same rotation about $P_0$, that is, $$ T(v) = R(v - P_0) + P_0 = Rv + (I-R)P_0, $$ corresponding to the homogeneous coordinates transformation matrix $$ \pmatrix{R & (I-R)P_0\\0 & 1}. $$
There are infinitely many choices for a rotation (about the center) from $v_1$ to $v_2$. However, I will suppose that we want the rotation by the smallest possible angle that achieves this result. I will also assume that $v_1$ and $v_2$ are linearly independent; the edge cases of $v_1 = \pm v_2$ can be handled separately.
I claim (without proof, for now) that the rotation that achieves this is the smaller of the two rotations along the axis orthogonal to both $v_1$ and $v_2$. In particular, it is the counterclockwise rotation by angle $$ \theta = \arccos\left(\frac{v_1 \cdot v_2}{\|v_1\|\,\|v_2\|}\right) $$ about the axis $v_1 \times v_2$. Let $u = (x_u,y_u,z_u)$ denote the unit vector parallel to $v_1 \times v_2$, $u = \frac{v_1\times v_2}{\|v_1 \times v_2\|}$. With Rodrigues' rotation formula, we find that this rotation has the matrix $R = I + \sin\theta K + (1-\cos \theta) K^2$, where $K$ denotes the "cross-product matrix" for the vector $u$, $$ K = \pmatrix{ 0 & -z_u & y_u\\ z_u&0&x_u\\ -y_u&x_u&0}. $$ From there, the $z$-$x$-$Z$ extrinsic Euler angles for the rotation can be obtained using the formulas given here: \begin{align} \phi &= \operatorname{arctan2}\left(R_{31}, R_{32}\right),\\ \theta &= \arccos\left(R_{33}\right),\\ \psi &= -\operatorname{arctan2}\left(R_{13}, R_{23}\right). \end{align}
As it turns out, we can neatly rewrite the formula for $R$ using the fact that $u = \frac 1{\sin\theta}v_1 \times v_2$. I claim (again, without proof for now) that if $v_1,v_2$ are normalized to have length $1$, then $$ R = (v_1^\top v_2) I + (v_2v_1^\top - v_1v_2^\top) + (1 + v_1^\top v_2)(v_1 \times v_2)(v_1 \times v_2)^\top. $$ Another nice reformulation: $$ M = \pmatrix{v_1 & v_2 & v_1\times v_2}, \\ R = (v_1^\top v_2) I + M\pmatrix{0&-1&0\\1&0&0\\0&0&1 + v_1^\top v_2}M^\top. $$