[MLton] [Fwd: Port my ray tracer]

Daniel C. Wang danwang@CS.Princeton.EDU
Thu, 02 Jun 2005 04:14:24 -0400


For those of you who don't troll comp.lang.functional.
MLton doing better than FORTRAN!

-------- Original Message --------
Subject: Port my ray tracer
Date: Thu, 02 Jun 2005 02:25:45 +0100
From: Jon Harrop <usenet@jdh30.plus.com>
Organization: Flying Frog Consultancy Ltd.
Newsgroups: comp.lang.functional


I've written a mini ray tracer for the computer language shootout. The
original version was a 66-line OCaml program which I ported to C++:

   http://www.ffconsultancy.com/free/ray_tracer/comparison.html

I have since ported the program to C (icc and gcc), SML (sml-nj and mlton)
and Java (j2se and gcj) and Simon Geard kindly ported it to Fortran 90.

I'd like to see it ported to other languages as well, particularly Haskell
and I'm keen to hear any constructive criticisms of my existing
implementations.

It is interesting to note that the SML is 50% longer than the OCaml. When
compiled with Mlton it is also significantly faster. Here are my timings
(1.2GHz Athlon t-bird) for 128x128 resolution, 4x4 oversampling and 6
levels of spheres:

Mlton    1.316s  mlton ray.sml
IFC      1.361s  ifort -O3 -u -static-libcxa -o raytracer raytracer.f90
C++      1.605s  g++-3.4 -O3 -funroll-all-loops -ffast-math ray.cpp -o ray
ocamlopt 1.932s  ocamlopt -inline 100 ray.ml -o ray
SML-NJ   2.245s  sml ray.sml
G95      3.351s  g95 -O3 -ffast-math ray.f90 -o ray
C        5.971s  gcc-3.4 -lm -std=c99 -O3 -ffast-math ray.c -o ray
Java     6.492s  javac ray.java
GCJ     20.316s  gcj-3.4 --main=ray -Wall -O3 -ffast-math ray.java -o ray
ocamlc  41.047s  ocamlc ray.ml -o ray

Note: the OCaml and C++ implementations timed here are not those on the site
- they are more optimised and longer (keeping within 100 LOC).

The compile times are also interesting:

ocamlc   0.099s
IFC      0.333s
SML-NJ   0.625s
C        0.655s
ocamlopt 0.714s
GCJ      0.723s
G95      0.941s
Java     1.362s
Mlton    8.267s
C++      8.676s

Here is the SML (for Mlton):

(* The Great Computer Language Shootout
    http://shootout.alioth.debian.org/
    contributed by Jon Harrop, 2005
    Compile: mlton ray.sml *)

fun real n = Real.fromInt n
(* Mlton does not provide a "for loop" construct. *)
fun for (s, e, f) = if s=e then () else (f (real s); for (s+1, e, f))

val delta = 0.00000001 (* FIXME: Where is mach eps in the std libs? *)
val infinity = 1.0 / 0.0 (* FIXME: Where is infinity? *)
val pi = 4.0 * Real.Math.atan 1.0 (* FIXME: Where is pi? *)

(* 3D vector
    SML typically uses tuples rather than records (as in the OCaml
    implementation). *)
type vec = real * real * real
(* SML allows operator precedences to be declared when defining infix
    operators. *)
infix 2 *| fun s *| (x, y, z) = (s*x+0.0, s*y, s*z)
infix 1 +| fun (x1, y1, z1) +| (x2, y2, z2) =
                (x1+x2+0.0, y1+y2+0.0, z1+z2+0.0)
infix 1 -| fun (x1, y1, z1) -| (x2, y2, z2) =
                (x1-x2+0.0, y1-y2+0.0, z1-z2+0.0)
fun dot (x1, y1, z1) (x2, y2, z2) = x1*x2 + y1*y2 + z1*z2+0.0
fun unitise r = (1.0 / Real.Math.sqrt (dot r r)) *| r

(* Node in the scene tree *)
datatype scene =
          Sphere of vec * real
        | Group of vec * real * scene list

(* Find the first intersection of the given ray with this sphere *)
fun ray_sphere (orig, dir) center radius =
     let
         val v = center -| orig
         val b = dot v dir
         val disc = b * b - dot v v + radius * radius
     in
         if disc < 0.0 then infinity else
         let val disc = Real.Math.sqrt disc in
             let val t2 = b + disc in
             if t2 < 0.0 then infinity else
             (fn t1 => if t1 > 0.0 then t1 else t2) (b - disc)
             end end end

(* Accumulate the first intersection of the given ray with this sphere *)
(* This function is used both for primary and shadow rays. *)
fun intersect (orig, dir) scene =
     let fun of_scene (scene, (l, n)) =
             case scene of
                 Sphere (center, radius) =>
                 let val l' = ray_sphere (orig, dir) center radius in
                     if l' >= l then (l, n) else
                     (l', unitise (orig +| l' *| dir -| center))
                 end
               | Group (center, radius, scenes) =>
                 let val l' = ray_sphere (orig, dir) center radius in
                     if l' >= l then (l, n) else
                     foldl of_scene (l, n) scenes
                 end in
         of_scene (scene, (infinity, (0.0, 0.0, 0.0)))
     end

(* Trace a single ray into the scene *)
fun ray_trace light (orig, dir) scene =
     let val (lambda, n) = intersect (orig, dir) scene in
         if lambda >= infinity then 0.0 else
         let val g = 0.0 - dot n light in
             (* If we are on the shadowed side of a sphere then don't
                bother casting a shadow ray as we know it will intersect
                the same sphere. *)
             if g <= 0.0 then 0.0 else
             let
                 val orig = orig +| lambda *| dir +| delta *| n
                 val dir = (0.0, 0.0, 0.0) -| light
             in
                 let val (l, _) = intersect (orig, dir) scene in
                     if l >= infinity then g else 0.0
                 end end end end

(* Recursively build the scene tree *)
fun create level r (x, y, z) =
     let val obj = Sphere ((x, y, z), r) in
         if level = 1 then obj else
         let val r' = 3.0 * r / Real.Math.sqrt 12.0 in
             let fun aux x' z' =
                     create (level-1) (0.5 * r) (x-x', y+r', z+z') in
                 let val objs = [aux (~r') (~r'), aux r' (~r'),
                                 aux (~r') r', aux r' r', obj] in
                     Group ((x, y, z), 3.0 * r, objs)
                     end end end end

(* Build a scene and trace many rays into it, outputting a PGM image *)
val () =
     let
         val level = 6 (* Number of levels of spheres *)
         val ss = 4 (* Oversampling *)
         (* Resolution *)
         val n = case CommandLine.arguments () of
                     [s] => (case Int.fromString s of
                                 SOME n => n | _ => 256)
                   | _ => 256
         val scene = create level 1.0 (0.0, ~1.0, 0.0) (* Scene tree *)
     in
         (fn s => print ("P5\n"^s^" "^s^"\n255\n")) (Int.toString n);
         for (0, n, fn y => for (0, n, fn x =>
             let val g = ref 0.0 in
                 for (0, ss, fn dx => for (0, ss, fn dy =>
                     let val n = real n
                         val x = x + dx / real ss - n / 2.0
                         val y = n - 1.0 - y + dy / real ss - n / 2.0 in
                         let val eye = unitise (~1.0, ~3.0, 2.0)
                             val ray = ((0.0, 0.0, ~4.0),
                                        unitise (x, y, n)) in
                             g := !g + ray_trace eye ray scene
                         end
                     end));
                 let val g = 0.5 + 255.0 * !g / real (ss*ss) in
                     print (String.str(Char.chr(Real.trunc g)))
                 end
             end))
     end