Signed Distance Functions are a pretty simple concept. Basically for each point in the world, you return a distance to the nearest surface, negative distances are inside geometry. One neat application of them to represent scene geometry, and ray trace into the SDF. Why would anyone do such a thing? Well, it turns out that when marching along a ray looking for intersections (which is obviously not the only way to trace rays), it’s jolly useful to know a minimum bound on when you might expect to encounter a surface. This way you can keep taking large jumps along the ray, using the SDF to figure out how large, until you finally reach some threshold distance. Here’s the end result of the code I’m presenting in this blog post.

Note the cool twisted pillar, and the nice soft ambient occlusion around contact points. Both of these features are simple to do when tracing SDFs, but much harder with other methods.

So basically we represent our scene as a SDF, but in order to have varying materials we need this function to not only return the distance to the closest surface, but the material of that surface too, so we end up with the following Haskell type for our scenes:

type Scene = Vec3 -> (Double, Material)

In order to define a simple scene we just need to write a function which returns the distance to the surfaces, for example here’s a unit sphere:

sphere pt = (mag pt - 1.0, defaultMaterial )

In other words, we take the magnitude of the position (i.e. distance from the origin) and subtract 1.0 since this is our radius. Simple! To build larger scenes from simple scenes such as this, we need a way to combine two scenes. This is done by just taking the minimum of the distances, and picking the corresponding material.

mergeScenes :: Scene -> Scene -> Scene
mergeScenes scene1 scene2 pt 
    | d1 < d2 = res1
    | otherwise = res2
    where   res1@(d1,_) = scene1 pt
            res2@(d2,_) = scene2 pt

Since our scenes are just functions from position to distance, it’s very easy to manipulate the geometry by warping the inputs to the SDF – for example, in order to build the twisted column shown in the screenshot above we first build a regular column, and then twist it by rotating the input by varying amounts (related to the y component):

column pt@(Vec3 x y z) = (mag pt - mag closestPt, defaultMaterial )
    where   clamp x = max (-0.5) (min 0.5 x)
            closestPt = Vec3 (clamp x) y (clamp z)              
twistedColumn pt@(Vec3 _ y _) = column (rotateY y pt)

There are many other possibilities for these kinds of modifications, for example we could add some random noise to the distance function to get bumpy surfaces etc.

Once we have a scene to ray trace against, we simply march through the SDF until we reach a distance small enough to consider “on” the surface:

type Ray = (Vec3, Vec3) -- origin, direction

rayTrace :: Scene -> Double -> Ray -> Maybe Color
rayTrace s end (pos,dir)
    | end <= 0 = Nothing
    | dist < 0.0001 = Just $ shade s material dir n pos defaultLight `scale` getAO s pos n
    | otherwise = rayTrace s (end-dist) (pos + dir `scale` dist, dir)
    where   (dist, material) = s pos                           
            n = getNormal dist s pos 

The first equation of rayTrace simply checks if we’ve reached the far plane, in which case we return Nothing, since there was no hit. The second guard checks if we’ve reached a surface, and if so shades the point using our light and the normal is computed by looking at the distances in the neighbourhood of the point. The whole thing is modulated by our ambient occlusion factor (which may not be strictly correct, I know).

This ambient occlusion factor is computed by looking at a small number of samples “above” our surface point and comparing the distance to these samples with the distances in the SDF for those samples. If the sampled distances are smaller than the distance to the surface, then that means there’s some amount of “occlusion” for that sample (since there’s evidently some object other than the surface nearby). We weight these occlusion estimates and sum them to get the final ambient occlusion. This is pretty hacky with several magic numbers to tweak, but it works really well, giving smooth 3D ambient occlusion with very little cost.

Here’s the full code listing. I’ve implemented a small vector library with the stuff we need in order to keep it reasonably self-contained. You will need some way of actually displaying/saving the final image, though - the code below uses the ppm package (”cabal install ppm” should do it). This code is not intended to be an example of “nice” Haskell code – it’s pretty much just a copy of the code I ended up with after experimenting with SDFs, so there’s lots of hard coded values and other limitations (e.g. just one light source) which are, erm, ”left as an exercise to the reader”.

{-# LANGUAGE  ParallelListComp #-}
import Codec.Image.PPM hiding ( Color )
import System.IO
import Data.List
import Control.Parallel.Strategies

type Material = (Vec3, Double, Double) -- Color, Specular power, gloss
type Scene = Vec3 -> (Double, Material)
type Light = (Vec3, Vec3, Double)           
type Color = Vec3
type Position = Vec3
type Direction = Vec3
type Ray = (Vec3, Vec3)

main = do
    let rays = getRays 256 (pi/2)
    let colors = parMap rnf (map (maybe (0,0,0) toRGB24 . rayTrace scene 100)) rays
    writePPM "output.ppm" colors

colorize :: Color -> Scene -> Scene
colorize c s pt = let (d, (_,p,g)) = s pt in (d, ( c,p,g ))

red, green :: Scene -> Scene
red = colorize (Vec3 1 0 0)
green = colorize (Vec3 0 1 0)

True `xor` a = not a
False `xor` a = a         

checker x y = checker1D x `xor` checker1D y
    where checker1D x = floor x `mod` 2 == 0

-- the scenes
sphere pt = (mag pt - 1.0, (Vec3 1 1 1, 20, 0.5) )
plane (Vec3 x y z) = (y + 2, if checker x z
                                then (Vec3 1 1 0.5, 10,0.4)
                                else (Vec3 0.3 0.3 0.1, 10,1) )
plane2 (Vec3 _ y _) = (1.5-y, (Vec3 1 1 1, 10, 1))

column pt@(Vec3 x y z) = (mag pt - mag closestPt, (Vec3 0.4 0.7 1, 2, 0.5) )
    where   clamp x = max (-0.5) (min 0.5 x)
            closestPt = Vec3 (clamp x) y (clamp z)              
twistedColumn pt@(Vec3 _ y _) = column (rotateY y pt)

mergeScenes :: Scene -> Scene -> Scene
mergeScenes scene1 scene2 pt 
    | d1 < d2 = res1
    | otherwise = res2
    where   res1@(d1,_) = scene1 pt
            res2@(d2,_) = scene2 pt
scene :: Scene
scene = finalScene . (+ Vec3 0 0 5)
    where finalScene = foldl1' mergeScenes [
            green sphere,
            plane2,
            plane,
            twistedColumn . (+ Vec3 2.2 0 0),
            red sphere . (+Vec3 (-1.8) 1 0) ]

getNormal :: Double -> Scene -> Position -> Direction
getNormal d s pt = normalize (Vec3 x y z)
    where   eps = 0.0001
            x = fst ( s (pt + Vec3 eps 0 0 )) - d
            y = fst ( s (pt + Vec3 0 eps 0 )) - d
            z = fst ( s (pt + Vec3 0 0 eps )) – d

shade :: Scene -> Material -> Direction -> Direction -> Position -> Light -> Color
shade s (materialColor, specPower, gloss) eye n pt (color, pos, att)
    = combinedColor `scale` (attenuation*lambert*shadow)
    where   attenuation = 1/((att * mag lightDir )^2)
            lightDir = pos-pt
            lightDirN = normalize lightDir
            lambert = max 0 (n `dot` lightDirN )
            combinedColor = color * materialColor + spec
            r = eye - n `scale` (n `dot` eye)*2
            spec = color `scale` ( gloss * (max 0 (r `dot` lightDirN))**specPower)
            shadow  | lambert == 0 = 0
                    | otherwise = case rayTrace s (mag lightDir) (pt + n `scale` 0.0001, lightDirN) of
                                    Nothing -> 1.0
                                    Just _ -> 0.4                   
defaultLight = (Vec3 15 15 14, Vec3 2 1 0, 0.5)

getAO :: Scene -> Position -> Direction -> Double
getAO s pt n = clamp $ 1.0-k*sum (take count [(x-y)*w | w <- iterate (*falloff) 1 | x <- dists | y <- sampledDists])
    where   dists = [0, delta ..] :: [Double]
            sampledDists = map (\d -> fst $ s ( pt + n `scale` d ) ) dists :: [Double]
            clamp x = max 0 (min 1 x)           
            k = 2.75
            count = 5
            delta = 0.1
            falloff = 0.57            

rayTrace :: Scene -> Double -> Ray -> Maybe Color
rayTrace s end (pos,dir)
    | end <= 0 = Nothing
    | dist < 0.0001 = Just $ shade s material dir n pos defaultLight `scale` getAO s pos n
    | otherwise = rayTrace s (end-dist) (pos + dir `scale` dist, dir)
    where   (dist, material) = s pos                           
            n = getNormal dist s pos 

toRGB24 (Vec3 r g b) = (exposure r, exposure g, exposure b)                           
    where exposure f = max 0 (min 255 (round ( 255 * (1-exp (-f)) )))

writePPM fname img = do
    let imgData = ppm_p6 img
    withBinaryFile fname WriteMode (\h -> hPutStr h imgData )

getRays sz horizFOV = map computeRow pixelCoords1D
    where   pixelCoords1D = take sz [-1, -1 + 2/fromIntegral sz .. ]
            z = -(1/atan (horizFOV / 2))
            computeRow y = [ (Vec3 0 0 0, normalize (Vec3 x (-y) z) ) | x <- pixelCoords1D ]  

-- quick and dirty vector maths
data Vec3 = Vec3 !Double !Double !Double deriving (Show,Eq)

instance Num Vec3 where
    (Vec3 x y z) + (Vec3 x' y' z') = Vec3 (x+x') (y+y') (z+z')
    negate (Vec3 x y z) = Vec3 (-x) (-y) (-z)
    x - y = x + negate y
    abs (Vec3 x y z) = Vec3 (abs x) (abs y) (abs z)
    (Vec3 x y z) * (Vec3 x' y' z') = Vec3 (x*x') (y*y') (z*z')
    fromInteger x = let x' = fromInteger x in Vec3 x' x' x'
    signum = error "signum on Vec3 is not implemented"

(Vec3 x y z) `dot` (Vec3 x' y' z') = x*x' + y*y' + z*z'
(Vec3 x y z) `scale` a = Vec3 (a*x) (a*y) (a*z)
mag v = sqrt (v `dot` v)
normalize v = v `scale` (1/mag v)
rotateY theta (Vec3 x y z) = Vec3 (x*costheta - z*sintheta) y (z*costheta + x*sintheta)
    where   costheta = cos theta
            sintheta = sin theta

Side note: the outer loop for our ray tracer uses “parMap” for almost perfect linear scaling (I get 1.9x on my dual core machine, and I suspect a lot of the overhead is just the serial business of actually writing the image to disk). It’s interesting to note that this was added pretty much as an afterthought; while playing with the code and adding more and more stuff the runtime got slower and slower, so as a quick’n’dirty speedup I changed a “map” to “parMap rnf” and it ran twice as fast. This wasn’t a carefully considered optimization pass, it was more like “This is slow, I’ll just hack in some parallelization real quick”. I dare say I’ll never utter that sentence when writing C – even if it had a “parMap” equivalent, I’d be very careful about using it since there’s no way of knowing that there won’t be any hidden side effects somewhere.