Tuesday, 25 August 2009

Some Ray Tracing Functions in Haskell

I'm still trying to find my way in Haskell, and I've always thought the best way is to implement some real code that does something. I thought I'd port across the little Clojure ray tracer example that I ported from ANSI Common Lisp.

The first version I wrote handily ignored type definitions, used tuples everywhere. I didn't write out the type signatures of my functions with the net result that they often didn't do what I thought. Lesson learnt. Haskell is a strongly typed language so to work with it, I should use types! I used the following types which are hopefully self explanatory.

data Point = Point { x :: Float
, y :: Float
, z :: Float
} deriving (Show)

data Sphere = Sphere { color :: Float
, radius :: Float
, centre :: Point
} deriving (Show)

data ObjectHit = ObjectHit { object :: Sphere
, location :: Point
} deriving (Show)

data Brightness = Brightness { value :: Float } deriving (Show)

deriving Show is is used to say that this type definition is printable and should be printed out in the default way. The names of the types inside can be used as accessor functions (e.g. x p gives you the x co-ordinate of Point p).

Firstly we need some basics functions to compute various number and vector properties.

square :: (Num a) => a -> a
square x = x * x

magnitude :: Point -> Float
magnitude p = sqrt ((square (x p)) + (square (y p)) + (square (z p)))

unitVector :: Point -> Point
unitVector p = let d = magnitude p
in Point ((x p)/d) ((y p)/d) ((z p)/d)

pointSubtract :: Point -> Point -> Point
pointSubtract p1 p2 = Point (x p1-x p2) (y p1-y p2) (z p1-z p2)

distance :: Point -> Point -> Float
distance p1 p2 = magnitude (pointSubtract p1 p2)

sphereNormal :: Sphere -> Point -> Point
sphereNormal s p = unitVector (pointSubtract (centre s) p)

lambert :: Sphere -> Point -> Point -> Float
lambert s i r = let n = sphereNormal s i
in max 0 ((x r * x n) + (y r * y n) + (z r * z n))

The definitions of the next set of functions is slightly more interesting.

minroot :: Float -> Float -> Float -> Maybe Float
minroot a b c
| a == 0 = Just ((- c) / b)
| otherwise = let disc = (square b) - (4 * a * c)
in if (disc > 0)
then Just (min (((-b) + sqrt disc) / (2 * a)) (((-b) - sqrt disc) / (2 * a)))
else Nothing

sphereIntersect :: Sphere -> Point -> Point -> Maybe ObjectHit
sphereIntersect s pt r = let c = centre s
n = minroot (square (x r) + square (y r) + square (z r))
(2 * ((x r * (x pt - x c)) + (y r * (y pt - y c)) + (z r * (z pt - z c))))
((square (x pt - x c)) + (square (y pt - y c)) + (square (z pt - z c)) - (square (radius s)))
in if (isNothing n)
then Nothing
else Just (ObjectHit
((x pt) + (fromJust n) * (x r))
((y pt) + (fromJust n) * (y r))
((z pt) + (fromJust n) * (z r))))

spheresHit :: [Sphere] -> Point -> Point -> [ObjectHit]
spheresHit sw pt r = mapMaybe (\x -> sphereIntersect x pt r) sw

nearestHit :: [Sphere] -> Point -> Point -> Maybe ObjectHit
nearestHit sp pt r = let hitSpheres = spheresHit sp pt r
case hitSpheres of
[] -> Nothing
x -> Just (head (sortBy
(\h1 h2 -> (compare (distance (location h1) pt) (distance (location h2) pt)))

Maybe is a type that might be null (Nothing in Haskell). The type system enforces that you handle both cases. This is useful for solving the quadratic equation (minroot) because we can indicate that an equation has no solution without having to resort to either exceptions or picking a sentinel value and hoping it never occurs. Similarly, spheresIntersect returns the intersection of a ray and a sphere, but that intersection might never occur, hence we can return Nothing. mapMaybe only performs the map if the element exists and throws out the results of anything with Nothing as a value.

Apparently Maybe is a monad too. I'm deliberately trying to avoid understanding the deep meaning of monad and trying to use it instead. There's already too many articles about monads!

Finally we can perform the actual ray tracing with:

sendRay :: [Sphere] -> Point -> Point -> Brightness
sendRay world src ray = let hit = nearestHit world src ray
in if (isNothing hit)
then (Brightness 0)
else let sp = object (fromJust hit) in
(Brightness ((color sp) * (lambert sp src ray)))

colorAt :: [Sphere] -> Point -> Float -> Float -> Brightness
colorAt world eye x y = let ray = unitVector (pointSubtract (Point x y 0) eye)
in (Brightness (255 * value (sendRay world eye ray)))

Now we get to the slight downside of Haskell (at least for a beginner like me). I have my code, and I'm confident it is correct - how do I write this image out as a PNG? Hackage lists lots of graphics packages, but they are all quite heavy going to a newbie to understand. At the moment I'm trying to write out a PPM file (the approach taken here), but it feels painful. Clojure doesn't suffer from this problem because Java has a library for everything (and I know Java!). I'm sure this is a temporary roadblock though, so time for me to try and find a suitable package on Hackage and understand enough bits and pieces to be able to use it!