[Shootout-list] Ray tracer

[Jon Harrop]

--Boundary-00=_2mMcCfbj7Z9UnU8
Content-Type: text/plain;
  charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

On Thursday 28 April 2005 11:38, Sebastien Loisel wrote:
> > 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.
>
> This program is very good. To please everyone perhaps we can try to
> simplify it a smidgen further. If use only translations to create a
> sphere flakes (as opposed to translations+rotations), do you suppose
> that would still be an interesting test? And then, would we need the
> matrix infrastructure?

Thanks for the tip. Now its only 62 LOC:

let delta = sqrt epsilon_float and pi = 4. *. atan 1.
let ( *| ) s r = [| s *. r.(0); s *. r.(1); s *. r.(2)|]
let ( +| ) a b = [| a.(0) +. b.(0); a.(1) +. b.(1); a.(2) +. b.(2)|]
let ( -| ) a b = [| a.(0) -. b.(0); a.(1) -. b.(1); a.(2) -. b.(2)|]
let dot a b = a.(0) *. b.(0) +. a.(1) *. b.(1) +. a.(2) *. b.(2)
let unitise r = (1. /. sqrt (dot r r)) *| r
type sphere = { center: float array; radius: float }
type ray = { orig: float array; dir: float array }
type obj = Sphere of sphere * float | Group of sphere * obj list
let ray_sphere ray sphere =
  let v = sphere.center -| ray.orig in
  let b = dot v ray.dir in
  let disc = b *. b -. dot v 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, [|0.; 0.; 0.|], 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 = [|0.; 0.; 0.|] -| 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 [|-1.; -3.; 2.|] in
  let scene =
    let rec aux level r (x, y, z) =
      let sphere = { center = [|x;y;z|]; radius = r } in
      let obj = Sphere (sphere, 1.) in
      if level = 1 then obj else begin
 let aux l (x', y', z') = aux (level-1) (0.5 *. r)
   (x -. x', y +. y', z +. z') :: l in
 let objs = let r' = r /. sqrt 3. in List.fold_left aux [obj]
   [-.r', r', -.r'; r', r', -.r'; -.r', r', r'; r', r', r'] in
 Group ({sphere with radius = 2. *. r}, objs)
      end in
    aux level 1. (0., 0., 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 -. float h in
 {orig = [|0.;2.;-5.|]; dir = unitise [|x; y; float (max w h)|] } in
      Printf.printf "%d " (int_of_float (256. *.ray_trace 1. light ray 
scene));
    done;
    Printf.printf "\n";
  done;

Unless I've failed to use OCaml's built-in ray_trace function, I don't think 
it's going to get a lot smaller... ;-)

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

--Boundary-00=_2mMcCfbj7Z9UnU8
Content-Type: text/x-csrc;
  charset="iso-8859-1";
  name="ray8.ml"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
	filename="ray8.ml"

let delta = sqrt epsilon_float and pi = 4. *. atan 1.
let ( *| ) s r = [| s *. r.(0); s *. r.(1); s *. r.(2)|]
let ( +| ) a b = [| a.(0) +. b.(0); a.(1) +. b.(1); a.(2) +. b.(2)|]
let ( -| ) a b = [| a.(0) -. b.(0); a.(1) -. b.(1); a.(2) -. b.(2)|]
let dot a b = a.(0) *. b.(0) +. a.(1) *. b.(1) +. a.(2) *. b.(2)
let unitise r = (1. /. sqrt (dot r r)) *| r
type sphere = { center: float array; radius: float }
type ray = { orig: float array; dir: float array }
type obj = Sphere of sphere * float | Group of sphere * obj list
let ray_sphere ray sphere =
  let v = sphere.center -| ray.orig in
  let b = dot v ray.dir in
  let disc = b *. b -. dot v 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, [|0.; 0.; 0.|], 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 = [|0.; 0.; 0.|] -| 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 [|-1.; -3.; 2.|] in
  let scene =
    let rec aux level r (x, y, z) =
      let sphere = { center = [|x;y;z|]; radius = r } in
      let obj = Sphere (sphere, 1.) in
      if level = 1 then obj else begin
	let aux l (x', y', z') = aux (level-1) (0.5 *. r)
	  (x -. x', y +. y', z +. z') :: l in
	let objs = let r' = r /. sqrt 3. in List.fold_left aux [obj]
	  [-.r', r', -.r'; r', r', -.r'; -.r', r', r'; r', r', r'] in
	Group ({sphere with radius = 2. *. r}, objs)
      end in
    aux level 1. (0., 0., 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 -. float h in
	{orig = [|0.;2.;-5.|]; dir = unitise [|x; y; float (max w h)|] } in
      Printf.printf "%d " (int_of_float (256. *.ray_trace 1. light ray scene));
    done;
    Printf.printf "\n";
  done;

--Boundary-00=_2mMcCfbj7Z9UnU8--