Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It only takes a minute to sign up.
Sign up to join this community
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.
How do I rotate the coordinate system to obtain a new $uvw$ system such that the $w$-axis is parallel to the vector $n$?
How do I obtain the coordinates of $p$ in the $uvw$-system?
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:
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:
The new coordinate of $p$ becomes $$ (p_u,p_v,p_w)=(p\cdot\hat{u},p\cdot\hat{v},p\cdot\hat{w}) $$
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.
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...
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.:)
