Поверните трехмерную систему координат так, чтобы ось Z была параллельна заданному вектору.

На основе ответа @Muphrid и с использованием формулы, приведенной по адресу:

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

Я резюмирую шаги для дальнейшего использования:

Позволять$\hat{n} = n/\|n\|$,$\theta=\cos^{-1}(\hat{k}\cdot \hat{n})$,$b=\hat{k}\times \hat{n}$, и$\hat{b}=(\hat{b}_x,\hat{b}_y,\hat{b}_z)=b/\|b\|$. Затем определите:

• $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$,

и

$$Q=\begin{bmatrix} q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2-q_0q_3) & 2(q_1q_3+q_0q_2)\\ 2(q_2q_1+q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3-q_0q_1\\ 2(q_3q_1-q_0q_2) & 2(q_3q_2+q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2 \end{bmatrix}.$$

Тогда мы имеем:

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

Новая координата$p$становится

$$(p_u,p_v,p_w)=(p\cdot\hat{u},p\cdot\hat{v},p\cdot\hat{w})$$

Есть два эквивалентных способа сделать это.

Позволять$\hat x, \hat y, \hat z$быть вашими базисными векторами. первый подход заключается в построении карты вращения, отображающей$\hat z \mapsto \hat z' = \hat n$. Затем, изображения$\hat x', \hat y', \hat z'$равны новым координатным базисным векторам$\hat u, \hat v, \hat w$. Вы можете извлечь компоненты из$p$используя скалярные произведения (т.е.$p_u = p \cdot \hat u$и так далее) в обычном порядке. Геометрически речь идет о вращении осей, чтобы найти новые базисные векторы.

Другой вариант — построить обратную карту, которая отображает$\hat n \mapsto \hat z$. Затем$xyz$-компоненты любого вектора изображения$uvw$-компоненты соответствующего входного вектора. Геометрически вы вращаете все входные данные, сохраняя фиксированную систему координат, поэтому вы можете извлекать новые компоненты, используя старую основу координат. Вы никогда не найдете новые базисные векторы таким образом, но они вам абсолютно не нужны. Таким образом, образ$p$при таком вращении будет$p_u \hat x + p_v \hat y + p_w \hat z$.

---

Таким образом, большой вопрос заключается в том, как вычислить карту вращения. Делать это с матрицами может быть довольно обременительно. Я предлагаю решение с использованием кватернионов. Найдите ось, перпендикулярную$\hat z$и$n$, вероятно, используя перекрестное произведение. Назовите эту ось$\hat b$. Найдите угол$\theta$между$\hat z$и$n$, вероятно, используя$\hat z \cdot \hat n = \cos \theta$. Запишите все векторы, используя чисто мнимые кватернионы, и вы сможете создать карту вращения.

$$R(a) = e^{\hat b \theta/2} a e^{-\hat b \theta/2}$$

Если вы нашли$\hat b$используя$\hat z \times \hat n$, то это соответствует первому описанному мною случаю, вращению$\hat z \mapsto \hat n$. Если вы выбрали$-\hat b$вместо этого это вращается$\hat n \mapsto \hat z$, что соответствует второму случаю.

Кватернионы (или их старшие братья, роторы в алгебре Клиффорда) очень полезны для вычисления вращений отдельных векторов или вращений, которые не соответствуют использованию вращений только относительно осей координат.

Спасибо, что разъяснили ваш предыдущий ответ в ответ на мой комментарий, @HåkonHægland. Я очень внимательно изучил его, потому что хочу/намереваюсь закодировать его в своей программе C "stringart" (технически, "symmography"), которая делает такие вещи. .

Таким образом, код проективной геометрии уже вроде как работает, но является случайным/уродливым/и т. д., и я намеревался разработать небольшую библиотеку для этих манипуляций и попытаться написать ее как можно элегантнее. Итак, я хочу сначала полностью и правильно понять подход (кватернионы, углы Эйлера, что у вас есть) и соответствующую математику.

Изменить. Я собрал первую часть этой маленькой библиотеки на основе ваших алгоритмов выше (и на странице ответов. Google, на которую вы ссылаетесь), которая находится под gpl'ed и которую вы можете получить с http: //www.forkosh .com/quatrot.c

Имеет встроенный тестовый драйвер

cc -DQRTESTDRIVE quatrot.c -lm -o quatrot

и, похоже, работает именно так, как вы рекламировали. Блок комментариев в тестовом драйвере main() содержит простые инструкции по использованию. Я еще не на 100% доволен функциональной декомпозицией, т.е. она должна быть одновременно простой в использовании для программиста и интуитивно понятной для математика. Но главное, что это работает (насколько я проверял, кажется). Всего 460 строк, поэтому я не буду приводить все это здесь (во всяком случае, проще взять по ссылке выше). Две функции, в основном основанные на answer.google, это...

/* ==========================================================================
 * 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() --- */

Это тестовый драйвер main(), который содержит код, основанный на вашем дополнительном обсуждении.$x,y,z\to u,v,w$, а затем проецирование заданной точки$p$чтобы получить компоненты по новым осям (фрагмент)...

  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) --- */

Фрагмент выше может быть немного труден для чтения вне контекста (полный контекст см. в функции main()), но там есть такие вещи, как...

  /* --- 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?*/

что в значительной степени дословно, так сказать, то, что вы написали выше. По крайней мере, я на это надеюсь. :)