/* 
(Transliteration of
<http://canonical.org/~kragen/sw/aspmisc/my-very-first-raytracer.html>,
which has the images.)

![image](http://canonical.org/~kragen/sw/aspmisc/scene.jpg)

I ran across [Sebastian Sylvan's post on raytracing using signed
distance
functions](http://sebastiansylvan.wordpress.com/2009/07/04/ray-tracing-signed-distance-functions/)
the other day, and it occurred to me that I'd never actually written a
raytracer.

So here's a simple raytracer (not using SDFs, just the usual basic
approach, taking dot products and stuff). The [scene
description](http://canonical.org/~kragen/sw/aspmisc/scene.txt) lists
five spheres:

    5
    x=1.2 y=0 z=3     r=.5 g=.25 b=0
    x=0 y=0 z=4       r=1  g=.5  b=0  refl=0.8 rad=1.2
    x=0 y=1.1 z=4     r=1  g=.5  b=1  refl=0.1 rad=1
    x=-.5 y=-.5 z=2.5 r=.2 g=.2  b=.2 refl=1 (silver) rad=.2
    y=-1  r=1 g=1  texture=1 refl=0.1

You can see that the small (rad=0.2) silver (refl=1) one is reflected in
the large reflective one.

It generates this image in about a third of a second. The [source is
here.](http://canonical.org/~kragen/sw/aspmisc/raytracer.c).

Isn't it *precious?*

It's a tiny program, four pages of code. Some notable **things it
doesn't contain**: matrices, transcendental functions, heap allocation,
multiple kinds of geometry, or rotation. (It does invoke `sqrt`,
`floor`, and `pow`, and use GCC's variable-sized arrays.). You could
probably run it on an Arduino, especially if you rewrote it to use
fixed-point math.

It *does* contain a tiny parser for a key-value-pair-based regular
language in which objects inherit from previous objects (32 lines), a
TrueColor PPM-format output file generator (12 lines), and a couple of
volumetric generative textures (11 lines). Toys, but fun ones!

More images below:

* * * * *

![image](http://canonical.org/~kragen/sw/aspmisc/gray.jpg)

I was talking with Brandon Moore about writing raytracers, and I was
struck by the calm beauty of the test image from [his test raytracer in
Haskell](https://github.com/bmmoore/ray/blob/master/Ray.hs). This was as
close as I could come with mine, because mine doesn't support changing
the direction of light or the background "texture", since both are quick
one-line hacks. (You could say the same of his, but they're *different*
one-line hacks.) Worse, mine doesn't support turning off specular
highlights.

    5
    z=15 r=1 g=1 b=1 rad=10   A big dull sphere for camera background.
    z=0  rad=40               And another for behind the camera.
    z=1  rad=.33 refl=.9      A reflective sphere in front.
    x=.33 rad=.17  r=0 g=0 b=0  Another to the side.
    x=0 y=-.67 z=0 rad=.33  r=.5 g=.5 b=.5  And another we see only in 
reflection.

* * * * *

![image](http://canonical.org/~kragen/sw/aspmisc/bigfield.jpg)

I was looking at [RTRT3, Mark VandeWettering's tiny realtime raytracer
from
2000](http://brainwagon.org/2011/04/12/code-from-the-past-a-tiny-real-time-raytracer-from-2000/),
and I was curious how performance compared between the two programs. So
I added this code to RTRT3 and found that it was getting 6 or 7 frames
per second on my netbook, compared to 96fps on his 2011 workstation:

    #include <sys/time.h>

    struct timeval old_tv;

    static void
    fps()
    {
        struct timeval new_tv, diff;
        double fps;
        gettimeofday(&new_tv, NULL);
        timersub(&new_tv, &old_tv, &diff);
        old_tv = new_tv;
        if (diff.tv_sec) return; // your system is too SLOW
        fps = (1000*1000.0)/diff.tv_usec;
        printf("fps: %.1lf\n", fps);
    }

But to do an apples-to-apples comparison, I needed a scene of similar
complexity: a sphere bouncing over an infinite textured plain. He also
implemented shadows, which I'm still scared of. Turns out that My Very
First Raytracer got about 4fps to RTRT3's 7fps, despite rendering a
smaller canvas, but perhaps losing some time in generating the output
file.

    2
    y=-1.1  r=1 g=.5 b=.25  refl=.5
    r=1 g=1 b=1  z=0 y=100000 rad=99999.7 refl=0  texture=1

The complex moiré aliasing in the reflected texture is *totally
tubular*.

Mark's code makes me want to hook My Very First Raytracer up to GLUT
too, in order to make it interactive.

* * * * *

![image](http://canonical.org/~kragen/sw/aspmisc/bbtext.jpg)

Another simple volumetric procedural texture, because checkerboards are
boring.

The source code for this texture is just:

        else if (obj->ma.texture == 'b') {
          value = (((int)(point.x*128) ^ (int)(point.y*64) ^ (int)(point.z/8)) 
               & 255) / 255.0;
        }

and the scene is composed as follows:

    8
    x=-1 y=-1  texture=b  r=1 g=1 b=1  refl=.5
    y=1 b=.5
    x=1 g=.5
    y=-1 r=.5
    z=4 g=1
    y=1 b=1
    x=-1 r=0
    z=4.5 x=0 y=0 texture=0 refl=1 r=.2 g=.2 b=.2 rad=.2

*/

// My First Ray Tracer
// One object type: a sphere.  Scalar reflectivity.  No textures.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

enum { debug_distance = 0, debug_reflection = 0 };



/* Vectors. */

typedef float sc;               // scalar
typedef struct { sc x, y, z; } vec;

static sc
dot(vec aa, vec bb) { return aa.x*bb.x + aa.y*bb.y + aa.z*bb.z; }

static sc
magsq(vec vv) { return dot(vv, vv); }

static vec
scale(vec vv, sc c) { vec rv = { vv.x*c, vv.y*c, vv.z*c }; return rv; }

static vec
normalize(vec vv) { return scale(vv, 1/sqrt(dot(vv, vv))); }

static vec
add(vec aa, vec bb) { vec rv = { aa.x+bb.x, aa.y+bb.y, aa.z+bb.z }; return rv; }

static vec
sub(vec aa, vec bb) { return add(aa, scale(bb, -1)); }


/* Ray-tracing proper. */

typedef vec color;              // So as to reuse dot(vv,vv) and scale
typedef struct { color co; sc reflectivity; char texture; } material;
typedef struct { vec cp; material ma; sc r; } sphere;
typedef struct { sphere *spheres; int nn; } world;
typedef struct { vec start; vec dir; } ray; // dir is normalized!
typedef int public_static_void, bool;


/* If there are intersections, store the parametric values t for which
 * rr.start + t * rr.dir is an intersection into `intersections` and 
 * return true; otherwise return false. */
static bool
find_intersections(ray rr, sphere ss, sc *intersections)
{
  vec center_rel = sub(rr.start, ss.cp);
  // Quadratic coefficients of parametric intersection equation.  a == 1.
  sc b = 2*dot(center_rel, rr.dir), c = magsq(center_rel) - ss.r*ss.r;
  sc discrim = b*b - 4*c;
  if (discrim < 0) return 0;
  intersections[0] = (-b - sqrt(discrim))/2;
  intersections[1] = (-b + sqrt(discrim))/2;
  return 1;
}

static ray
reflect(ray rr, vec start, vec normal)
{
  // Project ray direction onto normal
  vec proj = scale(normal, dot(rr.dir, normal));
  // Subtract that off twice to move the ray to the other side of surface
  vec reflected_dir = sub(rr.dir, scale(proj, 2));
  ray rv = { add(start, scale(reflected_dir, 0.001)), reflected_dir };
  return rv;
}

static sc
sc_max(sc a, sc b) { return a > b ? a : b; }

static color
clamp(color co)
{
  color rv = { sc_max(co.x, 0), sc_max(co.y, 0), sc_max(co.z, 0) };
  return rv;
}

static color
grey(sc value) { color rv = { value, value, value }; return rv; }

static color trace(world here, ray rr, sc importance);

static color
surface_color(world here, sphere *obj, ray rr, vec point, sc importance)
{
  vec light = { .58, -.58, -.58 };
  vec normal = normalize(sub(point, obj->cp));
  sc diffuse_intensity = dot(light, normal);
  sc refl_importance = importance * obj->ma.reflectivity;
  ray refl = reflect(rr, point, normal);
  color diffuse_color = obj->ma.co;

  if (obj->ma.texture) {
    sc value = 1.0;
    if (obj->ma.texture == '1') {
      sc d = dot(point, light), df = d*32 - floor(d*32);
      value = df * (1-df) * 4;
    } else if (obj->ma.texture == 'b') {
      value = (((int)(point.x*128) ^ (int)(point.y*64) ^ (int)(point.z/8)) 
               & 255) / 255.0;
    }
    diffuse_color = scale(diffuse_color, value);
  }

  color rv = scale(scale(diffuse_color,
                         diffuse_intensity < 0 ? 0 : diffuse_intensity),
                   1 - obj->ma.reflectivity);
  sc specular_cos = dot(refl.dir, light);

  if (specular_cos > 0.9)
    rv = add(rv, scale(grey(pow(specular_cos, 64)), obj->ma.reflectivity));

  if (refl_importance > 0.01) {
    if (debug_reflection) {
      color rv2 = {1, 0, 0};
      return rv2;
    }
    rv = add(rv, clamp(scale(trace(here, refl, refl_importance),
                             obj->ma.reflectivity)));
  }

  return rv;
}

static color
trace(world here, ray rr, sc importance)
{
  int ii;
  vec crap = { 0, 0, -0.5 };
  sc intersections[2];
  sc nearest_t = 1/.0;
  sphere *nearest_object = NULL;

  for (ii = 0; ii < here.nn; ii++) {
    if (find_intersections(rr, here.spheres[ii], intersections)) {
      int jj;
      for (jj = 0; jj < 2; jj++) {
        if (intersections[jj] < 0 || intersections[jj] >= nearest_t) continue;

        nearest_t = intersections[jj];
        nearest_object = &here.spheres[ii];
      }
    }
  }

  if (debug_distance && nearest_object) return grey(nearest_t-floor(nearest_t));

  if (nearest_object) {
    return surface_color(here, nearest_object, rr,
                         add(rr.start, scale(rr.dir, nearest_t)),
                         importance);
  }

  return normalize(add(crap, rr.dir));
}

static color
pixel_color(world here, int ww, int hh, int xx, int yy)
{
  vec pv = { (double)xx/ww - 0.5, (double)yy/hh - 0.5, 1 };
  ray rr = { {0}, normalize(pv) };
  return trace(here, rr, 1.0);
}


/* PPM P6 file format; see <http://netpbm.sourceforge.net/doc/ppm.html> */

static void
output_header(int ww, int hh)
{ printf("P6\n%d %d\n255\n", ww, hh); }

static unsigned char
byte(double dd) { return dd > 1 ? 255 : dd < 0 ? 0 : dd * 255 + 0.5; }

static void
encode_color(color co)
{ putchar(byte(co.x)); putchar(byte(co.y)); putchar(byte(co.z)); }


/* Input file format */

static char *
kv(char *line, char *name)      // name is up to 61 bytes
{
  char buf[64] = " ", *v = strstr(line, strcat(strncat(buf, name, 61), "="));
  return (v ? v + strlen(buf) : NULL);
}

static void
parse_input(world here)
{
  char line[512] = " ";
  // Defaults:
  sphere s = { .ma = { .co = {0, .5, 0} }, .r = 1, .cp = {0, 0, 5} };

  int ii;
  for (ii = 0; ii < here.nn; ii++) {
    fgets(line + 1, sizeof(line) - 1, stdin);

    #define field(name, where)                                          \
      do { char *v = kv(line, #name); if (v) s.where = atof(v); } while(0)

    field(x, cp.x); field(y, cp.y); field(z, cp.z); field(rad, r);
    field(r, ma.co.x); field(g, ma.co.y); field(b, ma.co.z);
    field(refl, ma.reflectivity);
    { char *v = kv(line, "texture"); if (v) s.ma.texture = *v; }

    #undef field

    here.spheres[ii] = s;
  }                              
}


/* High-level flow */

static void
render_pixel(world here, int ww, int hh, int xx, int yy)
{ encode_color(pixel_color(here, ww, hh, xx, yy)); }

static void
render(int nobjects, int ww, int hh)
{
  int ii, jj;
  sphere spheres[nobjects];
  world here = { spheres, nobjects };
  parse_input(here);

  output_header(ww, hh);
  for (ii = 0; ii < hh; ii++)
    for (jj = 0; jj < ww; jj++)
      render_pixel(here, ww, hh, jj, ii);
}

static int
usage(char **argv)
{
  fprintf(stderr, "Usage: %s width height < scene.txt > scene.ppm\n", argv[0]);
  return -1;
}


public_static_void
main(int argc, char **argv)
{
  int nobjects;
  if (argc < 3) return usage(argv);
  if (!scanf("%d\n", &nobjects)) usage(argv);
  render(nobjects, atoi(argv[1]), atoi(argv[2]));
  return 0;
}


-- 
To unsubscribe: http://lists.canonical.org/mailman/listinfo/kragen-hacks

Reply via email to