This is the first part of the ray tracing series of tutorials. We've seen in the introduction what a raytracer was and how it was different from other solutions. Now let's concentrate on what it would requires to implement one in C/C++.

Our raytracer will have the following functionality. It will not try to be real time, and so will not take shortcuts because of performance only. Its code will try to stay as straightforward as possible, that means we'll try not to introduce complex algorithms. If we can avoid it. We'll make sure the basic concepts explained here will be visible in the code itself. The raytracer will work as a command line executable, that will take a scene file such as this one and output an image such as this one :

What does the raytracer do ?

First start with a simple pseudo language description of the basic ray tracing algorithm.

for each pixel of the screen
{
  Final color = 0;
  Ray = { starting point, direction };
  Repeat 
  {
    for each object in the scene
    {
      determine closest ray object/intersection;
    }  
    if intersection exists 
    {
      for each light in the scene
      {
        if the light is not in shadow of another object
        {
          add this light contribution to computed color;
        }      
      }
    }
    Final color = Final color + computed color * previous reflection factor;
    reflection factor = reflection factor * surface reflection property;
    increment depth;
  } until reflection factor is 0 or maximum depth is reached;
}  

For this first part, we'll limit ourselves to spheres as objects (a center point and a radius), and point lights as light sources (a center point and an intensity).

Scanning the screen

We have to cast a ray, at each one of our pixel to determine its color. This translates as the following C code :

 for (int y = 0; y < myScene.sizey; ++y) {
 for (int x = 0; x < myScene.sizex; ++x) {
 // ...
 }
 }

Well so far no difficulty. Once we computed the final color we'll simply write it to our TGA file as this :

  imageFile.put(min(blue*255.0f,255.0f)).
            put(min(green*255.0f, 255.0f)).
            put(min(red*255.0f, 255.0f));

Okay, a few comments on this. First it should be noted that all the colors inside our code will be treated as three components floating point numbers. We'll see the importance of that later. But the TGA file itself is storing its data using a fixed precision 8 bits per component color. So we'll need to do a conversion. In this first page, a naive "saturation operator", that is clamp any value above 1.0f, then convert the result into an integer number between 0 and 255. We'll also see later why the clamping is considered "naive". Please bear with us in the meantime.

Now casting rays for good

We've got our starting point. It's somewhere on our projection screen. We arbitrarily decide to make all our rays point in a single direction (towards the positive Z coordinate).

 ray viewRay = { {float(x), float(y), -1000.0f}, { 0.0f, 0.0f, 1.0f}};

The single direction is important because it defines our projection. Every ray will be perpendicular to the screen, that means that we'll do the equivalent of splatting our scene along parallel lines toward the screen. It's called perpendicular projection (also known as orthographic projection). There are other types of projection (conic projection being another big one, and we'll introduce it later).

Types of projection change the way the ray are related to each other. But it should be noted that it doesn't change anything to the behavior of a single ray. So you'll see that even later when we change our code to handle conic perspective, the rest of the code, intersection and lighting will NOT change the least bit. That's a nice property of raytracing compared to other methods, for example rasterization hardware relies heavily on the property of orthographic/conic perspective.

We set the starting point at an arbitrary position that we hope will enclose the whole visible scene. Ideally there should be no clipping, so the starting point will depend on the scene. But it should not be too far, as too big a floating point number will cause floating point precision errors in the rendering.

When we'll go "conic", such a problem of determination of the starting point will disappear because the positions behind the oberver are then ILLEGAL.

The closest intersection

The ray is defined by its starting point and its direction (normalized) and is then parameterized by a arithmetic distance called "t". It denotes the distance from any point of the ray to the starting point of the ray.

We only move in one direction on the ray from the smaller t to the bigger t. We'll iterate through each object, determine if there is an intersection point, we'll find the corresponding t parameter of this intersection and we'll find the closest one by taking the smallest t that is positive (that means it's in front of our starting point not behind).

Sphere-Ray intersection

Since we're only dealing with spheres our code only contains a sphere-ray intersection function. It is one of the fastest and simplest intersection code that exists. That's why we start with spheres here. The math themselves are not that interesting we'll introduce our parameter inside an intersection equation. That will give us a second degree equation with one unknown t.

We compute the delta of the equation and if the delta is less than zero there is no solution, so no intersection and if it's greater than zero there is at most two solutions. We take the closest one. If you want more details to how we got to this please consult our Answers section (What is the intersection of a ray and a sphere ?). Here is the code of that function :

    vecteur dist = s.pos - r.start; 
    double B = r.dir * dist;
    double D = B*B - dist * dist + s.size * s.size; 
    if (D < 0.0f) 
        return false; 
    double t0 = B - sqrt(D); 
    double t1 = B + sqrt(D);
    bool retvalue = false;  
    if ((t0 > 0.1f) && (t0 < t)) 
    {
        t = (float)t0;
        retvalue = true; 
    } 
    if ((t1 > 0.1f) && (t1 < t)) 
    {
        t = (float)t1; 
        retvalue = true; 
    }

We take the closest one, but it has to be further than the starting point, that is why we compare t0 and t1 to 0.1f. Why not 0.0f ? Simply because there is a high risk that given our limited precision (in float), that after a reflection from a point we find that our ray, intersects the same object around the same point when it shouldn't with an infinite precision. By taking our starting point at a reasonable distance but close enough to not cause "gaps" we can avoid some artifacts.

Lighting our intersection point

The lighting in one point is equal to the sum of the contribution of each individual light source. The problem here is way oversimplified by providing a finite list of well defined point lights.

We rapidly cull any light that is not visible from our side of the surface. That is if the exiting normal to our sphere and the light direction, point towards opposite direction. This is achieved by this simple code :

  // sign of dot product tells if they're pointing in the same direction
  if ( n * dist <= 0.0f )
    continue;

Then we determine if the point is in the shadow of another object for this particular light.

The code for this is similar to the intersection code that we had previously. We will cast a ray from our intersection point towards the light source (it's possible because it is a point light, that is a single point in space). But this time we don't care to find what the ray is intersecting or at what distance. As soon as we find an intersection point between this new ray and any object in the scene we stop the iteration because we know that the object is in the light shadow.

Lambert

Lambert was a European mathematician of the 18th century. Lambert cosine's law is a relation between the angle an incoming light ray has with a surface, the light intensity of this ray and the amount of light energy that is received by unit of surface. In short if the angle is low (grazing light), the received energy is low, and if the angle is high (zenith) the energy is high. That law will for example in an unrelated topic explain why the temperature at the equator is higher than the temperature at the earth poles.

How is that important to us ? If we are given a perfectly diffusing surface, that's the only thing that will matter to compute the lighting. The property of diffusion is that the energy that is received in one direction is absorbed and partially reemitted in all directions. A perfect diffuse surface will emit the same amount of energy in all directions. That important property is why the perceived lighting of a surface that is perfectly diffuse will only depend on the angle of the incident light with that surface : it is totally independent of the viewing angle. Such a material is also called a lambertian, because its lighting property is only determined by the lambert's law described above (that's the only type of material that we support in this first chapter. For more complex models see the following chapters). On the illustration below, you can see the incoming light and the redistributed energy, it is constant in all direction that's why it's represented in a perfect circle.

When we hit such a surface in our raytracer code, we will compute the cosine of the angle theta that the incoming ray does with the surface (via the normal) :

 float lambert = (lightRay.dir * n) * coef;

Then we multiply that lambertian coeficient with the diffuse color property of the surface, that will give us the perceived lighting for the current viewing ray.

Reflection

What would be ray tracing without that good old perfect (though irrealist) reflection. In addition to the diffusion color attribute described above, we'll give our surface a reflection coefficient. If this coefficient is 0, then that surface is not reflective at all and we can stop our iteration for this pixel and go to the next one. If this coefficient is greater than 0, then part of the lighting for that surface will be computed by launching a new ray from a new starting point at the incident point on that surface. The direction of the new ray, is the "reflected" direction, reflection of the incident view ray by the surface normal.

 
 float reflet = 2.0f * (viewRay.dir * n);
 viewRay.start = newStart;
 viewRay.dir = viewRay.dir - reflet * n;
 

We could potentially end up in a situation where each new ray would hit a surface and be reflected and so on forever. In order to avoid infinite loop, we arbitrarily limit the number of iterations by pixel (for example, no more than ten iterations in the current code). Of course when we hit that limitation the lighting will be wrong, but as long as we have no perfect reflection (less than a hundred percent) each iteration will contribute less than the previous one to the overall lighting. So hopefully cutting after the tenth iteration will not cause too many visual issues. We could alter that if we are to render a special scene that still requires more (a hall of mirror type of scene for example).

Digression : Computing successive reflections seems to be a natural fit for an recursive algorithm "from the book", but in practice and in our case, we preferred to go for an iterative algorithm. We save on the C stack usage but it will also require us to do some trade off in flexibility later.

The lighting models that we use here (Lambert + perfect reflection) are way too limited to render the complex natural materials. Other types of models (Phong, Blinn, Anisotropic, BRDF theory...) need to be introduced to account for more. We'll come back to it later.

In order to compile the source code for this first page, you don't need any particular additional library. You just need a C++ compiler coming with the standard C++ library. Tested with the GCC 3.4.4 and Visual Studio .net/2005/2008. Decompress the .rar file with Winrar.

Get the source code of the raytracer in C++. Source code for the first page only.

To page 2 : "Specularity, supersampling and post processing".