home *** CD-ROM | disk | FTP | other *** search
/ 3D Color Clip Art / 3D_Color_Clip_Art.bin / 3d-progs / rad386 / utils.c < prev    next >
C/C++ Source or Header  |  1996-11-15  |  13KB  |  502 lines
  1. /*
  2.     Utility Routines Used By the Previewer.
  3.     V3Dot, V3Cross, V3Normalize, V3Length,ViewTransform,MatMul taken from GraphicsGems.
  4. -----
  5. Sumant
  6. 15 May 1991
  7. */
  8. #include <stdio.h>
  9. #include <math.h>
  10. #include "preview.h"
  11. extern int view_type;
  12. extern int oneview_flag;
  13. extern float eyepoint_distance;
  14. static Matrix ModelMatrix={
  15.     {1,0,0,0},
  16.     {0,1,0,0},
  17.     {0,0,1,0},
  18.     {0,0,0,1}
  19. };
  20. static double V3SquaredLength(a)
  21. Vector3 *a;
  22. {
  23.         return((a->x * a->x)+(a->y * a->y)+(a->z * a->z));
  24.         }
  25.  
  26. static double V3Dot(a, b)
  27. Vector3 *a, *b;
  28. {
  29.         return((a->x*b->x)+(a->y*b->y)+(a->z*b->z));
  30.         }
  31. static Vector3 *V3Cross(a, b, c)
  32. Vector3 *a, *b, *c;
  33. {
  34.         c->x = (a->y*b->z) - (a->z*b->y);
  35.         c->y = (a->z*b->x) - (a->x*b->z);
  36.         c->z = (a->x*b->y) - (a->y*b->x);
  37.         return(c);
  38.         }
  39. static double V3Length(a)
  40. Vector3 *a;
  41. {
  42.         return(sqrt(V3SquaredLength(a)));
  43.         }
  44.  
  45. Vector3 *V3Normalize(v,name)
  46. Vector3 *v;
  47. char *name;
  48. {
  49. double len = V3Length(v);
  50.         if (len != 0.0) { v->x /= len;  v->y /= len; v->z /= len; }
  51.         else fprintf(stderr,"V3Normalize : Zero Length Vector in %s.",name);
  52.         return(v);
  53.         }
  54. static void dump_matrix(T)
  55. Matrix T;
  56. {
  57.     int i,j;
  58.     for (i=0;i<4;i++){
  59.         for(j=0;j<4;j++)
  60.             fprintf(stderr,"%g\t",T[i][j]);
  61.         fprintf(stderr,"\n");
  62.     }
  63. }
  64. static void ViewTransform( VRP, EV, UP, T)
  65. Point3 VRP, EV, UP;
  66. Matrix T;
  67. {
  68.     Vector3    U, V, N;
  69.     float    dot;
  70.  
  71.     /*
  72.      *  N = EV  and normalize  N
  73.      */
  74.     N = EV;
  75.     V3Normalize(&N,"ViewTransform : Normal Calc.");
  76.  
  77.     /*
  78.      * V = UP
  79.      * Make vector  V  orthogonal to  N  and normalize  V
  80.      */
  81.     V=UP;
  82.     dot = V3Dot(&V,&N);
  83.     V.x -= dot * N.x; V.y -= dot * N.y; V.z -= dot * N.z;
  84.     V3Normalize(&V,"ViewTransform : UP Vector Calc.");
  85.  
  86.  
  87.     /*
  88.      * Compute vector  U = V x N  (cross product)
  89.      */
  90.     V3Cross(&V,&N,&U);
  91.  
  92.     /*
  93.      * Write the vectors U, V, and N as the first three rows of the
  94.      *       first, second, and third columns of  ViewMatrix, respectively
  95.      */
  96.     T[0][0] = U.x;        /* column 1 , vector U */
  97.     T[1][0] = U.y;
  98.     T[2][0] = U.z;
  99.     T[0][1] = V.x;        /* column 2 , vector V */
  100.     T[1][1] = V.y;
  101.     T[2][1] = V.z;
  102.     T[0][2] = N.x;        /* column 3 , vector N */
  103.     T[1][2] = N.y;
  104.     T[2][2] = N.z;
  105.  
  106.     /*
  107.      * Compute the fourth row of  ViewMatrix  to include the translation of
  108.      *       VRP  to the origin
  109.      */
  110.     T[3][0] = - U.x * VRP.x - U.y * VRP.y - U.z * VRP.z;
  111.     T[3][1] = - V.x * VRP.x - V.y * VRP.y - V.z * VRP.z;
  112.     T[3][2] = - N.x * VRP.x - N.y * VRP.y - N.z * VRP.z;
  113.  
  114.     /* Fourth Column is constant */
  115.     T[0][3] = T[1][3] = T[2][3] = 0.0; T[3][3] = 1.0;
  116.     return;
  117. }
  118. /* multiply together matrices c = ab */
  119. /* note that c must not point to either of the input matrices */
  120. static void MatMul(a, b, c)
  121. Matrix a, b, c;
  122. {
  123. int i, j, k;
  124.         for (i=0; i<4; i++) {
  125.                 for (j=0; j<4; j++) {
  126.                         c[i][j] = 0;
  127.                         for (k=0; k<4; k++) c[i][j] +=
  128.                                 a[i][k] * b[k][j];
  129.                 }
  130.          }
  131. }
  132.  
  133.  
  134. static void loadmatrix(mat)
  135. Matrix mat;
  136. {
  137.     int i;
  138.     float *s=(float *)mat;
  139.     float *t=(float *)ModelMatrix;
  140.     for (i = 0; i < 16; i++) *t++ = *s++;
  141. }
  142.  
  143. void identmatrix(a)
  144.     Matrix      a;
  145. {
  146.     int i;
  147.     register float    *p=(float *)a;
  148.  
  149.     for (i=0; i<16; i++) *p++ = 0.0;
  150.  
  151.     a[0][0] = a[1][1] = a[2][2] = a[3][3] = 1.0;
  152. }
  153. static void multvector(v, a, b)
  154. Vector    v, a;
  155. Matrix    b;
  156. {
  157.     int i,j;
  158.     for(i=0;i<4;i++) for(j=0,v[i]=0.;j<4;j++) v[i] += a[j]*b[j][i];
  159. }
  160.  
  161. #define MIN(a,b)        (((a)<(b))?(a):(b))
  162. #define MAX(a,b)        (((a)>(b))?(a):(b))
  163. #define check {\
  164.     multvector(v,a,m);\
  165.     if(v[3] ==0) error("Division By Zero","get_min_max");\
  166.     x = v[0]/v[3]; y = v[1]/v[3];\
  167.     *min_x=MIN(*min_x,x); *max_x=MAX(*max_x,x);\
  168.     *min_y=MIN(*min_y,y); *max_y=MAX(*max_y,y);\
  169. }
  170. static void get_min_max(m,left, right, front, back, bottom,top,min_x,min_y,max_x,max_y)
  171. Matrix m;
  172. float left, right, front, back, bottom,top;
  173. float *min_x,*min_y,*max_x,*max_y;
  174. {
  175.     Vector a,v;
  176.     float x,y;
  177.  
  178.     *min_x = *min_y = Large; 
  179.     *max_x = *max_y = -Large;
  180.  
  181.     a[0]=left;a[1]=front;a[2]=bottom;a[3]=1.0; check
  182.     a[0]=right;a[1]=front;a[2]=bottom;a[3]=1.0; check
  183.     a[0]=left;a[1]=back;a[2]=bottom;a[3]=1.0; check
  184.     a[0]=right;a[1]=back;a[2]=bottom;a[3]=1.0; check
  185.     a[0]=left;a[1]=front;a[2]=top;a[3]=1.0; check
  186.     a[0]=right;a[1]=front;a[2]=top;a[3]=1.0; check
  187.     a[0]=left;a[1]=back;a[2]=top;a[3]=1.0; check
  188.     a[0]=right;a[1]=back;a[2]=top;a[3]=1.0; check
  189. }
  190. #undef check
  191.  
  192. static Matrix ViewMatrix[Maxviews];
  193. static int Vx1[MaxviewsX2],Vy1[MaxviewsX2],Vx2[MaxviewsX2],Vy2[MaxviewsX2];
  194. static float Wx1[Maxviews],Wy1[Maxviews],Wx2[Maxviews],Wy2[Maxviews];
  195.  
  196. void set_screen_transformation(left, right, front, back, bottom,top)
  197. float left, right, front, back, bottom, top;
  198. {
  199.     extern int window_x_size,window_y_size;
  200.     Point3 VRP,EP,UP;
  201.     Vector box[8];
  202.     float min_x,min_y,max_x,max_y;
  203.     float deltax,deltay,deltaz,maxdelta;
  204.     if ((deltax=(right - left))== 0.0)
  205.         error("Left clipping plane same as right one.",
  206.             "utils.c");
  207.  
  208.     if ((deltaz=(top - bottom))== 0.0)
  209.         error("Bottom clipping plane same as top one.",
  210.             "utils.c");
  211.  
  212.     if ((deltay=(back - front))== 0.0)
  213.         error("Front clipping plane same as back one.",
  214.             "utils.c");
  215.     maxdelta = deltax;
  216.     if (deltay > maxdelta) maxdelta = deltay;
  217.     if (deltaz > maxdelta) maxdelta = deltaz;
  218.     left -= (maxdelta-deltax)/2.;
  219.     right += (maxdelta-deltax)/2.;
  220.     bottom -= (maxdelta-deltaz)/2.;
  221.     top += (maxdelta-deltaz)/2.;
  222.     front -= (maxdelta-deltay)/2.;
  223.     back += (maxdelta-deltay)/2.;
  224.     identmatrix(ModelMatrix);
  225.     VRP.x = 0.5*(left+right); VRP.y = 0.5*(front+back); VRP.z = 0.5*(bottom+top);
  226.     /* Elevation */
  227.     UP.x = UP.y = 0.0; UP.z = 1.0;
  228.     EP.y = -1.0; EP.x = EP.z = 0.0;
  229.     ViewTransform(VRP,EP,UP,ViewMatrix[0]);
  230.     get_min_max(ViewMatrix[0],left, right, front, back, bottom, top,
  231.         &min_x,&min_y,&max_x,&max_y);
  232.         Wx1[0] = min_x; Wy1[0] = min_y; Wx2[0] = max_x; Wy2[0] = max_y;
  233.     Vx1[0] = 1; Vy1[0] = window_y_size/2-2; Vx2[0] = window_x_size/2-2; Vy2[0] = 1;
  234.         /* Side */
  235.     EP.x =1.0; EP.z = EP.y = 0.0;
  236.     ViewTransform(VRP,EP,UP,ViewMatrix[1]);
  237.     get_min_max(ViewMatrix[1],left, right, front, back, bottom, top,
  238.         &min_x,&min_y,&max_x,&max_y);
  239.         Wx1[1] = min_x; Wy1[1] = min_y; Wx2[1] = max_x; Wy2[1] = max_y;
  240.     Vx1[1] = window_x_size/2+1; Vy1[1] = window_y_size/2-2;
  241.     Vx2[1] = window_x_size-2; Vy2[1] = 1;
  242.         /* Plan */
  243.     UP.z = UP.x = 0.0; UP.y = 1.0;
  244.     EP.z = 1.0; EP.y = EP.x = 0.0;
  245.     ViewTransform(VRP,EP,UP,ViewMatrix[2]);
  246.     get_min_max(ViewMatrix[2],left, right, front, back, bottom, top,
  247.         &min_x,&min_y,&max_x,&max_y);
  248.         Wx1[2] = min_x; Wy1[2] = min_y; Wx2[2] = max_x; Wy2[2] = max_y;
  249.     Vx1[2] = 1; Vy1[2] = window_y_size-2;
  250.     Vx2[2] = window_x_size/2-2; Vy2[2] = window_y_size/2+1;
  251.         /* Isometry */
  252.     UP.x = UP.y = 0.0; UP.z = 1.0;
  253.     EP.z = 1.0; EP.y = 1.0;EP.x = 1.0;
  254.     ViewTransform(VRP,EP,UP,ViewMatrix[3]);
  255.     ViewTransform(VRP,EP,UP,ViewMatrix[4]);
  256.     get_min_max(ViewMatrix[3],left, right, front, back, bottom, top,
  257.         &min_x,&min_y,&max_x,&max_y);
  258.         Wx1[3] = Wx1[4] = min_x; Wy1[3] = Wy1[4] = min_y;
  259.     Wx2[3] = Wx2[4] = max_x; Wy2[3] = Wy2[4] = max_y;
  260.     Vx1[3] = Vx1[4] = window_x_size/2+1;
  261.         Vy1[3] = Vy1[4] = window_y_size-2;
  262.     Vx2[3] = Vx2[4] = window_x_size-2;
  263.         Vy2[3] = Vy2[4] = window_y_size/2+1;
  264.  
  265.     Vx1[5] = Vx1[6] = Vx1[7] = Vx1[8] = Vx1[9] = 1;
  266.     Vy1[5] = Vy1[6] = Vy1[7] = Vy1[8] = Vy1[9] = window_y_size-2;
  267.     Vx2[5] = Vx2[6] = Vx2[7] = Vx2[8] = Vx2[9] = window_x_size-2;
  268.     Vy2[5] = Vy2[6] = Vy2[7] = Vy2[8] = Vy2[9] = 1;
  269. }
  270.  
  271. compute_perspective_matrix(left, right, front, back, bottom, top)
  272. float left, right, front, back, bottom, top;
  273. {
  274.     extern int window_x_size,window_y_size;
  275.     Matrix m,m1;
  276.     Point3 EP,VRP,UP;
  277.     float min_x,min_y,max_x,max_y;
  278.  
  279.     fprintf(stderr,"Screen Extents are <%g %g %g> to <%g %g %g>.\n",
  280.         bound_min_x,bound_min_y,bound_min_z,
  281.         bound_max_x,bound_max_y,bound_max_z);
  282.     fprintf(stderr,"Give the view params:\n");
  283.     fprintf(stderr,"Eye Point : ");
  284.     scanf("%lf%lf%lf",&(EP.x),&(EP.y),&(EP.z));
  285.     fprintf(stderr,"View Reference Point : ");
  286.     scanf("%lf%lf%lf",&(VRP.x),&(VRP.y),&(VRP.z));
  287.     fprintf(stderr,"Up Point : ");
  288.     scanf("%lf%lf%lf",&(UP.x),&(UP.y),&(UP.z));
  289.  
  290.     ViewTransform(VRP,EP,UP,ViewMatrix[4]);
  291.     /* Project on a view plane 1 unit away */
  292.     identmatrix(m);
  293.     m[2][2]=m[3][3]=0.; m[2][3] = 1.;
  294.     MatMul(m1,m,ViewMatrix[4]);
  295.     get_min_max(ViewMatrix[4],left, right, front, back, bottom, top,
  296.         &min_x,&min_y,&max_x,&max_y);
  297.            Wx1[4] = min_x; Wy1[4] = min_y; Wx2[4] = max_x; Wy2[4] = max_y;
  298.     Vx1[4] = window_x_size/2+1; Vy1[4] = window_y_size-2;
  299.     Vx2[4] = window_x_size-2; Vy2[4] = window_y_size/2+1;
  300. }
  301.  
  302. static int get_screen_xy(n,deltan,x,y,z,sx,sy)
  303. int n,deltan;
  304. float x,y,z;
  305. int *sx,*sy;
  306. {
  307.     Vector v,v1,vv;
  308.     v[0]=x; v[1]=y; v[2]=z; v[3]=1.0;
  309.  
  310.     /* Model Transform the Point */
  311.     multvector(v1,v,ModelMatrix);
  312.     /* View Transform the point */
  313.     multvector(vv,v1,ViewMatrix[n]);
  314.     if (vv[3] == 0.0){
  315.         warn("get_screen_xy : Homogeneous Coordinate.", "utils.c");
  316.         return(1);
  317.     }
  318.     /* Window to View port Map */
  319.     *sx = Vx1[n+deltan] + (vv[0]/vv[3] - Wx1[n])*
  320.                 (Vx2[n+deltan]-Vx1[n+deltan])/(Wx2[n]-Wx1[n]);
  321.     *sy = Vy1[n+deltan] + (vv[1]/vv[3] - Wy1[n])*
  322.                 (Vy2[n+deltan]-Vy1[n+deltan])/(Wy2[n]-Wy1[n]);
  323.     return(0);
  324. }
  325. static int previous_screen_x[4]={0,320,0,320},
  326.        previous_screen_y[4]={319,319,639,639};
  327.  
  328. #define Code(x,y,c) {\
  329.         c = 0x00;\
  330.         if (x < min_x) c |= left;\
  331.         else if (x > max_x) c |= right;\
  332.         if (y > min_y) c |= bottom;\
  333.         else if (y < max_y) c |= top;\
  334.     }
  335. void clip_and_draw(x1,y1,x2,y2,min_x,min_y,max_x,max_y)
  336. int x1,y1,x2,y2,min_x,max_y,max_x,min_y;
  337. {
  338.     unsigned char left=0x01,right=0x02,bottom=0x04,top=0x08;
  339.     int x,y;
  340.     unsigned char c,c1,c2;
  341.     Code(x1,y1,c1); Code(x2,y2,c2);
  342.     while (c1 || c2){
  343.         if (c1 & c2) return;
  344.         c = c1; if (!c) c = c2;
  345.         if (c & left){
  346.             y = y1 + ((float)(y2 - y1)*(min_x-x1))/(x2-x1);
  347.             x = min_x;
  348.         }
  349.         else if (c & right){
  350.             y = y1 + ((float)(y2 - y1)*(max_x-x1))/(x2-x1);
  351.             x = max_x;
  352.         }    
  353.         else if (c & bottom){
  354.             x = x1 + ((float)(x2 - x1)*(min_y-y1))/(y2-y1);
  355.             y = min_y;
  356.         }    
  357.         else if (c & top){
  358.             x = x1 + ((float)(x2 - x1)*(max_y-y1))/(y2-y1);
  359.             y = max_y;
  360.         }    
  361.         if (c == c1) {
  362.             x1 = x; y1 = y; Code(x,y,c1);
  363.         }
  364.         else{
  365.             x2 = x; y2 = y; Code(x,y,c2);
  366.         }
  367.     }
  368.     Xdraw(x1,y1,x2,y2);
  369. }
  370. #undef Code
  371. void draw(x,y,z)
  372. float x,y,z;
  373. {
  374.     int i;
  375.     int sx,sy;
  376.     if (view_type == ONE){
  377.         if (get_screen_xy(oneview_flag,Maxviews,x,y,z,&sx,&sy))
  378.             return;
  379.         Xdraw(previous_screen_x[0],previous_screen_y[0],sx,sy);
  380.         previous_screen_x[0]=sx; previous_screen_y[0]=sy;
  381.     }
  382.     else for(i=0;i<4;i++){
  383.         if (get_screen_xy(i,0,x,y,z,&sx,&sy))
  384.             continue;
  385.         clip_and_draw(
  386.             previous_screen_x[i],previous_screen_y[i],sx,sy,
  387.             Vx1[i],Vy1[i],Vx2[i],Vy2[i]);
  388.         previous_screen_x[i]=sx; previous_screen_y[i]=sy;
  389.     }
  390. }
  391. void move(x,y,z)
  392. float x,y,z;
  393. {
  394.     int i;
  395.     int sx,sy;
  396.     if (view_type == ONE){
  397.         if (get_screen_xy(oneview_flag,Maxviews,x,y,z,&sx,&sy))
  398.             return;
  399.         previous_screen_x[0]=sx; previous_screen_y[0]=sy;
  400.     }
  401.     else for(i=0;i<4;i++){
  402.         if (get_screen_xy(i,0,x,y,z,&sx,&sy))
  403.             continue;
  404.         previous_screen_x[i]=sx; previous_screen_y[i]=sy;
  405.     }
  406. }
  407. #define MaxStackSize 100
  408. static float stack[MaxStackSize];
  409. static int topofstack=0;
  410. void pushmatrix()
  411. {
  412.     if ((topofstack+16)< MaxStackSize){
  413.         int i;
  414.         float *m=(float *)ModelMatrix;
  415.         for(i=0;i<16;i++) stack[topofstack++] = *m++;
  416.     }
  417.     else error("pushmatrix : Stack overflow","utils.c");
  418. }
  419. void popmatrix()
  420. {
  421.     if (topofstack < 16) error("popmatrix : Empty Stack","utils.c");
  422.     else {
  423.         int i;
  424.         float *m=(float *)ModelMatrix;
  425.         for(i=0,m+=15;i<16;i++) *m--=stack[--topofstack];
  426.     }
  427. }
  428.  
  429. void translate(x,y,z)
  430. float x,y,z;
  431. {
  432.     Matrix m,m1;
  433.     identmatrix(m);
  434.     m[3][0] = x; m[3][1] = y; m[3][2] = z;
  435.     MatMul(m,ModelMatrix,m1);
  436.     loadmatrix(m1);
  437. }
  438.  
  439. void rotate(r,axis)
  440. float r;
  441. char axis;
  442. {
  443.     float sintheta=sin((double)(DTOR * r));
  444.     float costheta=cos((double)(DTOR * r));
  445.     float tmp;
  446.     Matrix m,m1;
  447.     identmatrix(m);
  448.     switch(axis){
  449.     case 'x':
  450.     case 'X': 
  451.           m[1][1]=m[2][2]= costheta;
  452.           m[1][2]= sintheta;
  453.           m[2][1]= -sintheta;
  454.         break;
  455.     case 'y':
  456.     case 'Y': 
  457.           m[0][0] = m[2][2] = costheta;
  458.           m[0][2] = -sintheta;
  459.           m[2][0] = sintheta;
  460.         break;
  461.     case 'z':
  462.     case 'Z': 
  463.           m[0][0] = m[1][1] = costheta;
  464.           m[0][1] = sintheta;
  465.           m[1][0] = -sintheta;
  466.         break;
  467.     default: error("rotate : Illegal axis.","utils.c"); break;
  468.     }
  469.     MatMul(m,ModelMatrix,m1);
  470.     loadmatrix(m1);
  471. }
  472.  
  473. void multmatrix(m)
  474. Matrix m;
  475. {
  476.     Matrix m1;
  477.     MatMul(m,ModelMatrix,m1);
  478.     loadmatrix(m1);
  479. }
  480.  
  481. void circle(radius)
  482. float radius;
  483. {
  484.     extern int circle_precision;
  485.     int i;
  486.     float   cx, cy, dx, dy;
  487.     float   angle=2.0 * PI /circle_precision;
  488.         float cosine = cos((double)angle);
  489.         float sine = sin((double)angle);
  490.         cx = radius;
  491.         cy = 0.0;
  492.     move(cx,cy,0.0);
  493.         for (i = 0; i < (circle_precision - 1); i++) {
  494.                 dx = cx;
  495.                 dy = cy;
  496.                 cx = dx * cosine - dy * sine;
  497.                 cy = dx * sine + dy * cosine;
  498.                 draw(cx, cy,0.0);
  499.     }
  500.     draw(radius,0.0,0.0);
  501. }
  502.