[Publicado como una respuesta separada ya que mi anterior ya ha sido aceptada].
A veces me resulta más fácil generar varias soluciones potenciales y luego reducirlas en lugar de intentar calcular directamente la solución correcta. Las pruebas que se realizan para descartar posibles soluciones pueden ser más sencillas que el análisis de casos que se requiere para acotar las cosas desde el principio. Este es el enfoque que adoptaré para este problema.
Asumo que ya has determinado que hay que hacer un ajuste en el eje de rotación, así que sólo me ocuparé de encontrar esa rotación ajustada. La tarea consiste en encontrar un eje (y un ángulo) de rotación tal que el vector de rotación sólo roce el cono azul de bloqueo y el ángulo de rotación sea mínimo entre todos esos ejes. De forma equivalente, buscamos un plano de rotación tal que el arco circular trazado en ese plano por el vector de rotación intersecte el cono exactamente en un punto. Para que esto ocurra, el vector debe ser paralelo a la generatriz del cono en ese punto. WLOG, que el vértice del cono esté en el origen y su eje el $z$ -eje. Una condición equivalente a la anterior es que el plano de rotación sea tangente a la intersección del cono con uno de los planos $z=\pm\|\vec v_1\|\cos\alpha_B$ un círculo de radio $\|\vec v_1\|\sin\alpha_B$ . Si consideramos sólo el plano de la circunferencia, el problema se reduce al bien estudiado de encontrar las intersecciones de las tangentes a esta circunferencia a través del punto común de intersección $p$ de los planos de rotación y el plano del círculo. Una forma sencilla de calcular estos puntos es encontrar la intersección del círculo con $p$ de la línea polar. Al elevar estas intersecciones a 3D se obtienen hasta cuatro puntos de roce candidatos, de los cuales se puede seleccionar el que da lugar a la mínima rotación.
Estos problemas de intersección pueden resolverse de varias maneras. Voy a describir un procedimiento que calcula los puntos directamente. Este procedimiento salta de un lado a otro entre espacios bidimensionales y tridimensionales y también cambia entre coordenadas cartesianas y homogéneas, así que un poco de notación para ayudar a mantener las cosas claras: letras minúsculas en negrita como $\mathbf p$ denotan vectores columna de coordenadas homogéneas en $\mathbb{RP}^2$ y las mayúsculas $\mathbf P$ denota lo mismo para $\mathbb{RP}^3$ . En ambos casos, una tilde ( $\tilde{\mathbf p}$ , $\tilde{\mathbf P}$ ) indica coordenadas no homogéneas (cartesianas). También me permitiré un pequeño abuso de notación al mostrar dos tuplas homogéneas como iguales cuando representan el mismo objeto, cuando en realidad sólo son equivalentes. Para reducir un poco el desorden, dejemos que $v=\|\vec v_1\|=\|\vec v_2\|$ .
Tenemos $\tilde{\mathbf V}_1=[x_1,y_1,z_1]^T$ y $\tilde{\mathbf V}_2=[x_2,y_2,z_2]^T$ . El plano del círculo superior es $$\mathbf\Pi_+=[0,0,1,-v\cos\alpha_B]^T.$$ La intersección de este plano con la línea que pasa por $\mathbf V_1$ y $\mathbf V_2$ se encuentra fácilmente a través de la matriz de Plücker $\mathscr L(\mathbf V_1,\mathbf V_2)=\mathbf V_1\mathbf V_2^T-\mathbf V_2\mathbf V_1^T$ de la línea: es simplemente $\mathbf P_+=\mathscr L(\mathbf V_1,\mathbf V_2)\mathbf\Pi_+$ . Este será un punto en el infinito si $z_1=z_2$ pero el resto del cálculo también se ocupa de este caso.
Proyecto $\mathbf P_+$ en el $x$ - $y$ plano dejando caer su tercera coordenada. Más concretamente, calcula $\mathbf p_+=\mathscr H\mathbf P_+$ , donde $$\mathscr H=\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&0&1\end{bmatrix}.$$ The circle has matrix $C=\operatorname{diag}(1,1,-v^2\sin^2\alpha_B)$, and the polar line to $\mathbf p_+$ is $\mathbf l_+=C\mathbf p_+$. Putting this cascade all together and simplifying produces the point $$\mathbf p_+=\mathscr H\mathscr L(\mathbf V_1,\mathbf V_2)\mathbf\Pi_+ = \begin{bmatrix} x_1z_2-x_2z_1-(x_1-x_2)v\cos\alpha_B \\ y_1z_2-y_2z_1-(y_1-y_2)v\cos\alpha_B \\z_2-z_1\end{bmatrix}$$ y la línea $$\mathbf l_+ = C\mathbf p_+ = \begin{bmatrix} x_1z_2-z_1x_2-(x_1-x_2)v\cos\alpha_B \\ y_1z_2-z_1y_2-(y_1-y_2)v\cos\alpha_B \\ (z_1-z_2)v^2\sin^2\alpha_B \end{bmatrix}.\tag1$$ Se presenta aquí sin pruebas† una fórmula general para la(s) intersección(es) de la recta $[\lambda,\mu,\nu]^T$ con un círculo de radio $r$ centrado en el origen: $$\begin{bmatrix}-\lambda\nu\pm\mu\sqrt{r^2(\lambda^2+\mu^2)-\nu^2} \\ -\mu\nu\mp\lambda\sqrt{r^2(\lambda^2+\mu^2)-\nu^2} \\ \lambda^2+\mu^2 \end{bmatrix}.\tag2$$ Si prefieres trabajar con la pértiga $\mathbf p_+$ en cambio, la fórmula de las intersecciones de las tangentes por un punto $[x,y,w]^T$ y un círculo de radio $r$ es bastante similar: $$\begin{bmatrix} r^2xw\pm ry\sqrt{x^2+y^2-r^2w^2} \\ r^2yw\mp rx\sqrt{x^2+y^2-r^2w^2} \\ x^2+y^2 \end{bmatrix}.$$ (La similitud de estas dos fórmulas no es casual. De hecho, están relacionadas por la relación polo/polar $[\lambda,\mu,\nu]^T=[x,y,-r^2w]^T$ del punto y la línea).
Levantar cada punto tangente $\mathbf g=[x_g,y_g,w_g]^T$ volver a 3-D. Es mejor convertir a coordenadas cartesianas ahora porque el siguiente paso utilizará exclusivamente coordenadas cartesianas, por lo que esta transformación es $$[x_g,y_g,w_g]^T\mapsto\left[\frac{x_g}{w_g},\frac{y_g}{w_g},v\cos\alpha_B\right]^T.$$ (Si $w_g=0$ uno de los puntos extremos de la rotación es interior al cono, por lo que no es posible la solución). Repite el algoritmo de intersección para el plano $z=-v\cos\alpha$ . Si utilizas las fórmulas anteriores, una forma un poco furtiva de conseguirlo es darle la vuelta al signo de $v$ .
Ahora tienes un conjunto de hasta cuatro puntos de pastoreo $\tilde{\mathbf G}_i$ a partir de la cual se debe seleccionar la rotación mínima. A medida que el plano de rotación se aleja del eje ideal $\tilde{\mathbf A}_0=\tilde{\mathbf V}_1\times\tilde{\mathbf V}_2$ el ángulo de rotación aumenta, pero el radio del arco de rotación disminuye. Afortunadamente, el ángulo gana, por lo que tanto el ángulo de rotación como la longitud del arco son funciones crecientes del ángulo entre el eje de rotación y $\tilde{\mathbf A}_0$ . Queremos que la orientación del eje de rotación capture el hecho de que $\tilde{\mathbf G}_i$ también se encuentra en el arco de rotación, por lo que se establece $\tilde{\mathbf A}_i=(\tilde{\mathbf G}_i-\tilde{\mathbf V}_1)\times(\tilde{\mathbf V}_2-\tilde{\mathbf G}_i)$ normalizado. Si es cero, significa que el punto de roce es uno de los dos puntos finales de la rotación. Ya sabes que la rotación con ese punto de rozamiento no es buena -por eso estás haciendo todo esto en primer lugar- así que descarta ese punto. Efectivamente, la rotación tendrá que dar "la vuelta larga". Por último, calcula $\tilde{\mathbf A}_0^T\tilde{\mathbf A}_i$ para cada uno de los puntos supervivientes y seleccionar el que maximiza este valor.
El cálculo anterior puede adaptarse fácilmente a un cono en posición general, aunque no estoy seguro de que merezca la pena la complejidad extra: Los planos $\mathbf\Pi_+$ y $\mathbf\Pi_-$ sólo tiene que utilizar $\tilde{\mathbf V}_B$ normalizado, en lugar de $[0,0,1]^T$ pero todavía tendrá que trasladar implícitamente el vértice al origen para calcular los ejes de rotación, y ambos $\mathscr H$ y el impulso de vuelta a las 3-D implicará una rotación. Me parece mucho más fácil transformar $\mathbf V_1$ y $\mathbf V_2$ y luego transformar de nuevo una vez que haya encontrado el nuevo eje de rotación.
† Derivé estas fórmulas utilizando el procedimiento para calcular la intersección de una recta y una cónica descrito en la obra de Richter-Gebert Perspectivas de la geometría proyectiva una variante de la cual se describe aquí . El algoritmo del libro es un cálculo directo que no requiere resolver ninguna ecuación.
Por ejemplo, dejemos que $\alpha=\pi/6$ , $\vec v_1=(1,1,2)$ y $\vec v_2=(-2,-1,-1)$ . Introduciendo estos valores en la fórmula (1) se obtiene $\mathbf l_+\approx[-3.36396, -3.24264, 4.5]^T$ y la fórmula (2) da como resultado $\mathbf g_1\approx[3.675, 26.4836, 21.831]^T$ , $\mathbf g_2\approx[26.6007, 2.70018, 21.831]^T$ . Repitiendo esto para el plano inferior se obtiene $\mathbf g_3\approx[22.6048, -139.23, 115.169]^T$ y $\mathbf g_4\approx[-106.88, 92.0459, 115.169]^T$ . Convirtiendo a cartesiana y elevando a 3D se obtienen los cuatro puntos de roce $$\tilde{\mathbf G}_1\approx[0.168339, 1.21312, 2.12132]^T \\ \tilde{\mathbf G}_2\approx[1.21848, 0.123686, 2.12132]^T \\ \tilde{\mathbf G}_3\approx[0.196275, -1.20892, -2.12132]^T \\ \tilde{\mathbf G}_4\approx[-0.928031, 0.799224, -2.12132]^T.$$ Los productos de puntos de los ejes correspondientes son $3.22549$ , $0.31171$ , $-5.79138$ y $-5.38299$ respectivamente, por lo que la mejor es la primera. Los resultados de este cálculo se ilustran a continuación.