3
$\begingroup$

Assume I have a Cartesian $xyz$ coordinate system with a point $p=(p_x,p_y,p_z)$. In addition a vector $n=(n_x,n_y,n_z)$ is given.

  1. How do I rotate the coordinate system to obtain a new $uvw$ system such that the $w$-axis is parallel to the vector $n$?

  2. How do I obtain the coordinates of $p$ in the $uvw$-system?

$\endgroup$
5
$\begingroup$

There are two equivalent ways of doing this.

Let $\hat x, \hat y, \hat z$ be your basis vectors. the first approach is to construct a rotation map that maps $\hat z \mapsto \hat z' = \hat n$. Then, the images $\hat x', \hat y', \hat z'$ are equal to the new coordinate basis vectors $\hat u, \hat v, \hat w$. You can extract the components of $p$ using dot products (i.e. $p_u = p \cdot \hat u$, and so on) in the usual manner. Geometrically, this is about rotating the axes to find the new basis vectors.

The other option is to construct the inverse map, which maps $\hat n \mapsto \hat z$. Then, the $xyz$-components of any image vector are the $uvw$-components of the corresponding input vector. Geometrically, you're rotating all inputs while keeping the coordinate system fixed, so you can extract new components while using the old coordinate basis. You never find the new basis vectors this way, but you don't absolutely need them. Thus, the image of $p$ under this rotation would be $p_u \hat x + p_v \hat y + p_w \hat z$.


The big question, then, is how to compute the rotation map. Doing this with matrices could be pretty onerous. I suggest a solution using quaternions. Find the axis perpendicular to $\hat z$ and $n$, probably using the cross product. Call this axis $\hat b$. Find the angle $\theta$ between $\hat z$ and $n$, probably using $\hat z \cdot \hat n = \cos \theta$. Write all the vectors using pure imaginary quaternions, and you can create a rotation map

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

If you found $\hat b$ by using $\hat z \times \hat n$, then this corresponds to the first case I described, rotating $\hat z \mapsto \hat n$. If you chose $-\hat b$ instead, then this rotates $\hat n \mapsto \hat z$, and it corresponds to the second case.

Quaternions (or their big brothers, rotors in a clifford algebra) are very useful for calculating rotations of single vectors, or rotations that do not conform to using rotations only with respect to coordinate axes.

$\endgroup$
  • $\begingroup$ Thanks for the interesting answer. But I am going to use this in a computer program and I am not familiar with quaternions.. $\endgroup$ – Håkon Hægland Oct 28 '13 at 14:43
  • $\begingroup$ Is that so? Quaternions are heavily used in computer programming, so there is often support for using them (so you don't usually have to implement them from scratch). What language are you using? $\endgroup$ – Muphrid Oct 28 '13 at 14:45
  • $\begingroup$ I am using Matlab. $\endgroup$ – Håkon Hægland Oct 28 '13 at 14:58
  • $\begingroup$ Good, Matlab seems to have decent quaternion support, and can compute the rotation matrix that corresponds to a given quaternion. See mathworks.com/discovery/quaternion.html $\endgroup$ – Muphrid Oct 28 '13 at 20:22
  • $\begingroup$ Thanks! I will try that. What is actually $a$ in the above formula? $\endgroup$ – Håkon Hægland Oct 28 '13 at 21:59
5
$\begingroup$

Based on the answer of @Muphrid and using the formula given at:

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

I summarize the steps for later reference:

Let $\hat{n} = n/\|n\|$, $\theta=\cos^{-1}(\hat{k}\cdot \hat{n})$, $b=\hat{k}\times \hat{n}$, and $\hat{b}=(\hat{b}_x,\hat{b}_y,\hat{b}_z)=b/\|b\|$. Then 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$,

and

$$ 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}. $$

We then have:

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

The new coordinate of $p$ becomes $$ (p_u,p_v,p_w)=(p\cdot\hat{u},p\cdot\hat{v},p\cdot\hat{w}) $$

$\endgroup$
  • $\begingroup$ Yeah, I guess matlab may not have a good quaternion to matrix conversion function built in, or it may be they assume you'll just use quats to do rotations without going to a matrix representation. $\endgroup$ – Muphrid Oct 29 '13 at 1:02
  • $\begingroup$ Excuse me for commenting so long after your post, but I came across this Q&A while almost posting a similar question just now. Your $\hat k\cdot n$ argument to $\cos^{-1}$ ... that should be $\hat n$, right? A unit vector in the $\vec n$-direction (otherwise $\cos^{-1}$ might get out-of-bounds args)? Thanks. $\endgroup$ – John Forkosh Sep 18 '16 at 5:25
  • 2
    $\begingroup$ @JohnForkosh Thanks for the comment! I have updated the answer. $\endgroup$ – Håkon Hægland Sep 18 '16 at 7:44
  • $\begingroup$ Thanks for the clarification, Håkon. I've added an "answer" below that just extends this comment becuase the image can't be put in a comment... $\endgroup$ – John Forkosh Sep 18 '16 at 9:49
1
$\begingroup$

Thanks for clarifying your preceding answer in reply to my comment, @HåkonHægland I was studying it so very closely because I want/intend to code it in my "stringart" (technically, "symmography") C program that does stuff like this...

enter image description here

So the projective geometry code already kind of works, but is ad-hoc/ugly/etc, and I've been meaning to develop a little library for these manipulations, and try to write it as elegantly as I can. So I want to get the approach (quaternions, euler angles, what have you) and the corresponding math completely and correctly figured out first.

Edit I've cobbled together a first cut of that little library, based on your algorithms above (and on the answers.google page you linked), which is gpl'ed, and which you can get from http://www.forkosh.com/quatrot.c

It has an embedded test driver

cc -DQRTESTDRIVE quatrot.c -lm -o quatrot

and seems to be working exactly as you advertised. The comment block in the test driver main() has simple usage instructions. I'm not 100% happy with the functional decomposition yet, i.e., should simultaneously be easy-to-use for the programmer while intuitively obvious for the mathematician. But the main thing is it works (seems to as far as I've tested it). The whole thing is 460 lines, so I won't reproduce it all here (easier to just get it from the link above, anyway). Two of the functions, mainly based on answers.google are...

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

It's the main() test driver that contains the code based on your additional discussion rotating $x,y,z\to u,v,w$, and then projecting given point $p$ to get components along the new axes (snippet)...

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

That snippet above may be a little difficult to read out-of-context (see the main() function for full context), but it's the stuff there like...

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

which is pretty much word-for-word, so to speak, what you wrote above. At least I hope it is.:)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.