home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 March / macformat-022.iso / Shareware City / Graphics / SPD / Sources / libvec.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-10-27  |  12.6 KB  |  540 lines  |  [TEXT/R*ch]
  1. /*
  2.  * libvec.c - a library of vector operations and a random number generator
  3.  *
  4.  * Author:  Eric Haines, 3D/Eye, Inc.
  5.  *
  6.  * Modified: 1 October 1992
  7.  *         Alexander R. Enzmann
  8.  *         Object generation routines split off to keep file size down
  9.  *
  10.  * Modified: 17 Jan 1993
  11.  *         Eduard [esp] Schwan
  12.  *         Removed unused local variables, added <stdlib.h> include for
  13.  *         exit() call
  14.  */
  15.  
  16. #include <stdio.h>
  17. #include <stdlib.h> /* exit() */
  18. #include <math.h>
  19. #include "libvec.h"
  20.  
  21. /*
  22.  * Normalize the vector (X,Y,Z) so that X*X + Y*Y + Z*Z = 1.
  23.  *
  24.  * The normalization divisor is returned.  If the divisor is zero, no
  25.  * normalization occurs.
  26.  *
  27.  */
  28. double
  29. lib_normalize_vector(cvec)
  30.     COORD3 cvec;
  31. {
  32.     double divisor;
  33.  
  34.     divisor = sqrt( (double)DOT_PRODUCT(cvec, cvec) );
  35.     if (divisor > 0.0) {
  36.     cvec[X] /= divisor;
  37.     cvec[Y] /= divisor;
  38.     cvec[Z] /= divisor;
  39.     }
  40.     return divisor;
  41. }
  42.  
  43.  
  44. /*
  45.  * Set all matrix elements to zero.
  46.  */
  47. void
  48. lib_zero_matrix(mx)
  49.     MATRIX mx;
  50. {
  51.     int i, j;
  52.  
  53.     for (i=0;i<4;++i)
  54.     for (j=0;j<4;++j)
  55.         mx[i][j] = 0.0;
  56. }
  57.  
  58. /*
  59.  * Create identity matrix.
  60.  */
  61. void
  62. lib_create_identity_matrix(mx)
  63.     MATRIX mx;
  64. {
  65.     int i;
  66.  
  67.     lib_zero_matrix(mx);
  68.     for (i=0;i<4;++i)
  69.     mx[i][i] = 1.0;
  70. }
  71.  
  72. /*
  73.  * Copy one matrix to another
  74.  */
  75. void
  76. lib_copy_matrix(mres, mx)
  77.     MATRIX mres, mx;
  78. {
  79.     int i, j;
  80.  
  81.     for (i=0;i<4;i++)
  82.     for (j=0;j<4;j++)
  83.         mres[i][j] = mx[i][j];
  84. }
  85.  
  86. /*
  87.  * Create a rotation matrix along the given axis by the given angle in radians.
  88.  */
  89. void
  90. lib_create_rotate_matrix(mx, axis, angle)
  91.     MATRIX mx;
  92.     int axis;
  93.     double angle;
  94. {
  95.     double cosine, sine;
  96.  
  97.     lib_zero_matrix(mx);
  98.     cosine = cos((double)angle);
  99.     sine = sin((double)angle);
  100.     switch (axis) {
  101.     case X_AXIS:
  102.         mx[0][0] = 1.0;
  103.         mx[1][1] = cosine;
  104.         mx[2][2] = cosine;
  105.         mx[1][2] = sine;
  106.         mx[2][1] = -sine;
  107.         break ;
  108.     case Y_AXIS:
  109.         mx[1][1] = 1.0;
  110.         mx[0][0] = cosine;
  111.         mx[2][2] = cosine;
  112.         mx[2][0] = sine;
  113.         mx[0][2] = -sine;
  114.         break;
  115.     case Z_AXIS:
  116.         mx[2][2] = 1.0;
  117.         mx[0][0] = cosine;
  118.         mx[1][1] = cosine;
  119.         mx[0][1] = sine;
  120.         mx[1][0] = -sine;
  121.         break;
  122.     }
  123.     mx[3][3] = 1.0;
  124. }
  125.  
  126. /*
  127.  * Create a rotation matrix along the given axis by the given angle in radians.
  128.  * The axis is a set of direction cosines.
  129.  */
  130. void
  131. lib_create_axis_rotate_matrix(mx, axis, angle)
  132.     MATRIX mx;
  133.     COORD3 axis;
  134.     double angle;
  135. {
  136.     double  cosine, one_minus_cosine, sine;
  137.  
  138.     cosine = cos((double)angle);
  139.     sine = sin((double)angle);
  140.     one_minus_cosine = 1.0 - cosine;
  141.  
  142.     mx[0][0] = SQR(axis[X]) + (1.0 - SQR(axis[X])) * cosine;
  143.     mx[0][1] = axis[X] * axis[Y] * one_minus_cosine + axis[Z] * sine;
  144.     mx[0][2] = axis[X] * axis[Z] * one_minus_cosine - axis[Y] * sine;
  145.     mx[0][3] = 0.0;
  146.  
  147.     mx[1][0] = axis[X] * axis[Y] * one_minus_cosine - axis[Z] * sine;
  148.     mx[1][1] = SQR(axis[Y]) + (1.0 - SQR(axis[Y])) * cosine;
  149.     mx[1][2] = axis[Y] * axis[Z] * one_minus_cosine + axis[X] * sine;
  150.     mx[1][3] = 0.0;
  151.  
  152.     mx[2][0] = axis[X] * axis[Z] * one_minus_cosine + axis[Y] * sine;
  153.     mx[2][1] = axis[Y] * axis[Z] * one_minus_cosine - axis[X] * sine;
  154.     mx[2][2] = SQR(axis[Z]) + (1.0 - SQR(axis[Z])) * cosine;
  155.     mx[2][3] = 0.0;
  156.  
  157.     mx[3][0] = 0.0;
  158.     mx[3][1] = 0.0;
  159.     mx[3][2] = 0.0;
  160.     mx[3][3] = 1.0;
  161. }
  162.  
  163. /* Create translation matrix */
  164. void
  165. lib_create_translate_matrix(mx, vec)
  166.     MATRIX mx;
  167.     COORD3 vec;
  168. {
  169.     lib_create_identity_matrix(mx);
  170.     mx[3][0] = vec[X];
  171.     mx[3][1] = vec[Y];
  172.     mx[3][2] = vec[Z];
  173. }
  174.  
  175. /* Create scaling matrix */
  176. void
  177. lib_create_scale_matrix(mx, vec)
  178.     MATRIX mx;
  179.     COORD3 vec;
  180. {
  181.     lib_zero_matrix(mx);
  182.     mx[0][0] = vec[X];
  183.     mx[1][1] = vec[Y];
  184.     mx[2][2] = vec[Z];
  185. }
  186.  
  187. /* Given a point and a direction, find the transform that brings a point
  188.    in a canonical coordinate system into a coordinate system defined by
  189.    that position and direction.    Both the forward and inverse transform
  190.    matrices are calculated. */
  191. void
  192. lib_create_canonical_matrix(trans, itrans, origin, up)
  193.     MATRIX trans, itrans;
  194.     COORD3 origin;
  195.     COORD3 up;
  196. {
  197.     MATRIX trans1, trans2;
  198.     COORD3 tmpv;
  199.     double ang;
  200.  
  201.     /* Translate "origin" to <0, 0, 0> */
  202.     SET_COORD3(tmpv, -origin[X], -origin[Y], -origin[Z]);
  203.     lib_create_translate_matrix(trans1, tmpv);
  204.  
  205.     /* Determine the axis to rotate about */
  206.     if (fabs(up[Z]) == 1.0)
  207.     SET_COORD3(tmpv, 1.0, 0.0, 0.0)
  208.     else
  209.     SET_COORD3(tmpv, -up[Y], up[X], 0.0)
  210.     lib_normalize_vector(tmpv);
  211.     ang = acos(up[Z]);
  212.  
  213.     /* Calculate the forward matrix */
  214.     lib_create_axis_rotate_matrix(trans2, tmpv, -ang);
  215.     lib_matrix_multiply(trans, trans1, trans2);
  216.  
  217.     /* Calculate the inverse transform */
  218.     lib_create_axis_rotate_matrix(trans1, tmpv, ang);
  219.     lib_create_translate_matrix(trans2, origin);
  220.     lib_matrix_multiply(itrans, trans1, trans2);
  221. }
  222.  
  223. /* Find the transformation that takes the current eye position and
  224.    orientation and puts it at {0, 0, 0} with the up vector aligned
  225.    with {0, 1, 0} */
  226. void
  227. lib_create_view_matrix(trans, from, at, up, xres, yres, angle, aspect)
  228.     MATRIX trans;
  229.     COORD3 from, at, up;
  230.     int xres, yres;
  231.     double angle, aspect;
  232. {
  233.     double x, y, z, d, theta;
  234.     double xscale, yscale, view_dist;
  235.     COORD3 Rt, right, up_dir;
  236.     COORD3 Va, Ve;
  237.     MATRIX Tv, Tt, Tx;
  238.  
  239.     /* Change the view angle into radians */
  240.     angle = PI * angle / 180.0;
  241.  
  242.     /* Translate point of interest to the origin */
  243.     SET_COORD3(Va, -at[X], -at[Y], -at[Z]);
  244.     lib_create_translate_matrix(Tv, Va);
  245.  
  246.     /* Move the eye by the same amount */
  247.     ADD3_COORD3(Ve, from, Va);
  248.  
  249.     /* Make sure this is a valid sort of setup */
  250.     if (Ve[X] == 0.0 && Ve[Y] == 0.0 && Ve[Z] == 0.0) {
  251.     fprintf(stderr, "Degenerate perspective transformation\n");
  252.     exit(1);
  253.     }
  254.  
  255.     /* Get the up vector to be perpendicular to the view vector */
  256.     SET_COORD3(Rt, -Ve[X], -Ve[Y], -Ve[Z]);
  257.     lib_normalize_vector(Rt);
  258.     COPY_COORD3(up_dir, up);
  259.     CROSS(right, up_dir, Rt);
  260.     lib_normalize_vector(right);
  261.     CROSS(up_dir, Rt, right);
  262.     lib_normalize_vector(up_dir);
  263.  
  264.     /* Determine the angle of rotation about the x-axis that takes
  265.       the eye into the x-z plane. */
  266.     y = Ve[Y];
  267.     if (fabs(y) > EPSILON) {
  268.     z = Ve[Z];
  269.     d = sqrt(y * y + z * z);
  270.     if (fabs(d) > EPSILON) {
  271.         if (y > 0)
  272.         theta = acos(z/d) - PI;
  273.         else
  274.         theta = -acos(z/d) - PI;
  275.         lib_create_rotate_matrix(Tt, X_AXIS, theta);
  276.         lib_matrix_multiply(Tx, Tv, Tt);
  277.         lib_copy_matrix(Tv, Tx);
  278.     }
  279.     }
  280.  
  281.     /* Determine the angle of rotation about the y-axis that takes
  282.     the eye onto the minus z-axis */
  283.     lib_transform_point(Ve, from, Tv);
  284.     x = Ve[X];
  285.     z = Ve[Z];
  286.     if (fabs(x) > EPSILON) {
  287.     d = sqrt(x * x + z * z);
  288.     if (fabs(d) > EPSILON) {
  289.         if (x > 0)
  290.         theta = PI - acos(z / d);
  291.         else
  292.         theta = PI + acos(z / d);
  293.         lib_create_rotate_matrix(Tt, Y_AXIS, theta);
  294.         lib_matrix_multiply(Tx, Tv, Tt);
  295.         lib_copy_matrix(Tv, Tx);
  296.     }
  297.     } else if (z > 0) {
  298.     lib_create_rotate_matrix(Tt, Y_AXIS, PI);
  299.     lib_matrix_multiply(Tx, Tv, Tt);
  300.     lib_copy_matrix(Tv, Tx);
  301.     }
  302.  
  303.     /* Determine the angle of rotation about the z-axis that takes
  304.     the up vector into alignment with the y-axis */
  305.     lib_transform_vector(Ve, up_dir, Tv);
  306.     y = Ve[Y];
  307.     if (y < 1.0) {
  308.     x = Ve[X];
  309.     d = sqrt(x * x + y * y);
  310.     if (fabs(d) > EPSILON) {
  311.         if (x > 0)
  312.         theta = acos(y / d);
  313.         else
  314.         theta = -acos(y / d);
  315.         lib_create_rotate_matrix(Tt, Z_AXIS, theta);
  316.         lib_matrix_multiply(Tx, Tv, Tt);
  317.         lib_copy_matrix(Tv, Tx);
  318.     }
  319.     }
  320.  
  321.     /* Finally, translate the eye position to the point {0, 0, 0} */
  322.     SUB3_COORD3(Ve, at, from);
  323.     view_dist = sqrt(DOT_PRODUCT(Ve, Ve));
  324.     lib_transform_point(Ve, from, Tv);
  325.     SET_COORD3(Va, 0, 0, -Ve[Z]);
  326.     lib_create_translate_matrix(Tt, Va);
  327.     lib_matrix_multiply(Tx, Tv, Tt);
  328.     lib_copy_matrix(Tv, Tx);
  329.  
  330.     /* Determine how much to scale things by */
  331.     yscale = (double)yres / (2.0 * view_dist * tan(angle/2.0));
  332.     xscale = yscale * (double)xres / ((double)yres * aspect);
  333.     SET_COORD3(Va, xscale, yscale, 1.0);
  334.     lib_create_scale_matrix(Tt, Va);
  335.     lib_matrix_multiply(Tx, Tv, Tt);
  336.     lib_copy_matrix(Tv, Tx);
  337.  
  338.     /* Now add in the perspective transformation */
  339.     lib_create_identity_matrix(Tt);
  340.     Tt[2][3] = 1.0 / view_dist;
  341.     lib_matrix_multiply(Tx, Tv, Tt);
  342.  
  343.     /* We now have the final transformation matrix */
  344.     lib_copy_matrix(trans, Tx);
  345. }
  346.  
  347. /*
  348.  * Multiply a vector by a matrix.
  349.  */
  350. void
  351. lib_transform_vector(vres, vec, mx)
  352.     COORD3 vres, vec;
  353.     MATRIX mx;
  354. {
  355.     vres[X] = vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0];
  356.     vres[Y] = vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1];
  357.     vres[Z] = vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2];
  358. }
  359.  
  360. /*
  361.  * Multiply a point by a matrix.
  362.  */
  363. void
  364. lib_transform_point(vres, vec, mx)
  365.     COORD3 vres, vec;
  366.     MATRIX mx;
  367. {
  368.     vres[X] = vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0] + mx[3][0];
  369.     vres[Y] = vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1] + mx[3][1];
  370.     vres[Z] = vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2] + mx[3][2];
  371. }
  372.  
  373. /*
  374.  * Multiply a 4 element vector by a matrix.
  375.  */
  376. void
  377. lib_transform_coord(vres, vec, mx)
  378.     COORD4 vres, vec;
  379.     MATRIX mx;
  380. {
  381.     vres[X] = vec[X]*mx[0][0]+vec[Y]*mx[1][0]+vec[Z]*mx[2][0]+vec[W]*mx[3][0];
  382.     vres[Y] = vec[X]*mx[0][1]+vec[Y]*mx[1][1]+vec[Z]*mx[2][1]+vec[W]*mx[3][1];
  383.     vres[Z] = vec[X]*mx[0][2]+vec[Y]*mx[1][2]+vec[Z]*mx[2][2]+vec[W]*mx[3][2];
  384.     vres[W] = vec[X]*mx[0][3]+vec[Y]*mx[1][3]+vec[Z]*mx[2][3]+vec[W]*mx[3][3];
  385. }
  386.  
  387. /*
  388.  * Compute transpose of matrix.
  389.  */
  390. void
  391. lib_transpose_matrix( mxres, mx )
  392.     MATRIX mxres, mx;
  393. {
  394.     int i, j;
  395.  
  396.     for (i=0;i<4;i++)
  397.     for (j=0;j<4;j++)
  398.         mxres[j][i] = mx[i][j];
  399. }
  400.  
  401. /*
  402.  * Multiply two 4x4 matrices.
  403.  */
  404. void
  405. lib_matrix_multiply(mxres, mx1, mx2)
  406.     MATRIX mxres, mx1, mx2;
  407. {
  408.     int i, j;
  409.  
  410.     for (i=0;i<4;i++)
  411.     for (j=0;j<4;j++)
  412.         mxres[i][j] = mx1[i][0]*mx2[0][j] + mx1[i][1]*mx2[1][j] +
  413.                   mx1[i][2]*mx2[2][j] + mx1[i][3]*mx2[3][j];
  414. }
  415.  
  416. /* Performs a 3D clip of a line segment from start to end against the
  417.     box defined by bounds.  The actual values of start and end are modified. */
  418. int
  419. lib_clip_to_box(start, end,  bounds)
  420.     COORD3 start, end;
  421.     double bounds[2][3];
  422. {
  423.     int i;
  424.     double d, t, dir, pos;
  425.     double tmin, tmax;
  426.     COORD3 D;
  427.  
  428.     SUB3_COORD3(D, end, start);
  429.     d = lib_normalize_vector(D);
  430.     tmin = 0.0;
  431.     tmax = d;
  432.     for (i=0;i<3;i++) {
  433.     pos = start[i];
  434.     dir = D[i];
  435.     if (dir < -EPSILON) {
  436.         t = (bounds[0][i] - pos) / dir;
  437.         if (t < tmin)
  438.         return 0;
  439.         if (t <= tmax)
  440.         tmax = t;
  441.         t = (bounds[1][i] - pos) / dir;
  442.         if (t >= tmin) {
  443.         if (t > tmax)
  444.             return 0;
  445.         tmin = t;
  446.         }
  447.     } else if (dir > EPSILON) {
  448.         t = (bounds[1][i] - pos) / dir;
  449.         if (t < tmin)
  450.         return 0;
  451.         if (t <= tmax)
  452.         tmax = t;
  453.         t = (bounds[0][i] - pos) / dir;
  454.         if (t >= tmin) {
  455.         if (t > tmax)
  456.             return 0;
  457.         tmin = t;
  458.         }
  459.     }
  460.     else if (pos < bounds[0][i] || pos > bounds[1][i])
  461.         return 0;
  462.     }
  463.     if (tmax < d) {
  464.     end[X] = start[X] + tmax * D[X];
  465.     end[Y] = start[Y] + tmax * D[Y];
  466.     end[Z] = start[Z] + tmax * D[Z];
  467.     }
  468.     if (tmin > 0.0) {
  469.     start[X] = start[X] + tmin * D[X];
  470.     start[Y] = start[Y] + tmin * D[Y];
  471.     start[Z] = start[Z] + tmin * D[Z];
  472.     }
  473.     return 1;
  474. }
  475.  
  476. /*
  477.  * Rotate a vector pointing towards the major-axis faces (i.e. the major-axis
  478.  * component of the vector is defined as the largest value) 90 degrees to
  479.  * another cube face.  Mod_face is a face number.
  480.  *
  481.  * If the routine is called six times, with mod_face=0..5, the vector will be
  482.  * rotated to each face of a cube.  Rotations are:
  483.  *   mod_face = 0 mod 3, +Z axis rotate
  484.  *   mod_face = 1 mod 3, +X axis rotate
  485.  *   mod_face = 2 mod 3, -Y axis rotate
  486.  */
  487. void
  488. lib_rotate_cube_face(vec, major_axis, mod_face)
  489.     COORD3 vec;
  490.     int major_axis, mod_face;
  491. {
  492.     double swap;
  493.  
  494.     mod_face = (mod_face+major_axis) % 3 ;
  495.     if (mod_face == 0) {
  496.     swap     = vec[X];
  497.     vec[X] = -vec[Y];
  498.     vec[Y] = swap;
  499.     } else if (mod_face == 1) {
  500.     swap     = vec[Y];
  501.     vec[Y] = -vec[Z];
  502.     vec[Z] = swap;
  503.     } else {
  504.     swap     = vec[X];
  505.     vec[X] = -vec[Z];
  506.     vec[Z] = swap;
  507.     }
  508. }
  509.  
  510. /*
  511.  * Portable gaussian random number generator (from "Numerical Recipes", GASDEV)
  512.  * Returns a uniform random deviate between 0.0 and 1.0.  'iseed' must be
  513.  * less than M1 to avoid repetition, and less than (2**31-C1)/A1 [= 300718]
  514.  * to avoid overflow.
  515.  */
  516. #define M1  134456
  517. #define IA1   8121
  518. #define IC1  28411
  519. #define RM1 1.0/M1
  520.  
  521. double
  522. lib_gauss_rand(iseed)
  523.     long iseed;
  524. {
  525.     double fac, v1, v2, r;
  526.     long     ix1, ix2;
  527.  
  528.     ix2 = iseed;
  529.     do {
  530.     ix1 = (IC1+ix2*IA1) % M1;
  531.     ix2 = (IC1+ix1*IA1) % M1;
  532.     v1 = ix1 * 2.0 * RM1 - 1.0;
  533.     v2 = ix2 * 2.0 * RM1 - 1.0;
  534.     r = v1*v1 + v2*v2;
  535.     } while ( r >= 1.0 );
  536.  
  537.     fac = sqrt((double)(-2.0 * log((double)r) / r));
  538.     return v1 * fac;
  539. }
  540.