-------------------------------------------------
-- Haskell Raytracer
-- 26 Dec 2000
-- adb@tardis.ed.ac.uk
-- A raytracer written in Haskell for scenes consisting of spheres and
-- planes with procedurally defined 3d textures, and a lighting model
-- consisting of diffuse illumination, shadows and reflections (to a
-- fixed depth). Lights are colored and can be point sources (which
-- cast shadows) or non-occludable directional sources (which don't).
-- Camera and viewing screen can be placed anywhere in the scene, and
-- parallel and perspective projections are provided.
-- To run this, install hugs from www.haskell.org and run it.
-- Type :load Raytrace.hs to load the code and type 'test' to
-- do the raytracing and produce the file c:/test.ppm.
-- Ray tracing details: This is a naive raytracer with no spatial
-- partitioning of the scene (but then again, there's only spheres and
-- planes), no adaptive depth control (fixed depth of 2 currently) or
-- first-hit speedups. Also, there's only one ray fired per screen
-- pixel (mmm, look at that aliasing).
-- Haskell details: This is pretty basic Haskell, and probably
-- contains lots of naive coding. It doesn't use qualified types,
-- monads or higher order types. If, like me, you have the SOE book
-- then everything in this program is explained by the end of chapter
-- 9. I haven't read any further yet! Oh, and apologies for the wide
-- code formatting .. emacs fullscreen at 1024x768.
-- Known bugs:
-- * Having lights inside spheres makes for (pretty looking) weirdness
-- I wrote this in a day to learn Haskell. I was pretty impressed by
-- the speed of development. Compared to C/C++/Java, you spend much
-- more time thinking about the problem rather than the implementation
-- language. You don't need to worry about memory management. Your
-- code isn't littered with detail about *how* to do things (ie. in
-- Haskell you'd use map in places where you'd have a for loop in C,
-- and you don't have to explicitly state the loop variable or the
-- number of iterations). The strict typechecking and the interpreted
-- nature of hugs means that most errors are flagged quickly. It
-- doesn't stop you making logic errors though, but it's quite easy to
-- test functions individual to work out what's going wrong (big win
-- for interpreters). It's pretty concise too (around 160 lines of
-- actual code). And the performance isn't all that bad.
-------------------------------------------------
-- 3d vector representation and manipulation
-- NB: We have a right handed coordinate system. If x increases to your right, and Y increases downwards then
-- you are looking in the direction of increasing Z.
type Vector3 = (Float, Float, Float)
add :: Vector3 -> Vector3 -> Vector3
add (x,y,z) (a,b,c) = (a+x, b+y, c+z)
sub :: Vector3 -> Vector3 -> Vector3
sub (a,b,c) (x,y,z) = (a-x, b-y, c-z)
squared_mag :: Vector3 -> Float
squared_mag (x,y,z) = (x*x + y*y + z*z)
mag :: Vector3 -> Float
mag v = sqrt (squared_mag v)
scalarmult :: Vector3 -> Float -> Vector3
scalarmult (x,y,z) c = (x*c, y*c, z*c)
dot :: Vector3 -> Vector3 -> Float
dot (x,y,z) (a,b,c) = x*a + b*y + c*z
cross :: Vector3 -> Vector3 -> Vector3
cross (a,b,c) (x,y,z) = (b*z + c*y, -(a*z + c*x), a*y + b*x)
normalize :: Vector3 -> Vector3
normalize v
| (mag v) /= 0 = scalarmult v (1 / mag v)
| otherwise = (0,0,0)
neg :: Vector3 -> Vector3
neg (x,y,z) = (-x,-y,-z)
-------------------------------------------------
-- 3d Ray: Used as a semi-infinite 3d line.
type Point3 = Vector3
type Direction3 = Vector3
type Time = Float
type Ray = (Point3, Direction3) -- base and direction
position_at_time :: Ray -> Time -> Point3
position_at_time (base, dir) t = base `add` (scalarmult dir t)
-------------------------------------------------
-- Generic maths functions
-- Calculate the roots of the equation a * x^2 + b * x + c = 0
roots :: Float -> Float -> Float -> [Float]
roots a b c = let discriminant = b*b - 4*a*c
in if (discriminant < 0.0) then []
else [ 0.5 * (-b + sqrt discriminant), 0.5 * (-b - sqrt discriminant) ]
-- Arg, I can't find haskells xor function!
xor :: Bool -> Bool -> Bool
xor True b = not b
xor False b = b
---------------------------------------------------------------------------------------
-- Color model: We store red, green and blue values as Floats between 0.0 and 1.0
-- However, the bounds are not enforced in any way, except by calling clamp
type Color = (Float, Float, Float)
-- Predefined colors for convience
red, green, blue, white, black, midgrey, nearlywhite :: Color
red = (1.0, 0.0, 0.0)
green = (0.0, 1.0, 0.0)
blue = (0.0, 0.0, 1.0)
white = (1.0, 1.0, 1.0)
black = (0.0, 0.0, 0.0)
midgrey = (0.5, 0.5, 0.5)
nearlywhite = (0.8,0.8,0.8)
scale_col :: Color -> Float -> Color
scale_col (r,g,b) k = (r*k, g*k, b*k)
add_col :: Color -> Color -> Color
add_col (r1,g1,b1) (r2,g2,b2) = (r1+r2, g1+g2, b1+b2)
clamp :: Color -> Color
clamp (r,g,b) = ( clampfloat r, clampfloat g, clampfloat b)
where clampfloat f = (max 0.0 (min 1.0 f))
combine_col :: Color -> Color -> Color
combine_col (r1,g1,b1) (r2,g2,b2) = (r1*r2, g1*g2, b1*b2)
---------------------------------------------------------------------------------------
-- Procedural textures: Various predefined 3d texture functions
flatred, shinyred, semishinygreen, shinywhite :: Point3 -> Material
flatred _ = (red, 0.0, 1.0)
shinyred _ = (red, 0.3, 0.9)
semishinygreen _ = (green, 0.5, 0.7)
shinywhite _ = (white, 0.3, 0.9)
-- alternate 20x20x20 black and white cubes
checked_matt :: Point3 -> Material
checked_matt (x,y,z) = let xeven = even (truncate (x / 20.0))
yeven = even (truncate (y / 20.0))
zeven = even (truncate (z / 20.0))
in if (xeven `xor` yeven `xor` zeven) then (white, 0.0, 1.0) else (black, 0.0, 1.0)
---------------------------------------------------------------------------------------
-- Materials: Each point on the surface of a shape has its own color, reflectivity (kr)
-- and diffuseness (kd)
type Reflectivity = Float
type Diffuseness = Float
type Material = (Color, Reflectivity, Diffuseness)
---------------------------------------------------------------------------------------
-- Shapes: A shape is something which we check for intersection with a ray, and get
-- information about the intersection point (material, normal and intersection position)
-- Currently we have a plane and a sphere
type Normal = Vector3
type Radius = Float
data Shape = Sphere Point3 Radius (Point3 -> Material)
| Plane Normal Float (Point3 -> Material)
-- Plane is defined by a normal (its a 2 sided plane though) and a distance.
-- The plane coincident with y=5 and normal (0,0,1) has distance -5.
---------------------------------------------------------------------------------------
-- Intersection calculations
-- An intersection is represented by the normal at the intersection point, the point of intersection,
-- the direction of the viewing ray coming in to the intersection and the material at the intersection
-- point.
-- These intersection equations are taken from www.education.siggraph.org/materials/HyperGraph
type Intersection = (Normal, Point3, Ray, Material)
-- When we calculate reflected rays, they start on the surface of the shape. Unfortunately,
-- our limited numerical precision can make them be under the surface and so the reflected
-- ray immediately (at t very close to 0) hits the surface of the shape which it's meant to
-- be going away from. So, we filter out any intersections which occur at t < epsilon.
epsilon :: Float
epsilon = 0.001
intersect :: Ray -> Shape -> [(Time, Intersection)]
intersect ray@(base, dir) (Sphere center rad materialfn) =
let a = squared_mag dir
b = 2 * ( dir `dot` (base `sub` center))
c = (squared_mag (base `sub` center)) - rad^2
times = filter (> epsilon) (roots a b c)
normal_at_time t = normalize ((position_at_time ray t) `sub` center)
intersection_at_time t = (normal_at_time t, position_at_time ray t, ray, materialfn (position_at_time ray t))
in map (\t -> (t,intersection_at_time t)) times
intersect ray@(base, dir) (Plane normal d materialfn) =
let vd = (normalize normal) `dot` dir
v0 = negate (((normalize normal) `dot` base) + d)
in if (vd == 0) then []
else let t = v0 / vd
hitpoint = position_at_time ray t
in if t > epsilon then [ (t, (if (vd > 0) then (neg normal) else normal, hitpoint, ray, materialfn hitpoint)) ]
else []
-- Extract the closest intersection (lowest time) from a list
closest :: [ (Time,Intersection) ] -> Intersection
closest xs = let select_nearest (t1,i1) (t2,i2) = if (t1 Point3 -> Bool
point_is_lit point lightpos = let path = lightpos `sub` point
time_at_light = mag path
ray = (point, normalize path)
hits = concat (map (intersect ray) shapes)
times = fst (unzip hits)
in if (null times) then True else (minimum times) > time_at_light
-- Helper to calculate the diffuse light at the surface normal, given
-- the light direction (from light source to surface)
diffuse_coeff :: Vector3 -> Vector3 -> Float
diffuse_coeff light_dir normal = max 0.0 (negate ((normalize light_dir) `dot` (normalize normal)))
local_light :: Intersection -> Light -> Color
-- Simple case of a non-shadowable directional light
local_light (normal,_,_,(materialcol,_,kd)) (Directional dir lightcol) =
let mixed_color = combine_col materialcol lightcol
diffuse = scale_col mixed_color ((diffuse_coeff dir normal) * kd)
in diffuse
-- Spotlight - shadowable
local_light (normal, hitpoint,_,(materialcol,_,kd)) (Spotlight lightpos lightcol) =
let mixed_color = combine_col materialcol lightcol
diffuse = scale_col mixed_color (kd * (diffuse_coeff (hitpoint `sub` lightpos) normal))
in if (point_is_lit hitpoint lightpos) then diffuse else black
---------------------------------------------------------------------------------------
-- Reflections (part of the global lighting model)
-- Ray trace the outgoing reflected ray from an intersection (depth is the level of recursion
-- which we're at in the ray tracing)
reflected_ray :: Integer -> Intersection -> Color
reflected_ray depth (normal, hitpoint,(_,in_ray_dir),(color,kr,_))
| kr == 0.0 = black
| otherwise = let k = 2 * ((normalize normal) `dot` (normalize (neg in_ray_dir)))
out_ray_dir = (scalarmult (normalize normal) k) `sub` (neg in_ray_dir)
reflected_col = raytrace (depth + 1) (hitpoint, out_ray_dir)
in scalarmult reflected_col kr
---------------------------------------------------------------------------------------
-- Image output: We can write a ppm (Portable Pixmap) file by converting a list of
-- colors (length is width * height) into a big string
make_pgm :: Integer -> Integer -> [ Color ] -> String
make_pgm width height xs = "P3\n" ++ show width ++ " " ++ show height ++ "\n255\n" ++ stringify(xs)
where stringify [] = ""
stringify ((r,g,b):xs) = show (round (r*255)) ++ " "
++ show (round (g*255)) ++ " "
++ show (round (b*255)) ++ " "
++ stringify xs
---------------------------------------------------------------------------------------
-- Viewing screen and camera functions: We define the camera position and the point which
-- we're looking at. We also define an 'up' vector and a distance to the screen. The
-- following functions generate a list of points (in raster order - a fact which is relied
-- upon when we write the image out as a ppm) which correspond to the 3d positions of the
-- pixels on our virtual screen.
-- Camera position, distance to screen, "Looking at" position, up vector
type View = (Point3, Float, Point3, Vector3)
pixel_grid :: View -> Float -> Float -> [ Point3 ]
pixel_grid (camerapos, viewdist, lookingat, viewup) width height =
let grid = [ (x, y, 0) | y <- [0..width-1], x <- [0..height-1] ]
centering_offset = (- width / 2.0, -height / 2.0 , 0)
pixel_offsets = map (add centering_offset) grid
viewdir = normalize (lookingat `sub` camerapos)
screen_center = camerapos `add` (scalarmult viewdir viewdist)
viewright = viewdir `cross` viewup
transform (x,y,_) = screen_center `add` (scalarmult viewright x) `add` (scalarmult (neg viewup) y)
in map transform pixel_offsets
-- Parallel projection function which creates rays parallel to the viewing screen
parallel_projection :: View -> Point3 -> Ray
parallel_projection (camerapos,_,lookingat,_) point = (point, normalize (lookingat `sub` camerapos))
-- Perspective projection which creates rays through (0,0,-distance) through the point
perspective_projection :: View -> Point3 -> Ray
perspective_projection (camerapos,_,_,_) point = (point, normalize (point `sub` camerapos))
---------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------
-- Main rendering functions
-- Calculate the overall color of a ray/shape intersection, taking into account
-- local lighting (diffuse only) and global lighting (reflections only, to a depth
-- of 2 bounces)
overall_lighting :: Integer -> Intersection -> Color
overall_lighting depth hit = let sum_colors = foldr add_col black
local_lighting = ambient_light `add_col` sum_colors (map (local_light hit) lights)
global_lighting = if (depth < 2) then (reflected_ray depth hit) else black
in clamp (local_lighting `add_col` global_lighting)
-- Trace a ray through the scene and work out what color it should be.
-- Takes a 'depth' argument which is 0 for top level viewing rays increases
-- by one for each level of recursive raytracing we do (as a result of reflections
-- or transmissions)
raytrace :: Integer -> Ray -> Color -- uses global 'shapes'
raytrace depth ray = let hits = concat (map (intersect ray) shapes)
in if (null hits) then background_color
else overall_lighting depth (closest hits)
---------------------------------------------------------------------------------------
-- Testing: We define a function 'test' which renders a fixed view and writes
-- out the result as c:/test.ppm
-- This writes out a pgm of our trivial rendering, given the screen width and height
render_to_pgm :: Float -> Float -> String
render_to_pgm width height = let view = ( (0,0,-100), 100, (0,0,100), (0,-1,0))
projection = perspective_projection view
ray_collection = map projection (pixel_grid view width height)
color_collection = map (raytrace 0) ray_collection
in make_pgm (round width) (round height) color_collection
main = do writeFile "test.ppm" (render_to_pgm 500 500)