home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / drdobbs / 1990 / 09 / lyke.asc < prev    next >
Text File  |  1990-07-25  |  10KB  |  376 lines
  1. _RAY TRACING_
  2. by Daniel Lyke
  3.  
  4. [LISTING ONE]
  5.  
  6. /* RAYTRACE.HPP */
  7. class RAY
  8. {
  9.    double dx, dy, dz; /* Direction vector */
  10.    double ox, oy, oz; /* Origin */
  11.    public:
  12.    RAY(double x, double y, double z, double vx, double vy, double vz);
  13.    friend class PLANE;
  14.    friend class SPHERE;
  15. };
  16. class PLANE
  17. {
  18.    double nx, ny, nz; /* Vector normal (perpendicular) to plane */
  19.    double px, py, pz; /* Point on plane */
  20.    public:
  21.    PLANE(double x, double y, double z, double vx, double vy, double vz);
  22.    double Intersect(RAY ray);
  23.    int Pattern(RAY ray, double time, int light);
  24. };
  25. class SPHERE
  26. {
  27.    double cx, cy, cz; /* Center of sphere */
  28.    double r2; /* Radius squared */
  29.    public:
  30.    double Intersect(RAY ray);
  31.    Reflect(RAY iray,double time, RAY &rray);
  32.    SPHERE(double x, double y, double z, double r);
  33. };
  34. class VECTOR
  35. {
  36.    public:
  37.    double dx, dy, dz; /* Three dimensional vector */
  38. };
  39.  
  40.  
  41.  
  42. [LISTING TWO] 
  43.  
  44. /* RAYTRACE.CPP */
  45.  
  46. #include <fg.h>
  47. #include <math.h>
  48. #include <stdio.h>
  49. #include <conio.h>
  50.  
  51. #include "raytrace.hpp"
  52.  
  53. #define WIDTH 640
  54. #define HEIGHT 350
  55.  
  56. inline void plot(int x,int y,int c)
  57. {
  58.    fg_drawdot(c,FG_MODE_SET,~0,x,y);
  59. }
  60.  
  61. int PlanePattern(unsigned int x, unsigned int y, int light)
  62. {
  63.    /* Put code for different plane patterns in here */
  64. //    return ((x + y) % 8) + 8 * light;
  65. //    return (x % 8) ^ (y % 8) + 8 * light;
  66.     return ((x * x + y * y) % 8) + 8 * light;
  67. } /*----- End: PlanePattern() -----*/
  68.  
  69. RAY::RAY(double x, double y, double z, double vx, double vy, double vz)
  70. {
  71.    this->ox = x;
  72.    this->oy = y;
  73.    this->oz = z;
  74.    this->dx = vx;
  75.    this->dy = vy;
  76.    this->dz = vz;
  77. } /*----- End: RAY::RAY() -----*/
  78.  
  79. SPHERE::SPHERE(double x, double y, double z, double r)
  80. {
  81.    this->cx = x;
  82.    this->cy = y;
  83.    this->cz = z;
  84.    this->r2 = r * r;
  85. } /*----- End: SPHERE::SPHERE ------*/
  86.  
  87. double SPHERE::Intersect(RAY ray)
  88. {
  89.    double a, b, c, t1, t2, t3, close, farther;
  90.    a = ray.dx * ray.dx + ray.dy * ray.dy + ray.dz * ray.dz;
  91.    close = farther = -1.0;
  92.    if(a)
  93.       {
  94.       b = 2.0 * ((ray.ox - this->cx) * ray.dx
  95.          + (ray.oy - this->cy) * ray.dy
  96.          + (ray.oz - this->cz) * ray.dz);
  97.       c = (ray.ox - this->cx) * (ray.ox - this->cx)
  98.          + (ray.oy - this->cy) * (ray.oy - this->cy)
  99.          + (ray.oz - this->cz) * (ray.oz - this->cz) - this->r2;
  100.       t1 = b * b - 4.0 * a * c;
  101.       if(t1 > 0)
  102.          {
  103.          t2 = sqrt(t1);
  104.          t3 = 2.0 * a;
  105.          close = -(b + t2) / t3;
  106.          farther = -(b - t2) / t3;
  107.          }
  108.       }
  109.    return (double)((close < farther) ? close : farther);
  110. } /*----- End: SPHERE::Intersect() -----*/
  111.  
  112. SPHERE::Reflect(RAY iray, double time, RAY &rray)
  113. {
  114.    VECTOR normal;                 /* Used for readability      */
  115.    double ndotn;                  /* Used for readability      */
  116.    double idotn;                  /* Used for readability      */
  117.    double idotn_div_ndotn_x2;     /* Used for optimization     */
  118.  
  119.    rray.ox = iray.dx * time + iray.ox; /* Find the point of     */
  120.    rray.oy = iray.dy * time + iray.oy; /* intersection between  */
  121.    rray.oz = iray.dz * time + iray.oz; /* iray and sphere.      */
  122.  
  123.    normal.dx = rray.ox - this->cx;      /* Find the ray normal  */ 
  124.    normal.dy = rray.oy - this->cy;      /* to the sphere at the */
  125.    normal.dz = rray.oz - this->cz;      /* intersection         */
  126.  
  127.    ndotn = (normal.dx * normal.dx +
  128.             normal.dy * normal.dy +
  129.             normal.dz * normal.dz);
  130.    idotn = (normal.dx * iray.dx +
  131.             normal.dy * iray.dy +
  132.             normal.dz * iray.dz);
  133.    idotn_div_ndotn_x2 = (2.0 * (idotn) / ndotn);
  134.  
  135.    rray.dx = iray.dx - idotn_div_ndotn_x2 * normal.dx;
  136.    rray.dy = iray.dy - idotn_div_ndotn_x2 * normal.dy;
  137.    rray.dz = iray.dz - idotn_div_ndotn_x2 * normal.dz;
  138. } /*----- End: SPHERE::Reflect() ------*/
  139.  
  140. PLANE::PLANE(double x, double y, double z, double vx, double vy, double vz)
  141. {
  142.    this->nx   = vx;
  143.    this->ny = vy;
  144.    this->nz = vz;
  145.  
  146.    this->px = x;
  147.    this->py = y;
  148.    this->pz = z;
  149. } /*----- End: PLANE::PLANE() -----*/
  150.  
  151. int PLANE::Pattern(RAY ray, double time, int light)
  152. {
  153.    PlanePattern((unsigned)(time * ray.dz + ray.oz),
  154.       (unsigned)(time * ray.dx + ray.ox),light);
  155. } /*----- End: PLANE::Pattern ------*/
  156.  
  157. double PLANE::Intersect(RAY ray)
  158. {
  159.    double p1, p2, p3;
  160.  
  161.    p1 = this->px * this->ny + this->py * this->ny + this->pz * this->nz;
  162.    p2 = ray.ox * this->nx + ray.oy * this->ny + ray.oz * this->nz;
  163.    p3 = ray.dx * this->nx + ray.dy * this->ny + ray.dz * this->nz;
  164.    return (double)((p1-p2)/p3);
  165. } /*----- End: PLANE::Intersect() -----*/
  166.  
  167. int trace(double x, double y)
  168. {
  169.    static PLANE plane(-8.0, 0.0, 0.0,   0.0, 1.0, 0.001);
  170.    static SPHERE sphere( 0.0, 0.0, 5.0, 1.0 );
  171.    RAY ray(0.0,0.0,0.0,(x - (double)WIDTH / 2.0) *.75,
  172.       y - (double)HEIGHT / 2.0,HEIGHT);
  173.    double time1, time2,;
  174.    time1 = sphere.Intersect(ray);
  175.    time2 = plane.Intersect(ray);
  176.    if(time1 > 0.0  && (time2 < 0.0 || time2 > time1)) /* Circle in fore */
  177.       {
  178.       sphere.Reflect(ray,time1,ray);
  179.       time2 = plane.Intersect(ray);
  180.       if(time2 > 0.0)
  181.          {
  182.          return plane.Pattern(ray,time2,0);
  183.          }
  184.       else
  185.          {
  186.          return 1;
  187.          }
  188.       }
  189.      else if(time2 > 0.0)
  190.       {
  191.       return plane.Pattern(ray,time2,1);
  192.       }
  193.    return 0;
  194. } /*----- End: trace() -----*/
  195.  
  196. draw()
  197. {
  198.    int x,y;
  199.    for(x = 0; x < WIDTH && !kbhit(); x ++)
  200.       {
  201.       for(y = 0; y < HEIGHT; y ++)
  202.          {
  203.          plot(x,y,trace((double)x,(double)y));
  204.          }
  205.       }
  206. } /*----- End: draw() -----*/
  207.  
  208. main()
  209. {
  210.    fg_init_egaecd();
  211.    draw();
  212.    getch();
  213.    fg_term();
  214. }
  215.  
  216.  
  217. [LISTING THREE]
  218.  
  219. #include <fg.h>
  220. #include <math.h>
  221. #include <dos.h>
  222.  
  223. #define WIDTH 640
  224. #define HEIGHT 350
  225. #define QUIT_OUT (!kbhit())
  226.  
  227. plot(x,y,c)
  228. int x,y,c;
  229. {
  230. fg_drawdot(c,FG_MODE_SET,~0,x,HEIGHT - y);
  231. }
  232.  
  233. typedef struct S_RAY
  234. {
  235.     double dx, dy, dz; /* Direction vector */
  236.     double ox, oy, oz; /* Origin */
  237. } RAY;
  238.  
  239. typedef struct S_PLANE
  240. {
  241.     double nx, ny, nz; /* Vector normal (perpendicular) to plane */
  242.     double px, py, pz; /* Point on plane */
  243. } PLANE;
  244.  
  245. typedef struct S_SPHERE
  246. {
  247.     double cx, cy, cz; /* Center of sphere */
  248.     double r2; /* Radius squared */
  249. } SPHERE;
  250.  
  251. typedef struct S_VECTOR
  252. {
  253.     double dx, dy, dz; /* Three dimensional vector */
  254. } VECTOR;
  255.  
  256. double sphere_intersect(RAY, SPHERE);
  257. double plane_intersect(RAY, PLANE);
  258. void reflect(VECTOR *, VECTOR *, VECTOR *);
  259. double sphere_intersect(RAY ray, SPHERE sphere)
  260. {
  261.     double a, b, c, t1, t2, t3, close, farther;
  262.     a = ray.dx * ray.dx + ray.dy * ray.dy + ray.dz * ray.dz;
  263.     close = farther = -1.0;
  264.     if(a)
  265.         {
  266.         b = 2.0 * ((ray.ox - sphere.cx) * ray.dx
  267.           + (ray.oy - sphere.cy) * ray.dy
  268.           + (ray.oz - sphere.cz) * ray.dz);
  269.         c = (ray.ox - sphere.cx) * (ray.ox - sphere.cx)
  270.           + (ray.oy - sphere.cy) * (ray.oy - sphere.cy)
  271.           + (ray.oz - sphere.cz) * (ray.oz - sphere.cz) - sphere.r2;
  272.         t1 = b * b - 4.0 * a * c;
  273.         if(t1 > 0)
  274.           {
  275.           t2 = sqrt(t1);
  276.           t3 = 2.0 * a;
  277.           close = -(b + t2) / t3;
  278.           farther = -(b - t2) / t3;
  279.           }
  280.         }
  281.     return (double)((close < farther) ? close : farther);
  282. } /*----- End: sphere_intersect() -----*/
  283.  
  284. double plane_intersect(RAY ray, PLANE plane)
  285. {
  286.     double p1, p2, p3;
  287.     p1 = plane.px * plane.ny + plane.py * plane.ny + plane.pz * plane.nz;
  288.     p2 = ray.ox * plane.nx + ray.oy * plane.ny + ray.oz * plane.nz;
  289.     p3 = ray.dx * plane.nx + ray.dy * plane.ny + ray.dz * plane.nz;
  290.     return (double)((p1-p2)/p3);
  291. } /*----- End: plane_intersect() -----*/
  292.  
  293. void reflect(VECTOR *normal, VECTOR *incident, VECTOR *r)
  294. {
  295.     double ndotn, idotn;
  296.     ndotn = (normal->dx * normal->dx +
  297.                 normal->dy * normal->dy +
  298.                 normal->dz * normal->dz);
  299.     idotn = (normal->dx * incident->dx +
  300.                 normal->dy * incident->dy +
  301.                 normal->dz * incident->dz);
  302.     r->dx = incident->dx - (2.0 * (idotn) / ndotn) * normal->dx;
  303.     r->dy = incident->dy - (2.0 * (idotn) / ndotn) * normal->dy;
  304.     r->dz = incident->dz - (2.0 * (idotn) / ndotn) * normal->dz;
  305. } /*----- End: reflect() -----*/
  306.  
  307. int plane_pattern(int x, int y, int light)
  308. {
  309.     return ((x + 16384) % 8) ^ ((y + 16384) % 8) + 8 * light;
  310. } /*----- End: plane_pattern() -----*/
  311.  
  312. int trace(int x, int y)
  313. {
  314.     static PLANE plane = { 0.0, 1.0, 0.001, -8.0, 0.0, 0.0};
  315.     static SPHERE sphere = { 0.0, 0.0, 5.0, 9.0 };
  316.     VECTOR v1, v2, v3;
  317.     RAY ray;
  318.     double t1, t2, time;
  319.     ray.ox = 0.0;         /* Set the ray origin to the eye */
  320.     ray.oy = 0.0;
  321.     ray.oz = 0.0;
  322.     ray.dz = 1.0;        /* Set the direction through the pixel    */
  323.     ray.dy = -((double)y - (double)HEIGHT / 2.0) / 100;
  324.     ray.dx = ((double)x - (double)WIDTH / 2.0) / 120;
  325.     t1 = sphere_intersect(ray,sphere);
  326.     t2 = plane_intersect(ray,plane);
  327.     if(t1 > 0.0 && (t2 < 0.0 || t2 > t1)) /* Circle in fore */
  328.       {
  329.       v1.dx = ray.dx; v1.dy = ray.dy; v1.dz = ray.dz;
  330.       v2.dx = ((ray.dx * t1 + ray.ox) - sphere.cx);
  331.       v2.dy = ((ray.dy * t1 + ray.oy) - sphere.cy);
  332.       v2.dz = ((ray.dz * t1 + ray.oz) - sphere.cz);
  333.       reflect(&v2,&v1, &v3);
  334.       ray.ox += ray.dx * t1; ray.oy += ray.dy * t1; ray.oz += ray.dz * t1;
  335.       ray.dx = v3.dx; ray.dy = v3.dy; ray.dz = v3.dz;
  336.       t2 = plane_intersect(ray,plane);
  337.       if(t2 > 0.0)
  338.         {
  339.                 return plane_pattern((int)(t2 * ray.dz + ray.oz),(int)(t2 * 
  340.                                                           ray.dx + ray.ox),0);
  341.         }
  342.       else
  343.         {
  344.         return 1;
  345.         }
  346.       }
  347.     else if(t2 > 0.0)
  348.       {
  349.       return plane_pattern((int)(t2 * ray.dz + ray.oz),(int)(t2 * 
  350.                                                           ray.dx + ray.ox),1);
  351.       }
  352.     return 0;
  353. } /*----- End: trace() -----*/
  354.  
  355. draw()
  356. {
  357.     int x,y;
  358.     for(x=0;x< WIDTH && QUIT_OUT; x++)
  359.         {
  360.         for(y = 0; y < HEIGHT; y++)
  361.             {
  362.             plot(x,y,trace(x,y));
  363.             }
  364.         }
  365. } /*----- End: draw() -----*/
  366.  
  367. main()
  368. {
  369.     fg_init_all();
  370.     draw();
  371.     while(QUIT_OUT);
  372.     getch();
  373.     fg_term();
  374. }
  375.  
  376.