home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / quatspin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-17  |  1.5 KB  |  60 lines
  1. /*
  2.  * Spherical linear interpolation of unit quaternions with spins
  3.  */
  4. #include <math.h>
  5.  
  6. #define EPSILON 1.0E-6             /* a tiny number */
  7. #define TRUE  1
  8. #define FALSE 0
  9.  
  10. typedef struct {            /* quaternion type */
  11.     double w, x, y, z;
  12. } Quaternion;
  13.  
  14. slerp(alpha, a, b, q, spin)
  15.     double alpha;            /* interpolation parameter (0 to 1) */
  16.     Quaternion *a, *b;        /* start and end unit quaternions */
  17.     Quaternion *q;            /* output interpolated quaternion */
  18.     int spin;            /* number of extra spin rotations */
  19. {
  20.     double beta;            /* complementary interp parameter */
  21.     double theta;            /* angle between A and B */
  22.     double sin_t, cos_t;        /* sine, cosine of theta */
  23.     double phi;            /* theta plus spins */
  24.     int bflip;            /* use negation of B? */
  25.  
  26.     /* cosine theta = dot product of A and B */
  27.     cos_t = a->x*b->x + a->y*b->y + a->z*b->z + a->w*b->w;
  28.  
  29.     /* if B is on opposite hemisphere from A, use -B instead */
  30.      if (cos_t < 0.0) {
  31.         cos_t = -cos_t;
  32.         bflip = TRUE;
  33.     } else
  34.         bflip = FALSE;
  35.  
  36.     /* if B is (within precision limits) the same as A,
  37.      * just linear interpolate between A and B.
  38.      * Can't do spins, since we don't know what direction to spin.
  39.      */
  40.     if (1.0 - cos_t < EPSILON) {
  41.         beta = 1.0 - alpha;
  42.      } else {                /* normal case */
  43.          theta = acos(cos_t);
  44.          phi = theta + spin * M_PI;
  45.          sin_t = sin(theta);
  46.          beta = sin(theta - alpha*phi) / sin_t;
  47.          alpha = sin(alpha*phi) / sin_t;
  48.      }
  49.  
  50.     if (bflip)
  51.         alpha = -alpha;
  52.  
  53.     /* interpolate */
  54.      q->x = beta*a->x + alpha*b->x;
  55.      q->y = beta*a->y + alpha*b->y;
  56.      q->z = beta*a->z + alpha*b->z;
  57.      q->w = beta*a->w + alpha*b->w;
  58. }
  59.  
  60.