Given a vector $v$, a rotation axis $r$ and a plane with normal $n$ (all three vectors are unit length), how can I obtain a rotation matrix or quaternion which rotates $v$ along $r$ so it lies on the plane? This can be expressed with the following equation:
$$qvq^{-1} \cdot n = 0$$
Where $q$ is a rotation quaternion along $r$. It can be written as follows in $(vector, scalar)$ form:
$$q = (r \sin \frac{\theta}{2}, \cos \frac{\theta}{2})$$
Now we have to find $\theta$.
It is acceptable to assume $\sin \theta \approx \theta$ in my specific case.
I quickly expanded the first equation (quite cumbersome to work with so it might be wrong) and I fell on a quadratic equation which makes sense because you can have two possible angles you can rotate about to make $v$ lie on the plane, or you can also have no solution in certain configurations.
I believe solving the quadratic is a possible solution to this problem despite being rather cumbersome. I wonder if there is another more elegant and stablished solution to this problem.
Thanks
Best Answer
Assuming that the plane passes through the origin, this problem comes down to finding the intersection of the plane with the nappe of the cone defined by $\mathbf v$ and $\mathbf r$. This intersection will consist of up to two rays. Finding unit vectors in the directions of these rays is easy, from which vectors you can get the corresponding rotation angles or otherwise construct the rotation using your favorite method.
The first thing to do is to check whether or not there is a non-trivial intersection. This can be done by comparing the angles between the rotation axis $\mathbf r$, and $\mathbf v$ and the plane represented by $\mathbf n$, respectively, or, more conveniently, their cosines. Since we’re working with unit vectors, there will be an intersection when $$(\mathbf v\cdot\mathbf r)^2\le1-(\mathbf n\cdot\mathbf r)^2.$$
You can solve the cone-plane intersection problem directly, but I find it easier to transform everything so that the cone is aligned with the $z$-axis. A rotation will do the trick, but I think it’s faster and easier to use a Householder reflection instead. So, compute the transformation matrix $$H=I-2{\mathbf a\mathbf a^T\over\mathbf a^T\mathbf a}$$ where $\mathbf a=\mathbf r-(0,0,1)^T$. This is a reflection in the plane that bisects the angle between $\mathbf v$ and $\mathbf r$. The advantage of using this transformation instead of a rotation is that $H=H^T=H^{-1}=H^{-T}$, so the same matrix can be used to transform back to the original coordinate system and to transform the normal vector $\mathbf n$. (That’s a covariant vector, so if points are transformed as $M\mathbf p$, it transforms as $M^{-T}\mathbf n$.) In practice, you might not need to construct this matrix, since you can directly compute $H\mathbf p=\mathbf p-2{\mathbf p\cdot\mathbf a\over\mathbf a\cdot\mathbf a}\mathbf a$ and only need to map a few vectors.
Set $\mathbf v'=H\mathbf v$ and $\mathbf n'=H\mathbf n$. The vectors that we’re looking for are the intersections of the transformed plane with the circle $x'^2+y'^2=1-(z_{\mathbf v}')^2$, $z'=z_{\mathbf v}'$ (i.e., the circle that is the intersection of the transformed cone nappe with the plane $z'=z_{\mathbf v}'$). Note, by the way, that $z_{\mathbf v}'=\mathbf v\cdot\mathbf r$. Following the method described here, reduce this to a two-dimensional problem by computing the intersection of the transformed target plane with $z'=z_{\mathbf v}'$ and project onto the $x'$-$y'$ plane. For the present problem, this amounts to multiplying the third coordinate of $\mathbf n'$ by $z_{\mathbf v}'$, i.e., $\mathbf l=\operatorname{diag}(1,1,z_{\mathbf v}')\,\mathbf n'$. Proceeding as in the linked answer, find the intersection of this line with the circle $x'^2+y'^2=1-(z_{\mathbf v}')^2$ and project back onto the $z'=z_{\mathbf v}'$ plane to get unit vectors $\mathbf w_1'$ and $\mathbf w_2'$ (which might be identical).
With these vectors in hand, you can proceed in several ways to derive the corresponding rotations. The rotation angle required to take $\mathbf v$ to $\mathbf w_i=H\mathbf w_i'$ is just the negation of the angle between the projections of $\mathbf v'$ and $\mathbf w_i'$ onto the $x'$-$y'$ plane (negated because $H$ is orientation-reversing), and once you have this angle, you can use well-known formulas to construct an appropriate rotation matrix or quaternion. You could also construct a rotation matrix $R'$ in the transformed frame, which will just be a rotation about the $z'$-axis, and then transform it back to the original coordinate system: $R=HR'H$. (This is another instance in which the identity $H=H^T=H^{-1}=H^{-T}$ comes in handy.)
For numerical stability, you might want to align the rotation axis with one of the other coordinate axes. The best choice in that regard is the standard basis vector $\mathbf e_i$ which maximizes $\|\mathbf r-\mathbf e_i\|$. Getting back to the above solution is a simple matter of permuting coordinates after applying $H$, but be sure you keep track of orientation changes if you’re going to use the computed rotation angles directly.