# Ray tracing: the basics

304

## Introduction

Ray tracing is all about vector algebra and analytic geometry. Most of the time you do just one thing: solving a problem of intersection of a ray and a surface. Oh, and inventing various tricks to make this solution work faster.

Imagine a point in the center of your display in terms of XY-coordinates, and located a bit beyond of your display it terms of Z coordinate. Now cast an imaginary ray from this point through each pixel of your display. Now see where did every ray bump into and set its pixel to the color of the surface where the ray has stopped. Repeating this for every pixel you'll get some image. This image is a product of ray tracing technique.

The most simple way to start doing ray tracing is WebGL GLSL and that's why:

1. You don't need to mess with compilers. Everything works right in the web browser;
2. GPU is much powerful than CPU, and ray tracing requires a lot of calculations. Also rays can be processed in parallel, and that's what GPU does really good;
3. GPU is so good so you can even get smooth real-time animation at good FPS;
4. Fragment program works as an isolated ray processor, no need to mess with distracting details, so you can focus on the ray tracing math essentially.

## Into the math

So, we have a starting ray point, and we have the ray itself. Point is, well, a point, and the ray is a point plus a direction vector. So the ray equation looks like this:

$$\vec r = \vec p + \vec a t,$$

or, in coordinate form:

$$\left\{ \begin{array}{l} x = p_x + a_x t \\ y = p_y + a_y t \\ z = p_z + a_z t \end{array} \right.$$

Where $$\vec p$$ is a starting point, $$\vec a$$ is a direction vector, and $$t$$ is a non-negative parameter. Let's demand that $$\vec a$$ must be normalized. In this case $$t$$ is equal to the distance from $$\vec p$$ to $$\vec r$$.

In general any surface can be represented by the equation:

$$F(\vec r) = 0,$$

where $$\vec r$$ is some point and $$F(\vec r)$$ is some function. This equation means that if $$\vec r$$ belongs to the surface then $$F(\vec r) = 0$$. If $$F(\vec r) \ne 0$$ then $$\vec r$$ doesn't belong to the surface. For example, sphere equation looks like this:

$$|\vec r - \vec r_0| - R = 0,$$

where $$\vec r_0$$ is the center the sphere and $$R$$ is sphere radius.

Or, in coordinate form:

$$\sqrt{(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2} - R = 0$$

So putting $$\vec r$$ of the ray equation into the $$F(\vec r) = 0$$ surface equation we get an equation $$f(t) = 0$$. Solving this equation for a lowest positive $$t$$ we get a point where the ray intersects the surface for the first time.

The next step is to calculate a color of this point. First thing we need is to find normal vector of the surface in that point. Normal vector's direction is equal to the direction of the gradient vector of $$F(\vec r)$$ at that point. Gradient vector is a vector of partial derivatives:

$$\nabla F(\vec r) = \begin{bmatrix}\frac{\partial F(\vec r)}{\partial x} \\ \frac{\partial F(\vec r)}{\partial y} \\ \frac{\partial F(\vec r)}{\partial z}\end{bmatrix}$$

So unit normal vector is:

$$\vec n = \frac{\nabla F(\vec r)}{|\nabla F(\vec r)|}$$

## Distance fields

General surface equation $$F(\vec r) = 0$$ isn't very convenient to work with. There are numerical methods which can be used to solve the equation $$f(t) = 0$$, but in general, in case of general form of $$f(t)$$, this is a tough problem. To overcome this there is a special class of surface functions which is much more convenient for numeric methods — a distance field.

Distance field is a special case of surface function. In general, when $$F(\vec r) \ne 0$$, the value of $$F(\vec r)$$ doesn't have special meaning at this point $$\vec r$$. Distance field, on the other hand, does. Distance field is a function which value at the point $$\vec r$$ equals to the shortest distance to the surface from this point $$\vec r$$.

For example, sphere equation mentioned above is a distance field:

$$F(\vec r) = |\vec r - \vec r_0| - R,$$

however the same sphere surface in the more familiar form —

$$F(\vec r) = |\vec r - \vec r_0|^2 - R^2$$

— is not a distance field, because $$|\vec r - \vec r_0|^2 - R^2$$ doesn't equal to the shortest distance to the sphere from this point $$\vec r$$.

In case of distance field, numeric method of finding the intersection of a ray and a surface becomes fast and simple:

1. Calculate the value of $$F$$ at the current point;
2. Advance from the current point in the ray direction for the value of $$F$$. Since this value is the shortest distance, such an advancement won't go beyond the surface;
3. Repeat until intersection is considered as found or non-existing.

For many cases this method converges pretty fast.

Also since a distance field is a special case of surface function, you compute normal vector of distance field in the same way as you do with general surface function — by computing a vector of partial derivatives.

## Getting back to the code

Alright, done with math introduction, let's write some code. First — prepare an HTML page template and set up WebGL context:

<!DOCTYPE html>
<html>
<body>
<canvas id="canvas" width="1280" height="720"></canvas>
<textarea id="fs">
#version 300 es
precision highp float;

uniform float time;

in vec2 screen;

out vec4 frag_color;

void main(void)
{
vec3 p = vec3(0.0, 0.0, -1.0);
vec3 a = normalized(vec3(screen, 1.0));
vec3 rgb = vec3(0.0);

frag_color = vec4(rgb, 1.0);
}
</textarea>
<textarea id="vs">
#version 300 es
precision highp float;

uniform float aspect;

in vec2 position;

out vec2 screen;

void main(void)
{
screen = position * vec2(aspect, 1.0);
gl_Position = vec4(position, 0.0, 1.0);
}
</textarea>
<script src="script.js"></script>
</body>
</html

The only thing you're interested in here is the fragment program. All of the rest — the vertex program, WebGL context setting up, etc, are absolutely unrelated, so I've moved it into script.js. You will find a link to the full source code at the end of this article.

Let's take a closer look at the fragment program main function:

void main(void)
{
vec3 p = vec3(0.0, 0.0, -1.0);
vec3 a = normalized(vec3(screen, 1.0));
vec3 rgb = vec3(0.0);

frag_color = vec4(rgb, 1.0);
}

p and a are vectors from the ray equation $$\vec r = \vec p + \vec a t$$. rgb is a resulting color of the pixel. And that's it! Now we just need to add some surface, trace the ray, find the intersection, and calculate the color.

Alright, since we've talked about a sphere, let's implement a sphere. First, let's write a sphere distant field function, just before the main:

float F(vec3 r)
{
vec3 r0 = vec3(0.0, 0.0, 3.0);
float R = 1.0;

return distance(r, r0) - R;
}

Now we need a function which trace a ray and find an intersection of a ray and the surface:

float trace(vec3 p, vec3 a)
{
float eps = 0.001;
float t = 0.0;

/* If intersection wasn't found after 256 iterations
then it is considered that there is no intersection */
for(int n = 0; n < 256; ++n)
{
float dist = F(p + a * t);

if(dist < 0.0) // Is point beyond the surface?
{
return -1.0;
}
else if(dist < eps) // Is point considered as intersection?
{
return t;
}
else
{
t += dist;
}
}

return -1.0;
}

We have distace field function, we have tracing algorithm. Now let't use the algorithm to find the intersection of the current ray (current pixel):

void main(void)
{
vec3 p = vec3(0.0, 0.0, -1.0);
vec3 a = normalize(vec3(screen, 1.0));
vec3 rgb = vec3(0.0);

float t = trace(p, a);

if(t >= 0.0) // Intersection is found if t >= 0.0
{
rgb = vec3(1.0);
}

frag_color = vec4(rgb, 1.0);
}

And that's it! White pixel if intersection is found, black otherwise. Let's put everything together:

#version 300 es
precision highp float;

uniform float time;

in vec2 screen;

out vec4 frag_color;

float F(vec3 r)
{
vec3 r0 = vec3(0.0, 0.0, 3.0);
float R = 1.0;

return distance(r, r0) - R;
}

float trace(vec3 p, vec3 a)
{
float eps = 0.001;
float t = 0.0;

/* If intersection wasn't found after 256 iterations
then it is considered that there is no intersection */
for(int n = 0; n < 256; ++n)
{
float dist = F(p + a * t);

if(dist < 0.0) // Is point beyond the surface?
{
return -1.0;
}
else if(dist < eps) // Is point considered as intersection?
{
return t;
}
else
{
t += dist;
}
}

return -1.0;
}

void main(void)
{
vec3 p = vec3(0.0, 0.0, -1.0);
vec3 a = normalize(vec3(screen, 1.0));
vec3 rgb = vec3(0.0);

float t = trace(p, a);

if(t >= 0.0) // Intersection is found if t >= 0.0
{
rgb = vec3(1.0);
}

frag_color = vec4(rgb, 1.0);
}

And see the result:

Not too exciting yet, however it works, and this is not just a circle, but a sphere produced by pure ray tracing algorithm. To make it look 3D we should add some lighting. And we can do that in a couple of easy steps.

First, let's implement normalized gradient function of F. Remember your days back in school how you calculated partial derivatives. Note that we don't need to divide everything by eps since we normalize the vector anyway:

vec3 dF(vec3 r)
{
float eps = 0.0001;
float f = F(r);

return normalize(vec3
(
F(r + vec3(eps, 0.0, 0.0)) - f,
F(r + vec3(0.0, eps, 0.0)) - f,
F(r + vec3(0.0, 0.0, eps)) - f
));
}

Next, let's add diffuse lighting. Diffuse color in a simple way is calculated like that:

$$C_d = C_m C_l \frac{\vec l \cdot \vec n}{A(|\vec r - \vec l|)},$$

where $$C_m$$ is material color, $$C_l$$ is light source color, $$\vec l \cdot \vec n$$ is a dot product of normalized light vector and surface normal vector, and $$A(|\vec r - \vec l|)$$ is an attenuation function which determines how fast illumination should decrease with distance growing between the point and the light source. Let's implement diffuse function in GLSL. To keep thing simple we'll use just one color for now:

vec3 diffuse(vec3 point, vec3 light, vec3 normal)
{
vec3 color = vec3(0.567, 1.0, 0.0);
return color * dot(normalize(light - point), normal) / attenuation(distance(light, point));
}

For the attenuation we'll use quadratic attenuation principle:

float attenuation(float dist)
{
return dist * dist / 16.0 + 1.0;
}

Now let's modify main function:

if(t >= 0.0) // Intersection is found if t >= 0.0
{
vec3 r = p + a * t; // Intersection point
vec3 n = dF(r); // Surface normal vector at the interserction point
vec3 light = vec3(-1.0, 1.0, 0.0); // Light position
rgb = diffuse(r, light, n);
}

And see the result:

Much better, huh?

## More spheres!

To make it even more interesing, let's clone the sphere into the endless grid of spheres. Cloning an object in ray tracing is an easy trick. Assume you have some linear function $$f(x)$$ and you want to loop it around some $$x_0$$ with raduis equal to $$R$$. This can be done as following:

$$g(x) = mod(f(x - x_0) - \frac{R}{2}, R) - \frac{R}{2}$$

So let's implement it in GLSL:

vec3 repeat(vec3 v, vec3 v0, vec3 radius)
{
}

And adjust our F distance field function:

float F(vec3 r)
{
vec3 r0 = vec3(1.5, 1.5, 3.0); // We move it away from the screen center
r = repeat(r, r0, vec3(4.0)); // And loop r coordinates around the sphere center
float R = 1.0;

return distance(r, r0) - R;
}

Such a simple trick, and just look at the result!

## Shininess effect

Let's add one more thing: a shininess (or a specular) effect. A shininess strength could be calculated like that:

$$S = k(\vec h \cdot \vec n)^P,$$

where $$k$$ is strength factor, $$P$$ is specular power which determines how focused is the reflection, $$\vec n$$ is surface normal vector and $$\vec h$$ is so-called halfway vector. Halfway vector is normalized vector which is bisector between direction from the surface point to the light source position and from the surface point to the observer's eye. Let's implement it in GLSL:

float specular(vec3 point, vec3 light, vec3 eye, vec3 normal)
{
vec3 halfway = normalize(normalize(light - point) + normalize(eye - point));
return pow(dot(halfway, normal), 32.0) / 2.0;
}

And modify our main function:

if(t >= 0.0) // Intersection is found if t >= 0.0
{
vec3 r = p + a * t; // Intersection point
vec3 n = dF(r); // Surface normal vector at the interserction point
vec3 light = vec3(0.0, 0.0, -1.0); // Light position
rgb = diffuse(r, light, n) + specular(r, light, p, n); // Add specular to diffuse
}

And here we go:

And the last thing for this part — let's add some animation. To do so we'll use time value as Z-coordinate of the sphere:

float F(vec3 r)
{
vec3 r0 = vec3(1.5, 1.5, -mod(time, 3.0));
r = repeat(r, r0, vec3(3.0));
float R = 1.0;

return distance(r, r0) - R;
}

If everything is correct then you should see:

Nice and smooth. And no polygons at all.

Next time we'll implement more visual effects. As for now, take a break and play with the code to see how it works.

Full source code for this article: source.zip. Unzip the archive and open page.html in your web browser.

Rate this post: