home *** CD-ROM | disk | FTP | other *** search
/ PC Graphics Unleashed / PC Graphics Unleashed.iso / ch18 / rad386 / radiosit / src / ray_sphe.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-21  |  2.9 KB  |  109 lines
  1. /*
  2.         Author: S.N.pattanaik
  3.         Date : 25 March 1991.
  4.  
  5.     NOTE : Ray Direction must be NORMALISED.
  6. */
  7.  
  8. #include <stdio.h>
  9. #include <math.h>
  10. #include "GraphicsGems.h"
  11. #include "data_structure.h"
  12.  
  13. double sphere_intersection(ray_start,ray_direction,
  14.                 sphere_radius, sphere_center)
  15. /*
  16.         This routine returns the parameter of intersection
  17.         of a ray defined as (ray_start+t.ray_direction) and
  18.         a sphere.
  19. NOTE : Assumption of the Ray Direction being NORMALISED has been made.
  20.                                              ----------
  21. */
  22. Point3 *ray_start;                    /* Input */
  23. Vector3 *ray_direction;                    /* Input */
  24. double sphere_radius;                    /* Input */
  25. Point3 *sphere_center;                    /* Input */
  26. {
  27.     int quadraticRoots();
  28.     double t = -1.0,roots[2];
  29.     double dx=ray_start->x-sphere_center->x;    
  30.     double dy=ray_start->y-sphere_center->y;    
  31.     double dz=ray_start->z-sphere_center->z;    
  32.     double B=2.0*(ray_direction->x*dx+
  33.         ray_direction->y*dy+ray_direction->z*dz);
  34.     double C=(dx*dx+dy*dy+dz*dz)-sphere_radius*sphere_radius;
  35.     int nroots = quadraticRoots((double)1.,B,C,roots);
  36. #define check(r,ray_start,ray_direction,lo,hi) 1
  37.     if (nroots)
  38.         set_t(0.,1.); /* Dummy parameters. */
  39. #undef check
  40. /*
  41.     double discriminant=B*B-4.0*C;
  42.     if (discriminant >= 0.0){
  43.         if (discriminant>0.0) discriminant = sqrt(discriminant);
  44.         t=(-B-discriminant)/2.0;
  45.         if ((discriminant>0.0) && (t <= 0.0)){
  46.             double t1;
  47.             t1 = (-B+discriminant)/2.0;
  48.             if (t > t1) t= t1;
  49.         }
  50.     }
  51. */
  52.     return(t);
  53. }
  54.  
  55. double ray_sphere(ray_start,ray_direction,
  56.         sphere_radius,
  57.         sphere_center,
  58.         bubble_flag,
  59.         pu,pv)
  60. /*
  61.     This routine finds out the intersection of ray with Sphere.
  62.     If the ray intersects the sphere then the routine returns 
  63.         the ray parameter t
  64.         and
  65.         surface parameters (u,v) of intersection.
  66.     t > 0 and
  67.     (u,v) are in the range 0 to 1 for a valid intersection.
  68.  
  69. Source: "Essential Ray Tracing Algorithms" by Eric Haines in
  70.         "An Introduction to RayTracing", page 48.
  71.  
  72. */    
  73. Point3 *ray_start;                    /* Input */
  74. Vector3 *ray_direction;                    /* Input */
  75. double sphere_radius;                    /* Input */
  76. Point3 *sphere_center;                    /* Input */
  77. int bubble_flag;                    /* Input */
  78. double *pu,*pv;                        /* Output */
  79. {
  80.     double t = sphere_intersection(ray_start,ray_direction,
  81.             sphere_radius,sphere_center);
  82.         if (t > EPSILON){
  83.                 Point3 R,N;
  84.         double u,v,phi;
  85.  
  86.         point_on_line(*ray_start,*ray_direction,t,R);
  87.         N.x = (R.x - sphere_center->x)/sphere_radius;
  88.         N.y = (R.y - sphere_center->y)/sphere_radius;
  89.         N.z = (R.z - sphere_center->z)/sphere_radius;
  90.  
  91.         if (N.z >= 1.) phi = 0.;
  92.         else if (N.z <= -1.) phi = PI;
  93.         else phi = acos(N.z);
  94.         v = phi / PI;
  95.         if ((v == 0.) || (v == 1.)) u = 0;
  96.         else{
  97.             double r = sin(phi);
  98.             circle_inverse_mapping(r,N.x,N.y,u)
  99.         }
  100. /*
  101.         *pu = (bubble_flag)?(1.-u):u;
  102. */
  103.         *pu = u;
  104.         *pv = v;
  105.         return(t);
  106.     }
  107.     return(-1.0);
  108. }
  109.