home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / GRAPHICS / rayshade.lzh / matrix.c < prev    next >
Text File  |  1990-05-08  |  10KB  |  386 lines
  1. /*
  2.  * matrix.c
  3.  *
  4.  * Copyright (C) 1989, Craig E. Kolb
  5.  *
  6.  * This software may be freely copied, modified, and redistributed,
  7.  * provided that this copyright notice is preserved on all copies.
  8.  *
  9.  * There is no warranty or other guarantee of fitness for this software,
  10.  * it is provided solely .  Bug reports or fixes may be sent
  11.  * to the author, who may or may not act on them as he desires.
  12.  *
  13.  * You may not include this software in a program or other software product
  14.  * without supplying the source, or without informing the end-user that the
  15.  * source is available for no extra charge.
  16.  *
  17.  * If you modify this software, you should include a notice giving the
  18.  * name of the person performing the modification, the date of modification,
  19.  * and the reason for such modification.
  20.  *
  21.  * $Id: matrix.c,v 3.0.1.1 89/11/18 14:05:11 craig Exp Locker: craig $
  22.  *
  23.  * $Log:    matrix.c,v $
  24.  * Revision 3.0.1.1  89/11/18  14:05:11  craig
  25.  * patch1: Renamed rotate(), translate() and scale() to avoid problems
  26.  * patch1: with external libraries.
  27.  * 
  28.  * Revision 3.0  89/10/27  02:05:56  craig
  29.  * Baseline for first official release.
  30.  * 
  31.  */
  32. #include <math.h>
  33. #include <stdio.h>
  34. #include "typedefs.h"
  35. #include "constants.h"
  36. #include "funcdefs.h"
  37.  
  38. /*
  39.  * Matrices are indexed row-first; that is:
  40.  * matrix[ROW][COLUMN]
  41.  */
  42. /*
  43.  * Allocate new structure which holds both object-to-world and
  44.  * world-to-object space transformation structures.  It probably
  45.  * should hold pointers to these structures.
  46.  */
  47. Trans *
  48. new_trans(t, it)
  49. TransInfo *t, *it;
  50. {
  51.     Trans *res;
  52.  
  53.     res = (Trans *)Malloc(sizeof(Trans));
  54.     res->obj2world = *t;
  55.     res->world2obj = *it;
  56.     return res;
  57. }
  58.  
  59. /*
  60.  * Allocate new transformation structure.
  61.  */
  62. TransInfo *
  63. new_transinfo()
  64. {
  65.     TransInfo *res;
  66.  
  67.     res = (TransInfo *)Malloc(sizeof(TransInfo));
  68.     init_trans(res);
  69.     return res;
  70. }
  71.  
  72. /*
  73.  * Set up matrix for counter-clockwise rotation about vector (x, y, z) by
  74.  * theta radians.
  75.  */
  76. rotate_trans(trans, x, y, z, theta)
  77. TransInfo *trans;
  78. double x, y, z, theta;
  79. {
  80.     double n1, n2, n3, sintheta, costheta;
  81.     Vector vector;
  82.  
  83.     vector.x = x;
  84.     vector.y = y;
  85.     vector.z = z;
  86.     (void)normalize(&vector);
  87.  
  88.     sintheta = sin(theta);
  89.     costheta = cos(theta);
  90.  
  91.     n1 = vector.x; n2 = vector.y; n3 = vector.z;
  92.     trans->matrix[0][0] = (double)(n1*n1 + (1. - n1*n1)*costheta);
  93.     trans->matrix[1][0] = (double)(n1*n2*(1 - costheta) + n3*sintheta);
  94.     trans->matrix[2][0] = (double)(n1*n3*(1 - costheta) - n2*sintheta);
  95.     trans->matrix[0][1] = (double)(n1*n2*(1 - costheta) - n3*sintheta);
  96.     trans->matrix[1][1] = (double)(n2*n2 + (1 - n2*n2)*costheta);
  97.     trans->matrix[2][1] = (double)(n2*n3*(1 - costheta) + n1*sintheta);
  98.     trans->matrix[0][2] = (double)(n1*n3*(1 - costheta) + n2*sintheta);
  99.     trans->matrix[1][2] = (double)(n2*n3*(1 - costheta) - n1*sintheta);
  100.     trans->matrix[2][2] = (double)(n3*n3 + (1 - n3*n3)*costheta);
  101.     trans->translate.x = trans->translate.y = trans->translate.z = 0.;
  102. }
  103.  
  104. /*
  105.  * Apply translation by (vec) to trans.
  106.  */
  107. RS_translate(trans, vec)
  108. TransInfo *trans;
  109. Vector *vec;
  110. {
  111.     trans->translate.x += vec->x;
  112.     trans->translate.y += vec->y;
  113.     trans->translate.z += vec->z;
  114. }
  115.  
  116. /*
  117.  * Apply rotation about (dir) to matrix.
  118.  */
  119. RS_rotate(trans, dir, radians)
  120. TransInfo *trans;
  121. double radians;
  122. Vector *dir;
  123. {
  124.     TransInfo ttmp;
  125.  
  126.     rotate_trans(&ttmp, dir->x, dir->y, dir->z, radians);
  127.     mmult(trans, &ttmp, trans);
  128. }
  129.  
  130. /*
  131.  * Apply scale of (x, y, z) to trans.
  132.  */
  133. RS_scale(trans, x, y, z)
  134. TransInfo *trans;
  135. double x, y, z;
  136. {
  137.     trans->matrix[0][0] *= x;
  138.     trans->matrix[1][0] *= x;
  139.     trans->matrix[2][0] *= x;
  140.     trans->matrix[0][1] *= y;
  141.     trans->matrix[1][1] *= y;
  142.     trans->matrix[2][1] *= y;
  143.     trans->matrix[0][2] *= z;
  144.     trans->matrix[1][2] *= z;
  145.     trans->matrix[2][2] *= z;
  146.  
  147.     trans->translate.x *= x;
  148.     trans->translate.y *= y;
  149.     trans->translate.z *= z;
  150.  
  151. }
  152.  
  153. /*
  154.  * Multiply m1 and m2, copy result into "res".
  155.  */
  156. mmult(t1, t2, res)
  157. TransInfo *t1, *t2, *res;
  158. {
  159.     register int i, j, k;
  160.     TransInfo tmp;
  161.  
  162.     for(i = 0; i < 3; i++) {
  163.         for(j = 0; j < 3;j++) {
  164.             tmp.matrix[i][j] = 0.;
  165.             for(k = 0; k < 3;k++)
  166.                 tmp.matrix[i][j] += t1->matrix[i][k] *
  167.                             t2->matrix[k][j];
  168.  
  169.         }
  170.     }
  171.     tmp.translate.x = t1->translate.x * t2->matrix[0][0] +
  172.                 t1->translate.y * t2->matrix[1][0] +
  173.                 t1->translate.z * t2->matrix[2][0] +
  174.                     t2->translate.x;
  175.     tmp.translate.y = t1->translate.x * t2->matrix[0][1] +
  176.                 t1->translate.y * t2->matrix[1][1] +
  177.                 t1->translate.z * t2->matrix[2][1] +
  178.                     t2->translate.y;
  179.     tmp.translate.z = t1->translate.x * t2->matrix[0][2] +
  180.                 t1->translate.y * t2->matrix[1][2] +
  181.                 t1->translate.z * t2->matrix[2][2] +
  182.                     t2->translate.z;
  183.     copy_trans(res, &tmp);
  184. }
  185.  
  186. /*
  187.  * Copy a given transformation structure.
  188.  */
  189. copy_trans(into, from)
  190. TransInfo *into, *from;
  191. {
  192.     register int i;
  193.  
  194.     for(i = 0;i < 3;i++) {
  195.         into->matrix[i][0] = from->matrix[i][0];
  196.         into->matrix[i][1] = from->matrix[i][1];
  197.         into->matrix[i][2] = from->matrix[i][2];
  198.     }
  199.  
  200.     into->translate = from->translate;
  201. }
  202.  
  203. /*
  204.  * Initialize transformation structure.
  205.  */
  206. init_trans(trans)
  207. TransInfo *trans;
  208. {
  209.     register int i, j;
  210.  
  211.     for(i = 0;i < 3;i++)
  212.         for(j = 0;j < 3;j++)
  213.             trans->matrix[i][j] = (double)(i == j);
  214.  
  215.     trans->translate.x = trans->translate.y = trans->translate.z = 0.;
  216. }
  217.  
  218. /*
  219.  * Calculate inverse of the given transformation structure.
  220.  */
  221. invert_trans(inverse, trans)
  222. TransInfo *inverse, *trans;
  223. {
  224.     TransInfo ttmp;
  225.     int i;
  226.     double d;
  227.     extern int yylineno;
  228.  
  229.     ttmp.matrix[0][0] = trans->matrix[1][1]*trans->matrix[2][2] -
  230.                 trans->matrix[1][2]*trans->matrix[2][1];
  231.     ttmp.matrix[1][0] = trans->matrix[1][0]*trans->matrix[2][2] -
  232.                 trans->matrix[1][2]*trans->matrix[2][0];
  233.     ttmp.matrix[2][0] = trans->matrix[1][0]*trans->matrix[2][1] -
  234.                 trans->matrix[1][1]*trans->matrix[2][0];
  235.  
  236.     ttmp.matrix[0][1] = trans->matrix[0][1]*trans->matrix[2][2] -
  237.                 trans->matrix[0][2]*trans->matrix[2][1];
  238.     ttmp.matrix[1][1] = trans->matrix[0][0]*trans->matrix[2][2] -
  239.                 trans->matrix[0][2]*trans->matrix[2][0];
  240.     ttmp.matrix[2][1] = trans->matrix[0][0]*trans->matrix[2][1] -
  241.                 trans->matrix[0][1]*trans->matrix[2][0];
  242.  
  243.     ttmp.matrix[0][2] = trans->matrix[0][1]*trans->matrix[1][2] -
  244.                 trans->matrix[0][2]*trans->matrix[1][1];
  245.     ttmp.matrix[1][2] = trans->matrix[0][0]*trans->matrix[1][2] -
  246.                 trans->matrix[0][2]*trans->matrix[1][0];
  247.     ttmp.matrix[2][2] = trans->matrix[0][0]*trans->matrix[1][1] -
  248.                 trans->matrix[0][1]*trans->matrix[1][0];
  249.  
  250.     d = trans->matrix[0][0]*ttmp.matrix[0][0] -
  251.         trans->matrix[0][1]*ttmp.matrix[1][0] +
  252.         trans->matrix[0][2]*ttmp.matrix[2][0];
  253.  
  254.     if (abs(d) < EPSILON*EPSILON) {
  255.         fprintf(stderr,"Singular matrix (line %d):\n",yylineno);
  256.         for (i = 0; i < 3; i++)
  257.             fprintf(stderr,"%g %g %g\n", trans->matrix[i][0],
  258.                 trans->matrix[i][1], trans->matrix[i][2]);
  259.         exit(0);
  260.     }
  261.  
  262.     ttmp.matrix[0][0] /= d;
  263.     ttmp.matrix[0][2] /= d;
  264.     ttmp.matrix[1][1] /= d;
  265.     ttmp.matrix[2][0] /= d;
  266.     ttmp.matrix[2][2] /= d;
  267.  
  268.     d = -d;
  269.  
  270.     ttmp.matrix[0][1] /= d;
  271.     ttmp.matrix[1][0] /= d;
  272.     ttmp.matrix[1][2] /= d;
  273.     ttmp.matrix[2][1] /= d;
  274.  
  275.     ttmp.translate.x = -(ttmp.matrix[0][0]*trans->translate.x +
  276.                  ttmp.matrix[1][0]*trans->translate.y +
  277.                  ttmp.matrix[2][0]*trans->translate.z);
  278.     ttmp.translate.y = -(ttmp.matrix[0][1]*trans->translate.x +
  279.                  ttmp.matrix[1][1]*trans->translate.y +
  280.                  ttmp.matrix[2][1]*trans->translate.z);
  281.     ttmp.translate.z = -(ttmp.matrix[0][2]*trans->translate.x +
  282.                  ttmp.matrix[1][2]*trans->translate.y +
  283.                  ttmp.matrix[2][2]*trans->translate.z);
  284.  
  285.     copy_trans(inverse, &ttmp);
  286. }
  287.  
  288. /*
  289.  * Apply a transformation to a point (translation affects the point).
  290.  */
  291. transform_point(vec, trans)
  292. Vector *vec;
  293. TransInfo *trans;
  294. {
  295.     Vector tmp;
  296.  
  297.     tmp.x = vec->x * trans->matrix[0][0] + vec->y * trans->matrix[1][0] +
  298.             vec->z * trans->matrix[2][0] + trans->translate.x;
  299.     tmp.y = vec->x * trans->matrix[0][1] + vec->y * trans->matrix[1][1] +
  300.             vec->z * trans->matrix[2][1] + trans->translate.y;
  301.     tmp.z = vec->x * trans->matrix[0][2] + vec->y * trans->matrix[1][2] +
  302.             vec->z * trans->matrix[2][2] + trans->translate.z;
  303.     *vec = tmp;
  304. }
  305.  
  306. /*
  307.  * Apply an explicit transformation to the given transformation structure.
  308.  * 'c1x' is the X (0th) component of the first column, and so on.
  309.  */
  310. explicit_trans(trans, c1x,c1y,c1z, c2x,c2y,c2z, c3x,c3y,c3z, tx, ty, tz, res)
  311. TransInfo *trans, *res;
  312. double c1x, c1y, c1z, c2x, c2y, c2z, c3x, c3y, c3z, tx, ty, tz;
  313. {
  314.     TransInfo tmp;
  315.  
  316.     tmp.matrix[0][0] = c1x;
  317.     tmp.matrix[1][0] = c1y;
  318.     tmp.matrix[2][0] = c1z;
  319.  
  320.     tmp.matrix[0][1] = c2x;
  321.     tmp.matrix[1][1] = c2y;
  322.     tmp.matrix[2][1] = c2z;
  323.  
  324.     tmp.matrix[0][2] = c3x;
  325.     tmp.matrix[1][2] = c3y;
  326.     tmp.matrix[2][2] = c3z;
  327.  
  328.     tmp.translate.x = tx;
  329.     tmp.translate.y = ty;
  330.     tmp.translate.z = tz;
  331.  
  332.     mmult(trans, &tmp, res);
  333. }
  334.  
  335. /*
  336.  * Apply transformation to a vector (translations have no effect).
  337.  */
  338. transform_vector(vec, trans)
  339. Vector *vec;
  340. TransInfo *trans;
  341. {
  342.     Vector tmp;
  343.  
  344.     tmp.x = vec->x*trans->matrix[0][0] +
  345.         vec->y*trans->matrix[1][0] + vec->z*trans->matrix[2][0];
  346.     tmp.y = vec->x*trans->matrix[0][1] +
  347.         vec->y*trans->matrix[1][1] + vec->z*trans->matrix[2][1];
  348.     tmp.z = vec->x*trans->matrix[0][2] +
  349.         vec->y*trans->matrix[1][2] + vec->z*trans->matrix[2][2];
  350.  
  351.     *vec = tmp;
  352. }
  353.  
  354. /*
  355.  * Transform "ray" by transforming the origin point and direction vector.
  356.  */
  357. double
  358. TransformRay(ray, trans)
  359. Ray *ray;
  360. TransInfo *trans;
  361. {
  362.     transform_point(&ray->pos, trans);
  363.     transform_vector(&ray->dir, trans);
  364.     return normalize(&ray->dir);
  365. }
  366.  
  367. /*
  368.  * Transform normal -- multiply by the transpose of the given
  369.  * matrix (which is the inverse of the desired transformation).
  370.  */
  371. TransformNormal(norm, it)
  372. Vector *norm;
  373. TransInfo *it;
  374. {
  375.     Vector onorm;
  376.  
  377.     onorm = *norm;
  378.  
  379.     norm->x = onorm.x*it->matrix[0][0] + onorm.y*it->matrix[0][1] +
  380.                 onorm.z*it->matrix[0][2];
  381.     norm->y = onorm.x*it->matrix[1][0] + onorm.y*it->matrix[1][1] +
  382.                 onorm.z*it->matrix[1][2];
  383.     norm->z = onorm.x*it->matrix[2][0] + onorm.y*it->matrix[2][1] +
  384.                 onorm.z*it->matrix[2][2];
  385. }
  386.