Haskell Raytracer

Lately I've been working on writing a ray tracer in Haskell.

Mathematically it's fairly complicated. In general what you do is substitute the formula for a ray into the formula for an object, then find the roots of the resulting equation. So for example if you have a sphere equation:

where Sr is the sphere radius,
and the ray equation:


where P is the starting point, and D is the direction.
Then you substitute the ray equation for the XYZ portions:

and then simplify:

Which is in quadratic form, so there's possibly two solutions. Note that if is normalized then so it can be omitted from the calculation. Once you have the roots T then you just pick the smallest that is in front of the ray starting point by at least epsilon, and convert that to the XYZ point.

sphereIntersect :: Vector -> Double -> Ray -> Maybe Hit
sphereIntersect loc radius ray@(pnt,dir) =
  let pl = pnt - loc
      b = 2*((x dir)*(x pl) + (y dir)*(y pl) + (z dir)*(z pl))
      c = (x pl)*(x pl) + (y pl)*(y pl) + (z pl)*(z pl) - radius*radius
      q = filter (>=epsilon) $ quadRootsIgnoreA b c
  in if (q == [])
      then Nothing
      else let ti = minimum q
               iPoint = evaluateRay ray ti
               iNorm = (iPoint-loc)/(toVector radius)
           in Just (iPoint, iNorm)

Of course you can't end with just spheres, there's planes, boxes, triangles, flat discs, cones, etc.

Transforms are complex too. Say for example that you'd like to declare a right cone frustum as:

cone {basePnt, baseRad, topPnt, topRad}

then you need to rotate and scale the ray along with the cone so that you can test against a standardized cone where you know the formula for the surface:
So you end up with pairs of 4x4 matrices describing the transform and it's inverse. Interestingly, once you get a list of transforms to be applied (translate, rotate, scale) they can be composed together into a single transform by matrix multiplication.

The intention is to get this heavy matrix calculation to run once per object, and then use the result transform in testing multiple ray intersects. In Object Oriented languages like C++ this would be done at construction of the object, in Haskell I'd like to partially apply an intersect function, and then use the resulting function on a ray... so: objectFn :: Ray -> Maybe Hit
All that is to be tested and implemented still.

Currently I'm working on torus, a torus-ray intersect results in a quartic equation:

So that requires a good quartic solver...