Imagine the following situation. You have 2 coordinate systems that are coincident at the beginning (one global and one auxiliary) and 4 points in 3D space.
The four points will always be defined such that they can describe the right hand Cartesian coordinate system.
Now let's say that the 4th point describes the origin where I need to translate my auxiliary coordinate system, and the other 3 points would define the rotation. Those 3 points define rotation such that when the auxiliary coordinate system is in its final position, the 1st point is on the positive part of x-axis, the 2nd point is on the positive part of y-axis, and the 3rd point is on the positive part of z-axis.
Finally, suppose that auxiliary CS is already translated to 4th point. Then my question is: how can I get the rotation value in quaternions (preferably quaternions) that would tell me how to rotate the auxiliary coordinate system w.r.t global CS, so the axes are coincident with the before mentioned points?
Thanks in advance.
Best Answer
Like Aretino commented, use a matrix. Then, you can trivially calculate the quaternion that corresponds to that matrix.
Let's say $\vec{p}_0$ is the origin for the rotated coordinate system, $\vec{p}_1 $ is on the positive $x$ axis in the rotated coordinate system, $\vec{p}_2$ is on the positive $y$ axis in the rotated coordinate system, and $\vec{p}_3$ is on the positive $z$ axis in the rotated coordinate system.
First, calculate the new unit axis vectors: $$\begin{aligned} \displaystyle \hat{\epsilon}_1 &= \frac{\vec{p}_1 - \vec{p}_0}{\left\lVert \vec{p}_1 - \vec{p}_0 \right\rVert} = \left [ \begin{matrix} x_1 \\ y_1 \\ z_1 \end{matrix} \right ] \\ \displaystyle \hat{\epsilon}_2 &= \frac{\vec{p}_2 - \vec{p}_0}{\left\lVert \vec{p}_2 - \vec{p}_0 \right\rVert} = \left [ \begin{matrix} x_2 \\ y_2 \\ z_2 \end{matrix} \right ] \\ \displaystyle \hat{\epsilon}_3 &= \frac{\vec{p}_3 - \vec{p}_0}{\left\lVert \vec{p}_3 - \vec{p}_0 \right\rVert} = \left [ \begin{matrix} x_3 \\ y_3 \\ z_3 \end{matrix} \right ] \\ \end{aligned} \tag{1}\label{EQ1}$$ Note that while the above ensures that $\left\lVert\hat{\epsilon}_1\right\rVert = 1$, $\left\lVert\hat{\epsilon}_2\right\rVert = 1$, and $\left\lVert\hat{\epsilon}_3\right\rVert = 1$, it does not ensure that the three unit vectors are orthogonal, i.e. that $\hat{\epsilon}_1\cdot\hat{\epsilon}_2 = 0$, $\hat{\epsilon}_1\cdot\hat{\epsilon}_3 = 0$, and $\hat{\epsilon}_2\cdot\hat{\epsilon}_3 = 0$. You can use the Gram-Schmidt process on the vectors relative to $\vec{p}_0$, before normalizaton to unit length, if the points might not be exactly on the axes. That ensures the matrix $\mathbf{R}$ does not have significant error, which is important for the next step (conversion to quaternion).
The matrix $\mathbf{R}$ that rotates to that coordinate system is $$\mathbf{R} = \left [ \begin{matrix} \hat{\epsilon}_1 & \hat{\epsilon}_2 & \hat{\epsilon}_3 \end{matrix} \right ] = \left [ \begin{matrix} x_1 & x_2 & x_3 \\ y_1 & y_2 & y_3 \\ z_1 & z_2 & z_3 \\ \end{matrix} \right ] \tag{2}\label{EQ2}$$ Side note: because $\mathbf{R}$ is a pure rotation matrix, it is orthonormal, and therefore its inverse is its transpose, i.e. $$\mathbf{R}^{-1} = \mathbf{R}^T = \left [ \begin{matrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ \end{matrix} \right ]$$ Also note that we assume matrix-vector post-multiplication, i.e. that we use $\mathbf{M}\vec{v}$ to describe the rotation of vector $\vec{v}$ by matrix $\mathbf{M}$.
We can recover the unit quaternion corresponding to the pure rotation matrix $\mathbf{R}$ using one of the schemes described in the Wikipedia Rotation matrix article.
The idea is very simple, and is based on the fact that unit quaternion $\mathbf{q} = w + x\mathbf{i} + y\mathbf{j} + z\mathbf{k}$ corresponds to rotation matrix $\mathbf{R}$, $$\mathbf{R} = \left [ \begin{matrix} w^2 + x^2 - y^2 - z^2 & 2 x y - 2 w z & 2 x z + 2 w y \\ 2 x y + 2 w z & w^2 - x^2 + y^2 - z^2 & 2 y z - 2 w x \\ 2 x z - 2 w y & 2 y z + 2 w x & w^2 - x^2 - y^2 + z^2 \\ \end{matrix} \right ] \tag{3}\label{EQ3}$$ When the matrix $\mathbf{R}$ entries are known, we only need some of the entries to fully recover the corresponding unit quaternion $\mathbf{q}$ components $w$, $x$, $y$, and $z$.
When we are using numerical values, especially in a computer program, we want to minimize the effects of rounding or limited precision. This is called numerical stability. Because of this, a function or algorithm recovering the quaternion from the rotation matrix contains branches; typically using one approach if $\operatorname{Tr}\mathbf{R} = x_1 + y_2 + z_3$ is positive, and otherwise using one of three approaches, depending on which of $\lvert x_1 \rvert$, $\lvert y_2 \rvert$, or $\lvert z_3 \rvert$ is greatest. (To simplify a bit, we are trying to avoid division by any value close to zero, both implicit and explicit. We assume all the matrix entries have roughly the same absolute magnitude of error, which means values close to zero have the largest relative error; thus dividing by any matrix entry close to zero yields large relative error -- and we want to avoid/minimize that.)
If you want a purely mathematical approach, you can construct a 4×4 matrix $\mathbf{M}$ using the components of matrix $\mathbf{R}$ (as described in the Wikipedia article). Its eigenvectors correspond to the recovered rotation quaternion $\mathbf{Q}$ (that you should normalize to unit length), and the eigenvalues correspond roughly to the numerical accuracy or reliability of each approach; so picking the eigenvector corresponding to the largest eigenvalue yields the most sensible approximation for the rotation quaternion.