[Shootout-list] Ray tracer

[Jon Harrop]

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

On Thursday 28 April 2005 20:41, Brent Fulgham wrote:
> ...Maybe 120 lines is a better number...

Ye of little faith. ;-)

I have now ported the ray tracer to a 97-line C++. I've changed the task 
slightly to do supersampling in order to increase the amount of CPU work done 
per byte of output (to make it a better benchmark). So there is a new OCaml 
implementation too (which is slightly faster than the last).

The C++ version uses natural representations. Specifically, the OCaml used a 
variant type and pattern matched functions over this variant type:

type obj = Sphere of vec * float * float | Group of vec * float * obj array

let intersect ray scene =
  let rec of_scene ((l, _, _) as first) = function
      Sphere (center, radius, material) -> ...
    | Group (center, radius, scenes) -> ...

In contrast, the C++ version uses an inheritance hierarchy instead of the 
variant type, and implements each expression corresponding to each pattern in 
the appropriate derived class:

struct Scene { // Abstract base class representing a scene
  virtual void intersect(double &, Vec &, Ray) = 0;
  virtual void del() {};
};

struct Sphere : public Scene { // Derived class representing a sphere
  Vec center;
  double radius;
  ...
  void intersect(double &lambda, Vec &normal, Ray ray) { ... }
}

struct Group : public Scene { // Derived class representing a group of scenes
  Sphere bound;
  vector<Scene *> objs;
  void intersect(double &lambda, Vec &normal, Ray ray) { ... }
};

Note that the "intersect" function was a purely functional expression in OCaml 
but acts by way of side-effects in C++ (altering the "lambda" and "normal" 
variables passed to it by reference).

I've done some timings and the results are very wierd. On my AMD64, the C++ 
and OCaml implementations perform similarly (17.281s in C++ and 19.428s in 
OCaml). On my 32-bit Athlon I was expecting OCaml to be half the speed of the 
C++ as ocamlopt is much better at optimising floating point arithmetic for 
AMD64 than for x86. However, the OCaml is actually more than twice as fast 
(92.805s in C++ vs 44.862s in OCaml).

The C++ could be optimised in various ways but this would bloat the code and 
would be fiddly to do.

Given that this can be done in <100 lines of C++, I'd imagine it can also be 
done in Java and C#. For one thing, a few lines are dedicated to memory 
deallocation.

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

--Boundary-00=_qmZcC/AKsXUTR7l
Content-Type: text/x-c++src;
  charset="iso-8859-1";
  name="ray.cpp"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
	filename="ray.cpp"

#include <vector>
#include <iostream>
#include <limits>
#include <cmath>
using namespace std;
numeric_limits<double> dbl;
double delta = sqrt(dbl.epsilon()), infinity = dbl.infinity(), pi = M_PI;
struct Vec { // 3D vector
  double x, y, z;
  Vec(double x2, double y2, double z2) : x(x2), y(y2), z(z2) {}
};
Vec operator+(Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y, a.z + b.z); }
Vec operator-(Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y, a.z - b.z); }
Vec operator*(double a, Vec b) { return Vec(a * b.x, a * b.y, a * b.z); }
double dot(Vec a, Vec b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
Vec unitise(Vec a) { return (1 / sqrt(dot(a, a))) * a; }
struct Ray { Vec orig, dir; Ray(Vec o, Vec d) : orig(o), dir(d) {} };
struct Scene { // Abstract base class representing a scene
  virtual void intersect(double &, Vec &, Ray) = 0;
  virtual void del() {};
};
struct Sphere : public Scene { // Derived class representing a sphere
  Vec center;
  double radius;
  Sphere(Vec c, double r) : center(c), radius(r) {}
  double ray_sphere(Ray ray) { // Intersection of a ray with a sphere
    Vec v = center - ray.orig;
    double b = dot(v, ray.dir), disc = b*b - dot(v, v) + radius * radius;
    if (disc < 0) return infinity;
    double d = sqrt(disc), t2 = b + d;
    if (t2 < 0) return infinity;
    double t1 = b - d;
    return (t1 > 0 ? t1 : t2);
  }
  void intersect(double &lambda, Vec &normal, Ray ray) {
    double l = ray_sphere(ray);
    if (l >= lambda) return;
    lambda = l;
    normal = unitise(ray.orig + l * ray.dir - center);
  }
};
struct Group : public Scene { // Derived class representing a group of scenes
  Sphere bound;
  vector<Scene *> objs;
  Group(Sphere b) : bound(b), objs(0) {}
  void del() {
    for (vector<Scene *>::iterator it=objs.begin(); it!=objs.end(); ++it)
      delete *it;
  }
  void intersect(double &lambda, Vec &normal, Ray ray) {
    double l = bound.ray_sphere(ray);
    if (l >= lambda) return;
    for (vector<Scene *>::iterator it=objs.begin(); it!=objs.end(); ++it)
      (*it)->intersect(lambda, normal, ray);
  }
};
Vec light = unitise(Vec(-1, -3, 2));
double ray_trace(double weight, Ray ray, Scene *scene) { // Trace a ray
  double lambda = infinity;
  Vec normal(0, 0, 0);
  scene->intersect(lambda, normal, ray);
  if (lambda == infinity) return 0;
  Vec o = ray.orig + lambda * ray.dir + delta * normal;
  double g = max(0., -dot(normal, light)), l = infinity;
  scene->intersect(l, normal, Ray(o, Vec(0, 0, 0) - light));
  return (l == infinity ? g : 0);
}
Scene *create(int level, double r, double x, double y, double z) {
  Sphere *sphere = new Sphere(Vec(x, y, z), r);
  if (level == 1) return sphere;
  Group group = Group(Sphere(Vec(x, y, z), 3*r));
  group.objs.push_back(sphere);
  double rn = 3*r/sqrt(12.);
  for (int dz=-1; dz<=1; dz+=2)
    for (int dx=-1; dx<=1; dx+=2)
      group.objs.push_back(create(level-1, r/2, x - dx*rn, y + rn, z - dz*rn));
  return (new Group(group));
}
int main(int argc, char *argv[]) {
  int level = (argc==2 ? atoi(argv[1]) : 6), w = 512, h = 512, ss = 4;
  Scene *scene=create(level, 1, 0, -1, 0); // Build the scene
  cout << "P2\n" << w << " " << h << "\n256\n";
  for (int y=h-1; y>=0; --y) {
    for (int x=0; x<w; ++x) {
      double g=0;
      for (int dx=0; dx<ss; ++dx)
	for (int dy=0; dy<ss; ++dy) {
	  Vec d(x+double(dx)/ss-w/2, y+double(dy)/ss-h/2, max(w, h));
	  g += ray_trace(1, Ray(Vec(0, 0, -4), unitise(d)), scene);
	}
      cout << int(256. * g / (ss*ss)) << " ";
    }
    cout << endl;
  }
  scene->del(); // Deallocate the scene
  return 0;
}

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

let delta = sqrt epsilon_float and pi = 4. *. atan 1.
type vec = {x:float; y:float; z:float}
let ( *| ) s r = {x=s *. r.x; y=s *. r.y; z=s *. r.z}
let ( +| ) a b = {x=a.x +. b.x; y=a.y +. b.y; z=a.z +. b.z}
let ( -| ) a b = {x=a.x -. b.x; y=a.y -. b.y; z=a.z -. b.z}
let dot a b = a.x *. b.x +. a.y *. b.y +. a.z *. b.z
let unitise r = (1. /. sqrt (dot r r)) *| r
type sphere = { center: vec; radius: float } and ray = { orig: vec; dir: vec }
type obj = Sphere of vec * float * float | Group of vec * float * obj array
let ray_sphere ray center radius =
  let v = center -| ray.orig in
  let b = dot v ray.dir in
  let disc = b *. b -. dot v v +. radius *. 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 (center, radius, material) ->
	let l' = ray_sphere ray center radius in
	if l' >= l then first else (* No nearer *)
	  let normal = unitise (ray.orig +| l' *| ray.dir -| center) in
	  l', normal, material (* Replace with nearer intersection *)
    | Group (center, radius, scenes) ->
	let l' = ray_sphere ray center radius in (* Cull if possible *)
	if l' >= l then first else Array.fold_left of_scene first scenes in
  of_scene (infinity, {x=0.; y=0.; z=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 = {x=0.; y=0.; z=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 {x= -1.; y= -3.; z=2.} in
  let scene =
    let rec aux level r (x, y, z) =
      let obj = Sphere ({x=x; y=y; z=z}, r, 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' = 3. *. r /. sqrt 12. in
	Array.fold_left aux [obj]
	  [|-.r', r', -.r'; r', r', -.r'; -.r', r', r'; r', r', r'|] in
	Group ({x=x; y=y; z=z}, 3. *. r, Array.of_list objs)
      end in
    aux level 1. (0., -1., 0.) in
  let w, h = 512, 512 and ss = 4 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 g = ref 0. in
      for dx = 0 to ss - 1 do
	for dy = 0 to ss - 1 do
	  let ray =
	    {orig = {x=0.; y=0.; z= -4.};
	     dir = unitise {x = float (x - w / 2) +. float dx /. float ss;
			    y = float (y - h / 2) +. float dy /. float ss;
			    z=float (max w h)} } in
	  g := !g +. ray_trace 1. light ray scene;
	done;
      done;
      Printf.printf "%d " (int_of_float (256. *. !g /. float (ss * ss)));
    done;
    Printf.printf "\n";
  done;

--Boundary-00=_qmZcC/AKsXUTR7l--