[Shootout-list] Ray tracer

[Jon Harrop]

On Wednesday 27 April 2005 22:55, Greg Buchholz wrote:
> > This is a nice program but if I read this right it is brute force. It
> > maps isect to objects; the point of the other program is to do
> > something less than map isect to objects. Unless the lazyness is doing
> > something I'm not getting?
>
>     Nope, brute force.  Its merely a proof-of-concept that a ray tracer
> can be done in <100 LOC.  We could always make it a "same thing" test and
> let the masses produce all the optimizations their little hearts desire.

99 LOC! That's incredible... for a Haskell programmer. ;-)

I've boiled my program down to 94 LOC but managed to keep the hierarchical
spherical bounding volumes, so it can still render a sphere-flake containing
600,000 spheres in under 8.5 seconds on my cheapy laptop.

In contrast, your Haskell program takes 6 minutes and renders significantly
fewer objects. However, my OCaml code doesn't support checkered planes and
reflections but it does do shadows.

Here's the code:

let delta = sqrt epsilon_float and pi = 4. *. atan 1.
let sin t = sin(pi *. t /. 180.) and cos t = cos(pi *. t /. 180.)
let vec3 x y z = [| x; y; z; 1. |] and zero = [|0.; 0.; 0.; 1.|]
let ( *| ) s r = [| s *. r.(0); s *. r.(1); s *. r.(2); 1. |]
let ( +| ) a b = [| a.(0) +. b.(0); a.(1) +. b.(1); a.(2) +. b.(2); 1. |]
let ( -| ) a b = [| a.(0) -. b.(0); a.(1) -. b.(1); a.(2) -. b.(2); 1. |]
let dot a b = a.(0) *. b.(0) +. a.(1) *. b.(1) +. a.(2) *. b.(2)
let length2 r = dot r r and length r = sqrt(dot r r)
let unitise r = (1. /. length r) *| r
let init f = Array.init 4 (fun i -> Array.init 4 (f i))
let trans r = init (fun i j -> if i=j then 1. else if j=3 then r.(i) else 0.)
let identity = trans zero
let sum4 f = f 0 +. f 1 +. f 2 +. f 3
let mat_mul a b = init (fun i j -> sum4 (fun k -> a.(i).(k) *. b.(k).(j)))
let transform m r = Array.init 4 (fun i -> dot m.(i) r +. m.(i).(3))
let rot_x t = let cos, sin = cos t, sin t in [|[|    1.;    0.;    0.; 0. |];
            [|    0.;   cos;   sin; 0. |];
            [|    0.; -.sin;   cos; 0. |];
            [|    0.;    0.;    0.; 1. |]|]
let rot_y t = let cos, sin = cos t, sin t in [|[|   cos;    0.;   sin; 0. |];
            [|    0.;    1.;    0.; 0. |];
            [| -.sin;    0.;   cos; 0. |];
            [|    0.;    0.;    0.; 1. |]|]
let rot_z t = let cos, sin = cos t, sin t in [|[|   cos;   sin;    0.; 0. |];
            [| -.sin;   cos;    0.; 0. |];
            [|    0.;    0.;    1.; 0. |];
            [|    0.;    0.;    0.; 1. |]|]
type sphere = { center: float array; radius: float }
type obj = Sphere of sphere * float | Group of sphere * obj list
type ray = { orig: float array; dir: float array }
let ray_sphere ray sphere =
  let v = sphere.center -| ray.orig in
  let b = dot v ray.dir in
  let disc = b *. b -. length2 v +. sphere.radius *. sphere.radius in
  if disc < 0. then infinity else
    let disc = sqrt disc in
    (fun t2 -> if t2 < 0. then infinity else
       ((fun t1 -> if t1 > 0. then t1 else t2) (b -. disc))) (b +. disc)
let intersect ray scene =
  let rec of_scene ((l, _, _) as first) = function
      Sphere (s, material) ->
 let l' = ray_sphere ray s in
 if l' >= l then first else (* No nearer *)
   let normal = unitise (ray.orig +| l' *| ray.dir -| s.center) in
   l', normal, material (* Replace with nearer intersection *)

    | Group (bound, scenes) ->

 let l' = ray_sphere ray bound in (* Cull if possible *)
 if l' >= l then first else List.fold_left of_scene first scenes in
  of_scene (infinity, zero, 0.) scene
let rec ray_trace weight light ray scene = match intersect ray scene with
    lambda, n, color ->
      if lambda = infinity then 0. else
 let o = ray.orig +| lambda *| ray.dir +| delta *| n in
 (match intersect { orig = o; dir = zero -| light } scene with
    l, _, _ when l = infinity -> max 0. (-. dot n light)

  | _ -> 0.) *. color

let () =
  let level = match Sys.argv with [| _; l |] -> int_of_string l | _ -> 6 in
  let light = unitise (vec3 (-1.) (-3.) 2.) and s = 1. /. 3. in
  let scene =
    let rec aux level r m =
      let sphere = { center = transform m zero; radius = r } in
      let obj = Sphere (sphere, 1.) in
      if level = 1 then obj else begin
 let m' = trans (vec3 0. ((1. +. s) *. r) 0.) in
 let rec aux2 f z i objs =
   if i=0 then objs else
     let m = List.fold_left mat_mul m [rot_y (float (f i)); z; m'] in
     aux2 f z (i-1) (aux (level-1) (s*.r) m :: objs) in
 let objs = aux2 (fun i -> 60*i) (rot_z 115.) 6
   (aux2 (fun i -> 30+120*i) (rot_z 45.) 3 []) in
 let aux2 r = function Sphere (sph', _) | Group (sph', _) ->
   max r (length (sph'.center -| sphere.center) +. sph'.radius) in
 let r' = (* Bounding radius *)
   List.fold_left
     (fun r -> function Sphere (sph', _) | Group (sph', _) ->
        max r (length (sph'.center -| sphere.center) +. sph'.radius))
     r objs in
 Group ({sphere with radius = List.fold_left aux2 r' objs}, obj::objs)
      end in
    aux level 1. (trans (vec3 0. 1. 0.)) in
  let w, h = 768, 768 in
  Printf.printf "P2\n%d %d\n256\n" w h;
  for y = h - 1 downto 0 do
    for x = 0 to w - 1 do
      let ray =
 let x, y = float x -. 0.5 *. float w, float y -. 0.5 *. float h in
 let d = unitise (vec3 x y (float (max w h))) in
 {orig = vec3 0. 2. (-4.); dir = transform (rot_x (-14.)) d} in
      Printf.printf "%d " (int_of_float (256. *.ray_trace 1. light ray
scene));
    done;
    Printf.printf "\n";
  done;

compile it with:

ocamlopt ray.ml -o ray

then run it with:

./ray >image.pgm

-- 
Dr Jon D Harrop, Flying Frog Consultancy Ltd.
Objective CAML for Scientists
http://www.ffconsultancy.com/products/ocaml_for_scientists