home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / pov22f.zip / source / matrices.c < prev    next >
C/C++ Source or Header  |  1993-07-28  |  14KB  |  513 lines
  1. /****************************************************************************
  2. *                matrices.c
  3. *
  4. *  This module contains code to manipulate 4x4 matrices.
  5. *
  6. *  from Persistence of Vision Raytracer
  7. *  Copyright 1993 Persistence of Vision Team
  8. *---------------------------------------------------------------------------
  9. *  NOTICE: This source code file is provided so that users may experiment
  10. *  with enhancements to POV-Ray and to port the software to platforms other 
  11. *  than those supported by the POV-Ray Team.  There are strict rules under
  12. *  which you are permitted to use this file.  The rules are in the file
  13. *  named POVLEGAL.DOC which should be distributed with this file. If 
  14. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  15. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  16. *  Forum.  The latest version of POV-Ray may be found there as well.
  17. *
  18. * This program is based on the popular DKB raytracer version 2.12.
  19. * DKBTrace was originally written by David K. Buck.
  20. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  21. *
  22. *****************************************************************************/
  23.  
  24. #include "frame.h"
  25. #include "vector.h"
  26. #include "povproto.h"
  27.  
  28. void MZero (result)
  29. MATRIX *result;
  30.   {
  31.   /* Initialize the matrix to the following values:
  32.    0.0   0.0   0.0   0.0
  33.    0.0   0.0   0.0   0.0
  34.    0.0   0.0   0.0   0.0
  35.    0.0   0.0   0.0   0.0
  36. */
  37.   register int i, j;
  38.  
  39.   for (i = 0 ; i < 4 ; i++)
  40.     for (j = 0 ; j < 4 ; j++)
  41.     (*result)[i][j] = 0.0;
  42.   }
  43.  
  44. void MIdentity (result)
  45. MATRIX *result;
  46.   {
  47.   /* Initialize the matrix to the following values:
  48.    1.0   0.0   0.0   0.0
  49.    0.0   1.0   0.0   0.0
  50.    0.0   0.0   1.0   0.0
  51.    0.0   0.0   0.0   1.0
  52. */
  53.   register int i, j;
  54.  
  55.   for (i = 0 ; i < 4 ; i++)
  56.     for (j = 0 ; j < 4 ; j++)
  57.     if (i == j)
  58.       (*result)[i][j] = 1.0;
  59.     else
  60.       (*result)[i][j] = 0.0;
  61.   }
  62.  
  63. void MTimes (result, matrix1, matrix2)
  64. MATRIX *result, *matrix1, *matrix2;
  65.   {
  66.   register int i, j, k;
  67.   MATRIX temp_matrix;
  68.  
  69.   for (i = 0 ; i < 4 ; i++)
  70.     for (j = 0 ; j < 4 ; j++) 
  71.     {
  72.     temp_matrix[i][j] = 0.0;
  73.     for (k = 0 ; k < 4 ; k++)
  74.       temp_matrix[i][j] += (*matrix1)[i][k] * (*matrix2)[k][j];
  75.     }
  76.  
  77.   for (i = 0 ; i < 4 ; i++)
  78.     for (j = 0 ; j < 4 ; j++)
  79.     (*result)[i][j] = temp_matrix[i][j];
  80.   }
  81.  
  82. /*  AAC - These are not used, so they are commented out to save code space...
  83.  
  84. void MAdd (result, matrix1, matrix2)
  85.    MATRIX *result, *matrix1, *matrix2;
  86.    {
  87.    register int i, j;
  88.  
  89.    for (i = 0 ; i < 4 ; i++)
  90.       for (j = 0 ; j < 4 ; j++)
  91.          (*result)[i][j] = (*matrix1)[i][j] + (*matrix2)[i][j];
  92.    }
  93.  
  94. void MSub (result, matrix1, matrix2)
  95.    MATRIX *result, *matrix1, *matrix2;
  96.    {
  97.    register int i, j;
  98.  
  99.    for (i = 0 ; i < 4 ; i++)
  100.       for (j = 0 ; j < 4 ; j++)
  101.          (*result)[i][j] = (*matrix1)[i][j] - (*matrix2)[i][j];
  102.    }
  103.  
  104. void MScale (result, matrix1, amount)
  105. MATRIX *result, *matrix1;
  106. DBL amount;
  107. {
  108.    register int i, j;
  109.  
  110.    for (i = 0 ; i < 4 ; i++)
  111.       for (j = 0 ; j < 4 ; j++)
  112.      if (amount == 1.0)
  113.         (*result)[i][j] = (*matrix1)[i][j]; * just copy *
  114.      else
  115.             (*result)[i][j] = (*matrix1)[i][j] * amount;
  116.    return;
  117. }
  118. ... up to here! */
  119.  
  120. void MTranspose (result, matrix1)
  121. MATRIX *result, *matrix1;
  122.   {
  123.   register int i, j;
  124.   MATRIX temp_matrix;
  125.  
  126.   for (i = 0 ; i < 4 ; i++)
  127.     for (j = 0 ; j < 4 ; j++)
  128.     temp_matrix[i][j] = (*matrix1)[j][i];
  129.  
  130.   for (i = 0 ; i < 4 ; i++)
  131.     for (j = 0 ; j < 4 ; j++)
  132.     (*result)[i][j] = temp_matrix[i][j];
  133.   }
  134.  
  135.  
  136. void MTransPoint (result, vector, transform)
  137. VECTOR *result, *vector;
  138. TRANSFORM *transform;
  139.   {
  140.   register int i;
  141.   DBL answer_array[4];
  142.   MATRIX *matrix;
  143.  
  144.   matrix = (MATRIX *) transform -> matrix;
  145.  
  146.   for (i = 0 ; i < 4 ; i++)
  147.     answer_array[i] = vector -> x * (*matrix)[0][i]
  148.     + vector -> y * (*matrix)[1][i]
  149.     + vector -> z * (*matrix)[2][i]
  150.     + (*matrix)[3][i];
  151.  
  152.   result -> x  = answer_array[0];
  153.   result -> y  = answer_array[1];
  154.   result -> z  = answer_array[2];
  155.   }
  156.  
  157. void MInvTransPoint (result, vector, transform)
  158. VECTOR *result, *vector;
  159. TRANSFORM *transform;
  160.   {
  161.   register int i;
  162.   DBL answer_array[4];
  163.   MATRIX *matrix;
  164.  
  165.   matrix = (MATRIX *) transform -> inverse;
  166.  
  167.   for (i = 0 ; i < 4 ; i++)
  168.     answer_array[i] = vector -> x * (*matrix)[0][i]
  169.     + vector -> y * (*matrix)[1][i]
  170.     + vector -> z * (*matrix)[2][i]
  171.     + (*matrix)[3][i];
  172.  
  173.   result -> x  = answer_array[0];
  174.   result -> y  = answer_array[1];
  175.   result -> z  = answer_array[2];
  176.   }
  177.  
  178. void MTransDirection (result, vector, transform)
  179. VECTOR *result, *vector;
  180. TRANSFORM *transform;
  181.   {
  182.   register int i;
  183.   DBL answer_array[4];
  184.   MATRIX *matrix;
  185.  
  186.   matrix = (MATRIX *) transform -> matrix;
  187.  
  188.   for (i = 0 ; i < 4 ; i++)
  189.     answer_array[i] = vector -> x * (*matrix)[0][i]
  190.     + vector -> y * (*matrix)[1][i]
  191.     + vector -> z * (*matrix)[2][i];
  192.  
  193.   result -> x  = answer_array[0];
  194.   result -> y  = answer_array[1];
  195.   result -> z  = answer_array[2];
  196.   }
  197.  
  198. void MInvTransDirection (result, vector, transform)
  199. VECTOR *result, *vector;
  200. TRANSFORM *transform;
  201.   {
  202.   register int i;
  203.   DBL answer_array[4];
  204.   MATRIX *matrix;
  205.  
  206.   matrix = (MATRIX *) transform -> inverse;
  207.  
  208.   for (i = 0 ; i < 4 ; i++)
  209.     answer_array[i] = vector -> x * (*matrix)[0][i]
  210.     + vector -> y * (*matrix)[1][i]
  211.     + vector -> z * (*matrix)[2][i];
  212.  
  213.   result -> x  = answer_array[0];
  214.   result -> y  = answer_array[1];
  215.   result -> z  = answer_array[2];
  216.   }
  217.  
  218. void MTransNormal (result, vector, transform)
  219. VECTOR *result, *vector;
  220. TRANSFORM *transform;
  221.   {
  222.   register int i;
  223.   DBL answer_array[3];
  224.   MATRIX *matrix;
  225.  
  226.   matrix = (MATRIX *) transform -> inverse;
  227.  
  228.   for (i = 0 ; i < 3 ; i++)
  229.     answer_array[i] = vector -> x * (*matrix)[i][0]
  230.     + vector -> y * (*matrix)[i][1]
  231.     + vector -> z * (*matrix)[i][2];
  232.  
  233.   result -> x  = answer_array[0];
  234.   result -> y  = answer_array[1];
  235.   result -> z  = answer_array[2];
  236.   }
  237.  
  238. void MInvTransNormal (result, vector, transform)
  239. VECTOR *result, *vector;
  240. TRANSFORM *transform;
  241.   {
  242.   register int i;
  243.   DBL answer_array[3];
  244.   MATRIX *matrix;
  245.  
  246.   matrix = (MATRIX *) transform -> matrix;
  247.  
  248.   for (i = 0 ; i < 3 ; i++)
  249.     answer_array[i] = vector -> x * (*matrix)[i][0]
  250.     + vector -> y * (*matrix)[i][1]
  251.     + vector -> z * (*matrix)[i][2];
  252.  
  253.   result -> x  = answer_array[0];
  254.   result -> y  = answer_array[1];
  255.   result -> z  = answer_array[2];
  256.   }
  257.  
  258. void Compute_Scaling_Transform (result, vector)
  259. TRANSFORM *result;
  260. VECTOR *vector;
  261.   {
  262.   MIdentity ((MATRIX *)result -> matrix);
  263.   (result -> matrix)[0][0] = vector -> x;
  264.   (result -> matrix)[1][1] = vector -> y;
  265.   (result -> matrix)[2][2] = vector -> z;
  266.  
  267.   MIdentity ((MATRIX *)result -> inverse);
  268.   (result -> inverse)[0][0] = 1.0 / vector -> x;
  269.   (result -> inverse)[1][1]= 1.0 / vector -> y;
  270.   (result -> inverse)[2][2] = 1.0 / vector -> z;
  271.   }
  272.  
  273. /* AAC - This is not used, so it's commented out...
  274.  
  275. void Compute_Inversion_Transform (result)
  276.    TRANSFORM *result;
  277.    {
  278.    MIdentity ((MATRIX *)result -> matrix);
  279.    (result -> matrix)[0][0] = -1.0;
  280.    (result -> matrix)[1][1] = -1.0;
  281.    (result -> matrix)[2][2] = -1.0;
  282.    (result -> matrix)[3][3] = -1.0;
  283.  
  284.  
  285.    (result -> inverse)[0][0] = -1.0;
  286.    (result -> inverse)[1][1] = -1.0;
  287.    (result -> inverse)[2][2] = -1.0;
  288.    (result -> inverse)[3][3] = -1.0;
  289.    }
  290. ... up to here! */
  291.  
  292. void Compute_Translation_Transform (transform, vector)
  293. TRANSFORM *transform;
  294. VECTOR *vector;
  295.   {
  296.   MIdentity ((MATRIX *)transform -> matrix);
  297.   (transform -> matrix)[3][0] = vector -> x;
  298.   (transform -> matrix)[3][1] = vector -> y;
  299.   (transform -> matrix)[3][2] = vector -> z;
  300.  
  301.   MIdentity ((MATRIX *)transform -> inverse);
  302.   (transform -> inverse)[3][0] = 0.0 - vector -> x;
  303.   (transform -> inverse)[3][1] = 0.0 - vector -> y;
  304.   (transform -> inverse)[3][2] = 0.0 - vector -> z;
  305.   }
  306.  
  307. void Compute_Rotation_Transform (transform, vector)
  308. TRANSFORM *transform;
  309. VECTOR *vector;
  310.   {
  311.   MATRIX Matrix;
  312.   VECTOR Radian_Vector;
  313.   register DBL cosx, cosy, cosz, sinx, siny, sinz;
  314.  
  315.   VScale (Radian_Vector, *vector, M_PI/180.0);
  316.   MIdentity ((MATRIX *)transform -> matrix);
  317.   cosx = cos (Radian_Vector.x);
  318.   sinx = sin (Radian_Vector.x);
  319.   cosy = cos (Radian_Vector.y);
  320.   siny = sin (Radian_Vector.y);
  321.   cosz = cos (Radian_Vector.z);
  322.   sinz = sin (Radian_Vector.z);
  323.  
  324.   (transform -> matrix) [1][1] = cosx;
  325.   (transform -> matrix) [2][2] = cosx;
  326.   (transform -> matrix) [1][2] = sinx;
  327.   (transform -> matrix) [2][1] = 0.0 - sinx;
  328.   MTranspose ((MATRIX *)transform -> inverse, (MATRIX *)transform -> matrix);
  329.  
  330.   MIdentity ((MATRIX *)Matrix);
  331.   Matrix [0][0] = cosy;
  332.   Matrix [2][2] = cosy;
  333.   Matrix [0][2] = 0.0 - siny;
  334.   Matrix [2][0] = siny;
  335.   MTimes ((MATRIX *)transform -> matrix, (MATRIX *)transform -> matrix, (MATRIX *)Matrix);
  336.   MTranspose ((MATRIX *)Matrix, (MATRIX *)Matrix);
  337.   MTimes ((MATRIX *)transform -> inverse, (MATRIX *)Matrix, (MATRIX *)transform -> inverse);
  338.  
  339.   MIdentity ((MATRIX *)Matrix);
  340.   Matrix [0][0] = cosz;
  341.   Matrix [1][1] = cosz;
  342.   Matrix [0][1] = sinz;
  343.   Matrix [1][0] = 0.0 - sinz;
  344.   MTimes ((MATRIX *)transform -> matrix, (MATRIX *)transform -> matrix, (MATRIX *)Matrix);
  345.   MTranspose ((MATRIX *)Matrix, (MATRIX *)Matrix);
  346.   MTimes ((MATRIX *)transform -> inverse, (MATRIX *)Matrix, (MATRIX *)transform -> inverse);
  347.   }
  348.  
  349. /* AAC - This is not used so it's commented out...
  350.  
  351. void Compute_Look_At_Transform (result, Look_At, Up, Right)
  352.    TRANSFORM *result;
  353.    VECTOR *Look_At, *Up, *Right;
  354.    {
  355.    MIdentity ((MATRIX *)result -> inverse);
  356.    (result -> matrix)[0][0] = Right->x;
  357.    (result -> matrix)[0][1] = Right->y;
  358.    (result -> matrix)[0][2] = Right->z;
  359.    (result -> matrix)[1][0] = Up->x;
  360.    (result -> matrix)[1][1] = Up->y;
  361.    (result -> matrix)[1][2] = Up->z;
  362.    (result -> matrix)[2][0] = Look_At->x;
  363.    (result -> matrix)[2][1] = Look_At->y;
  364.    (result -> matrix)[2][2] = Look_At->z;
  365.  
  366.    MIdentity ((MATRIX *)result -> matrix);
  367.    MTranspose ((MATRIX *)result -> matrix, (MATRIX *)result -> inverse);   
  368.    }
  369.  
  370. ... up to here! */
  371.  
  372. void Compose_Transforms (Original_Transform, New_Transform)
  373. TRANSFORM *Original_Transform, *New_Transform;
  374.   {
  375.   MTimes ((MATRIX *)Original_Transform -> matrix,
  376.     (MATRIX *)Original_Transform -> matrix,
  377.     (MATRIX *)New_Transform -> matrix);
  378.  
  379.   MTimes ((MATRIX *)Original_Transform -> inverse,
  380.     (MATRIX *)New_Transform -> inverse,
  381.     (MATRIX *)Original_Transform -> inverse);
  382.   }
  383.  
  384. /* Rotation about an arbitrary axis - formula from:
  385.       "Computational Geometry for Design and Manufacture",
  386.       Faux & Pratt
  387.    Note that the angles for this transform are specified in radians.
  388. */
  389. void Compute_Axis_Rotation_Transform (transform, V, angle)
  390. TRANSFORM *transform;
  391. VECTOR *V;
  392. DBL angle;
  393.   {
  394.   DBL l, cosx, sinx;
  395.  
  396.   VLength(l, *V);
  397.   VInverseScaleEq(*V, l);
  398.  
  399.   MIdentity(&transform->matrix);
  400.   cosx = cos(angle);
  401.   sinx = sin(angle);
  402.   transform->matrix[0][0] = V->x * V->x + cosx * (1.0 - V->x * V->x);
  403.   transform->matrix[0][1] = V->x * V->y * (1.0 - cosx) + V->z * sinx;
  404.   transform->matrix[0][2] = V->x * V->z * (1.0 - cosx) - V->y * sinx;
  405.   transform->matrix[1][0] = V->x * V->y * (1.0 - cosx) - V->z * sinx;
  406.   transform->matrix[1][1] = V->y * V->y + cosx * (1.0 - V->y * V->y);
  407.   transform->matrix[1][2] = V->y * V->z * (1.0 - cosx) + V->x * sinx;
  408.   transform->matrix[2][0] = V->x * V->z * (1.0 - cosx) + V->y * sinx;
  409.   transform->matrix[2][1] = V->y * V->z * (1.0 - cosx) - V->x * sinx;
  410.   transform->matrix[2][2] = V->z * V->z + cosx * (1.0 - V->z * V->z);
  411.   MTranspose(&transform->inverse, &transform->matrix);   
  412.   }
  413.  
  414. /* Given a point and a direction and a radius, find the transform
  415.    that brings these into a canonical coordinate system */
  416. void
  417. Compute_Coordinate_Transform(trans, origin, up, radius, length)
  418. TRANSFORM *trans;
  419. VECTOR *origin;
  420. VECTOR *up;
  421. DBL radius;
  422. DBL length;
  423.   {
  424.   TRANSFORM trans2;
  425.   VECTOR tmpv;
  426.  
  427.   Make_Vector(&tmpv, radius, radius, length);
  428.   Compute_Scaling_Transform(trans, &tmpv);
  429.   if (fabs(up->z) == 1.0)
  430.     Make_Vector(&tmpv, 1.0, 0.0, 0.0)
  431. else
  432.   Make_Vector(&tmpv, -up->y, up->x, 0.0)
  433.     Compute_Axis_Rotation_Transform(&trans2, &tmpv, acos(up->z));
  434. Compose_Transforms(trans, &trans2);
  435. Compute_Translation_Transform(&trans2, origin);
  436. Compose_Transforms(trans, &trans2);
  437. }
  438.  
  439. TRANSFORM *Create_Transform()
  440. {
  441. TRANSFORM *New;
  442.  
  443. if ((New = (TRANSFORM *) malloc (sizeof (TRANSFORM))) == NULL)
  444. MAError ("transform");
  445.  
  446. MIdentity ((MATRIX *) &(New -> matrix[0][0]));
  447. MIdentity ((MATRIX *) &(New -> inverse[0][0]));
  448. return (New);
  449. }
  450.  
  451. TRANSFORM *Copy_Transform (Old)
  452. TRANSFORM *Old;
  453. {
  454. TRANSFORM *New;
  455. if (Old != NULL)
  456.   {
  457.   New  = Create_Transform ();
  458.   *New = *Old;
  459.   }
  460. else New = NULL;
  461. return (New);
  462. }
  463.  
  464. VECTOR *Create_Vector ()
  465. {
  466. VECTOR *New;
  467.  
  468. if ((New = (VECTOR *) malloc (sizeof (VECTOR))) == NULL)
  469. MAError ("vector");
  470.  
  471. Make_Vector (New, 0.0, 0.0, 0.0);
  472.  
  473. return (New);
  474. }
  475.  
  476. VECTOR *Copy_Vector (Old)
  477. VECTOR *Old;
  478. {
  479. VECTOR *New;
  480. if (Old != NULL)
  481.   {
  482.   New  = Create_Vector ();
  483.   *New = *Old;
  484.   }
  485. else New = NULL;
  486. return (New);
  487. }
  488.  
  489. DBL *Create_Float ()
  490. {
  491. DBL *New_Float;
  492.  
  493. if ((New_Float = (DBL *) malloc (sizeof (DBL))) == NULL)
  494. MAError ("float");
  495.  
  496. *New_Float = 0.0;
  497. return (New_Float);
  498. }
  499.  
  500. DBL *Copy_Float (Old)
  501. DBL *Old;
  502. {
  503. DBL *New;
  504. if (Old)
  505.   {
  506.   New  = Create_Float ();
  507.   *New = *Old;
  508.   }
  509. else New = NULL;
  510. return (New);
  511. }
  512.  
  513.