Gire el sistema de coordenadas 3D de modo que el eje z sea paralelo a un vector dado

Basado en la respuesta de @Muphrid y usando la fórmula dada en:

http://answers.google.com/answers/threadview/id/361441.html

Resumo los pasos para referencia posterior:

Dejar$\hat{n} = n/\|n\|$,$\theta=\cos^{-1}(\hat{k}\cdot \hat{n})$,$b=\hat{k}\times \hat{n}$, y$\hat{b}=(\hat{b}_x,\hat{b}_y,\hat{b}_z)=b/\|b\|$. Luego define:

• $q_0=\cos(\theta/2)$,
• $q_1=\sin(\theta/2)\hat{b}_x$,
• $q_2=\sin(\theta/2)\hat{b}_y$,
• $q_3=\sin(\theta/2)\hat{b}_z$,

y

Entonces tenemos:

• $\hat{u}=Q\hat{i}$,
• $\hat{v}=Q\hat{j}$,
• $\hat{w}=Q\hat{k}=\hat{n}$.

La nueva coordenada de$p$se convierte

Hay dos formas equivalentes de hacer esto.

Dejar$\hat x, \hat y, \hat z$ser sus vectores base. el primer enfoque es construir un mapa de rotación que mapee$\hat z \mapsto \hat z' = \hat n$. Luego, las imágenes$\hat x', \hat y', \hat z'$son iguales a los nuevos vectores de base de coordenadas$\hat u, \hat v, \hat w$. Puede extraer los componentes de$p$usando productos escalares (es decir,$p_u = p \cdot \hat u$, etc.) de la forma habitual. Geométricamente, se trata de rotar los ejes para encontrar los nuevos vectores base.

La otra opción es construir el mapa inverso, que mapea$\hat n \mapsto \hat z$. Entonces el$xyz$-los componentes de cualquier vector de imagen son los$uvw$-componentes del vector de entrada correspondiente. Geométricamente, está girando todas las entradas mientras mantiene fijo el sistema de coordenadas, por lo que puede extraer nuevos componentes mientras usa la base de coordenadas anterior. Nunca encuentras los nuevos vectores base de esta manera, pero no los necesitas absolutamente. Así, la imagen de$p$bajo esta rotación sería$p_u \hat x + p_v \hat y + p_w \hat z$.

La gran pregunta, entonces, es cómo calcular el mapa de rotación. Hacer esto con matrices podría ser bastante oneroso. Sugiero una solución usando cuaterniones. Encuentre el eje perpendicular a$\hat z$y$n$, probablemente usando el producto cruz. Llame a este eje$\hat b$. Encuentra el ángulo$\theta$entre$\hat z$y$n$, probablemente usando$\hat z \cdot \hat n = \cos \theta$. Escriba todos los vectores usando cuaterniones imaginarios puros y podrá crear un mapa de rotación

Si usted encontró$\hat b$mediante el uso$\hat z \times \hat n$, entonces esto corresponde al primer caso que describí, rotando$\hat z \mapsto \hat n$. si eliges$-\hat b$en cambio, entonces esto gira$\hat n \mapsto \hat z$, y corresponde al segundo caso.

Los cuaterniones (o sus hermanos mayores, los rotores en un álgebra de Clifford) son muy útiles para calcular rotaciones de vectores únicos o rotaciones que no se ajustan al uso de rotaciones solo con respecto a los ejes de coordenadas.

Gracias por aclarar su respuesta anterior en respuesta a mi comentario, @HåkonHægland Lo estaba estudiando muy de cerca porque quiero/tengo la intención de codificarlo en mi programa C "stringart" (técnicamente, "simografía") que hace cosas como esta ... .

Así que el código de geometría proyectiva ya funciona, pero es ad-hoc/feo/etc, y he tenido la intención de desarrollar una pequeña biblioteca para estas manipulaciones, y tratar de escribirlo tan elegantemente como pueda. Así que primero quiero obtener el enfoque (cuaterniones, ángulos de Euler, lo que sea) y las matemáticas correspondientes de forma completa y correcta.

Editar He improvisado un primer corte de esa pequeña biblioteca, basado en sus algoritmos anteriores (y en la página de respuestas.google que vinculó), que está gpl'ed, y que puede obtener de http: //www.forkosh .com/quatrot.c

Tiene un controlador de prueba incorporado.

cc -DQRTESTDRIVE quatrot.c -lm -o quatrot

y parece estar funcionando exactamente como lo anunciaste. El bloque de comentarios en el controlador de prueba main() tiene instrucciones de uso simples. Todavía no estoy 100% satisfecho con la descomposición funcional, es decir, debería ser al mismo tiempo fácil de usar para el programador e intuitivamente obvio para el matemático. Pero lo principal es que funciona (parece que lo he probado). Todo tiene 460 líneas, por lo que no lo reproduciré todo aquí (de todos modos, es más fácil obtenerlo del enlace de arriba). Dos de las funciones, principalmente basadas en answers.google son...

/* ==========================================================================
 * Function:    qrotate ( LINE axis, double theta )
 * Purpose:     returns quaternion corresponding to rotation
 *              through theta (**in radians**) around axis
 * --------------------------------------------------------------------------
 * Arguments:   axis (I)        LINE axis around which rotation by theta
 *                              is to occur
 *              theta (I)       double theta rotation **in radians**
 * --------------------------------------------------------------------------
 * Returns:     ( QUAT )        quaternion corresponding to rotation
 *                              through theta around axis
 * --------------------------------------------------------------------------
 * Notes:     o Rotation direction determined by right-hand screw rule
 *              (when subsequent qmultiply() is called with istranspose=0)
 * ======================================================================= */
/* --- entry point --- */
QUAT    qrotate ( LINE axis, double theta )
{
/* --- allocations and declarations --- */
QUAT    q = { cos(0.5*theta), 0.,0.,0. } ; /* constructed quaternion */
double  x = axis.pt2.x - axis.pt1.x,    /* length of x-component of axis */
        y = axis.pt2.y - axis.pt1.y,    /* length of y-component of axis */
        z = axis.pt2.z - axis.pt1.z;    /* length of z-component of axis */
double  r = sqrt((x*x)+(y*y)+(z*z));    /* length of axis */
double  qsin = sin(0.5*theta);          /* for q1,q2,q3 components */
/* --- construct quaternion and return it to caller */
if ( r >= 0.0000001 ) {                 /* error check */
  q.q1 = qsin*x/r;                      /* x-component */
  q.q2 = qsin*y/r;                      /* y-component */
  q.q3 = qsin*z/r; }                    /* z-component */
return ( q );
} /* --- end-of-function qrotate() --- */

/* ==========================================================================
 * Function:    qmatrix ( QUAT q )
 * Purpose:     returns 3x3 rotation matrix corresponding to quaternion q
 * --------------------------------------------------------------------------
 * Arguments:   q (I)           QUAT q for which a rotation matrix
 *                              is to be constructed
 * --------------------------------------------------------------------------
 * Returns:     ( double * )    3x3 rotation matrix, stored row-wise
 * --------------------------------------------------------------------------
 * Notes:     o The matrix constructed from input q = q0+q1*i+q2*j+q3*k is:
 *                    (q0²+q1²-q2²-q3²)    2(q1q2-q0q3)     2(q1q3+q0q2)
 *                Q  =   2(q2q1+q0q3)   (q0²-q1²+q2²-q3²)   2(q2q3-q0q1)
 *                       2(q3q1-q0q2)      2(q3q2+q0q1)  (q0²-q1²-q2²+q3²)
 *            o The returned matrix is stored row-wise, i.e., explicitly
 *                --------- first row ----------
 *                qmatrix[0] = (q0²+q1²-q2²-q3²)
 *                       [1] =    2(q1q2-q0q3)
 *                       [2] =    2(q1q3+q0q2)
 *                --------- second row ---------
 *                       [3] =    2(q2q1+q0q3)
 *                       [4] = (q0²-q1²+q2²-q3²)
 *                       [5] =    2(q2q3-q0q1)
 *                --------- third row ----------
 *                       [6] =    2(q3q1-q0q2)
 *                       [7] =    2(q3q2+q0q1)
 *                       [8] = (q0²-q1²-q2²+q3²)
 *            o qmatrix maintains a static buffer of 64 3x3 matrices
 *              returned to the caller one at a time. So you may issue
 *              64 qmatrix() calls and continue using all returned matrices.
 *              But the 65th call re-uses the memory used by the 1st call, etc.
 * ======================================================================= */
/* --- entry point --- */
double  *qmatrix ( QUAT q )
{
/* --- allocations and declarations --- */
static  double Qbuff[64][9];            /* up to 64 calls before wrap-around */
static  int    ibuff = (-1);            /* Qbuff[ibuff][] index 0<=ibuff<=63 */
double  *Q = NULL;                      /* returned ptr Q=Qbuff[ibuff] */
double  q0=q.q0, q1=q.q1, q2=q.q2, q3=q.q3; /* input quaternion components */
double  q02=q0*q0, q12=q1*q1, q22=q2*q2, q32=q3*q3; /* components squared */
/* --- first maintain Qbuff[ibuff][] buffer --- */
if ( ++ibuff > 63 ) ibuff=0;            /* wrap Qbuff[ibuff][] index */
Q = Qbuff[ibuff];                       /* ptr to current 3x3 buffer */
/* --- just do the algebra described in the above comments --- */
Q[0] = (q02+q12-q22-q32);
Q[1] =  2.*(q1*q2-q0*q3);
Q[2] =  2.*(q1*q3+q0*q2);
Q[3] =  2.*(q2*q1+q0*q3);
Q[4] = (q02-q12+q22-q32);
Q[5] =  2.*(q2*q3-q0*q1);
Q[6] =  2.*(q3*q1-q0*q2);
Q[7] =  2.*(q3*q2+q0*q1);
Q[8] = (q02-q12-q22+q32);
/* --- return constructed quaternion to caller */
return ( Q );
} /* --- end-of-function qmatrix() --- */

Es el controlador de prueba principal () que contiene el código basado en su discusión adicional rotando$x,y,z\to u,v,w$, y luego proyectando el punto dado$p$para obtener componentes a lo largo de los nuevos ejes (fragmento)...

  int   ipt=0, npts=4;                  /* testpoints[] index, #test points */
  POINT testpoints[] = {                /* test points for quatrot funcs */
        {000.,000.,000.}, {100.,000.,000.}, {000.,100.,000.}, {000.,000.,100.}
        } ;
  /* --- test data for simple rotations around given test axis --- */
  QUAT  q = qrotate(axis,pi*theta/180.); /* rotation quaternion */
  double *Q = qmatrix(q);               /* quat rotation matrix */
  /* --- test data to rotate z-axis to given axis and project points --- */
  double thetatoz = acos(dotprod(khat,unitvec(axis))); /* rotate axis to z */
  LINE  bhat = unitvec(crossprod(khat,unitvec(axis))); /* as per stackexch */
  QUAT  qtoz = qrotate(bhat,thetatoz);  /* quat to rotate z to axis */
  double *Qtoz = qmatrix(qtoz);         /* quat matrix to rotate z to axis */
  LINE  uhat = { {0.,0.,0.}, qmultiply(Qtoz,ihat.pt2,0) }, /*new x-axis*/
        vhat = { {0.,0.,0.}, qmultiply(Qtoz,jhat.pt2,0) }, /*new y-axis*/
        what = { {0.,0.,0.}, qmultiply(Qtoz,khat.pt2,0) }; /*z=testaxis?*/
  /* --- apply rotations/projections to testpoints[] --- */
  fprintf(msgfp,"      theta: %.3f\n",theta);
  prtline ("       axis:",&axis);
  for ( ipt=0; ipt<npts; ipt++ ) {      /* rotate each test point */
    /* --- select test point --- */
    POINT pt = testpoints[ipt];         /* current test point */
    /* --- rotate test point around axis in existing x,y,z coords --- */
    POINT ptrot = qmultiply(Q,pt,0);    /* rotate pt around axis by theta */
    POINT pttran= qmultiply(Q,pt,1);    /* transpose rotation */
    /* --- project test point to rotated axes, where given axis=z-axis --- */
    POINT pttoz = {dotprod(pt,uhat),dotprod(pt,vhat),dotprod(pt,what)} ;
    /* --- display test results --- */
    fprintf(msgfp,"testpoint#%d...\n",ipt+1);
    prtpoint("     testpoint: ",&pt);   /* current test point... */
    prtpoint("       rotated: ",&ptrot);  /* ...rotated around given axis */
    prtpoint(" transpose rot: ",&pttran); /* ...transpose rotation */
    prtpoint("   axis=z-axis: ",&pttoz);  /* ...coords when axis=z-axis */
    } /* --- end-of-for(ipt) --- */

Ese fragmento de arriba puede ser un poco difícil de leer fuera de contexto (vea la función main() para el contexto completo), pero son cosas como...

  /* --- test data to rotate z-axis to given axis and project points --- */
  double thetatoz = acos(dotprod(khat,unitvec(axis))); /* rotate axis to z */
  LINE  bhat = unitvec(crossprod(khat,unitvec(axis))); /* as per stackexch */
  QUAT  qtoz = qrotate(bhat,thetatoz);  /* quat to rotate z to axis */
  double *Qtoz = qmatrix(qtoz);         /* quat matrix to rotate z to axis */
  LINE  uhat = { {0.,0.,0.}, qmultiply(Qtoz,ihat.pt2,0) }, /*new x-axis*/
        vhat = { {0.,0.,0.}, qmultiply(Qtoz,jhat.pt2,0) }, /*new y-axis*/
        what = { {0.,0.,0.}, qmultiply(Qtoz,khat.pt2,0) }; /*z=testaxis?*/

que es más o menos palabra por palabra, por así decirlo, lo que escribiste anteriormente. Al menos espero que lo sea :)