home *** CD-ROM | disk | FTP | other *** search
/ Media Share 9 / MEDIASHARE_09.ISO / progmisc / frasrc18.zip / 3D.C next >
C/C++ Source or Header  |  1992-07-10  |  10KB  |  399 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(i) =        1     0     0      0
  39.               0   cosi  sini    0
  40.               0  -sini  cosi    0
  41.               0     0      0    1
  42.  
  43. Rotation about y axis i degrees:
  44. ROTY(i) =          cosi    0  -sini    0
  45.               0     1      0    0
  46.             sini    0   cosi    0
  47.               0     0      0    1
  48.  
  49. Rotation about z axis i degrees:
  50. ROTZ(i) =          cosi  sini    0     0
  51.            -sini  cosi    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 <float.h>
  60. #include <string.h>
  61. #include "fractint.h"
  62. #include "prototyp.h"
  63. extern int overflow;
  64. extern int bad_value;
  65.  
  66. /* initialize a matrix and set to identity matrix
  67.    (all 0's, 1's on diagonal) */
  68. void identity(MATRIX m)
  69. {
  70.    int i,j;
  71.    for(i=0;i<CMAX;i++)
  72.    for(j=0;j<RMAX;j++)
  73.       if(i==j)
  74.      m[j][i] = 1.0;
  75.       else
  76.       m[j][i] = 0.0;
  77. }
  78.  
  79. /* Multiply two matrices */
  80. void mat_mul(MATRIX mat1, MATRIX mat2, MATRIX mat3)
  81. {
  82.      /* result stored in MATRIX new to avoid problems
  83.     in case parameter mat3 == mat2 or mat 1 */
  84.      MATRIX new;
  85.      int i,j;
  86.      for(i=0;i<4;i++)
  87.      for(j=0;j<4;j++)
  88.     new[j][i] =  mat1[j][0]*mat2[0][i]+
  89.              mat1[j][1]*mat2[1][i]+
  90.              mat1[j][2]*mat2[2][i]+
  91.              mat1[j][3]*mat2[3][i];
  92.      memcpy(mat3,new,sizeof(new));
  93. }
  94.  
  95. /* multiply a matrix by a scalar */
  96. void scale (double sx, double sy, double sz, MATRIX m)
  97. {
  98.    MATRIX scale;
  99.    identity(scale);
  100.    scale[0][0] = sx;
  101.    scale[1][1] = sy;
  102.    scale[2][2] = sz;
  103.    mat_mul(m,scale,m);
  104. }
  105.  
  106. /* rotate about X axis    */
  107. void xrot (double theta, MATRIX m)
  108. {
  109.    MATRIX rot;
  110.    double sintheta,costheta;
  111.    sintheta = sin(theta);
  112.    costheta = cos(theta);
  113.    identity(rot);
  114.    rot[1][1] = costheta;
  115.    rot[1][2] = -sintheta;
  116.    rot[2][1] = sintheta;
  117.    rot[2][2] = costheta;
  118.    mat_mul(m,rot,m);
  119. }
  120.  
  121. /* rotate about Y axis    */
  122. void yrot (double theta, MATRIX m)
  123. {
  124.    MATRIX rot;
  125.    double sintheta,costheta;
  126.    sintheta = sin(theta);
  127.    costheta = cos(theta);
  128.    identity(rot);
  129.    rot[0][0] = costheta;
  130.    rot[0][2] = sintheta;
  131.    rot[2][0] = -sintheta;
  132.    rot[2][2] = costheta;
  133.    mat_mul(m,rot,m);
  134. }
  135.  
  136. /* rotate about Z axis    */
  137. void zrot (double theta, MATRIX m)
  138. {
  139.    MATRIX rot;
  140.    double sintheta,costheta;
  141.    sintheta = sin(theta);
  142.    costheta = cos(theta);
  143.    identity(rot);
  144.    rot[0][0] = costheta;
  145.    rot[0][1] = -sintheta;
  146.    rot[1][0] = sintheta;
  147.    rot[1][1] = costheta;
  148.    mat_mul(m,rot,m);
  149. }
  150.  
  151. /* translate  */
  152. void trans (double tx, double ty, double tz, MATRIX m)
  153. {
  154.    MATRIX trans;
  155.    identity(trans);
  156.    trans[3][0] = tx;
  157.    trans[3][1] = ty;
  158.    trans[3][2] = tz;
  159.    mat_mul(m,trans,m);
  160. }
  161.  
  162. /* cross product  - useful because cross is perpendicular to v and w */
  163. int cross_product (VECTOR v, VECTOR w, VECTOR cross)
  164. {
  165.    VECTOR tmp;
  166.    tmp[0] =  v[1]*w[2] - w[1]*v[2];
  167.    tmp[1] =  w[0]*v[2] - v[0]*w[2];
  168.    tmp[2] =  v[0]*w[1] - w[0]*v[1];
  169.    cross[0] = tmp[0];
  170.    cross[1] = tmp[1];
  171.    cross[2] = tmp[2];
  172.    return(0);
  173. }
  174.  
  175. /* cross product integer arguments (not fudged) */
  176. /*** pb, unused
  177. int icross_product (IVECTOR v, IVECTOR w, IVECTOR cross)
  178. {
  179.    IVECTOR tmp;
  180.    tmp[0] =  v[1]*w[2] - w[1]*v[2];
  181.    tmp[1] =  w[0]*v[2] - v[0]*w[2];
  182.    tmp[2] =  v[0]*w[1] - w[0]*v[1];
  183.    cross[0] = tmp[0];
  184.    cross[1] = tmp[1];
  185.    cross[2] = tmp[2];
  186.    return(0);
  187. }
  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.    memcpy(t,tmp,sizeof(tmp));
  227.    return(0);
  228. }
  229.  
  230. /* multiply vector s by matrix m, result in s */
  231. /* use with a function pointer in line3d.c */
  232. /* must coordinate calling conventions with */
  233. /* mult_vec_iit in general.asm */
  234. void mult_vec_c(s)
  235. VECTOR s;
  236. {
  237.    extern MATRIX m;
  238.    VECTOR tmp;
  239.    int i,j;
  240.    for(j=0;j<CMAX-1;j++)
  241.    {
  242.       tmp[j] = 0.0;
  243.       for(i=0;i<RMAX-1;i++)
  244.      tmp[j] += s[i]*m[i][j];
  245.       /* vector is really four dimensional with last component always 1 */
  246.       tmp[j] += m[3][j];
  247.    }
  248.    /* set target = tmp. Necessary to use tmp in case source = target */
  249.    memcpy(s,tmp,sizeof(tmp));
  250. }
  251.  
  252. /* perspective projection of vector v with respect to viewpont vector view */
  253. perspective(VECTOR v)
  254. {
  255.    extern VECTOR view;
  256.    double denom;
  257.    denom = view[2] - v[2];
  258.  
  259.    if(denom >= 0.0)
  260.    {
  261.       v[0] = bad_value;   /* clipping will catch these values */
  262.       v[1] = bad_value;   /* so they won't plot values BEHIND viewer */
  263.       v[2] = bad_value;
  264.       return(-1);
  265.    }
  266.    v[0] = (v[0]*view[2] - view[0]*v[2])/denom;
  267.    v[1] = (v[1]*view[2] - view[1]*v[2])/denom;
  268.  
  269.    /* calculation of z if needed later */
  270.    /* v[2] =  v[2]/denom;*/
  271.    return(0);
  272. }
  273.  
  274. /* long version of vmult and perspective combined for speed */
  275. longvmultpersp(s, m, t0, t, lview, bitshift)
  276. LVECTOR s;     /* source vector */
  277. LMATRIX m;     /* transformation matrix */
  278. LVECTOR t0;     /* after transformation, before persp */
  279. LVECTOR t;     /* target vector */
  280. LVECTOR lview;     /* perspective viewer coordinates */
  281. int bitshift;     /* fixed point conversion bitshift */
  282. {
  283.    LVECTOR tmp;
  284.    int i,j, k;
  285.    overflow = 0;
  286.    k = CMAX-1;            /* shorten the math if non-perspective and non-illum */
  287.    if (lview[2] == 0 && t0[0] == 0) k--;
  288.  
  289.    for(j=0;j<k;j++)
  290.    {
  291.       tmp[j] = 0;
  292.       for(i=0;i<RMAX-1;i++)
  293.      tmp[j] += multiply(s[i],m[i][j],bitshift);
  294.       /* vector is really four dimensional with last component always 1 */
  295.       tmp[j] += m[3][j];
  296.    }
  297.    if(t0[0]) /* first component of  t0 used as flag */
  298.    {
  299.       /* faster than for loop, if less general */
  300.       t0[0] = tmp[0];
  301.       t0[1] = tmp[1];
  302.       t0[2] = tmp[2];
  303.    }
  304.    if (lview[2] != 0)        /* perspective 3D */
  305.    {
  306.  
  307.       LVECTOR tmpview;
  308.       long denom;
  309.  
  310.       denom = lview[2] - tmp[2];
  311.       if (denom >= 0)        /* bail out if point is "behind" us */
  312.       {
  313.        t[0] = bad_value;
  314.        t[0] = t[0]<<bitshift;
  315.        t[1] = t[0];
  316.        t[2] = t[0];
  317.        return(-1);
  318.       }
  319.  
  320.       /* doing math in this order helps prevent overflow */
  321.       tmpview[0] = divide(lview[0],denom,bitshift);
  322.       tmpview[1] = divide(lview[1],denom,bitshift);
  323.       tmpview[2] = divide(lview[2],denom,bitshift);
  324.  
  325.       tmp[0] = multiply(tmp[0], tmpview[2], bitshift) -
  326.            multiply(tmpview[0], tmp[2], bitshift);
  327.  
  328.       tmp[1] = multiply(tmp[1], tmpview[2], bitshift) -
  329.            multiply(tmpview[1], tmp[2], bitshift);
  330.  
  331.       /* z coordinate if needed       */
  332.       /* tmp[2] = divide(lview[2],denom);  */
  333.    }
  334.  
  335.    /* set target = tmp. Necessary to use tmp in case source = target */
  336.    /* faster than for loop, if less general */
  337.    t[0] = tmp[0];
  338.    t[1] = tmp[1];
  339.    t[2] = tmp[2];
  340.    return(overflow);
  341. }
  342.  
  343. /* Long version of perspective. Because of use of fixed point math, there
  344.    is danger of overflow and underflow */
  345. longpersp(LVECTOR lv, LVECTOR lview, int bitshift)
  346. {
  347.    LVECTOR tmpview;
  348.    long denom;
  349.    overflow = 0;
  350.    denom = lview[2] - lv[2];
  351.    if (denom >= 0)        /* bail out if point is "behind" us */
  352.    {
  353.     lv[0] = bad_value;
  354.     lv[0] = lv[0]<<bitshift;
  355.     lv[1] = lv[0];
  356.     lv[2] = lv[0];
  357.     return(-1);
  358.    }
  359.  
  360.    /* doing math in this order helps prevent overflow */
  361.    tmpview[0] = divide(lview[0],denom,bitshift);
  362.    tmpview[1] = divide(lview[1],denom,bitshift);
  363.    tmpview[2] = divide(lview[2],denom,bitshift);
  364.  
  365.    lv[0] = multiply(lv[0], tmpview[2], bitshift) -
  366.        multiply(tmpview[0], lv[2], bitshift);
  367.  
  368.    lv[1] = multiply(lv[1], tmpview[2], bitshift) -
  369.        multiply(tmpview[1], lv[2], bitshift);
  370.  
  371.    /* z coordinate if needed           */
  372.    /* lv[2] = divide(lview[2],denom);  */
  373.    return(overflow);
  374. }
  375.  
  376. int longvmult(LVECTOR s,LMATRIX m,LVECTOR t,int bitshift)
  377. {
  378.    LVECTOR tmp;
  379.    int i,j, k;
  380.    overflow = 0;
  381.    k = CMAX-1;
  382.  
  383.    for(j=0;j<k;j++)
  384.    {
  385.       tmp[j] = 0;
  386.       for(i=0;i<RMAX-1;i++)
  387.      tmp[j] += multiply(s[i],m[i][j],bitshift);
  388.       /* vector is really four dimensional with last component always 1 */
  389.       tmp[j] += m[3][j];
  390.    }
  391.  
  392.    /* set target = tmp. Necessary to use tmp in case source = target */
  393.    /* faster than for loop, if less general */
  394.    t[0] = tmp[0];
  395.    t[1] = tmp[1];
  396.    t[2] = tmp[2];
  397.    return(overflow);
  398. }
  399.