home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / rad386 / radiosit / src / ray_poly.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-21  |  5.7 KB  |  185 lines
  1. /*
  2.     Author: S.N.pattanaik
  3.     Date : 24 March 1991.
  4.  
  5.     NOTE : The Quadrilateral must be CONVEX.
  6. */
  7.  
  8. #include <stdio.h>
  9. #include <math.h>
  10. #include "GraphicsGems.h"
  11. #include "data_structure.h"
  12.  
  13. double plane_intersect(ray_start,ray_direction,
  14.             plane_normal,plane_constant)
  15. /* 
  16.     This routine returns the parameter of intersection
  17.     of a ray defined as (ray_start+t.ray_direction) and
  18.     a polygon.
  19. NOTE : No assumption of the vectors being normalised has been made.
  20. */
  21. Point3 *ray_start;            /* Input */
  22. Vector3 *ray_direction,*plane_normal;    /* Input */
  23. double plane_constant;            /* Input */
  24. {
  25.         double t = -1.0,product;
  26.         if(fabs(product=V3Dot(ray_direction,plane_normal))>EPSILON){
  27.           if(product<0.0)
  28.         t = -(plane_constant+V3Dot(ray_start,plane_normal))/product;
  29.         }
  30.         return(t);
  31. }
  32.  
  33. double ray_polygon(ray_start,ray_direction,plane_normal, plane_constant,dominant_axis,n,vlist)
  34. /*
  35.     Adapted from  "An Introduction to Ray Tracing" by Andrew Glassner.
  36.  
  37.     Project the polygon points on a co-ordinate plane.
  38.            The projection does not preserve area but topology
  39.            is preserved.
  40.            Projection is simply throwing away that co-ordinate.
  41.            Translate the polygon so that the intersection point is at origin.
  42.            This simply means subtract the intersection point coordinate from
  43.            each vertex of the polygon.
  44.            Test each edge of the polygon against the +U axis. Count the total
  45.            number of crosses. If the number is odd then the point is inside
  46.            else outside.
  47.     The algo has been given a winding number flavour. - Check -
  48. */
  49. Point3 *ray_start;        /* Input */
  50. Vector3 *ray_direction;        /* Input */
  51. Vector3 *plane_normal;        /* Input */
  52. double plane_constant;        /* Input */
  53. char dominant_axis;        /* Input */
  54. int n;                /* Input */
  55. Point3 *vlist;            /* Input */
  56. {
  57. /*
  58.  * Quadrants are defined as:
  59.  *        |
  60.  *   1    |   0
  61.  *        |
  62.  * -------c--------
  63.  *        |
  64.  *   2    |   3
  65.  *        |
  66.  */
  67. #define quadrant(p, c) ((p.x<c.x) ? ((p.y<c.y) ? 2 : 1) : ((p.y<c.y) ? 3 : 0))
  68.  
  69.  
  70.     double t=plane_intersect(ray_start,ray_direction,
  71.             plane_normal,plane_constant);
  72.     if (t > EPSILON){ /* Check if point inside polygon or not */
  73.         int i,num_of_crossings=0;
  74.         Point2 C,last,cur;
  75.         char lastquad,curquad;
  76.  
  77.         if (dominant_axis == 0){ /* X Axis */
  78.             C.x = (*ray_start).y + t*(*ray_direction).y;
  79.             C.y = (*ray_start).z + t*(*ray_direction).z;
  80.         }
  81.         else if (dominant_axis == 1){ /* Y Axis */
  82.             C.x = (*ray_start).z + t*(*ray_direction).z;
  83.             C.y = (*ray_start).x + t*(*ray_direction).x;
  84.         }
  85.         else{ /* Z Axis */
  86.             C.x = (*ray_start).x + t*(*ray_direction).x;
  87.             C.y = (*ray_start).y + t*(*ray_direction).y;
  88.         }
  89.         project_point(last,vlist[n-1],dominant_axis);
  90.         lastquad = quadrant(last,C);
  91.         for(i=0; i< n;i++,last=cur){
  92.             project_point(cur,vlist[i],dominant_axis);
  93.             curquad = quadrant(cur,C);
  94.             if (curquad == lastquad) continue;
  95.             if (((lastquad+1) & 3)==curquad) num_of_crossings++;
  96.             else if (((curquad + 1) & 3) == lastquad) num_of_crossings--;
  97.             else{
  98.                         /*
  99.                          * Find where edge crosses
  100.                          * center's X axis.
  101.                          */
  102.                 double right,left;
  103.                 right = last.x - cur.x;
  104.                 left = (last.y - cur.y) * (C.x - last.x);
  105.                 if (left + last.y * right > right * C.y)
  106.                     num_of_crossings+=2;
  107.                 else
  108.                     num_of_crossings-=2;
  109.             }
  110.             lastquad = curquad;
  111.         }
  112.         if (num_of_crossings != 0) return(t);
  113.     }
  114.     return(-1.);
  115. }
  116.  
  117. double ray_quadrilateral(ray_start,ray_direction,
  118.             plane_normal, plane_constant,
  119.             Du0,Du1,Du2,Dux,Duy,
  120.             Dv0,Dv1,Dv2,Dvx,Dvy,
  121.             Qux,Quy,Qvx,Qvy,Na,Nb,Nc,
  122.             pu,pv)
  123. /*
  124.     This routine computes if
  125.         the ray intersects the CONVEX quadrilateral.
  126.         and
  127.         returns the ray parameter t and 
  128.         surface parameters (u,v) of intersection.
  129.     Because of the (u,v) requirement, the following method
  130.         uses inverse mapping of the point of intersection
  131.         of the ray with the plane to compute (u,v).
  132.         For the point to be inside the quadrilateral the
  133.         (u,v) will be in the range 0 to 1.
  134.         So t > 0 and
  135.         (u,v) are in the range 0 to 1 for a valid intersection.
  136.  
  137. Source: "Essential Ray Tracing Algorithms" by Eric Haines in
  138.     "An Introduction to RayTracing", page 59.
  139.  
  140. NOTE : Assume that the Quadrilateral is CONVEX.
  141. */
  142. Point3 *ray_start;        /* Input */
  143. Vector3 *ray_direction;        /* Input */
  144. Vector3 *plane_normal;        /* Input */
  145. double plane_constant;        /* Input */
  146. double  Du0,Du1,Du2,Dux,Duy,
  147.     Dv0,Dv1,Dv2,Dvx,Dvy;
  148. Vector3 *Qux,*Quy,*Qvx,*Qvy;
  149. Vector3 *Na,*Nb,*Nc;        /* Input. Precomputed according to the
  150.                    equations given in the Source. */
  151. double *pu, *pv;         /* Output */
  152. {
  153.     double t=plane_intersect(ray_start,ray_direction,
  154.             plane_normal,plane_constant);
  155.     if (t > EPSILON){
  156.         Point3 R;
  157.         R.x = (*ray_start).x + t*(*ray_direction).x;
  158.         R.y = (*ray_start).y + t*(*ray_direction).y;
  159.         R.z = (*ray_start).z + t*(*ray_direction).z;
  160.         if (fabs(Du2) < EPSILON)
  161.             *pu = (V3Dot(Nc,&R)-Du0)/(Du1-V3Dot(Na,&R));
  162.         else{
  163.             double Ka=Dux+V3Dot(Qux,&R),
  164.                 Kb=Duy+V3Dot(Quy,&R);
  165.             double discriminant=sqrt(Ka*Ka-Kb);
  166.             if (Ka < discriminant) *pu = Ka + discriminant;
  167.             else *pu = Ka - discriminant;
  168.         }
  169.         if ( (*pu < 0.0) || (*pu > 1.0)) return(-1.0);
  170.         if (fabs(Dv2) < EPSILON)
  171.             *pv = (V3Dot(Nb,&R)-Dv0)/(Dv1-V3Dot(Na,&R));
  172.         else{
  173.             double Ka=Dvx+V3Dot(Qvx,&R),
  174.                 Kb=Dvy+V3Dot(Qvy,&R);
  175.             double discriminant=sqrt(Ka*Ka-Kb);
  176.             if (discriminant < 0.0) return(-1.0);
  177.             if (Ka < discriminant) *pv = Ka + discriminant;
  178.             else *pv = Ka - discriminant;
  179.         }
  180.         if ( (*pv < 0.0) || (*pv > 1.0)) return(-1.0);
  181.         return(t);
  182.     }
  183.     return(-1.);
  184. }
  185.