home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / raymath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  10.7 KB  |  600 lines
  1. /*
  2.  * Raymath.c -- simple math routines.
  3.  * 
  4.  * (c) 1993, 1994 by Han-Wen Nienhuys <hanwen@stack.urc.tue.nl>
  5.  * 
  6.  * based on work by George Kyriazis, 1988
  7.  * 
  8.  * This program is free software; you can redistribute it and/or modify it
  9.  * under the terms of the GNU General Public License as published by the
  10.  * Free Software Foundation;
  11.  * 
  12.  * This program is distributed in the hope that it will be useful, but
  13.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  14.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15.  * General Public License for more details.
  16.  * 
  17.  * You should have received a copy of the GNU General Public License along
  18.  * with this program; if not, write to the Free Software Foundation, Inc.,
  19.  * 675 Mass Ave, Cambridge, MA 02139, USA.
  20.  */
  21.  
  22. /*
  23.  * math and vector operations
  24.  */
  25.  
  26. #include     <stdlib.h>
  27. #include    "ray.h"
  28. #include    "proto.h"
  29. #include    "extern.h"
  30. #include    "random.h"
  31.  
  32. /*
  33.  * floating point random number generator, from K&R chap 7.8.7 [deleted]
  34.  * now with table lookup.
  35.  */
  36.  
  37. PRIVATE         randidx;
  38.  
  39. PUBLIC double
  40. frand(void)
  41. {
  42.     if (++randidx > RANDOMTABSIZE)
  43.     randidx = rand() % RANDOMTABSIZE;
  44.  
  45.     return randtable[randidx];
  46. }
  47.  
  48.  
  49.  
  50. /*
  51.  * Random number between -1 and 1
  52.  */
  53. PUBLIC double
  54. rand1(void)
  55. {
  56.     return (frand() - 0.5) * 2;
  57. }
  58.  
  59.  
  60. /*
  61.  * approximate a gaussian random number generator between -1 and 1
  62.  */
  63. PUBLIC double
  64. Gauss_rand(void)
  65. {
  66.     double          t;
  67.  
  68.     /* get a random number between -1 and 1 */
  69.     t = rand1();
  70.  
  71.     /* and then square it.  Don't lose the sign! */
  72.     return (t * t * SGN(t));
  73. }
  74.  
  75.  
  76. /* check wether a vector is the nullvector */
  77. PUBLIC bool
  78. isvnull(vector a)
  79. {
  80.     if (quickveclen(a) < EPSILON)
  81.     return TRUE;
  82.     else
  83.     return FALSE;
  84. }
  85.  
  86. /*
  87.  * evaluate a polynomial given by coef, ord at the point x. Turbo C has a
  88.  * function poly() to do this, but it is not ANSI-C, so this one is
  89.  * included, just to be on the safe side.
  90.  * 
  91.  * It's taken from Vort.
  92.  */
  93.  
  94. PUBLIC double
  95. poly_eval(register double x, int ord, double *coef)
  96. {
  97.     register double *fp,
  98.                     f;
  99.  
  100.     fp = &coef[ord];
  101.     f = *fp;
  102.  
  103.     for (fp--; fp >= coef; fp--)
  104.     f = x * f + *fp;
  105.  
  106.     return f;
  107. }
  108.  
  109. /****************************************************************************
  110. Performs a safe arctangent calculation for the two provided numbers.  If the
  111. denominator is 0.o, atan2 is not called (it results in a floating exception
  112. on some machines).  Instead, + or - PI/2 is returned.
  113. ****************************************************************************/
  114.  
  115. PUBLIC double
  116. safe_arctangent(double a, double b)
  117. {
  118.     if (b != 0.0)
  119.     return atan2(a, b);
  120.     else if (a > 0.0)
  121.     return M_PI_2;
  122.     else
  123.     return -M_PI_2;
  124. }
  125.  
  126.  
  127.  
  128. /* These are not used, so they are commented out. */
  129. #ifdef UNDEFINED
  130. PUBLIC double
  131. quickcos(double x)
  132. {
  133.     double          val;
  134.  
  135.     val = 1 - x * x / 2.4684;
  136.  
  137.     return val;
  138. }
  139.  
  140. /* not used */
  141. PUBLIC double
  142. quickinvcos(double x)
  143. {
  144.     double          val;
  145.  
  146.     val = sqrt(2.4684 * (1 - x));
  147.  
  148.     return val;
  149. }
  150.  
  151. #endif
  152.  
  153. /* just a standard routine ... */
  154. PUBLIC void
  155. copy_vector(vector *dst, vector *src)
  156. {
  157.     assert(src != NULL && src != NULL);
  158.  
  159.     *dst = *src;
  160. }
  161.  
  162. /* transpose a matrix. WATCH IT! don't do transpose_matrix(A,A); */
  163. PUBLIC void
  164. transpose_matrix(matrix tr, matrix org)
  165. {
  166.     int             i,
  167.                     j;
  168.  
  169.     for (i = 0; i < 4; i++)
  170.     for (j = 0; j < 4; j++)
  171.         tr[i][j] = org[j][i];
  172.  
  173. }
  174.  
  175. PUBLIC void
  176. copy_matrix(matrix dst, matrix src)
  177. {
  178.     int             i,
  179.                     j;
  180.  
  181.     for (i = 0; i < 4; i++)
  182.     for (j = 0; j < 4; j++)
  183.         dst[i][j] = src[i][j];
  184. }
  185.  
  186. /*
  187.  * multiply matrix a with matrix, in that order. mmproduct(A,A,B) will
  188.  * give valid results
  189.  */
  190. PUBLIC void
  191. mmproduct(matrix c, matrix a, matrix b)
  192. {
  193.     int             i,
  194.                     j,
  195.                     k;
  196.     matrix          tmp;
  197.  
  198.     for (i = 0; i < 4; i++)
  199.     for (j = 0; j < 4; j++) {
  200.         tmp[i][j] = 0;
  201.         for (k = 0; k < 4; k++)
  202.         tmp[i][j] += a[i][k] * b[k][j];
  203.     }
  204.  
  205.     copy_matrix(c, tmp);
  206. }
  207.  
  208. /*
  209.  * matrix vector product: multiply  a 3x3 (or 4x4) matrix with
  210.  * columnvector x
  211.  */
  212. PUBLIC vector
  213. mvproduct(matrix A, vector x)
  214. {
  215.     vector          v;
  216.  
  217.     v.x = A[0][0] * x.x + A[0][1] * x.y + A[0][2] * x.z + A[0][3];
  218.     v.y = A[1][0] * x.x + A[1][1] * x.y + A[1][2] * x.z + A[1][3];
  219.     v.z = A[2][0] * x.x + A[2][1] * x.y + A[2][2] * x.z + A[2][3];
  220.  
  221.     return v;
  222. }
  223.  
  224. /* multiply with transpose of A */
  225. PUBLIC vector
  226. transform_normal(matrix A, vector n)
  227. {
  228.     vector          v;
  229.  
  230.     v.x = n.x * A[0][0] + n.y * A[1][0] + n.z * A[2][0];
  231.     v.y = n.x * A[0][1] + n.y * A[1][1] + n.z * A[2][1];
  232.     v.z = n.x * A[0][2] + n.y * A[1][2] + n.z * A[2][2];
  233.  
  234.     return v;
  235. }
  236.  
  237. /*
  238.  * matrix-vector multiply, but don't use translation.
  239.  */
  240. PUBLIC vector
  241. transform_direction(matrix A, vector x)
  242. {
  243.     vector          v;
  244.  
  245.     v.x = A[0][0] * x.x + A[0][1] * x.y + A[0][2] * x.z;
  246.     v.y = A[1][0] * x.x + A[1][1] * x.y + A[1][2] * x.z;
  247.     v.z = A[2][0] * x.x + A[2][1] * x.y + A[2][2] * x.z;
  248.  
  249.     return v;
  250. }
  251.  
  252. /*
  253.  * transform a ray using the object's inverse transform
  254.  */
  255. PUBLIC void
  256. transform_ray(struct ray *r, object *o)
  257. {
  258.  
  259.     if (o->inv_trans == NULL)
  260.     return;
  261.     r->dir = transform_direction(*o->inv_trans, r->dir);
  262.  
  263.     /* remember! r->dir is NOT normalized */
  264.  
  265.     r->pos = mvproduct(*o->inv_trans, r->pos);
  266. }
  267.  
  268.  
  269. /* initialize m to the null matrix */
  270. PUBLIC void
  271. null_matrix(matrix m)
  272. {
  273.     int             i,
  274.                     j;
  275.  
  276.     for (i = 0; i < 4; i++)
  277.     for (j = 0; j < 4; j++)
  278.         m[i][j] = 0;
  279. }
  280.  
  281. /* m := diag(1,1,..) */
  282. PUBLIC void
  283. unit_matrix(matrix m)
  284. {
  285.     int             i,
  286.                     j;
  287.  
  288.     for (i = 0; i < 4; i++)
  289.     for (j = 0; j < 4; j++)
  290.         if (i != j)
  291.         m[i][j] = 0;
  292.         else
  293.         m[i][j] = 1;
  294.  
  295. }
  296.  
  297. /* ?rotate() generates matrix for lefthanded rotation about ? axis. */
  298. PRIVATE void
  299. xrotate(matrix m, double angle)
  300. {
  301.     unit_matrix(m);
  302.     m[1][1] = m[2][2] = cos(angle);
  303.     m[1][2] = -sin(angle);
  304.     m[2][1] = -m[1][2];
  305. }
  306.  
  307. PRIVATE void
  308. yrotate(matrix m, double angle)
  309. {
  310.     unit_matrix(m);
  311.     m[0][0] = m[2][2] = cos(angle);
  312.     m[0][2] = sin(angle);
  313.     m[2][0] = -m[0][2];
  314. }
  315.  
  316. PRIVATE void
  317. zrotate(matrix m, double angle)
  318. {
  319.     unit_matrix(m);
  320.     m[0][0] = m[1][1] = cos(angle);
  321.     m[1][0] = sin(angle);
  322.     m[0][1] = -m[1][0];
  323. }
  324.  
  325. /*
  326.  * generate matrix for rotation around x (angle xa), y (angle ya), z
  327.  * (angle za), in that order.
  328.  */
  329. PUBLIC void
  330. rotate_matrix(matrix m, vector r)
  331. {
  332.     matrix          a,
  333.                     b,
  334.                     c,
  335.                     d;
  336.  
  337.     xrotate(a, r.x);
  338.     yrotate(b, r.y);
  339.     mmproduct(d, b, a);
  340.  
  341.     zrotate(c, r.z);
  342.     mmproduct(m, c, d);
  343. }
  344.  
  345.  
  346. /* generate a matrix which is the inverse of the rotation matrix of r */
  347. PUBLIC void
  348. inv_rotate_matrix(matrix m, vector r)
  349. {
  350.     matrix          rm;
  351.  
  352.     rotate_matrix(rm, r);
  353.     transpose_matrix(m, rm);
  354. }
  355.  
  356. /* generate a matrix for scaling with the vector s */
  357. PUBLIC void
  358. scale_matrix(matrix m, vector s)
  359. {
  360.     unit_matrix(m);
  361.     m[0][0] = s.x;
  362.     m[1][1] = s.y;
  363.     m[2][2] = s.z;
  364. }
  365.  
  366.  
  367. /* inverse of scale_matrix() */
  368. PUBLIC void
  369. inv_scale_matrix(matrix m, vector s)
  370. {
  371.     s.x = 1 / s.x;
  372.     s.y = 1 / s.y;
  373.     s.z = 1 / s.z;
  374.     scale_matrix(m, s);
  375. }
  376.  
  377. /* translate_matrix */
  378. PUBLIC void
  379. translate_matrix(matrix m, vector t)
  380. {
  381.     unit_matrix(m);
  382.     m[0][3] = t.x;
  383.     m[1][3] = t.y;
  384.     m[2][3] = t.z;
  385. }
  386.  
  387. /* inverse translation */
  388. PUBLIC void
  389. inv_translate_matrix(matrix m, vector t)
  390. {
  391.     vneg(t, t);
  392.     translate_matrix(m, t);
  393. }
  394.  
  395. /* rotate a vector, just a standard routine */
  396. PUBLIC void
  397. rotate_vector(vector *v, matrix rotmat)
  398. {
  399.     *v = mvproduct(rotmat, *v);
  400. }
  401.  
  402. /* conversion. Not really used. */
  403. PUBLIC double
  404. degtorad(double deg)
  405. {
  406.     return M_PI * deg / 180.0;
  407. }
  408.  
  409. PUBLIC double
  410. radtodeg(double rad)
  411. {
  412.     return 180.0 * rad / M_PI;
  413. }
  414.  
  415. /*
  416.  * convert a vector in degrees to a vector in radians, generate rotation
  417.  * matrix, put it in m
  418.  */
  419. PUBLIC void
  420. make_rotation_matrix(matrix m, vector v)
  421. {
  422.     v.x = degtorad(v.x);
  423.     v.y = degtorad(v.y);
  424.     v.z = degtorad(v.z);
  425.     rotate_matrix(m, v);
  426. }
  427.  
  428.  
  429. /*
  430.  * decompose an inverse transform into the orginal simple transformations
  431.  * Every combination of transforms should be writeable as
  432.  * 
  433.  * TRS
  434.  * 
  435.  * Translation, Rotation, Scale So any inverse transform is writeable as
  436.  * 
  437.  * S^-1 R^-1 T^-1
  438.  * 
  439.  */
  440.  
  441. /*
  442.  * this routine could be a lot faster, but this isn't critical, so I leave
  443.  * at this.
  444.  */
  445. PUBLIC void
  446. invert_trans(matrix out, matrix A)
  447. {
  448.     vector          row;
  449.     matrix          temp,
  450.                     rot,
  451.                     tm,
  452.                     sm;
  453.     vector          scal;
  454.  
  455.     unit_matrix(rot);
  456.     row.x = A[0][0];
  457.     row.y = A[0][1];
  458.     row.z = A[0][2];
  459.  
  460.     scal.x = 1 / veclen(row);
  461.     rot[0][0] = row.x * scal.x;
  462.     rot[1][0] = row.y * scal.x;
  463.     rot[2][0] = row.z * scal.x;
  464.  
  465.     row.x = A[1][0];
  466.     row.y = A[1][1];
  467.     row.z = A[1][2];
  468.  
  469.     scal.y = 1 / veclen(row);
  470.     rot[0][1] = row.x * scal.y;
  471.     rot[1][1] = row.y * scal.y;
  472.     rot[2][1] = row.z * scal.y;
  473.  
  474.     row.x = A[2][0];
  475.     row.y = A[2][1];
  476.     row.z = A[2][2];
  477.  
  478.     scal.z = 1 / veclen(row);
  479.     rot[0][2] = row.x * scal.z;
  480.     rot[1][2] = row.y * scal.z;
  481.     rot[2][2] = row.z * scal.z;
  482.     scale_matrix(sm, scal);
  483.  
  484.     mmproduct(temp, sm, A);
  485.     mmproduct(tm, rot, temp);
  486.     tm[0][3] = -tm[0][3];
  487.     tm[1][3] = -tm[1][3];
  488.     tm[2][3] = -tm[2][3];
  489.  
  490.     mmproduct(temp, rot, sm);
  491.     mmproduct(out, tm, temp);
  492. }
  493.  
  494.  
  495. /*
  496.  * The provided point is compared against the current minimum and maximum
  497.  * values in X, Y, and Z.  If a new maximum or minimum is found, the min
  498.  * and max variables are updated.
  499.  */
  500. /* from GGems III */
  501. PUBLIC void
  502. update_min_and_max(vector *min, vector *max, vector point)
  503. {
  504.     if (point.x < min->x)
  505.     min->x = point.x;
  506.     if (point.x > max->x)
  507.     max->x = point.x;
  508.     if (point.y < min->y)
  509.     min->y = point.y;
  510.     if (point.y > max->y)
  511.     max->y = point.y;
  512.     if (point.z < min->z)
  513.     min->z = point.z;
  514.     if (point.z > max->z)
  515.     max->z = point.z;
  516. }
  517.  
  518.  
  519. /*
  520.  * MEMORY MANAGEMENT
  521.  */
  522.  
  523.  
  524. PUBLIC matrix  *
  525. get_new_matrix(void)
  526. {
  527.     matrix         *p;
  528.     p = malloc(sizeof(matrix));
  529.  
  530.     CHECK_MEM(p, "transform");
  531.     unit_matrix(*p);
  532.     return p;
  533. }
  534.  
  535. PUBLIC void
  536. free_float(double *f)
  537. {
  538.     free((void *) f);
  539. }
  540.  
  541. PUBLIC void
  542. free_vector(vector *v)
  543. {
  544.     free((void *) v);
  545. }
  546.  
  547. PUBLIC vector  *
  548. get_new_vector(void)
  549. {
  550.     vector         *p;
  551.     p = malloc(sizeof(vector));
  552.  
  553.     CHECK_MEM(p, "vector");
  554.     setvector(*p, 0, 0, 0);
  555.     return p;
  556. }
  557.  
  558. PUBLIC double  *
  559. get_new_float(void)
  560. {
  561.     double         *p;
  562.  
  563.     p = ALLOC(double);
  564.  
  565.     CHECK_MEM(p, "float");
  566.  
  567.     *p = 0.0;
  568.     return p;
  569. }
  570.  
  571. /*
  572.  * DEBUGGING
  573.  */
  574.  
  575. PUBLIC void
  576. print_v(char *s, vector v)
  577. {
  578. #ifdef DEBUG
  579.     printf("%s: <%lf, %lf, %lf>\n", s, v.x, v.y, v.z);
  580. #endif
  581. }
  582.  
  583.  
  584. PUBLIC void
  585. print_matrix(matrix c)
  586. {
  587. #ifdef DEBUG
  588.     int             i,
  589.                     j;
  590.  
  591.     printf("\n");
  592.     for (i = 0; i < 4; i++) {
  593.     for (j = 0; j < 4; j++)
  594.         printf("%10lf ", c[i][j]);
  595.     printf("\n");
  596.     }
  597.  
  598. #endif
  599. }
  600.