home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / unix / volume18 / mtvraytrace / part02 / tri.c < prev   
Encoding:
C/C++ Source or Header  |  1989-03-26  |  5.9 KB  |  272 lines
  1. /***********************************************************************
  2.  * $Author: markv $
  3.  * $Revision: 1.1 $
  4.  * $Date: 88/11/12 16:23:13 $
  5.  * $Log:    tri.c,v $
  6.  * Revision 1.1  88/11/12  16:23:13  markv
  7.  * Initial revision
  8.  * 
  9.  * TRIs are triangular patches with normals defined at the vertices.
  10.  * When an intersection is found, it interpolates the normal to the 
  11.  * surface at that point.
  12.  *
  13.  * Algorithm is due to Jeff Arenburg, (arenburg@trwrb.uucp) and was 
  14.  * posted to USENET.
  15.  *
  16.  * Basically, for each triangle we calculate an inverse transformation
  17.  * matrix, and use it to determine the point of intersection in the plane
  18.  * of the triangle relative to the "base point" of the triangle.  We then
  19.  * figure its coordinates relative to that base point.  These base points
  20.  * are used to find the barycentric coordinates, and then in normal 
  21.  * interpolation...
  22.  ***********************************************************************/
  23.  
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include "defs.h"
  27. #include "extern.h"
  28.  
  29. typedef struct t_spheredata {
  30.     Vec    tri_P[3] ;
  31.     Vec    tri_N[3] ;
  32.     Vec    tri_bb[3] ;
  33. } TriData ;
  34.  
  35. int TriPrint ();
  36. int TriIntersect ();
  37. int TriNormal ();
  38.  
  39. ObjectProcs TriProcs = {
  40.     TriPrint,
  41.     TriIntersect,
  42.     TriNormal,
  43. } ;
  44.  
  45. int 
  46. TriPrint(obj)
  47.  Object *obj ;
  48. {
  49.     TriData *td ;
  50.     int i ;
  51.  
  52.     td = (TriData *) obj -> o_data ;
  53.     printf("pp 3\n") ;
  54.     for (i = 0 ; i < 3 ; i ++) {
  55.         VecPrint("\t", td -> tri_P[i]) ;
  56.         VecPrint("\t", td -> tri_N[i]) ;
  57.     }
  58. }
  59.  
  60. TriIntersect(obj, ray, hit)
  61.  Object * obj ;
  62.  Ray * ray ;
  63.  Isect * hit ;
  64. {
  65.     TriData *td ;
  66.     Flt n, d, dist ;
  67.     Flt r, s, t ;
  68.     Flt a, b ;
  69.     Vec P, Q ;
  70.  
  71.     td = (TriData *) obj -> o_data ;
  72.  
  73.     /*
  74.      * The matrix td -> tri_bb transforms vectors in the world 
  75.      * space into a space with the following properties.
  76.      *
  77.      * 1.  The sides of the triangle are coincident with the
  78.      *     x and y axis, and have unit length.
  79.      * 2.  The normal to the triangle is coincident with the 
  80.      *     z axis.
  81.      *
  82.      */
  83.  
  84.     /*
  85.      * d is the slope with respect to the z axis.  If d is zero, then
  86.      * the ray is parallel to the plane of the polygon, and we count 
  87.      * it as a miss...
  88.      */
  89.  
  90.     d = VecDot(ray -> D, td -> tri_bb[2]) ;
  91.     if (fabs(d) < rayeps)
  92.         return 0 ;
  93.  
  94.     /*
  95.      * Q is a vector from the eye to the triangles "origin" vertex.
  96.      * n is then set to be the distance of the tranformed eyepoint
  97.      * to the plane in the polygon.
  98.      * Together, n and d allow you to find the distance to the polygon, 
  99.      * which is merely n / d.
  100.      */
  101.  
  102.     VecSub(td -> tri_P[0], ray -> P, Q) ;
  103.  
  104.     n = VecDot(Q, td -> tri_bb[2]) ;
  105.  
  106.     dist = n / d ;
  107.  
  108.     if (dist < rayeps) 
  109.         return 0 ;
  110.     
  111.     /* 
  112.      * Q is the point we hit.  Find its position relative to the
  113.      * origin of the triangle.
  114.      */
  115.  
  116.     RayPoint(ray, dist, Q) ;
  117.     VecSub(Q, td -> tri_P[0], Q) ;
  118.  
  119.     a = VecDot(Q, td -> tri_bb[0]) ;
  120.     b = VecDot(Q, td -> tri_bb[1]) ;
  121.  
  122.     if (a < 0.0 || b < 0.0 || a + b > 1.0) 
  123.         return 0 ;
  124.     
  125.     r = 1.0 - a - b ;
  126.     s = a ;
  127.     t = b ;
  128.  
  129.     hit -> isect_t = dist ;
  130.     hit -> isect_enter = 0 ;
  131.     hit -> isect_prim = obj ;
  132.     hit -> isect_surf = obj -> o_surf ;
  133.  
  134.     VecZero(hit->isect_normal) ;
  135.     VecAddS(r, td -> tri_N[0], hit -> isect_normal, hit -> isect_normal) ;
  136.     VecAddS(s, td -> tri_N[1], hit -> isect_normal, hit -> isect_normal) ;
  137.     VecAddS(t, td -> tri_N[2], hit -> isect_normal, hit -> isect_normal) ;
  138.     VecNormalize(hit -> isect_normal) ;
  139.  
  140.     return(1) ;
  141. }
  142.  
  143. int
  144. TriNormal(obj, hit, P, N)
  145.  Object * obj ;
  146.  Isect * hit ;
  147.  Point P, N ;
  148. {
  149.     VecCopy(hit -> isect_normal, N) ;
  150. }
  151.  
  152. Object *
  153. MakeTri(point)
  154.  Vec *point ;
  155. {
  156.     Object * o ;
  157.     TriData * td ;
  158.     Vec     N, P, Q;
  159.     int i, j ;
  160.     Flt dmin, dmax, d ;
  161.     Vec B[3] ;
  162.  
  163.     o = (Object *) malloc (sizeof(Object)) ;
  164.     o -> o_type = T_TRI ;
  165.     o -> o_procs = & TriProcs ;
  166.     o -> o_surf = CurrentSurface ;
  167.  
  168.     td = (TriData *) malloc (sizeof(TriData)) ;
  169.  
  170.     /* 
  171.      * copy in the points....
  172.      */
  173.     VecCopy(point[0], td -> tri_P[0]) ;
  174.     VecCopy(point[2], td -> tri_P[1]) ;
  175.     VecCopy(point[4], td -> tri_P[2]) ;
  176.  
  177.     /*
  178.      * and the normals, then normalize them...
  179.      */
  180.     VecCopy(point[1], td -> tri_N[0]) ;
  181.     VecCopy(point[3], td -> tri_N[1]) ;
  182.     VecCopy(point[5], td -> tri_N[2]) ;
  183.     VecNormalize(td -> tri_N[0]) ;
  184.     VecNormalize(td -> tri_N[1]) ;
  185.     VecNormalize(td -> tri_N[2]) ;
  186.  
  187.     /*
  188.      * construct the inverse of the matrix...
  189.      * | P1 |
  190.      * | P2 |
  191.      * | N  |
  192.      * and store it in td -> tri_bb[]
  193.      */
  194.     
  195.     VecSub(td -> tri_P[1], td -> tri_P[0], B[0]) ;
  196.     VecSub(td -> tri_P[2], td -> tri_P[0], B[1]) ;
  197.     VecCross(B[0], B[1], B[2]) ;
  198.     VecNormalize(B[2]) ;
  199.  
  200.     InvertMatrix(B, td -> tri_bb) ;
  201.  
  202.     for (i = 0 ; i < NSLABS ; i ++) {
  203.         dmin = HUGE ;
  204.         dmax = - HUGE ;
  205.  
  206.         for (j = 0 ; j < 3 ; j ++) {
  207.             d = VecDot(Slab[i], td -> tri_P[j]) ;
  208.             if (d < dmin) dmin = d ;
  209.             if (d > dmax) dmax = d ;
  210.         }
  211.         o -> o_dmin[i] = dmin - rayeps ;
  212.         o -> o_dmax[i] = dmax + rayeps ;
  213.     }
  214.  
  215.     o -> o_data = (void *) td ;
  216.  
  217.     return o ;
  218. }
  219.  
  220. InvertMatrix(in, out)
  221.  Vec in[3] ;
  222.  Vec out[3] ;
  223. {
  224.     int i, j, k ;
  225.     Flt tmp, det, sum ;
  226.  
  227.     out[0][0] = (in[1][1] * in[2][2] - in[1][2] * in[2][1]) ;
  228.     out[1][0] = -(in[0][1] * in[2][2] - in[0][2] * in[2][1]) ;
  229.     out[2][0] = (in[0][1] * in[1][2] - in[0][2] * in[1][1]) ;
  230.  
  231.     out[0][1] = -(in[1][0] * in[2][2] - in[1][2] * in[2][0]) ;
  232.     out[1][1] = (in[0][0] * in[2][2] - in[0][2] * in[2][0]) ;
  233.     out[2][1] = -(in[0][0] * in[1][2] - in[0][2] * in[1][0]) ;
  234.  
  235.     out[0][2] = (in[1][0] * in[2][1] - in[1][1] * in[2][0]) ;
  236.     out[1][2] = -(in[0][0] * in[2][1] - in[0][1] * in[2][0]) ;
  237.     out[2][2] = (in[0][0] * in[1][1] - in[0][1] * in[1][0]) ;
  238.     
  239.     det = 
  240.     in[0][0] * in[1][1] * in[2][2] +
  241.     in[0][1] * in[1][2] * in[2][0] +
  242.     in[0][2] * in[1][0] * in[2][1] -
  243.     in[0][2] * in[1][1] * in[2][0] -
  244.     in[0][0] * in[1][2] * in[2][1] -
  245.     in[0][1] * in[1][0] * in[2][2] ;
  246.  
  247.     det = 1 / det ;
  248.  
  249.     for (i = 0 ; i < 3 ; i ++) {
  250.         for (j = 0 ; j < 3 ; j++) {
  251.             out[i][j] *= det ;
  252.         }
  253.     }
  254.     
  255. #ifdef DEBUG
  256.     for (i = 0 ; i < 3 ; i++) {
  257.         for (j = 0 ; j < 3 ; j++) {
  258.             sum = 0.0 ;
  259.             for (k = 0 ; k < 3 ; k++) {
  260.                 sum += in[i][k] * out[k][j] ;
  261.             }
  262.             if (fabs(sum) < rayeps) {
  263.                 sum = 0.0 ;
  264.             }
  265.             printf(" %g") ;
  266.         } 
  267.         printf("\n") ;
  268.     }
  269.     printf("\n") ;
  270. #endif /* DEBUG */
  271. }
  272.