Is this the most effective and efficient programming style for prototyping (a ray tracer) in Mathematica?

Question

I am trying out Mathematica as a prototyping tool. As an initial exercise I have put together a brute force ray tracer, with a view to using Mathematica's built in probability distribution functions to test out some fancy shaders. (My starting point was the F# ray tracer available http://code.msdn.microsoft.com/ParExtSamples.)

As this is a first attempt at using Mathematica, other than as a fancy calculator, I would welcome some guidance or critique on whether this coding style is going to be effective. (You will detect from the code that I am accustomed to an object-oriented paradigm.)

My own observations, in no particular order, are:

1. It's an enjoyable way to work, as a lot of the busy type definitions and additional syntax required by other languages is not needed here; the code is quite compact and the intention fairly readable (IMHO).
2. I have built the ray tracer from scratch, without using the built-in graphics primitives as that would defeat the learning purpose of the exercise and I would probably not be able to roll my own shaders. As a simple example, addition is not defined for RGBColor, for example.
3. Probably as a consequence, the tracer is very slow, even if I use all the 4 kernels that are available to me. I dare say that the existing code could be speeded up considerably by removing all the type / pattern matching, but the code would be much harder to follow, I think, and certainly much harder to debug. My approach seems to mean that I cannot take advantage of Compile optimization, so far as I can see. If I wanted a fast ray tracer I would write it in C++, but I wonder whether there are any easy optimizations that I have missed (and that don't involve mangling what I have to the extent that the prototyping benefits of Mathematica -- easy refactoring -- would be lost). However, I am surprised that even with 4 kernels, the cpu of this 4-core (8-core, if you count hyperthreading) never maxes out. For example, is it worth changing some of the := to =? Things would be even slower if I aliased the results by averaging 4 adjacent traces, for example. Real time ray tracing seems out of reach.
4. Ray tracing results in a lot of "corner cases" that are dealt with naturally with IEEE maths. Unfortunately Mathematica does produce +/-Infinity for +/-1/0, for example, so some of the code should really be extended to treat those cases properly.
5. It would be great if Mathematica had more built-in vector algebra so that I could write the equations defining the objects and the rays involved and getting Mathematica to calculate the ray intersection points that are at the heart of the ray tracer. As things stand, Reduce and Solve have not helped me to find better intersection algorithms, producing either nothing at all or something large an unintelligible, depending on how I posed the problem.

Anyway, here is what you get after 110s from raytrace[400, 400, basescene, 6]:

... and here is the code:

    (* Colour helpers *)
    black = {0., 0., 0.};
    darkgrey = {.25, .25, .25};
    grey = {.5, .5, .5};
    white = {1., 1., 1.};
    background = black;
    defaultcolor = black;

    brightness[{r_, g_, b_}] = Mean[{r, g, b}]; 

    scale[k_, c: {r_, g_, b_}] = k * c;

    zero = {0.,0.,0.};

    (* Mainly for reference; pattern matching normally used instead *)

    ray /: start[ray[s_, d_]] = s;
    ray /: dir[ray[s_, d_]] = d;

    camera /: pos[camera[p_, l_]] = p;
    camera /: lookat[camera[p_, l_]] = l;
    camera /: forward[camera[p_, l_]] = Normalize[l - p];
    camera /: down[camera[p_, l_]] = {0., -1., 0.};
    camera /: right[c : camera[p_, l_]] := 1.5 * Normalize[Cross[forward[c], down[c]]];
    camera /: up[c : camera[p_, l_]] := 1.5 * Normalize[Cross[forward[c], right[c]]];


    light /: pos[light[p_, c_]] = p;
    light /: color[light[p_, c_]] = c;

    scene /: things[scene[t_, l_, c_]] = t;
    scene /: lights[scene[t_, l_, c_]] = l;
    scene /: camera[scene[t_, l_, c_]] = c;

    surface /: diffuse[surface[d_, s_, re_, ro_]] = d;
    surface /: specular[surface[d_, s_, re_, ro_]] = s;
    surface /: reflect[surface[d_, s_, re_, ro_]] = re;
    surface /: roughness[surface[d_, s_, re_, ro_]] = ro;

    intersection /: thing[intersection[t_, r_, d_]] = t;
    intersection /: ray[intersection[t_, r_, d_]] = r;
    intersection /: dist[intersection[t_, r_, d_]] = d;
    miss = intersection[nothing, ray[zero, zero], Infinity];

    sceneobject /: surface[sceneobject[s_, i_, n_]] = s;
    sceneobject /: intersect[sceneobject[s_, i_, n_]] = i;
    sceneobject /: normal[sceneobject[s_, i_, n_]] = n;

    sphere /: center[sphere[c_, r_, s_]] = c;
    sphere /: radius[sphere[c_, r_, s_]] = r;
    sphere /: surface[sphere[c_, r_, s_]] = s;
    normal[sphere[center_, _, _], pos_] = Normalize[pos - center];

    plane /: normal[plane[n_, o_, s_]] = n;
    plane /: offset[plane[n_, o_, s_]] = o;
    plane /: surface[plane[n_, o_, s_]] = s;
    normal[plane[n_, _, _], _] = n;

    (* Axis-aligned bounding box *)
    (* TODO: not yet used; integrate into tracer *)
    box /: lowerb[box[l_, u_]] := l;
    box /: upperb[box[l_, u_]] := u;

    extendby[box[l_, u_], pt_] :=
        box[MapThread[Min, {l, pt}], MapThread[Max, {u, pt}]];
    size[box[l_, u_]] :=
        u - l;
    majoraxis[b : box[l_, u_]] :=
        Ordering[size[b], -1]; 


    (* TODO: This does not work for cases where dir has 0 compnent as Mathematic returns ComplexInfinity,
       not +/-Infinity for +/-1/0 *)
    intersectboxQ[b : box[l_, u_], r : ray[start_, dir_]] :=
        Module[ {tl = (l - start) / dir, tu = (u - start) / dir, tmin, tmax},

            (* Swap u[i] and l[i] if dir[i] < 0 to avoid 
               erroneous result because 0 == -0 *)
            tmin = Max[MapThread[Min, {tu, tl}]];
            tmax = Min[MapThread[Max, {tu, tl}]];
            Not[tmax < 0 && tmin > tmax];
        (* Use Not to cover some Infinity comparisons *)

        (*  Which[
            tmax < 0, False, (* Intersection at t = tmax, but it's behind us *)
            tmin > tmax, False, (* No intersection *)
            True, True (* Interesection at t = tmin *)
            ]*)
        ];



    intersect[s : sphere[center_, radius_, _], r : ray[start_, dir_], i : intersection[_, _, currentdist_]] :=
        Module[ {eo = center-start, v, dist, disc},
            v = eo.dir;
            dist = If[ v < 0.,
                       0.,
                       disc = radius * radius - (eo.eo - v * v);
                       If[ disc < 0.,
                           0.,
                           v - Sqrt[disc]
                       ]
                   ];
            If[ dist == 0. || dist > currentdist,
                i,
                intersection[s, r, dist]
            ]
        ];


    intersect[p : plane[norm_, offset_, _], r : ray[start_, dir_], i : intersection[_, _, currentdist_]] :=
        Module[ {denom =  norm . dir, candidatedist},
            If[ denom >= 0.,
                i,
                candidatedist = (norm . start + offset) / (-denom);
                If[ candidatedist > currentdist,
                    i,
                    intersection[p, r, candidatedist]
                ]
            ]
        ];


    testray[ray_, scene_] :=
        dist[Fold[intersect[#2, ray, #1]&, miss, things[scene]]];


    traceray[ray_, scene_, depth_, maxdepth_] :=
        shade[Fold[intersect[#2, ray, #1]&, miss, things[scene]], scene, depth, maxdepth];


    shade[miss, _, _, _] :=
        background;     
    shade[intersection[thing_, ray[start_, dir_], dist_], scene_, depth_, maxdepth_] :=
        Module[ {pos = dist * dir + start, n, reflectdir, naturalcolor, reflectedcolor},
            n = normal[thing, pos];
            reflectdir = dir - 2. * n . dir * n;
            naturalcolor = defaultcolor + getnaturalcolor[thing, pos, n, reflectdir, scene];
            reflectedcolor = If[ depth >= maxdepth,
                                 grey,
                                 getreflectioncolor[thing, pos + (0.001*reflectdir), n, 
                                     reflectdir, scene, depth, maxdepth]
                             ];
            naturalcolor + reflectedcolor
        ];

    getreflectioncolor[thing_, pos_, n_, rd_, scene_, depth_, maxdepth_] :=
        reflect[surface[thing]][pos] *
        traceray[ray[pos, rd], scene, depth + 1, maxdepth];

    getnaturalcolor[thing_, pos_, n_, rd_, scene_] :=
        Module[ {addlight, normraydir = Normalize[rd], howrough = roughness[surface[thing]]},
            SetAttributes[addlight, Listable];
            addlight[light[p_, c_]] :=
                Module[ {ldis = p - pos, livec, neatisect, isinshadow, illum, lcolor, spec, scolor},
                    livec = Normalize[ldis];
                    neatisect = testray[ray[pos, livec], scene];
                    isinshadow = neatisect <= Norm[ldis];
                    If[ isinshadow,
                        defaultcolor,
                        illum = livec . n;
                        lcolor = If[ illum > 0.,
                                     illum * c,
                                     defaultcolor
                                 ];
                        spec = livec . normraydir;
                        scolor = If[ spec > 0.,
                                     (spec ^ howrough) * c,
                                     defaultcolor
                                 ];
                        diffuse[surface[thing]][pos] * lcolor + specular[surface[thing]][pos] * scolor
                    ]
                ];
            defaultcolor + Total[addlight[lights[scene]]]
        ];



     raytrace[screenwidth_ : 64, screenheight_ : 64, scene_ : basescene, maxdepth_ : 1] :=
         Module[ {getpoint},
             getpoint[x_, y_, camera_] :=
                 With[ {recenterx = (x - screenwidth / 2.) / (2. * screenwidth), 
                 recentery = -(y - screenheight / 2.) / (2. * screenheight)},
                     Normalize[forward[camera] + recenterx * right[camera] + recentery * up[camera]]
                 ];
             Image[ParallelArray[traceray[ray[pos[camera[scene]], getpoint[#2-1, #1-1, camera[scene]]], 
                 scene, 0, maxdepth]&, {screenheight, screenwidth}]] // AbsoluteTiming
         ];


     (* Harness *)

    (* surface diffuse, specular reflect, roughness *)
    uniformsurface[diffuse_, specular_, reflect_, roughness_] = surface[diffuse&, specular&, reflect&, roughness];

    shiny = uniformsurface[white, grey, .7, 250.];
    matteshiny = uniformsurface[white, darkgrey, .7, 250.];

    checkerboard =
        surface[
        If[ OddQ[Floor[#[[3]]] + Floor[#[[1]]]],
            white,
            black
        ]&, 
        white &, 
        If[ OddQ[Floor[#[[3]]] + Floor[#[[1]]]],
            .1,
            .7
        ]&, 
        150.];



     basescene = scene[{  
         sphere[{0., 1., -.25}, 1., shiny], 
              sphere[{-.5, 1.3, 1.5}, 0.5, matteshiny],
         plane[{0., 1., 0.}, 0., checkerboard]},
          {light[{-2., 2.5, 0.}, {.5,.45,.41}], light[{2.,4.5,2.},{.99,.95,.8}]},
          camera[{2.75, 2.0, 3.75}, {-.6, .5, 0.}]];

---

To illustrate where I am going with this, here is an example of how one could generate a mesh for a more complex object (a polysphere):

    polyspherepoints[rad_Real, divs_Integer] :=
        With[ {u = -(Pi/2.), v = -Pi, 
          du = Pi/divs, dv = (2.*Pi)/divs},
            rad * 
             Flatten[Table[{Cos[du*i + u]*Cos[dv*j + v], Sin[du*i + u], 
                Cos[du*i + u]*Sin[dv*j + v]}, {j, 0, divs}, {i, 0, 
            divs}], 1]
        ];

    (* Put the polygon vertices in the right order *)
    ordervertices[{{a_, b_}, {c_, d_}}] :=
        {a, b, d, c};

    orderverticestotriangeles[{{a_, b_}, {c_, d_}}] :=
    {{a, b, d}, {a, c, d}}

    (* Generate a list of (polyspherepoint) vertice numbers, 
       partition them cyclically, and then into quads, and associate them
       with Polygons *)

    polyspheremeshtriangles[rad_Real, divs_Integer] :=
        Normal @ GraphicsComplex[polyspherepoints[rad, divs], 
          Map[Polygon, 
           Map[orderverticestotriangeles,
             Partition[Partition[Range[(divs+1)^2], divs+1], {2, 2}, 1, 1], {2}], 1]]; 

    polyspheremeshtriangles[rad_Real, divs_Integer] :=
        Normal @ GraphicsComplex[polyspherepoints[rad, divs], 
          Map[Polygon, 
           Map[orderverticestotriangeles,
             Partition[Partition[Range[(divs+1)^2], divs+1], {2, 2}, 1, 1], {2}], 1]];

(It would have been satisfying to use some of the geometric transform functions built into Mathematica to generate the vertices, but life was too short.)

And here is what Graphics3D @ polyspheremeshtriangles[1., 8] generates:

Answer 1

It's dangerous if you don't know what the argument is, but replacing sin, cos, sqrt with their Taylor series approximations when the argument is small may help (or writing your own Compile Taylor series versions). Another thing which comes to mind but may not be helpful in your case is Fast Inverse Square Root