RAY TRACING
Rendering 3-D solid objects is easier than you think
Daniel Lyke
Daniel is a programmer for Signal Data and he can be contacted at 5475 Hixson Pike, Suite I-202, Hixson, TN 37343.
I have always been fascinated by computer rendered three-dimensional (3-D) images, especially those that strive for realism. Various approaches exist for drawing 3-D images, ranging from algorithms for wire-frame drawings, where all objects can be seen through, and only the edges are shown, to techniques that create solid opaque objects, mirrored surfaces, shading for light sources, and shadows.
Of all the methods to render solids, ray tracing is probably the easiest to understand. Basically, you take a straight line (ray), drawn from your eye through a pixel on the screen, and see what objects it hits in the computer universe. You then find whichever object is closest to the origin of that line (your eye), and color the pixel on the screen the color of the object. If that object reflects or refracts, you just need to compute the new direction of the line and repeat the process from the new origin. Pretty simple, really. The hard part is representing the universe.
The Universe
Starting at the lowest level, a ray has an origin and a direction vector. The origin is a point in three-dimensional Cartesian space (x,y,z), and the direction is a change in each of the three axes. For the purposes of finding how far it is to an object, we'll multiply the vector by a "time" or distance.
For the purposes of this article, the objects I use are a sphere and a plane. The sphere is easy: It merely has an origin and a radius; although for simpler calculations the radius can be stored as the radius squared. The plane is a little more complicated.
Rather than express a plane as its equation in 3-D space (something such as "ax + by + cz = d"), I'll represent it as a vector normal (perpendicular) to the plane (which really means taking the "a," "b," and "c" constants in the expression above and using them as directions) and a point on the plane.
Before going further, I should mention that in the context of this article, the "." operator (not to be confused with ".", which is the structure element selection operator in C) means multiply the individual components of the specified vector together and sum them (for example, ray1 * ray2 results in ray1.x * ray2.x + ray1.y * ray2.y + ray1.z * ray2.z). Where it applies to items with multiple components, multiplication by scalars means that you multiply every component of the item by that scalar.
Finding the Intersection of a Ray and a Sphere
The sphere, in its simplest form, is expressed as the equation x2 + y2 + z2 = r2. This is why we store the radius as the radius squared. If we offset the ray's origin by the center of the sphere and use "time" to express where along the direction vector the ray intersects the sphere, we end up with this:
(time * ray_direction - ray_origin)<SUP>2</SUP> = sphere_radius<SUP>2</SUP>
which, when expanded, leads to:
a time<SUP>2</SUP> + b time + c = 0
where:
a = ray_direction * ray_direction b = 2(ray_origin - sphere_origin) * ray_ direction c = (ray_origin - sphere_origin) * (ray_origin - sphere_origin) - sphere_radius<SUP>2</SUP>
Delving into my old math textbooks, I come up with the equations for solving such a quadratic:
____________
\/ B<SUP>2</SUP> - 4 ac
time = -b +/- -------------
2a
We get two results because the ray can intersect the sphere at two points, when it enters the sphere and when it leaves it. There's a third possibility, which is that the argument to the square root is negative, in which case the ray never intersects the sphere at all. The function that does all of this returns the shortest positive time to intersection, or a negative number if the ray never intersects the sphere.
To find the intersection of a ray with a plane, follow this procedure:
a = plane_point * plane_normal
b = plane_normal * ray_origin
c = ray_direction * plane_normal
a - b
time = -----
c
Now we have the means to draw all of our objects. We merely take a ray through a screen pixel, find the time to intersect with all of our objects, take the color of the object that intersects in the shortest time, and color the pixel the color of that object. Life is nice in a non-reflecting universe. Unfortunately, interesting pictures aren't so easy. To do things like putting patterns on surfaces (such as the plane) or finding where to reflect objects that aren't flat (such as the sphere), you need to find the x, y, and z coordinates of the point of intersection by multiplying the direction of the ray by the time, and add it to the origin of the ray. To find the x coordinate use the equation:
x_intersect = ray_origin_x + ray_direction_x (time).
Repeat the sequence for the y and z coordinates, respectively.
Reflections
To determine reflection, use an incident vector (the one that I'm trying to reflect) and the normal vector (the one perpendicular to the surface that we're reflecting from). You can split the incident vector into two component vectors, the sum of which is the incident vector.
If one of these components is the normal vector, the other is the component of the vector perpendicular to the normal vector. So if you subtract this perpendicular vector from the normal vector, you end up with the reflected vector. The perpendicular vector is shown in Example 1.
Example 1: Perpendicular vector
incident . normal
perpendicular_vector = _________________ normal
normal . normal
If you subtract this vector perpendicular to the normal vector from the incident vector we get the normal vector and we can subtract two times the reflected vector from the incident vector to get the reflected vector, shown in Example 2.
Example 2: Reflected vector
incident . normal
reflected_vector = incident - 2 * _________________ normal
normal . normal
About the Program
The program presented here compiles and runs with Zortech C++ 2.x and Flash Graphics on an EGA in 640 x 350 mode. I've hardcoded it for a specific mode so that the screen sizes could be constants (and thereby nursed a few extra clock cycles from the system) and so that I could make the colors on the sphere darker versions of the reflected colors. Listing One (page 158) is the header file RAYTRACE.HPP and Listing Two (page 158) is the C++ source. Additionally, Listing Three (page 159) lists a C version of the program.
My first version ran on a CGA in four-color mode, so it should be relatively easy to port this to anything else. The machine-specific parts are in main( ) and plot( ), with some color-specific information in PlanePattern( ) and the defines WIDTH and HEIGHT being the screen width and height. Figure 1 illustrates the output of the program.
For ease of programming, I've put the eye at the origin (0,0,0). The screen is one unit in front of the eye: The delta (change) in the z component of the vector is 1. From there, calculate the position of each pixel relative to the center of the screen and put that into the x and y components of the vector.
trace( ) follows a ray through a pixel on the screen (the parameter) until it hits something, in this case through SPHERE:: Intersect( )and PLANE:: Intersect( ). Then it just checks the sphere and the plane to see which one has the shortest time to intersection. If both of the times are negative, the ray goes off into black sky and trace( ) returns 0.
If we hit the plane, we need to find out where, and draw some sort of tiling pattern so that the sphere has something interesting to reflect. PLANE:: Pattern( ) takes a ray and a time to intersection and returns the appropriate color.
If the ray hits the sphere, we need to find the point of intersection with the sphere and create a ray with that as an origin, which is reflected about the normal line from the center of the sphere through that point of intersection. This is all done in SPHERE:: Reflect( ). If the ray then hits the plane, PLANE:: Pattern() colors the spot appropriately. If the reflected ray does not hit the plane, then the spot is colored blue so that we can distinguish it from the black of the background.
Note that we call it with a trailing 0 PLANE:: Pattern() here so that these points are darker than the others.
Possible Enhancements
Before you get into the heavy work of new algorithms, some possibilities for enhancement exist. You could make the program draw the screen in low resolution, working upward fairly easily, so that the form of the image is apparent after a few of the points, and you can tell if the objects you're looking for actually show up in your viewport without rendering the entire images. Or add another pattern to the plane or another sphere (to make things simple, non-reflective).
You could also calculate several images, save them on a fast disk or in memory, and flash them quickly to the screen for animation. There's still quite a bit of precision loss at the far ends of the rays and reflections. I've got a couple of ideas on ways to solve this but I want to play with the theory a bit first. Also, there are several places where calculations could be optimized (and quite a few places where they can be done in parallel, given the right hardware).
Several additions suggest themselves immediately, none terribly complex. Most notable are bounded planes, light sources (for shadowing and shading), and refraction. Refraction is one of the easiest, just a modification of the reflection routine. Consult any physics textbook for the equations.
Light sources are also straightforward. Start by defining a point as a light source, and at each intersection to be displayed, trace a ray to the light source. If the ray hits nothing color the pixel as lighted, otherwise it's dark. Using angles relative to the incident ray, glossy (but nonreflective) surfaces can be handled with displays that have enough color to give the subtleties of shading. The only problem with this method is that it doesn't allow for the diffraction of light around the edges of objects (soft edged shadows have got to be faked) and prisms, lenses, and focusing mirrors don't work on the light.
To do more serious modeling, the algorithms here need to be expanded to work with bounded planes, so that more complex shapes are possible. One way to do this is to bound each plane as a convex shape with other planes. I'm still playing with ideas for clean data structures for all of this.
The basic ray tracing technique has many applications, including video special effects and examination of architectural lighting problems -- special effects to commercials to architects figuring lighting problems. I'd love to see what others are coming up with and welcome your comments.
_RAY TRACING_ by Daniel Lyke
[LISTING ONE]<a name="01ee_000c">
/* 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 */
};
<a name="01ee_000d"><a name="01ee_000d">
<a name="01ee_000e">[LISTING TWO]<a name="01ee_000e">
/* 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();
}
<a name="01ee_000f"><a name="01ee_000f">
<a name="01ee_0010">[LISTING THREE]<a name="01ee_0010">
#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();
}
