Rotations about an arbitrary axis

[Dan Christensen]

In article <SPENCER.90May22144823 at spline.eecs.umich.edu> spencer at eecs.umich.edu (Spencer W. Thomas) writes:
>In article <9005171647.AA14106 at baby.swmed.utexas.edu> rose at BABY.SWMED.UTEXAS.EDU (Rose Oguz) writes:
>[Problem of rotation about an arbitrary axis.]
>> Basically, I have two vectors coming
>> from the same point and one of the vectors needs to be rotated into the 
>> other.
>
>	You want to rotate the vector V1 into the vector V2.  We can
>do this in two steps: (1) rotate V1 to the X axis and then (2) rotate
>the X axis to V2.  (Why do this, you say?  Because each step is easy.)

But there is more than one rotation mapping a vector V1 to a vector V2.
The final result can be oriented in any way about V2.  The original
poster indicated that he wanted the one that rotates about V1xV2,
ie. the most "direct" rotation.  I don't think that your solution 
does this. 

For example, suppose in a right handed coordinate system the poster
wants to rotate the y axis to the z axis (by rotating about y cross z,
ie. the x axis).  This will map y -> z, z -> -y and x -> x.  I believe
that your solution will map y -> z, z -> -x and x -> -y which is
probably not what is wanted.  (I may have the details wrong, but I
think the idea is right.)

Here is some C code that will do what the poster wants.  To produce
the input vector use

axis = V1 x V2   (May have the order backwards.  Can't remember.)
len = length( axis )
axis = axis * arcsin( len ) / len

Since the cross product produces a vector whose length is the
sine of the angle between the two vectors, we must scale the
vector as above.

/* Takes a vector whose direction defines an axis of rotation and whose
 * length is the angle of rotation in radians and fills in the corresponding
 * (4x4) rotation matrix.  Obtained from a paper by Michael Pique */
void      vector_rot(drot, rotation)
double   *drot;
Matrix    rotation;
{
  double    s, c, t, a;
  int       i, j;

  a = sqrt(drot[0] * drot[0] + drot[1] * drot[1] + drot[2] * drot[2]);

  if (a != 0.0)
  {
    drot[0] /= a;
    drot[1] /= a;
    drot[2] /= a;
  }

  s = sin(a);
  c = cos(a);
  t = 1.0 - c;

  rotation[0][0] = t * drot[0] * drot[0] + c;
  rotation[0][1] = t * drot[0] * drot[1] + s * drot[2];
  rotation[0][2] = t * drot[0] * drot[2] - s * drot[1];
  rotation[0][3] = 0;
  rotation[1][0] = t * drot[0] * drot[1] - s * drot[2];
  rotation[1][1] = t * drot[1] * drot[1] + c;
  rotation[1][2] = t * drot[1] * drot[2] + s * drot[0];
  rotation[1][3] = 0;
  rotation[2][0] = t * drot[0] * drot[2] + s * drot[1];
  rotation[2][1] = t * drot[1] * drot[2] - s * drot[0];
  rotation[2][2] = t * drot[2] * drot[2] + c;
  rotation[2][3] = 0;
  rotation[3][0] = 0;
  rotation[3][1] = 0;
  rotation[3][2] = 0;
  rotation[3][3] = 1;
}

Dan Christensen
jdchrist at watcgl.uwaterloo.ca