Introduction

Raytracing is an ever popular form of rendering, mainly due to its accurate modelling of phenomena such as reflections, refractions, specular highlights and so on. Its an elegant algorithm, providing simple solutions to many of the problems that plague traditional rendering methods. It also happens to be extremely slow. In this page, I'll run through a the basics of the algorithm, plus a little information on how to implement the simpler effects.

The Algorithm

The algorithm is pretty straightforward. The idea is that we trace rays through a scene, adjusting colour information and so on as we go. We fire a ray for each pixel, trace it through the scene, and evaluate its colour. Light bounces off primitives, possibly spawning new rays, shadows are checked for, lighting calculations are performed. Everything is handled it one simple, slow algorithm. Clipping, projection to screen, lighting, shadows, the lot.

Rays

The basic equation of a ray is:

Origin is the point of origin for the ray, not the world. Direction is a normalized vector, pointing in the direction of the ray. t is a parameter, specifying a point along the path of the ray. Sometimes the equation is written as:

Here, t is in the range 0-1. The direction vector is not normalized. Obviously here we can use this to determine if a point is on a ray by finding t, and checking to see if 0 <= t <= 1.

Shooting a ray for a given pixel is simple. The origin is your camera position, the direction is built from the screens pixel, and a Z component. The Z component seems to behave like d in the perspective projection equation. Anyway, I use the following:

    origin.x = origin.y = 0
    origin.z = -256
    dir.x = pixel.x
    dir.y = pixel.y
    dir.z = 256
    Normalize dir

Note that pixel must not be in range 0..width, 0..height, but in the range -width/2..+width/2, -height/2..+height/2. Otherwise your routine will not work correctly.

Intersecting Spheres

A common primitive in raytraced scenes is the sphere. Its pretty easy to code, looks nice, useful for debugging. The equation for a sphere is

Where X-Y-Z are a point on the sphere, and A-B-C is the centre of the sphere. What we are doing here is basically specifying a point of constant distance from another point, the centre. Any point on a sphere is always distance R from the centre. Now on to intersection. Firstly, substitute in the ray equation:

We need to re-arrange to get the equation in terms of t. You can do the maths yourself if you want, and I would recommend that you do, but for the more mathematically challenged amongst us, the solution is:

I laid the formula out so you can see it is a quadratic. A quadratic is a 2nd degree polynomial, Ax^2 + Bx + C. Now, we can find t by using the equations:

Discriminant =

And this equation can be simplified even more, because some of the coefficients work out to be 1. So we have:

The discriminant can be used to find information about the roots of the equation. If we call the discriminant 'd', then:

• d > 0, we have 2 different real roots
• d = 0, 2 equal real roots
• d < 0, roots are complex

We're only interested in real roots, so we can bail out if d < 0. If d = 0, then the ray just grazes the sphere. If d > 0 though, we have 2 intersections. -b - sqrt(d) gives the nearer one, -b + sqrt(d) gives the further one. So when you put all this into code, you will have to:

• Calculate d
• Check to see if d is legal
• If not, exit
• If so, find closer of two solutions
• Return that to some procedure

Once we have determined the closest intersection, we need the normal, for calculation. Again, simple enough for a sphere:

(Intersection - Centre) / Radius

The length of the vector is constant for a point on the sphere, because all points on a sphere are the same distance from the centre. So calculate this, feed it to a lighting formula.

Polygons

Polygons are a little more tricky. First, you'll need the plane equation, which is:

Ax + By + Cz + D = 0

Then, subsitute in the ray equation for x-y-z:

A(o.x + dir.x*t) + B(o.y + dir.y*t) + C(o.z + dir.z*t) + D = 0

Factorise for t:

t*(A*dir.x + B*dir.y + C*dir.z) + (A*o.x + B*o.y + C*o.z + D) = 0

And finally re-arrange for t:

This will give you a point on the plane. A denominator of 0 means the ray and plane are parallel - so chuck it out. Also if t < 0 then its behind the origin of the ray - chuck it out! A simple way of checking to see if the point is contained in the polygon is to then project down to 2-d, and do a 2-d point-in-polygon check. Use the ABC co-efficients to determine which axis to project down, if A is greatest, use X, if B is greatest, use Y and so on.

There is another algorithm however, but this works fore triangles. The idea is to use barycentric co-ordinates. I found this from an article in the Ray Tracing News, so credit goes to that author. Barycentric co-ordinates can map one triangle onto another, via a 1x1x1 cube. They are specified as 3 co-efficients, abc, rst, whatever. a+b+c always equals 1. From this, we can derive the following triangle:

The idea here is that we split the triangle into another 3 triangles. If the total area of the triangle is A, then triangle 1-2-P has area rA, triangle P-2-3 has area sA, triangle 1-P-3 has area tA. So the summation of the area must equal A. If it doesn't, then we're not in the triangle.

So we need to calculate rst. This is done like so:

N = surface normal (you should know how to do this)

Recalling that of course x = cross product, * = dot product

So now you can find the intersection. rst can also be used to interpolate the vertex normals, to find the normal at intersection. This is done like so:

NormalX = r*vnorm[0].x + s*vnorm[1].x + t*vnorm[2].x

NormalY = r*vnorm[0].y + s*vnorm[1].y + t*vnorm[2].y

NormalZ = r*vnorm[0].z + s*vnorm[1].z + t*vnorm[2].z

vnorm = array of vertex normals

So back to the polygons - you could test for containment in n - 2 triangles, where n is the number of points in the polygon. Here's some example C code I knocked up late one night that works for this:

int TestTriangle(VECTOR *tri1, VECTOR *tri2, VECTOR *tri3, VECTOR *point)
{
	VECTOR normal, spana, spanb, vec;
	double len, len2;
	double r, s, t;

	spana.x = tri2->x - tri1->x;
	spana.y = tri2->y - tri1->y;
	spana.z = tri2->z - tri1->z;
	spanb.x = tri3->x - tri1->x;
	spanb.y = tri3->y - tri1->y;
	spanb.z = tri3->z - tri1->z;
	
	normal.x = spana.y * spanb.z - spana.z * spanb.y;
	normal.y = spana.z * spanb.x - spana.x * spanb.z;
	normal.z = spana.x * spanb.y - spana.y * spanb.x;
	len = sqrt(normal.x * normal.x + normal.y * normal.y + normal.z * normal.z);

	spana.x = point->x - tri1->x;
	spana.y = point->y - tri1->y;
	spana.z = point->z - tri1->z;
	spanb.x = tri3->x - tri1->x;
	spanb.y = tri3->y - tri1->y;
	spanb.z = tri3->z - tri1->z;
	vec.x = spana.y * spanb.z - spana.z * spanb.y;
	vec.y = spana.z * spanb.x - spana.x * spanb.z;
	vec.z = spana.x * spanb.y - spana.y * spanb.x;
	len2 = sqrt(vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
	s = (vec.x * normal.x + vec.y * normal.y + vec.z * normal.z) / len;

	spana.x = tri2->x - tri1->x;
	spana.y = tri2->y - tri1->y;
	spana.z = tri2->z - tri1->z;
	spanb.x = point->x - tri1->x;
	spanb.y = point->y - tri1->y;
	spanb.z = point->z - tri1->z;
	vec.x = spana.y * spanb.z - spana.z * spanb.y;
	vec.y = spana.z * spanb.x - spana.x * spanb.z;
	vec.z = spana.x * spanb.y - spana.y * spanb.x;
	len2 = sqrt(vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
	t = (vec.x * normal.x + vec.y * normal.y + vec.z * normal.z) / len;

	r = 1.0 - (s + t);

	printf("(%f %f %f)\n", r, s, t);
	if ((r + s + t) == 1.0)
		return 1;
	else
		return 0;
}

Its not the most efficient piece of code, but it should provide a working example of this.

Reflections

These are pretty straight forward to do. You need to mirror the rays direction about the point of intersection. And intersection normal. The formula for this is:

2*(N.D)*N - D

Where N is the normal, and D is -direction. Note it must be -direction. This particular little catch up took a while to figure out ...Set your rays origin to be the point of intersection, and set the direction to be the above formula. Trace the ray, and multiply the returned colour by some factor, the reflectivity level. Then just add it to the colour. Simple enough.

Refractions

This works in a similar way to reflection, but its a lot harder and slower. To get the direction of refraction, we need to use the following formula:

Where:

T is refracted direction

n is eta, index of refraction

I is incident light

N is normal

I won't give the full working of the formula; you'll wake up in the middle of the night screaming. If you want to look it up, its on pages 757/758 of Foley & van Damm. You then have to do a similar thing to your reflections. Also note that there is a sqrt in that formula. You may get values put into the sqrt that are < 0. As you may now, you cannot sqrt a value < 0. To get around this, I use the following (C++) code:

Vector incident = -dir, tr;
Real n = medium / nearest->refindex;
Real c1 = -incident*inormal;
Real c2;

c2 = 1 - n*n*(1 - c1*c1);
if(c2 > 0.0)
{
        c2 = sqrt(c2);
        tr = n*incident + (n*c1 - c2)*inormal;
        tr.Normalize();

        Ray trans(ipoint, tr);
        transcol = trans.TraceRay(depth + 1, nearest->refindex);
        finalcol += transcol*nearest->trans;
}

Shadows

Shadows are very easily incorporated into a raytracer. A shadow, is a point that is blocked from a light source. So, all you have to do, is track a ray from the point of intersection, to the light source. Very simple. If a 100% opaque object intersects the ray, then simply set the colour to black. If no intersection, then the point is in light, so carry on as normal. However, if the ray intersects 1 or more (semi)transparent objects, then you will need to return some fraction 0 <= f <= 1, and multiply the colour by 'f'. A way to get soft shadows is by using distributed raytracing. However, I haven't tried this myself yet, so I can't report my findings.

Code layout

Raytracer code is usually structured recursively. This allows us to build a ray tree without additional coding. You may want to try avoiding the recursion, but you'd probably be better off spending your time doing more enjoyable things, such as sandpapering your teeth. Recursion is an elegant solution for this; despite possible stack problems I believe its the best solution. Your basic layout will look something like this:

Colour TraceRay(Ray ray, int depth)
{
        Find nearest object

        If found an intersection
                Shade and texture that object
                Calculate reflected ray
                If Depth < MaxDepth
                        Reflected-Colour = TraceRay(reflected-ray, depth+1)
                        Add Reflected-Colour*reflectivity to Colour
                End

                Calculate refracted ray
                If Depth < MaxDepth
                        Refracted-Colour = TraceRay(Refracted-ray, depth+1)
                        Add Refracted-Colour*refractivity to Colour
                End
        Else
                Colour = Black
        End

        Return Colour
}

This provides a clear simple solution; we use the MaxDepth parameter to bail out if its getting out of hand. Obviously, you put in a condition, and only spawn child rays if the object is reflective and if the object is refractive. Now begins the task of coding, and one hell of a load of optimizing.

[Ref: CGPP, Advanced Animation + Rendering]