home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / radiance / simplerd.lha / simplerad / FinalFTP / WalkT / geo.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-07-20  |  7.7 KB  |  277 lines
  1. /**********************************************************************/
  2. /* geo.c :                                                            */
  3. /*                                                                    */
  4. /* Geometric operations                                               */
  5. /*                                                                    */
  6. /* Copyright (C) 1992, Bernard Kwok                                   */
  7. /* All rights reserved.                                               */
  8. /* Revision 1.0                                                       */
  9. /* May, 1992                                                          *
  10. /**********************************************************************/
  11. #include <stdio.h>
  12. #ifdef IRIS
  13. #include <gl/gl.h>
  14. #endif /* IRIS */
  15. #include <math.h>
  16. #include "misc.h"
  17. #include "geo.h"
  18.  
  19. #ifdef IRIS
  20. /*
  21.  * ResetWindMat() - Reset matrices
  22.  */
  23. void ResetWindMat()
  24. {
  25.   ResetMat(Mat);
  26.   loadmatrix(IdentityMat);
  27. }
  28. #endif IRIS
  29.  
  30. /**********************************************************************/
  31. /* find dot product of vectors v0 and v1 */
  32. /**********************************************************************/
  33. double dot(v0,v1)
  34.      Vector *v0, *v1;
  35. {  return (v0->x * v1->x) + (v0->y * v1->y) + (v0->z * v1->z); }
  36.  
  37. /**********************************************************************/
  38. /* Find cross product of vectors v0 and v1 */
  39. /**********************************************************************/
  40. Vector cross(v0, v1)
  41.      Vector *v0, *v1;
  42. {
  43.   Vector v2;
  44.   v2.x = (v0->y * v1->z) - (v0->z * v1->y);
  45.   v2.y = (v0->z * v1->x) - (v0->x * v1->z);
  46.   v2.z = (v0->x * v1->y) - (v0->y * v1->x);
  47.   return v2;
  48. }
  49.  
  50. /**********************************************************************/
  51. /* Normalized vector v0 and return length, return negative length if */
  52. /* close to zero, or less than zero */
  53. /**********************************************************************/
  54. double norm(v0)
  55.      Vector *v0;
  56. {
  57.   double len;
  58.  
  59.   if ((len = dot(v0,v0)) <= 0.0) return FALSE;
  60.   len = sqrt((double)len);
  61.   if (len < MIN_VECTOR_LENGTH) return FALSE;
  62.   v0->x = v0->x / len;
  63.   v0->y = v0->y / len;
  64.   v0->z = v0->z / len;
  65.   return TRUE;
  66. }
  67.  
  68. /**********************************************************************/
  69. /* Return the distance between two points */
  70. /**********************************************************************/
  71. double vdist(a, b)
  72.      Vector *a, *b;
  73. {
  74.   double dx, dy, dz;
  75.  
  76.   dx = a->x - b->x; dy = a->y - b->y;  dz = a->z - b->z;
  77.   return(sqrt((dx*dx)+(dy*dy)+(dz*dz)));
  78. }
  79.  
  80. /**********************************************************************/
  81. /* Create vector from 2 points, and return length */
  82. /**********************************************************************/
  83. double geo_line(p0, p1, v)
  84.      point3 *p0, *p1;
  85.      Line *v;
  86. {
  87.   v->start = *p0;
  88.   v->dir.x = p1->x - p0->x;
  89.   v->dir.y = p1->y - p0->y;
  90.   v->dir.z = p1->z - p0->z;
  91.   return norm(&v->dir);
  92. }
  93.  
  94. /**********************************************************************/
  95. /* Add 2 vectors */
  96. /**********************************************************************/
  97. Vector *vadd(v0, v1)
  98.      Vector *v0, *v1;
  99. {
  100.   static Vector *add;
  101.   
  102.   add->x = v0->x + v1->x;
  103.   add->y = v0->y + v1->y;
  104.   add->z = v0->z + v1->z;
  105.   return add;
  106. }
  107.  
  108. /**********************************************************************/
  109. /* Subtract vector v1 from v0 */
  110. /**********************************************************************/
  111. Vector *vsub(v0, v1)
  112.      Vector *v0, *v1;
  113. {
  114.   static Vector sub;
  115.   
  116.   sub.x = v0->x - v1->x;
  117.   sub.y = v0->y - v1->y;
  118.   sub.z = v0->z - v1->z;
  119.   return ⊂
  120. }
  121.  
  122. /**********************************************************************/
  123. /* Vector linear combination: a*v0 + b*v1 */
  124. /**********************************************************************/
  125. Vector *vcomb(v0,v1,a,b)
  126.      Vector *v0, *v1;
  127.      double a,b;
  128. {
  129.   static Vector comb;
  130.   
  131.   comb.x = a * v0->x + b * v1->x;
  132.   comb.y = a * v0->y + b * v1->y;
  133.   comb.z = a * v0->z + b * v1->z;
  134.   return &comb;
  135. }
  136.  
  137. /**********************************************************************/
  138. /* This function is similar to pow() but works right only for */
  139. /* positive integer b (but a lot faster) */
  140. /**********************************************************************/
  141. double Power(a,b)
  142.      double a,b;
  143. {
  144.   register int power;
  145.   double result;
  146.   
  147.   result= 1.0;
  148.   power= (int) b;
  149.   while(power > 0) {
  150.     if(power & 1)
  151.       result *= a;
  152.     a *= a;
  153.     power >>= 1;
  154.   }
  155.   return(result);
  156. }
  157.  
  158. /**********************************************************************/
  159. /* Calculate normal for a polygon from its vertices */
  160. /**********************************************************************/
  161. Vector polynorm(vert,size)
  162.      register Vector vert[];
  163.      register int size;
  164. {
  165.   Vector n;
  166.   register int i,j;
  167.   
  168.   n.x = n.y = n.z= 0.0;
  169.   for(i=0;i<size;i++) {
  170.     j= (i+1)%size;
  171.     n.x += (vert[i].y-vert[j].y) * (vert[i].z+vert[j].z);
  172.     n.y += (vert[i].z-vert[j].z) * (vert[i].x+vert[j].x);
  173.     n.z += (vert[i].x-vert[j].x) * (vert[i].y+vert[j].y);
  174.   }
  175.   if(norm(&(n)) < 0.) {
  176.     n.x = n.y = 0.;
  177.     n.z= 1.;
  178.   }
  179.   return(n);
  180. }
  181.  
  182. /**********************************************************************/
  183. /* Matrix * vector */
  184. /**********************************************************************/
  185. Vector vtransform(v,M)
  186.      Vector v;
  187.      register Matrix M;
  188. {
  189.   Vector vPrime;
  190.   
  191.   vPrime.x= v.x * M[0][0] + v.y * M[1][0] + v.z * M[2][0] + M[3][0];
  192.   vPrime.y= v.x * M[0][1] + v.y * M[1][1] + v.z * M[2][1] + M[3][1];
  193.   vPrime.z= v.x * M[0][2] + v.y * M[1][2] + v.z * M[2][2] + M[3][2];
  194.   return(vPrime);
  195. }
  196.  
  197. /**********************************************************************/
  198. /* Invert matrix a -> b (using determinant method) */
  199. /**********************************************************************/
  200. int invert_matrix(a,b)
  201.      Matrix a,b;
  202. {
  203.   double det;
  204.  
  205.   b[0][0]= a[1][1]*a[2][2] - a[1][2]*a[2][1];
  206.   b[1][0]= a[1][0]*a[2][2] - a[1][2]*a[2][0];
  207.   b[2][0]= a[1][0]*a[2][1] - a[1][1]*a[2][0];
  208.   
  209.   b[0][1]= a[0][1]*a[2][2] - a[0][2]*a[2][1];
  210.   b[1][1]= a[0][0]*a[2][2] - a[0][2]*a[2][0];
  211.   b[2][1]= a[0][0]*a[2][1] - a[0][1]*a[2][0];
  212.   
  213.   b[0][2]= a[0][1]*a[1][2] - a[0][2]*a[1][1];
  214.   b[1][2]= a[0][0]*a[1][2] - a[0][2]*a[1][0];
  215.   b[2][2]= a[0][0]*a[1][1] - a[0][1]*a[1][0];
  216.   
  217.   det= a[0][0]*b[0][0] - a[0][1]*b[1][0] + a[0][2]*b[2][0];
  218.   
  219.   /* matrix is singular */
  220.   if(det < DAMN_SMALL && det > -DAMN_SMALL) return(ERROR);
  221.  
  222.   b[0][0] /= det;
  223.   b[0][1] /= -det;
  224.   b[0][2] /= det;
  225.   b[1][0] /= -det;
  226.   b[1][1] /= det;
  227.   b[1][2] /= -det;
  228.   b[2][0] /= det;
  229.   b[2][1] /= -det;
  230.   b[2][2] /= det;
  231.   
  232.   b[3][0]= -(b[0][0]*a[3][0] + b[1][0]*a[3][1] + b[2][0]*a[3][2]);
  233.   b[3][1]= -(b[0][1]*a[3][0] + b[1][1]*a[3][1] + b[2][1]*a[3][2]);
  234.   b[3][2]= -(b[0][2]*a[3][0] + b[1][2]*a[3][1] + b[2][2]*a[3][2]);
  235.   b[3][3]= 1.;
  236.   b[0][3]= 0.;
  237.   b[1][3]= 0.;
  238.   b[2][3]= 0.;
  239.   return(OK);
  240. }
  241.  
  242. /*====================================================================*/
  243. /* Output routines */
  244. /*====================================================================*/
  245.  
  246. /**********************************************************************/
  247. /* Print a vector */
  248. /**********************************************************************/
  249. void print_vector(label,v)
  250.      char *label;
  251.      Vector *v;
  252. {
  253.   printf("\tVector: %s = %e,%e,%e\n", label, v->x, v->y, v->z);
  254. }
  255.  
  256. /**********************************************************************/
  257. /* Print a matrix */
  258. /**********************************************************************/
  259. /* originally printmatrix from trans.c, label parm added !! */
  260. void print_matrix(label, a)
  261.      char *label;
  262.      register Matrix a;
  263. {
  264.   int i,j;
  265.  
  266.   printf("\tMatrix: %s\n", label);
  267.   for(i=0;i<4;i++) {
  268.     printf("\t");
  269.     for(j=0;j<4;j++)
  270.       printf("%10.6f  ",a[i][j]);
  271.     printf("\n");
  272.   }
  273.   printf("\n");
  274. }
  275.  
  276.  
  277.