Basics of Raytracing – Lighting, Reflections and Isosurfaces

This article contains the formulas you can use to add lighting, reflections and isosurfaces to your raytracer. It is assumed that you already have a simple raytracer running, and that you are familiar with some graphics programming terminology.

Throughout this article, I assume that all normal and direction vectors are of unit length, meaning (x * x + y * y + z * z) = 1. The ⋅ symbol means dot product. Cross products aren't used in this article.

Phong Illumination Model

The phong illumination model (also phong shading) can be used to calculate an approximation of the color at a given point, given some material and light properties. It is by no means physically correct, but it's fast, and it is by far the most popular model.

First you'll want to calculate the intensity of the diffuse reflection:

diffuse_intensity = max(0, light.direction ⋅ -point.normal)

Then comes the intensity of the specular reflection, the shiny bit:

half_vector = (light.direction + ray.direction).normalize() specular_intensity = pow(max(0, half_vector ⋅ -point.normal), material.shininess)

To calculate the total output color, you can use this formula:

output_color = light.ambient * material.ambient + light.diffuse * material.diffuse * diffuse_intensity + light.specular * material.specular * specular_intensity

The normalize() function converts a vector to unit length. The ray.direction variable is the direction of the ray that was traced to get to the point we're currently illuminating.

In the picture at the top of this article, the balls have diffuse colors between (0.1, 0.3, 0.7) and (0.3, 0.7, 0.9). The specular color is 0.5 times the diffuse color, and the shininess is 50. You can also see some transparency and reflection, which I will explain later.

Intersection Between a Ray and a Sphere

To find the intersection point between two geometric primitives (A and B), you need to create an equation A = B (where A and B are the equations defining the surface or interior of the objects), and solve it for a variable, for example time (t). If you do this for a ray and a sphere, solved using the ray's time, you get the following formula:

B = 2 * (ray.direction ⋅ (ray.origin - sphere.origin)) C = (ray.origin - sphere.origin).square() - sphere.radius * sphere.radius d = B * B - 4 * C if d < 0: No Collision t = (-B - sqrt(d)) * 0.5 if t < 0: t = (-B + sqrt(d)) * 0.5 if d < 0: No collision point = ray.origin + ray.direction * t normal = point - sphere.origin / sphere.radius

The square() function returns the square of the magnitude of a vector. When performing comparisons in this algorithm, you might want to use a small value like 1.0e-4 in place of 0 (zero). The reason is that you don't want to collide with a sphere that is at or near the origin of the ray, for example when you trace the reflection ray from the surface of a sphere.

Calculating the Reflection Vector

When you have traced a ray that ended up on a reflective surface, you want to start a new trace from that position in a new direction. This new direction is the reverse of the original direction, rotated around the normal at the intersection point (point), by 180 degrees. Here is the formula:

reflection.direction = ray.direction - point.normal * (ray.direction ⋅ point.normal * 2)

Intersection Between a Ray and an Isosurface

The following picture shows a visual representation of a 3D array consisting of 128x128x128 floating point values, where transitions from values below 0.5 to values above 0.5 are drawn as a surface.

Here is a 2D illustration of how isosurfaces work. The isoline consists of the points where the interpolated value is exactly 0.5.

To raytrace an isosurface, you go through the following steps:

1. Find the first point at which the ray is inside the volume
2. Find the value at this point, by interpolating the 8 closest values from the volume
3. If the value is below 0.5, make a small move in the direction of the ray, and jump back to step 2
4. We are now at the point where the ray hits the surface. Find the normal by calculating the derived value in each direction

Intersection Between a Ray and an Axis Aligned Box

far = HUGE_VAL near = -HUGE_VAL for i = 0 to 2: if ray.direction(i) == 0: # Ray is parallell to two sides of the box; # Must start between these to be inside if ray.origin(i) < box.min(i) or ray.origin(i) > box.max(i): No collision else: t1 = (box.min(i) - ray.origin(i)) / ray.direction(i) t2 = (box.max(i) - ray.origin(i)) / ray.direction(i) if t2 > t1: swap(t1, t2) if t2 > near: near = t2 if t1 < far: far = t1 if far < 0: No collision if near > far: No collision if near < 0: # Ray starts inside box t = 0 else: t = near

Again, instead of comparing to 0, you should compare to a small value, like 1.0e-4. Meaning instead of X == 0 you use X > -epsilon and X < epsilon.

Finding the Normal at a Point on an Isosurface

Say you have the function sample(x, y, z) which returns the value at the given position. To calculate the normal at a point [x, y, z], use the following formula:

normal_unnormalized = vector((sample(x - 1, y, z) - sample(x + 1, y, z)), (sample(x, y - 1, z) - sample(x, y + 1, z)), (sample(x, y, z - 1) - sample(x, y, z + 1))) normal = normal_unnormalized.normalize()

Here's another isosurface picture for your viewing pleasure (click to enlarge):