_RAY TRACING_
by Daniel Lyke

[LISTING ONE]

/* RAYTRACE.HPP */
class RAY
{
   double dx, dy, dz; /* Direction vector */
   double ox, oy, oz; /* Origin */
   public:
   RAY(double x, double y, double z, double vx, double vy, double vz);
   friend class PLANE;
   friend class SPHERE;
};
class PLANE
{
   double nx, ny, nz; /* Vector normal (perpendicular) to plane */
   double px, py, pz; /* Point on plane */
   public:
   PLANE(double x, double y, double z, double vx, double vy, double vz);
   double Intersect(RAY ray);
   int Pattern(RAY ray, double time, int light);
};
class SPHERE
{
   double cx, cy, cz; /* Center of sphere */
   double r2; /* Radius squared */
   public:
   double Intersect(RAY ray);
   Reflect(RAY iray,double time, RAY &rray);
   SPHERE(double x, double y, double z, double r);
};
class VECTOR
{
   public:
   double dx, dy, dz; /* Three dimensional vector */
};



[LISTING TWO] 

/* RAYTRACE.CPP */

#include <fg.h>
#include <math.h>
#include <stdio.h>
#include <conio.h>

#include "raytrace.hpp"

#define WIDTH 640
#define HEIGHT 350

inline void plot(int x,int y,int c)
{
   fg_drawdot(c,FG_MODE_SET,~0,x,y);
}

int PlanePattern(unsigned int x, unsigned int y, int light)
{
   /* Put code for different plane patterns in here */
//	return ((x + y) % 8) + 8 * light;
//	return (x % 8) ^ (y % 8) + 8 * light;
	return ((x * x + y * y) % 8) + 8 * light;
} /*----- End: PlanePattern() -----*/

RAY::RAY(double x, double y, double z, double vx, double vy, double vz)
{
   this->ox = x;
   this->oy = y;
   this->oz = z;
   this->dx = vx;
   this->dy = vy;
   this->dz = vz;
} /*----- End: RAY::RAY() -----*/

SPHERE::SPHERE(double x, double y, double z, double r)
{
   this->cx = x;
   this->cy = y;
   this->cz = z;
   this->r2 = r * r;
} /*----- End: SPHERE::SPHERE ------*/

double SPHERE::Intersect(RAY ray)
{
   double a, b, c, t1, t2, t3, close, farther;
   a = ray.dx * ray.dx + ray.dy * ray.dy + ray.dz * ray.dz;
   close = farther = -1.0;
   if(a)
      {
      b = 2.0 * ((ray.ox - this->cx) * ray.dx
         + (ray.oy - this->cy) * ray.dy
         + (ray.oz - this->cz) * ray.dz);
      c = (ray.ox - this->cx) * (ray.ox - this->cx)
         + (ray.oy - this->cy) * (ray.oy - this->cy)
         + (ray.oz - this->cz) * (ray.oz - this->cz) - this->r2;
      t1 = b * b - 4.0 * a * c;
      if(t1 > 0)
         {
         t2 = sqrt(t1);
         t3 = 2.0 * a;
         close = -(b + t2) / t3;
         farther = -(b - t2) / t3;
         }
      }
   return (double)((close < farther) ? close : farther);
} /*----- End: SPHERE::Intersect() -----*/

SPHERE::Reflect(RAY iray, double time, RAY &rray)
{
   VECTOR normal;                 /* Used for readability      */
   double ndotn;                  /* Used for readability      */
   double idotn;                  /* Used for readability      */
   double idotn_div_ndotn_x2;     /* Used for optimization     */

   rray.ox = iray.dx * time + iray.ox; /* Find the point of     */
   rray.oy = iray.dy * time + iray.oy; /* intersection between  */
   rray.oz = iray.dz * time + iray.oz; /* iray and sphere.      */

   normal.dx = rray.ox - this->cx;      /* Find the ray normal  */ 
   normal.dy = rray.oy - this->cy;      /* to the sphere at the */
   normal.dz = rray.oz - this->cz;      /* intersection         */

   ndotn = (normal.dx * normal.dx +
            normal.dy * normal.dy +
            normal.dz * normal.dz);
   idotn = (normal.dx * iray.dx +
            normal.dy * iray.dy +
            normal.dz * iray.dz);
   idotn_div_ndotn_x2 = (2.0 * (idotn) / ndotn);

   rray.dx = iray.dx - idotn_div_ndotn_x2 * normal.dx;
   rray.dy = iray.dy - idotn_div_ndotn_x2 * normal.dy;
   rray.dz = iray.dz - idotn_div_ndotn_x2 * normal.dz;
} /*----- End: SPHERE::Reflect() ------*/

PLANE::PLANE(double x, double y, double z, double vx, double vy, double vz)
{
   this->nx   = vx;
   this->ny = vy;
   this->nz = vz;

   this->px = x;
   this->py = y;
   this->pz = z;
} /*----- End: PLANE::PLANE() -----*/

int PLANE::Pattern(RAY ray, double time, int light)
{
   PlanePattern((unsigned)(time * ray.dz + ray.oz),
      (unsigned)(time * ray.dx + ray.ox),light);
} /*----- End: PLANE::Pattern ------*/

double PLANE::Intersect(RAY ray)
{
   double p1, p2, p3;

   p1 = this->px * this->ny + this->py * this->ny + this->pz * this->nz;
   p2 = ray.ox * this->nx + ray.oy * this->ny + ray.oz * this->nz;
   p3 = ray.dx * this->nx + ray.dy * this->ny + ray.dz * this->nz;
   return (double)((p1-p2)/p3);
} /*----- End: PLANE::Intersect() -----*/

int trace(double x, double y)
{
   static PLANE plane(-8.0, 0.0, 0.0,   0.0, 1.0, 0.001);
   static SPHERE sphere( 0.0, 0.0, 5.0, 1.0 );
   RAY ray(0.0,0.0,0.0,(x - (double)WIDTH / 2.0) *.75,
      y - (double)HEIGHT / 2.0,HEIGHT);
   double time1, time2,;
   time1 = sphere.Intersect(ray);
   time2 = plane.Intersect(ray);
   if(time1 > 0.0  && (time2 < 0.0 || time2 > time1)) /* Circle in fore */
      {
      sphere.Reflect(ray,time1,ray);
      time2 = plane.Intersect(ray);
      if(time2 > 0.0)
         {
         return plane.Pattern(ray,time2,0);
         }
      else
         {
         return 1;
         }
      }
     else if(time2 > 0.0)
      {
      return plane.Pattern(ray,time2,1);
      }
   return 0;
} /*----- End: trace() -----*/

draw()
{
   int x,y;
   for(x = 0; x < WIDTH && !kbhit(); x ++)
      {
      for(y = 0; y < HEIGHT; y ++)
         {
         plot(x,y,trace((double)x,(double)y));
         }
      }
} /*----- End: draw() -----*/

main()
{
   fg_init_egaecd();
   draw();
   getch();
   fg_term();
}


[LISTING THREE]

#include <fg.h>
#include <math.h>
#include <dos.h>

#define WIDTH 640
#define HEIGHT 350
#define QUIT_OUT (!kbhit())

plot(x,y,c)
int x,y,c;
{
fg_drawdot(c,FG_MODE_SET,~0,x,HEIGHT - y);
}

typedef struct S_RAY
{
	double dx, dy, dz; /* Direction vector */
	double ox, oy, oz; /* Origin */
} RAY;

typedef struct S_PLANE
{
	double nx, ny, nz; /* Vector normal (perpendicular) to plane */
	double px, py, pz; /* Point on plane */
} PLANE;

typedef struct S_SPHERE
{
	double cx, cy, cz; /* Center of sphere */
	double r2; /* Radius squared */
} SPHERE;

typedef struct S_VECTOR
{
	double dx, dy, dz; /* Three dimensional vector */
} VECTOR;

double sphere_intersect(RAY, SPHERE);
double plane_intersect(RAY, PLANE);
void reflect(VECTOR *, VECTOR *, VECTOR *);
double sphere_intersect(RAY ray, SPHERE sphere)
{
	double a, b, c, t1, t2, t3, close, farther;
	a = ray.dx * ray.dx + ray.dy * ray.dy + ray.dz * ray.dz;
	close = farther = -1.0;
	if(a)
		{
		b = 2.0 * ((ray.ox - sphere.cx) * ray.dx
		  + (ray.oy - sphere.cy) * ray.dy
		  + (ray.oz - sphere.cz) * ray.dz);
		c = (ray.ox - sphere.cx) * (ray.ox - sphere.cx)
		  + (ray.oy - sphere.cy) * (ray.oy - sphere.cy)
		  + (ray.oz - sphere.cz) * (ray.oz - sphere.cz) - sphere.r2;
		t1 = b * b - 4.0 * a * c;
		if(t1 > 0)
		  {
		  t2 = sqrt(t1);
		  t3 = 2.0 * a;
		  close = -(b + t2) / t3;
		  farther = -(b - t2) / t3;
		  }
		}
	return (double)((close < farther) ? close : farther);
} /*----- End: sphere_intersect() -----*/

double plane_intersect(RAY ray, PLANE plane)
{
	double p1, p2, p3;
	p1 = plane.px * plane.ny + plane.py * plane.ny + plane.pz * plane.nz;
	p2 = ray.ox * plane.nx + ray.oy * plane.ny + ray.oz * plane.nz;
	p3 = ray.dx * plane.nx + ray.dy * plane.ny + ray.dz * plane.nz;
	return (double)((p1-p2)/p3);
} /*----- End: plane_intersect() -----*/

void reflect(VECTOR *normal, VECTOR *incident, VECTOR *r)
{
	double ndotn, idotn;
	ndotn = (normal->dx * normal->dx +
				normal->dy * normal->dy +
				normal->dz * normal->dz);
	idotn = (normal->dx * incident->dx +
				normal->dy * incident->dy +
				normal->dz * incident->dz);
	r->dx = incident->dx - (2.0 * (idotn) / ndotn) * normal->dx;
	r->dy = incident->dy - (2.0 * (idotn) / ndotn) * normal->dy;
	r->dz = incident->dz - (2.0 * (idotn) / ndotn) * normal->dz;
} /*----- End: reflect() -----*/

int plane_pattern(int x, int y, int light)
{
	return ((x + 16384) % 8) ^ ((y + 16384) % 8) + 8 * light;
} /*----- End: plane_pattern() -----*/

int trace(int x, int y)
{
	static PLANE plane = { 0.0, 1.0, 0.001, -8.0, 0.0, 0.0};
	static SPHERE sphere = { 0.0, 0.0, 5.0, 9.0 };
	VECTOR v1, v2, v3;
	RAY ray;
	double t1, t2, time;
	ray.ox = 0.0;		 /* Set the ray origin to the eye */
	ray.oy = 0.0;
	ray.oz = 0.0;
	ray.dz = 1.0;		/* Set the direction through the pixel	*/
	ray.dy = -((double)y - (double)HEIGHT / 2.0) / 100;
	ray.dx = ((double)x - (double)WIDTH / 2.0) / 120;
	t1 = sphere_intersect(ray,sphere);
	t2 = plane_intersect(ray,plane);
	if(t1 > 0.0 && (t2 < 0.0 || t2 > t1)) /* Circle in fore */
	  {
	  v1.dx = ray.dx; v1.dy = ray.dy; v1.dz = ray.dz;
	  v2.dx = ((ray.dx * t1 + ray.ox) - sphere.cx);
	  v2.dy = ((ray.dy * t1 + ray.oy) - sphere.cy);
	  v2.dz = ((ray.dz * t1 + ray.oz) - sphere.cz);
	  reflect(&v2,&v1, &v3);
	  ray.ox += ray.dx * t1; ray.oy += ray.dy * t1; ray.oz += ray.dz * t1;
	  ray.dx = v3.dx; ray.dy = v3.dy; ray.dz = v3.dz;
	  t2 = plane_intersect(ray,plane);
	  if(t2 > 0.0)
		{
                return plane_pattern((int)(t2 * ray.dz + ray.oz),(int)(t2 * 
                                                          ray.dx + ray.ox),0);
		}
	  else
		{
		return 1;
		}
	  }
	else if(t2 > 0.0)
	  {
	  return plane_pattern((int)(t2 * ray.dz + ray.oz),(int)(t2 * 
                                                          ray.dx + ray.ox),1);
	  }
	return 0;
} /*----- End: trace() -----*/

draw()
{
	int x,y;
	for(x=0;x< WIDTH && QUIT_OUT; x++)
		{
		for(y = 0; y < HEIGHT; y++)
			{
			plot(x,y,trace(x,y));
			}
		}
} /*----- End: draw() -----*/

main()
{
	fg_init_all();
	draw();
	while(QUIT_OUT);
	getch();
	fg_term();
}