home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch5_4 / sweep.cpp < prev   
Encoding:
C/C++ Source or Header  |  1995-03-04  |  3.7 KB  |  136 lines
  1. /*************************************************
  2.  *  SWEEP.CPP
  3.  *  Andreas Leipelt, "Ray Tracing a Swept Sphere"
  4.  *  from "Graphics Gems", Academic Press
  5.  *
  6.  *  This file contains the code to handle a swept sphere in
  7.  *  ray tracing
  8.  */
  9.  
  10. #include <math.h>
  11. #include "poly.h"
  12.  
  13. #define rayeps   1E-8  // tolerance for intersection test
  14.  
  15. // refer to  Andrew Woo, "Fast Ray-Box Intersection",
  16. // "Graphics Gems I"
  17. extern char HitBoundingBox(double*,double*,double*,double*);
  18.  
  19. // class of the swept sphere primitive
  20. class swept_sph {
  21.   polynomial m[3]; // center of the sphere
  22.   polynomial r;    // radius of the sphere
  23.   polynomial r2;   // r2 = r*r
  24.   double a, b;     // the interval [a;b], where  m  and  r  live
  25.   double minB[3],  // lower left corner of the bounding box
  26.          maxB[3];  // upper right corner of the bounding box
  27.   double param;    // parameter of last intersection, used for member
  28.                    // 'normal'
  29.   public:
  30.  
  31.   swept_sph() {}
  32.   swept_sph(polynomial*,polynomial,double,double);
  33.   int  intersect(double*,double*,double*);
  34.   void normal(double*,double*);
  35.   int  inside(double*);
  36. };
  37.  
  38. // constructor of the swept_sph-class
  39. swept_sph::swept_sph(polynomial *M, polynomial R, double A, double B)
  40. // M : trajectory of the center of the moving sphere.
  41. //     An array of polynomials, which is interpreted as a
  42. //     vector valued polynomial.
  43. // R : varying radius of the moving sphere. The radius is assumed
  44. //     to be non-negative.
  45. {
  46.   for (int i=0; i < 3; i++) m[i] = M[i];
  47.   r = R;
  48.   r2 = r*r;
  49.   a = A;  b = B;
  50.   // Calculate the axis aligned bounding box
  51.   for (i=0; i < 3; i++) {
  52.     minB[i] = (m[i] - r).min(a, b);
  53.     maxB[i] = (m[i] + r).max(a, b);
  54.   }
  55. }
  56.  
  57. int swept_sph::intersect(double *origin, double *dir, double *l)
  58. // origin : origin of the ray
  59. // dir    : unit direction of the ray
  60. // t      : intersection parameter of the ray
  61. {
  62.   polynomial p, q, dp, dq, s;
  63.   double save[3];
  64.   double roots[MAX_DEGREE];
  65.   double p_val, q_val, D, test;
  66.  
  67.   if (!HitBoundingBox(minB, maxB, origin, dir)) return 0;
  68.   // save the constant term of the trajectory
  69.   for (int i=0; i < 3; i++) {
  70.     save[i] = m[i].coef[0];
  71.     m[i].coef[0] -= origin[i];
  72.   }
  73.   p = dir[0]*m[0] + dir[1]*m[1] + dir[2]*m[2];
  74.   q = m[0]*m[0] + m[1]*m[1] + m[2]*m[2] - r2;
  75.   dp = p.derivative();
  76.   dq = q.derivative();
  77.   s = dq*dq + 4.0*dp*(dp*q - p*dq);
  78.   int n = s.roots_between(a, b, roots);
  79.   roots[n++] = a;
  80.   roots[n]   = b;
  81.   *l = 1E20;
  82.   // test all possible values
  83.   for (i=0; i <= n; i++) {
  84.     // calculate the real solutions of the equation
  85.     // l = p_val +- sqrt(p_val*p_val - q_val)
  86.     p_val = p.eval(roots[i]);
  87.     q_val = q.eval(roots[i]);
  88.     D = p_val*p_val - q_val;
  89.     if (D >= 0.0) {
  90.       // check, if the candidate  roots[i]  leads to a better
  91.       // intersection value  l
  92.       D = sqrt(D);
  93.       test = p_val - D;
  94.       if (test < rayeps) test = p_val + D;
  95.       if ((test >= rayeps) && (test < *l)) {
  96.     param = roots[i];
  97.     *l = test;
  98.       }
  99.     }
  100.   }
  101.   // restore the constant term of the trajectory
  102.   for (i=0; i < 3; i++) m[i].coef[0] = save[i];
  103.   if (*l < 1E20) return 1;
  104.   else return 0;
  105. }
  106.  
  107. void swept_sph::normal(double *IP, double* Nrm)
  108. // IP  : intersection point
  109. // Nrm : normal at IP
  110. {
  111.   double R = r.eval(param);
  112.   // if the radius is zero, return an arbitrary normal.
  113.   if (R < polyeps) {
  114.     Nrm[0] = Nrm[1] = 0.0;
  115.     Nrm[2] = 1.0;
  116.     return;
  117.   }
  118.   for (int i=0; i < 3; i++) Nrm[i] = (IP[i] - m[i].eval(param))/R;
  119. }
  120.  
  121. // returns 1, if the point P lies inside.
  122. int swept_sph::inside(double *P)
  123. {
  124.   double save[3];
  125.   int is_inside;
  126.  
  127.   for (int i=0; i < 3; i++) {
  128.     save[i] = m[i].coef[0];
  129.     m[i].coef[0] -= P[i];
  130.   };
  131.   is_inside =
  132.     ((m[0]*m[0]+m[1]*m[1]+m[2]*m[2]-r2).min(a, b) < rayeps);
  133.   for (i=0; i < 3; i++) m[i].coef[0] = save[i];
  134.   return is_inside;
  135. }
  136.