Kittler's Raytracer -- Code Critique

Author: Friedrich Kittler
Year: early 1980s-2000s

Hardware: Pentium IV,
x86 family of processors
x87 family of floating point coprocessors
Languages: C and Assembly

/*
v. 3.57
31.07.11
PTRACE.PAS (c't 1 / 93,167ff.) Extended
Super ellipsoid, rotational body and procedural textures from povray3
DOS version no longer supported

compile:
Normal: xgraf xsuptrace ray.s matrices.s
// Option -DNEWFRESNEL: simplified Fresnel light from povray3
Option -DJITTER: fuzzy (and time-consuming) shadows
Option -DGAMMA: gamma correction
SVGALIB or DGA: picture.i ".equ XWINDOW, 1" change!

RUNTIME:
<xsuptrace 1>: reproducible noise for runtime tests

FEATURES:
All reflective surfaces can be switched to ReflectionMapping: what
then appears on its surface is a picture to be loaded. This feature
But users have to demand standard interfaces first.

Prompts (':') for scalars and vectors can be defined either with <w [eiter]>
or <n [on]> acknowledge or answer with a new input

3 constant objects: sky, hell, ground (which only allows plein-air images)
Any number of variable objects (limited only by memory)
As variable objects by default 2 balls and 1 other object each predefined
ned. But standard objects can be deleted again.
If more than the standard number is required by an object, inquire
Prompt the new coordinates.
All objects editable (constant only by reassignment of color to normal)
procedure, variable also by location and size).
Some exotic objects are scalable and rotatable. This will be expanded.
Any number of lamps, the first 2 predefined.
Surfaces 1. global, 2. individually assignable after assignment to objects.
LINUX: any size * .24f images can be loaded as 2D textures.

Objects as a ring list for faster intershade () (see Foley, p.
Lamps as a simply linked list.
Left-hand coordinate system: left <x <right, front <y <behind,
below <z <above.

NEW:
01.04.97: AttCol () - acceleration away - transparencies would be worse
07.04.98: individual Fresnel coefficients with individual dullness for
opaque, but metallic surfaces (so these coefficients TransC and
to overwrite the dullness transparently). That's physically correct
and, mirabile dictu, better than POVRAY3. Color charts (ColTabT) global
edierbar
Nov. 22, 1998: ReflV () after Glassner, Image Synthesis, p. 726, again in
Light () calculated per lamp
24.12.98: Object slice selectable; Unzipped DOS - *. 24f files are loadable
08.03.99: Experimental support for Ohta / Maekawa algorithm
01.09.01: DOS no longer supported
13.08.04: Init_Ohta () different, untested
24.10.08: Steel new to /usr/ich/laptop/xsuptrace.c

BUGS:
Spheres of stei, sup, sor determined very empirically
2DMapping on Steiner, Agnesi, SuperQuadrik and Drehkoerper only as a reflection
Map implements: the transformation (x, y, z) -> (u, v) would be hard
MapProc () and Init_Fresnel () are not prepared for multiple calls from xsuptrace
Editing boxes and pyramids still uneven
lambda and thin globally, not variable per surface
No transparency shading as in CALCNEW.C
For floating objects, NULL pointer errors are inevitable; you
Change a midpoint coordinate to small amounts
SOR orb is detected in Edit_Sor () instead of Gravity ()
Refractive indices do not depend on wavelengths of light
08.04.99: intershade () returns L-> Shad correctly, but the object cohesion
breaks again and again, slows down more
30.12.00: xgraf (gcc with optimization) can falsify the SOR curve if
in SorInput the difference quotient Dy / Dx (ie the slope between
two x fixed points) becomes too large.
09.09.03: Changes to transformation matrices only apply to gcc -g. dark
31.07.11: QuColProc () dare

HINTS:
Between internal data structures and user displays, complex conversion
take place; So do not patch global data!
*/

#define  SUPTRACE
#define  COLTABSIZE   127       // unter DOS ggf. kleiner     
#define  PIII                   // bei schlechterer CPU dringend aendern!

// CompileTime: Globale Zaehler, bei neuen Objekten oder Oberflaechen erhoehen

#define SurfCnt       27        
#define FormCnt       16

#include <time.h>
#define  SPALTEN 640
#define  ZEILEN  480
#define  RAY
#include "xdefines.h"   // SVGALIB: defines.h 
#include "ray.h"        // SVGALIB: #include bild16m.c || lock16m.c 

extern float Infinit,halb,Epsilon;
VCT3  lambda = { 52.36, 56.6, 68.313 };
float duenne = 0.35, LampDiv = 1.;

// RunTime: Zaehler je Sitzung 
int ObjCnt  = 0;
int RingCnt = 0;
int LampCnt = 0;
int FormCnts[FormCnt] = {0};
int TestFormCnt = FormCnt-1;

#ifdef JITTER
float circle[19][2] = {
  {0.75, 0.433}, { 0.0,0.866},   {-0.75,0.433},{-0.75,-0.433 },{0.0,-0.866 },
  {0.75,-0.433}, { 0.0,0.0  },   { 0.75,0.0  },{0.375, 0.65  },{-0.375,0.65},
  {-0.75,0.0  }, {-0.375,-0.65}, {0.375,-0.65},{0.375,0.216  },{0.0,0.433  },
  {-0.375,0.217},{-0.375,-0.216},{0.0,-0.433 },{0.375,-0.217}};
float jitter_offset = 0.75/((float)RAND_MAX);
#endif

// global, um Debuggen und Fehler() zu erlauben
static uint x,y,RDepth = 0;
// global wegen ray.s
VCTd4  dquart;          
VCT4   quart;
VCT6   bicub;

#ifdef OHTA
typedef struct { VCT3 dist; float costheta; } dirT;
dirT    **direction;
int     lastring;
#endif

int   normal[FormCnt] = {0,1,4,7,2,11,10,7,6,11,3,14,12,5,2,23};
char  ObjStr[FormCnt][32] = {
  "den Himmel          ",
  "die Hoelle          ",
  "den Boden           ",
  "eine Kugel          ",
  "ein Moebiusband     ",
  "eine Steinerflaeche ",
  "einen Torus         ",
  "ein Ellipsoid       ",
  "einen Kegel         ",
  "eine Pyramide       ",
  "einen Quader        ",
  "eine Agnesi         ",             
  "ein Superellipsoid  ",
  "Rotationskoerper    ",
  "eine Scheibe        ",
  "eine Lemniskate     "
  };

FrameT Frame = {
  { 2.0, -8.0, 0.992},  // eye          
  { 0.1,  0.1, 0.1 },   // Ambient   
    0.96, 17,  5, 
    0.01, 0.0,          // AttenEps unbenutzt 
  { 0.50, 0.5, 0.6}, 0};// Fog_Color, Fog 

LampT *Lampen = NULL;
Prims *Boden,*Hoelle,*Himmel,*Ring,*Last;
int   arg = TRUE;

// VORDEFINIERTE OBJEKTE 

kugelT k1   =  {
  { 2.5, -2.7, 1.0 }, 1.0 };

kugelT k2   =  {
  {-0.2, -4.0, 1.3 }, 1.15 };

steiT  stei =  {
  {-2.0, -5.0, 1.0},
  { 1.0,  1.0, 1.0}};

moeT   moeb =  {
  {{ 2.4, 0.5, 4.0 }, 1.4 }, 
  0.25, 1 }; 

torT   tor  =  {
  { 2.9,-5.5, 0.6}, 0.3, 0.1 };

ellT   ell  =  {
  {3.8, -3.9, 1.4 },
  {0.3, -0.2, 0.6 },
  {0,0,0,0,0,0,0,0,0,0}};

coneT cone =  {
  {{ 1.4,-4.5,1.6 }, 0.29146 },
  0.3, 1.4, 0,1,20 }; 

pyrT pyr = {
  {{ 3,4,
  { 1.0, 1.0, 1.0 },    // LINKS - normal
    {{-0.9,-1.3, 2.2 }, // links 
    { 0.4,-2.5, 2.6 },  // Mitte 
    { 0.1,-0.5, 4.2 }}},// oben  
  { 3,4,
    { 1.0, 1.0, 1.0 },  // RECHTS 
    {{ 0.4,-2.5, 2.6 }, // Mitte  
     { 1.3, 0.2, 2.2 }, // rechts 
     { 0.1,-0.5, 4.2 }}},// oben   
  { 3,4,
    { 1.0, 1.0, 1.0 },   // UNTEN  
    {{ 1.3, 0.2, 2.2 },  // rechts 
    { 0.4,-2.5, 2.6 },   // Mitte  
    {-0.9,-1.3, 2.2 }}}, // links  
  { 3,4,
    { 1.0, 1.0, 1.0 },   // HINTEN 
    {{ 0.1,-0.5, 4.2 },  // oben   
    { 1.3, 0.2, 2.2 },   // rechts 
    {-0.9,-1.3, 2.2 }}}},// links  
  0, 0 };                // Polygon,dummy 

agnT agn = {{{-7.3,0.5,3.5},3.}, {-20, 0, 0}, 1.1 };

BoxInput B1 = { -1,1,3.0,1.5,-1,1 };

boxT box1; 

superT sup = {{ -1.5,-1.5,2.4 }, {2.0/0.7,0.7/0.6,2.0/0.6}, 
  {{-1,-1,-1},{1,1,1}},{1,1,1}};

sorInput sorI = {
  7,TRUE,
  {{ 1.80000,-0.100000},
    { 0.118143,0.020000},
  { 0.620253,0.540084},
    { 0.210970,0.827004},
    { 0.194093,0.962025},
  { 0.286920,1.035000},
    { 0.488354,1.150000}}};

discT scheibe = {{{ -1, 2, 2.7 }, 1.8 }, {0.4,-0.8,0.55}, 1, 0.7 };

lemniT lemnis = {{{ 4.1, -0.8, 1.85 }, 1.3 }, { 90., 0., 0. }};

// Vordefinierte Lampen 

LampT L1 = {
  {  4.0, -3.5, 4.5 },
  {  0.95, 0.8, 0.8 },
  {  1.0,  1.0, 1.0 },   
    0.6, NULL, NULL };
LampT L2 = {
  { -2.0, -3.5, 1.5 },  
  {  0.9,  0.9, 1.0 },
  {  1.0,  1.0, 1.0 },
    0.3, NULL, NULL };

// TEXTUREN 

uchar Lattice[LatCnt][LatCnt][LatCnt];

typedef struct { int Nr; VCT3 c; } ColTabT;
typedef VCT3 ColList[COLTABSIZE+1];

ColList HellTab   = {{ 0,0,0 }};
ColList FlameTab  = {{ 0,0,0 }};
ColList MarbleTab = {{ 0,0,0 }};

ColTabT MarbleCols[4]  = {
  {  0, { 0.41,0.25,0.15}},
  { 36, { 0.20,0.43,0.55}},
  { 66, { 0.91,0.93,0.88}},
  {127, { 0.64,0.66,0.55}}};

ColTabT HellCols[4]  = {
  {  0, { 0.7, 0.2, 0.1}},
  { 30, { 0.5, 0.8, 0.2}},
  { 60, { 0.9, 0.7, 1.0}},
  {127, { 0.7, 0.4, 0.0}}}; 

ColTabT FlameCols[6] = {
  {  0, { 1.00,1.00,1.00}},
  { 59, { 0.99,0.77,0.42}},
  { 63, { 0.75,0.80,0.51}},
  { 88, { 0.90,0.59,0.34}}, 
  {108, { 0.69,0.38,0.30}},
  {127, { 0.12,0.19,0.54}}}; 

#include "surfaces.h"

float  h,v;

/*
#ifdef DEBUG
VCT3   c;
ulong  nr_of_rays = 0; // 1199073 
ulong  nr_of_refl = 0; //  574234 
ulong  nr_of_refr = 0; //  438160 
ulong  nr_of_coh  = 0; //   66893
ulong  nr_of_rest = 0; //   42195
#endif
*/

void 
Translate_Object (Prims *Obj, VCT3 *t)
{
  TransformT trans;
  Compute_Translation_Transform (&trans, t);
  Compose_Transforms (&Obj->transform, &trans);
}

void 
Scale_Object (Prims *Obj, VCT3 *s)
{
  TransformT trans;
  Compute_Scaling_Transform (&trans, s);
  Compose_Transforms (&Obj->transform, &trans);
}

void 
Rotate_Object (Prims *Obj, VCT3 *r)
{
  TransformT trans;
  Compute_Rotation_Transform (&trans, r);
  Compose_Transforms (&Obj->transform, &trans);
}

void 
Fehler(void)
{
  blau();
  printf("%c%s",7,"Nullvektor!\n");
  printf("Spalte: %d, Zeile: %d, Tiefe: %d\n",x,y,RDepth);
  exit(1);
}

// RAYTRACER-ZUSATZ-FUNKTIONEN 

void 
FogProc (VCT3 *p,VCT3 *c)
{
  float factor;

  if ((factor = exp(-vdist(p,&Frame.eye)/Frame.Fog_Distance)) < 0.2)
    factor = 0.2;
  factor += 0.1 * turbulence(p);
  c->x = (c->x - Frame.Fog_Color.x) * factor + Frame.Fog_Color.x;
  c->y = (c->y - Frame.Fog_Color.y) * factor + Frame.Fog_Color.y;
  c->z = (c->z - Frame.Fog_Color.z) * factor + Frame.Fog_Color.z;
}

void    
Gravity(Prims *Obj) // ermittelt Schwerpunkt und von daher Umkugel 
{
  float  r,s,t;
  ellT   *ep;
  torT   *tp;
  moeT   *mp;  
  coneT  *cp;
  pyrT   *pp;
  steiT  *sp;
  agnT   *ap;
  boxT   *bp;
  superT *up;
  discT  *dp;
  lemniT *lp;
  int   i,j,ganz;             
  VCT3  d[2];

  switch(Obj->form)
  {
  case moebius:
    mp = Obj->UU.moeptr;
    r  = halb * mp->hoch * mp->hoch * (1 - M_SQRT_2);
    Obj->umkugel = mp->ax;
    Obj->umkugel.m.x += r;
    Obj->umkugel.m.y += r;
    Obj->umkugel.rad += mp->hoch * 2; 
    break;  
  case steiner:
    sp = Obj->UU.steiptr;
    Obj->umkugel.m = sp->m;
    r = fabs(sp->ax.x);
    s = fabs(sp->ax.y);
    Obj->umkugel.rad = t = fabs(sp->ax.z);
    if (r > s)
      {
      if (r > t)
        Obj->umkugel.rad = r;
      }
    else
      {
      if (s > t)
        Obj->umkugel.rad = s;
      }
    break;
  case torus:
    tp = Obj->UU.torptr;
    Obj->umkugel.m = tp->m;
    Obj->umkugel.rad = tp->a + tp->b;
    tp->a *= tp->a;
    tp->b *= tp->b;
    break;
  case ellipsoid:
    ep = Obj->UU.ellptr;
    ep->mat.a = r = 1./(ep->ax.x * ep->ax.x);
    ep->mat.e = s = 1./(ep->ax.y * ep->ax.y);
    ep->mat.h = t = 1./(ep->ax.z * ep->ax.z); // Hyperboloid = -h
    ep->mat.d = -ep->m.x*r;
    ep->mat.g = -ep->m.y*s;
    ep->mat.i = -ep->m.z*t;     // Hyperboloid = -i 
    ep->mat.j = (ep->m.x*ep->m.x*r + ep->m.y*ep->m.y*s
      + ep->m.z*ep->m.z*t)-1;
    Obj->umkugel.m = ep->m;
    r = fabs(ep->ax.x);
    s = fabs(ep->ax.y);
    Obj->umkugel.rad = t = fabs(ep->ax.z);
    if (r > s)
      {
      if (r > t)
        Obj->umkugel.rad = r;
      }
    else
      {
      if (s > t)
        Obj->umkugel.rad = s;
      }
    break;
  case kegel:
    cp = Obj->UU.coneptr;
    r  = cp->apex.rad * cp->apex.rad;
    Obj->umkugel.m = cp->apex.m;
    Obj->umkugel.m.z = halb*(r+1) * (cp->high+cp->low) - r * cp->apex.m.z;
    Obj->umkugel.rad = r * (cp->apex.m.z-cp->high) * (cp->apex.m.z-cp->high) 
      + (cp->high-Obj->umkugel.m.z) * (cp->high-Obj->umkugel.m.z);
    Obj->umkugel.rad = sqrt(Obj->umkugel.rad);
    break;
  case pyramide: 
    pp = Obj->UU.pyrptr;
    pp->planes = 4;
    for (i = 0; i < 4; i++)
      {
      pp->flaechen[i].vertnum = 3;
      pp->flaechen[i].polynum = 4;
      }
    ganz = pp->flaechen[0].vertnum * pp->planes;
    Obj->umkugel.rad = 0;
    Obj->umkugel.m.x = Obj->umkugel.m.y = Obj->umkugel.m.z = 0;
    for (i = 0; i < pp->planes; i++)
      for (j = 0; j < pp->flaechen[i].vertnum; j++)
        vaddeq(&Obj->umkugel.m,&pp->flaechen[i].vtx[j]);
    vscaleeq(&Obj->umkugel.m,1./(float)ganz);
    for (i = 0; i < pp->planes; i++)
      for (j = 0; j < pp->flaechen[i].vertnum; j++)
        if ((r = vdist(&Obj->umkugel.m,&pp->flaechen[i].vtx[j])) 
          > Obj->umkugel.rad)
         Obj->umkugel.rad = r;
    for (i = 0; i < pp->planes; i++)
      {
      vsub(d,&pp->flaechen[i].vtx[1],&pp->flaechen[i].vtx[0]);
      vsub(&d[1],&pp->flaechen[i].vtx[2],&pp->flaechen[i].vtx[1]);
      vcross(&pp->flaechen[i].normal,d,&d[1]);
      normalize(&pp->flaechen[i].normal);
      vsubnorm(d,&pp->flaechen[i].vtx[0],&Obj->umkugel.m);
      // gegen falsche Reihenfolge von Polygonecken-Eingaben
      negnorm(d,&pp->flaechen[i].normal);
      }
    break;
  case box:
    bp = Obj->UU.boxptr;
    bp->planes = 6;
    for (i = 0; i < 6; i++)
      {
      bp->flaechen[i].vertnum = 4;
      bp->flaechen[i].polynum = 6;    
      } 
    ganz = bp->flaechen[0].vertnum * bp->planes;
    Obj->umkugel.rad = 0;
    Obj->umkugel.m.x = Obj->umkugel.m.y = Obj->umkugel.m.z = 0;
    for (i = 0; i < bp->planes; i++)
       for (j = 0; j < bp->flaechen[i].vertnum; j++)
        vaddeq(&Obj->umkugel.m,&bp->flaechen[i].vtx[j]);
    vscaleeq(&Obj->umkugel.m,1./(float)ganz);
    for (i = 0; i < bp->planes; i++)
      for (j = 0; j < bp->flaechen[i].vertnum; j++)
        if ((r = vdist(&Obj->umkugel.m,&bp->flaechen[i].vtx[j]))
          > Obj->umkugel.rad)
          Obj->umkugel.rad = r;
    for (i = 0; i < bp->planes; i++)
      {
      vsub(d,&bp->flaechen[i].vtx[1],&bp->flaechen[i].vtx[0]);
      vsub(&d[1],&bp->flaechen[i].vtx[2],&bp->flaechen[i].vtx[1]);
      vcross(&bp->flaechen[i].normal,d,&d[1]);
      normalize(&bp->flaechen[i].normal);
      vsubnorm(d,&bp->flaechen[i].vtx[0],&Obj->umkugel.m);
      // gegen falsche Reihenfolge von Polygonecken-Eingaben 
      negnorm(d,&bp->flaechen[i].normal);
      }
    break;
  case agnesi:
    ap = Obj->UU.agnptr;
    Obj->umkugel.m = ap->ax.m;
    Obj->umkugel.rad = ap->scale*1.2;
    Create_Transform (&Obj->transform);
    Rotate_Object (Obj,&ap->rotate);
    d[0].x = d[0].y = d[0].z = ap->scale;
    Scale_Object(Obj,d);
    Translate_Object (Obj,&ap->ax.m);
    break;
  case super:
    up = Obj->UU.supptr;
    Obj->umkugel.m = up->m;
    if (up->scal.x > up->scal.y)
      {
      if (up->scal.x > up->scal.z)
        Obj->umkugel.rad = up->scal.x;
      else
        Obj->umkugel.rad = up->scal.z;
      }  
    else
      {
      if (up->scal.y > up->scal.z)
        Obj->umkugel.rad = up->scal.y;
      else
        Obj->umkugel.rad = up->scal.z;
      }
    Obj->umkugel.rad *= 1.7;   // groesser als theoretisches sqrt(2) 
    Create_Transform (&Obj->transform);
    Scale_Object (Obj,&up->scal);
    Translate_Object (Obj, &up->m); 
    break;
  case disc:
    dp = Obj->UU.discptr;
    Obj->umkugel = dp->ax;
    normalize (&dp->normal);
    dp->d = -vdot(&dp->ax.m,&dp->normal);   // v.Mangold, I,430  
    break;
  case lemni:
    lp = Obj->UU.lemniptr;
    Obj->umkugel = lp->ax;
    Create_Transform (&Obj->transform);
    Rotate_Object(Obj,&lp->rotate);       
    Translate_Object(Obj,&lp->ax.m);
    break;   
  default: ;
  }
}

// Spezielle Schnittpunkt-Routinen 

#define DEPTH_TOLERANCE 1.0e-3  // double 1.0E-4  
#define ZERO_TOLERANCE  1.0e-4  // double 1.0E-10 
#define SGNX(x) (((x) >= 0.) ? 1 : -1)
#define MAX_ITERATIONS  20
#define PLANECOUNT      9

static float planes[PLANECOUNT][4] =
   {{1, 1, 0, 0}, {1,-1, 0, 0},
    {1, 0, 1, 0}, {1, 0,-1, 0},
    {0, 1, 1, 0}, {0, 1,-1, 0},
    {1, 0, 0, 0}, {0, 1, 0, 0},
    {0, 0, 1, 0}};

char 
intersect_box (VCT3 *P, VCT3 *D, float *dmin, float *dmax)
{
  float tmin,tmax;

  // Links/rechts 
  if (fabs (D->x) > Epsilon)
    {
    if (D->x > Epsilon)
      {
      if ((*dmax = (1-P->x) / D->x) < Epsilon)
        return FALSE;
      *dmin = (-1-P->x) / D->x;
      }
    else
      {
      if ((*dmax = (-1-P->x) / D->x) < Epsilon)
        return FALSE;
      *dmin = (1-P->x) / D->x;
      }
    if (*dmin > *dmax) 
      return FALSE;
    }
  else
    {
    if (fabs(P->x) > 1)
      return FALSE;
    *dmin = -Infinit;
    *dmax =  Infinit;
    }
  tmin = tmax = 0.;
  // Oben/unten 
  if (fabs (D->y) > Epsilon)
    {
    if (D->y > Epsilon)
      {
      tmin = (-1-P->y) / D->y;
      tmax = (1-P->y) / D->y;
      }
    else
      {
      tmax = (-1-P->y) / D->y;
      tmin = (1-P->y) / D->y;
      }
    if (tmax < *dmax)
      {
      if (tmax < Epsilon) 
        return FALSE;
      if (tmin > *dmin)
        {
        if (tmin > tmax) 
          return FALSE;
        *dmin = tmin;
        }
      else
        {
        if (*dmin > tmax)
          return FALSE;
        }
      *dmax = tmax;
      }
    else
      {
      if (tmin > *dmin)
        {
        if (tmin > *dmax) 
          return FALSE;
        *dmin = tmin;
        }
      }
    }
  else
    {
    if (fabs(P->y) > 1)
      return FALSE;
    }
  // Vorn/hinten 
  if (fabs (D->z) > Epsilon)
    {
    if (D->z > Epsilon)
      {
      tmin = (-1-P->z) / D->z;
      tmax = (1-P->z) / D->z;
      }
    else
      {
      tmax = (-1-P->z) / D->z;
      tmin = (1-P->z) / D->z;
      }
    if (tmax < *dmax)
      {
      if (tmax < Epsilon) 
        return FALSE;
      if (tmin > *dmin)
        {
        if (tmin > tmax) 
          return FALSE;
        *dmin = tmin;
        }
      else
        {
        if (*dmin > tmax) 
          return FALSE;
        }
      *dmax = tmax;
      }
    else
      {
      if (tmin > *dmin)
        {
        if (tmin > *dmax) 
          return FALSE;
        *dmin = tmin;
        }
      }
    }
  else
    {
    if (fabs(P->z) > 1)
      return FALSE;
    }
  *dmin = tmin;
  *dmax = tmax;
  return TRUE;
}

So here's the code for the raytracer written by media philosopher Chapter 6 of Critical Code Studies. I'd like to invite you to an exploration of this portion. I can also make the full file available for exploration. I can also post the original German note. This version has been translated by the Hammermann family.

A few starter questions:
* What parts seem to be written by Kittler? And what parts by others?
* How does this raytracer code complement his theoretical writing?