home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / graf / fract3.zip / 3D.C next >
Text File  |  1989-10-24  |  11KB  |  378 lines
  1. /*
  2. A word about the 3D library. Even though this library supports 
  3. three dimensions, the matrices are 4x4 for the following reason. 
  4. With normal 3 dimensional vectors, translation is an ADDITION, 
  5. and rotation is a MULTIPLICATION. A vector {x,y,z} is represented 
  6. as a 4-tuple {x,y,z,1}. It is then possible to define a 4x4 
  7. matrix such that multiplying the vector by the matrix translates 
  8. the vector. This allows combinations of translation and rotation 
  9. to be obtained in a single matrix by multiplying a translation 
  10. matrix and a rotation matrix together. Note that in the code, 
  11. vectors have three components; since the fourth component is 
  12. always 1, that value is not included in the vector variable to 
  13. save space, but the routines make use of the fourth component 
  14. (see vec_mult()). Similarly, the fourth column of EVERY matrix is
  15. always
  16.          0
  17.          0
  18.          0
  19.          1
  20. but currently the C version of a matrix includes this even though 
  21. it could be left out of the data structure and assumed in the 
  22. routines. Vectors are ROW vectors, and are always multiplied with 
  23. matrices FROM THE LEFT (e.g. vector*matrix). Also note the order 
  24. of indices of a matrix is matrix[row][column], and in usual C 
  25. fashion, numbering starts with 0.
  26.  
  27. TRANSLATION MATRIX =  1     0     0     0
  28.                       0     1     0     0
  29.                       0     0     1     0
  30.                       Tx    Ty    Tz    1
  31.  
  32. SCALE MATRIX =        Sx    0     0     0
  33.                       0     Sy    0     0
  34.                       0     0     Sz    0
  35.                       0     0     0     1
  36.  
  37. Rotation about x axis i degrees:
  38. ROTX(Θ) =             1     0     0     0
  39.                       0   cosΘ  sinΘ    0
  40.                       0  -sinΘ  cosΘ    0
  41.                       0     0     0     1
  42.  
  43. Rotation about y axis i degrees:
  44. ROTY(Θ) =           cosΘ    0  -sinΘ    0
  45.                       0     1     0     0
  46.                     sinΘ    0   cosΘ    0
  47.                       0     0     0     1
  48.  
  49. Rotation about z axis i degrees:
  50. ROTZ(Θ) =           cosΘ  sinΘ    0     0
  51.                    -sinΘ  cosΘ    0     0
  52.                       0     0     1     0 
  53.                       0     0     0     1
  54.  
  55.                       --  Tim Wegner  April 22, 1989
  56. */
  57.  
  58. #include <stdio.h>
  59. #include <math.h>
  60. #include <float.h>
  61. #include "fractint.h"
  62. extern int bad_value;
  63.  
  64. /* initialize a matrix and set to identity matrix 
  65.    (all 0's, 1's on diagonal) */
  66. void identity(MATRIX m)
  67. {  
  68.    int i,j;
  69.    for(i=0;i<CMAX;i++)
  70.    for(j=0;j<RMAX;j++)
  71.       if(i==j)
  72.          m[j][i] = 1.0;
  73.       else
  74.           m[j][i] = 0.0;
  75. }
  76.  
  77. /* Multiply two matrices */
  78. void mat_mul(MATRIX mat1, MATRIX mat2, MATRIX mat3)
  79. {
  80.      /* result stored in MATRIX new to avoid problems
  81.         in case parameter mat3 == mat2 or mat 1 */
  82.      MATRIX new;
  83.      int i,j;
  84.      for(i=0;i<4;i++)
  85.      for(j=0;j<4;j++)
  86.         new[j][i] =  mat1[j][0]*mat2[0][i]+
  87.                      mat1[j][1]*mat2[1][i]+
  88.                      mat1[j][2]*mat2[2][i]+
  89.                      mat1[j][3]*mat2[3][i];
  90.                      
  91.      /* should replace this with memcpy */
  92.      for(i=0;i<4;i++)
  93.      for(j=0;j<4;j++)
  94.         mat3[j][i] =  new[j][i];
  95.  
  96. }
  97.  
  98. /* multiply a matrix by a scalar */
  99. void scale (double sx, double sy, double sz, MATRIX m)
  100. {
  101.    MATRIX scale;
  102.    identity(scale);
  103.    scale[0][0] = sx;
  104.    scale[1][1] = sy;
  105.    scale[2][2] = sz;
  106.    mat_mul(m,scale,m);
  107. }
  108.  
  109. /* rotate about X axis  */
  110. void xrot (double theta, MATRIX m)
  111. {
  112.    MATRIX rot;
  113.    double sintheta,costheta;
  114.    sintheta = sin(theta);
  115.    costheta = cos(theta);
  116.    identity(rot);
  117.    rot[1][1] = costheta;
  118.    rot[1][2] = -sintheta;
  119.    rot[2][1] = sintheta;
  120.    rot[2][2] = costheta;
  121.    mat_mul(m,rot,m);
  122. }
  123.  
  124. /* rotate about Y axis  */
  125. void yrot (double theta, MATRIX m)
  126. {
  127.    MATRIX rot;
  128.    double sintheta,costheta;
  129.    sintheta = sin(theta);
  130.    costheta = cos(theta);
  131.    identity(rot);
  132.    rot[0][0] = costheta;
  133.    rot[0][2] = sintheta;
  134.    rot[2][0] = -sintheta;
  135.    rot[2][2] = costheta;
  136.    mat_mul(m,rot,m);
  137. }
  138.  
  139. /* rotate about Z axis  */
  140. void zrot (double theta, MATRIX m)
  141. {
  142.    MATRIX rot;
  143.    double sintheta,costheta;
  144.    sintheta = sin(theta);
  145.    costheta = cos(theta);
  146.    identity(rot);
  147.    rot[0][0] = costheta;
  148.    rot[0][1] = -sintheta;
  149.    rot[1][0] = sintheta;
  150.    rot[1][1] = costheta;
  151.    mat_mul(m,rot,m);
  152. }
  153.  
  154. /* translate  */
  155. void trans (double tx, double ty, double tz, MATRIX m)
  156. {
  157.    MATRIX trans;
  158.    identity(trans);
  159.    trans[3][0] = tx;
  160.    trans[3][1] = ty;
  161.    trans[3][2] = tz;
  162.    mat_mul(m,trans,m);
  163. }
  164.  
  165. /* cross product  - useful because cross is perpendicular to v and w */
  166. int cross_product (VECTOR v, VECTOR w, VECTOR cross)
  167. {
  168.    VECTOR tmp;
  169.    tmp[0] =  v[1]*w[2] - w[1]*v[2];
  170.    tmp[1] =  w[0]*v[2] - v[0]*w[2];
  171.    tmp[2] =  v[0]*w[1] - w[0]*v[1];
  172.    cross[0] = tmp[0];
  173.    cross[1] = tmp[1];
  174.    cross[2] = tmp[2];
  175.    return(0);
  176. }
  177. /* cross product integer arguments (not fudged) */
  178. int icross_product (IVECTOR v, IVECTOR w, IVECTOR cross)
  179. {
  180.    IVECTOR tmp;
  181.    tmp[0] =  v[1]*w[2] - w[1]*v[2];
  182.    tmp[1] =  w[0]*v[2] - v[0]*w[2];
  183.    tmp[2] =  v[0]*w[1] - w[0]*v[1];
  184.    cross[0] = tmp[0];
  185.    cross[1] = tmp[1];
  186.    cross[2] = tmp[2];
  187.    return(0);
  188. }
  189.  
  190. /* normalize a vector to length 1 */
  191. normalize_vector(VECTOR v)
  192. {
  193.     double vlength;
  194.     vlength = dot_product(v,v);
  195.  
  196.     /* bailout if zero vlength */
  197.     if(vlength < FLT_MIN || vlength > FLT_MAX)
  198.        return(-1);
  199.     vlength = sqrt(vlength);
  200.     if(vlength < FLT_MIN)
  201.        return(-1);
  202.        
  203.     v[0] /= vlength;
  204.     v[1] /= vlength;
  205.     v[2] /= vlength;
  206.     return(0);
  207. }
  208.  
  209. /* multiply source vector s by matrix m, result in target t */
  210. /* used to apply transformations to a vector */
  211. int vmult(s,m,t)
  212. VECTOR s,t;
  213. MATRIX m;
  214. {
  215.    VECTOR tmp;
  216.    int i,j;
  217.    for(j=0;j<CMAX-1;j++)
  218.    {
  219.       tmp[j] = 0.0;
  220.       for(i=0;i<RMAX-1;i++)
  221.          tmp[j] += s[i]*m[i][j];
  222.       /* vector is really four dimensional with last component always 1 */   
  223.       tmp[j] += m[3][j];
  224.    }
  225.    /* set target = tmp. Necessary to use tmp in case source = target */
  226.    for(i=0;i<DIM;i++)
  227.       t[i] = tmp[i];
  228.    return(0);
  229. }
  230.  
  231. /* perspective projection of vector v with respect to viewpont vector view */
  232. perspective(VECTOR v)
  233. {
  234.    extern VECTOR view;
  235.    double denom;
  236.    denom = view[2] - v[2];
  237.  
  238.    if(denom >= 0.0)
  239.    {
  240.       v[0] = bad_value;   /* clipping will catch these values */
  241.       v[1] = bad_value;   /* so they won't plot values BEHIND viewer */
  242.       v[2] = bad_value;
  243.       return(-1);
  244.    }   
  245.    v[0] = (v[0]*view[2] - view[0]*v[2])/denom;
  246.    v[1] = (v[1]*view[2] - view[1]*v[2])/denom;
  247.  
  248.    /* calculation of z if needed later */
  249.    /* v[2] =  v[2]/denom;*/
  250.    return(0);
  251. }
  252.  
  253. /* long version of vmult and perspective combined for speed */
  254. longvmultpersp(s, m, t0, t, lview, bitshift)
  255. LVECTOR s;       /* source vector */
  256. LMATRIX m;       /* transformation matrix */
  257. LVECTOR t0;      /* after transformation, before persp */
  258. LVECTOR t;       /* target vector */
  259. LVECTOR lview;   /* perspective viewer coordinates */
  260. int bitshift;    /* fixed point conversion bitshift */
  261. {
  262.    LVECTOR tmp;
  263.    int i,j, k;
  264.  
  265.    k = CMAX-1;            /* shorten the math if non-perspective and non-illum */
  266.    if (lview[2] == 0 && t0[0] == 0) k--;
  267.  
  268.    for(j=0;j<k;j++)
  269.    {
  270.       tmp[j] = 0;
  271.       for(i=0;i<RMAX-1;i++)
  272.          tmp[j] += multiply(s[i],m[i][j],bitshift);
  273.       /* vector is really four dimensional with last component always 1 */   
  274.       tmp[j] += m[3][j];
  275.    }
  276.    if(t0[0]) /* first component of  t0 used as flag */
  277.    {
  278.       /* faster than for loop, if less general */
  279.       t0[0] = tmp[0];
  280.       t0[1] = tmp[1];
  281.       t0[2] = tmp[2];
  282.    }
  283.    if (lview[2] != 0)        /* perspective 3D */
  284.    {
  285.  
  286.       LVECTOR tmpview;
  287.       long denom;
  288.       
  289.       denom = lview[2] - tmp[2];
  290.       if (denom >= 0)         /* bail out if point is "behind" us */
  291.       {
  292.            tmp[0] = bad_value;
  293.            tmp[0] = tmp[0]<<bitshift;
  294.            tmp[1] = tmp[0];
  295.            tmp[2] = tmp[0];
  296.            return(-1);
  297.       }
  298.       
  299.       /* doing math in this order helps prevent overflow */ 
  300.       tmpview[0] = divide(lview[0],denom,bitshift);
  301.       tmpview[1] = divide(lview[1],denom,bitshift);
  302.       tmpview[2] = divide(lview[2],denom,bitshift);
  303.       
  304.       tmp[0] = multiply(tmp[0], tmpview[2], bitshift) -
  305.                multiply(tmpview[0], tmp[2], bitshift);
  306.    
  307.       tmp[1] = multiply(tmp[1], tmpview[2], bitshift) -
  308.                multiply(tmpview[1], tmp[2], bitshift);
  309.    
  310.       /* z coordinate if needed           */
  311.       /* tmp[2] = divide(lview[2],denom);  */
  312.    }
  313.  
  314.    /* set target = tmp. Necessary to use tmp in case source = target */
  315.    /* faster than for loop, if less general */
  316.    t[0] = tmp[0];
  317.    t[1] = tmp[1];
  318.    t[2] = tmp[2];
  319.    return(0);
  320. }
  321.  
  322. /* Long version of perspective. Because of use of fixed point math, there
  323.    is danger of overflow and underflow */
  324. longpersp(LVECTOR lv, LVECTOR lview, int bitshift)
  325. {
  326.    LVECTOR tmpview;
  327.    long denom;
  328.    
  329.    denom = lview[2] - lv[2];
  330.    if (denom == 0)         /* bail out if point is "behind" us */
  331.    {
  332.         lv[0] = bad_value;
  333.         lv[0] = lv[0]<<bitshift;
  334.         lv[1] = lv[0];
  335.         lv[2] = lv[0];
  336.         return(-1);
  337.    }
  338.  
  339.    /* doing math in this order helps prevent overflow */ 
  340.    tmpview[0] = divide(lview[0],denom,bitshift);
  341.    tmpview[1] = divide(lview[1],denom,bitshift);
  342.    tmpview[2] = divide(lview[2],denom,bitshift);
  343.    
  344.    lv[0] = multiply(lv[0], tmpview[2], bitshift) -
  345.            multiply(tmpview[0], lv[2], bitshift);
  346.  
  347.    lv[1] = multiply(lv[1], tmpview[2], bitshift) -
  348.            multiply(tmpview[1], lv[2], bitshift);
  349.  
  350.    /* z coordinate if needed           */
  351.    /* lv[2] = divide(lview[2],denom);  */
  352.    return(0);
  353. }
  354.  
  355. int longvmult(LVECTOR s,LMATRIX m,LVECTOR t,int bitshift)
  356. {
  357.    LVECTOR tmp;
  358.    int i,j, k;
  359.  
  360.    k = CMAX-1;
  361.  
  362.    for(j=0;j<k;j++)
  363.    {
  364.       tmp[j] = 0;
  365.       for(i=0;i<RMAX-1;i++)
  366.          tmp[j] += multiply(s[i],m[i][j],bitshift);
  367.       /* vector is really four dimensional with last component always 1 */   
  368.       tmp[j] += m[3][j];
  369.    }
  370.  
  371.    /* set target = tmp. Necessary to use tmp in case source = target */
  372.    /* faster than for loop, if less general */
  373.    t[0] = tmp[0];
  374.    t[1] = tmp[1];
  375.    t[2] = tmp[2];
  376.    return(0);
  377. }
  378.