[Shootout-list] Ray tracer
Sebastien Loisel
Sebastien Loisel <sloisel@gmail.com>
Thu, 28 Apr 2005 10:47:08 +0200
This is good! However, my email has word-wrapped it; can you stick it
somewhere on the web or otherwise attach it?
Sebastien Loisel
On 4/28/05, Jon Harrop <jon@ffconsultancy.com> wrote:
> 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 doin=
g
> > > 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 a=
nd
> > let the masses produce all the optimizations their little hearts desire=
.
>=20
> 99 LOC! That's incredible... for a Haskell programmer. ;-)
>=20
> I've boiled my program down to 94 LOC but managed to keep the hierarchica=
l
> spherical bounding volumes, so it can still render a sphere-flake contain=
ing
> 600,000 spheres in under 8.5 seconds on my cheapy laptop.
>=20
> In contrast, your Haskell program takes 6 minutes and renders significant=
ly
> fewer objects. However, my OCaml code doesn't support checkered planes an=
d
> reflections but it does do shadows.
>=20
> Here's the code:
>=20
> let delta =3D sqrt epsilon_float and pi =3D 4. *. atan 1.
> let sin t =3D sin(pi *. t /. 180.) and cos t =3D cos(pi *. t /. 180.)
> let vec3 x y z =3D [| x; y; z; 1. |] and zero =3D [|0.; 0.; 0.; 1.|]
> let ( *| ) s r =3D [| s *. r.(0); s *. r.(1); s *. r.(2); 1. |]
> let ( +| ) a b =3D [| a.(0) +. b.(0); a.(1) +. b.(1); a.(2) +. b.(2); 1. =
|]
> let ( -| ) a b =3D [| a.(0) -. b.(0); a.(1) -. b.(1); a.(2) -. b.(2); 1. =
|]
> let dot a b =3D a.(0) *. b.(0) +. a.(1) *. b.(1) +. a.(2) *. b.(2)
> let length2 r =3D dot r r and length r =3D sqrt(dot r r)
> let unitise r =3D (1. /. length r) *| r
> let init f =3D Array.init 4 (fun i -> Array.init 4 (f i))
> let trans r =3D init (fun i j -> if i=3Dj then 1. else if j=3D3 then r.(i=
) else 0.)
> let identity =3D trans zero
> let sum4 f =3D f 0 +. f 1 +. f 2 +. f 3
> let mat_mul a b =3D init (fun i j -> sum4 (fun k -> a.(i).(k) *. b.(k).(j=
)))
> let transform m r =3D Array.init 4 (fun i -> dot m.(i) r +. m.(i).(3))
> let rot_x t =3D let cos, sin =3D 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 =3D let cos, sin =3D 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 =3D let cos, sin =3D cos t, sin t in [|[| cos; sin; 0.=
; 0. |];
> [| -.sin; cos; 0.; 0. |];
> [| 0.; 0.; 1.; 0. |];
> [| 0.; 0.; 0.; 1. |]|]
> type sphere =3D { center: float array; radius: float }
> type obj =3D Sphere of sphere * float | Group of sphere * obj list
> type ray =3D { orig: float array; dir: float array }
> let ray_sphere ray sphere =3D
> let v =3D sphere.center -| ray.orig in
> let b =3D dot v ray.dir in
> let disc =3D b *. b -. length2 v +. sphere.radius *. sphere.radius in
> if disc < 0. then infinity else
> let disc =3D 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 =3D
> let rec of_scene ((l, _, _) as first) =3D function
> Sphere (s, material) ->
> let l' =3D ray_sphere ray s in
> if l' >=3D l then first else (* No nearer *)
> let normal =3D unitise (ray.orig +| l' *| ray.dir -| s.center) in
> l', normal, material (* Replace with nearer intersection *)
>=20
> | Group (bound, scenes) ->
>=20
> let l' =3D ray_sphere ray bound in (* Cull if possible *)
> if l' >=3D 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 =3D match intersect ray scene wi=
th
> lambda, n, color ->
> if lambda =3D infinity then 0. else
> let o =3D ray.orig +| lambda *| ray.dir +| delta *| n in
> (match intersect { orig =3D o; dir =3D zero -| light } scene with
> l, _, _ when l =3D infinity -> max 0. (-. dot n light)
>=20
> | _ -> 0.) *. color
>=20
> let () =3D
> let level =3D match Sys.argv with [| _; l |] -> int_of_string l | _ -> =
6 in
> let light =3D unitise (vec3 (-1.) (-3.) 2.) and s =3D 1. /. 3. in
> let scene =3D
> let rec aux level r m =3D
> let sphere =3D { center =3D transform m zero; radius =3D r } in
> let obj =3D Sphere (sphere, 1.) in
> if level =3D 1 then obj else begin
> let m' =3D trans (vec3 0. ((1. +. s) *. r) 0.) in
> let rec aux2 f z i objs =3D
> if i=3D0 then objs else
> let m =3D 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 =3D aux2 (fun i -> 60*i) (rot_z 115.) 6
> (aux2 (fun i -> 30+120*i) (rot_z 45.) 3 []) in
> let aux2 r =3D function Sphere (sph', _) | Group (sph', _) ->
> max r (length (sph'.center -| sphere.center) +. sph'.radius) in
> let r' =3D (* 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 =3D List.fold_left aux2 r' objs}, obj::objs)
> end in
> aux level 1. (trans (vec3 0. 1. 0.)) in
> let w, h =3D 768, 768 in
> Printf.printf "P2\n%d %d\n256\n" w h;
> for y =3D h - 1 downto 0 do
> for x =3D 0 to w - 1 do
> let ray =3D
> let x, y =3D float x -. 0.5 *. float w, float y -. 0.5 *. float h in
> let d =3D unitise (vec3 x y (float (max w h))) in
> {orig =3D vec3 0. 2. (-4.); dir =3D transform (rot_x (-14.)) d} in
> Printf.printf "%d " (int_of_float (256. *.ray_trace 1. light ray
> scene));
> done;
> Printf.printf "\n";
> done;
>=20
> compile it with:
>=20
> ocamlopt ray.ml -o ray
>=20
> then run it with:
>=20
> ./ray >image.pgm
>=20
> --
> Dr Jon D Harrop, Flying Frog Consultancy Ltd.
> Objective CAML for Scientists
> http://www.ffconsultancy.com/products/ocaml_for_scientists
>=20
> _______________________________________________
> Shootout-list mailing list
> Shootout-list@lists.alioth.debian.org
> http://lists.alioth.debian.org/mailman/listinfo/shootout-list
>