★ /prog/ Challenge Summer Edition ★

THE CHALLENGE:
I) Write three functions, union :: [(Float,Float)] -> [(Float,Float)], difference :: [(Float,Float)] -> [(Float,Float)], intersection :: [(Float,Float)] -> [(Float,Float)]
Which perform the respective set operations on their input, a list of closed intervals of the form (min, max). Your functions need not support sentinel values such as NaN or infinity. The number of intervals returned should be minimized where possible.

II) Use your functions from I to implement a simple ^CSG _raytracer. Your raytracer should support the following features:

• Spheres
• Diffuse surfaces
• Specular surfaces

Upon completion, post your source code here, along with one (1) example rendering, showcasing ^your raytracer. The example render should include an example of each set operation performed on two spheres, for a total of six spheres. A single directional light source should be used to light the scene.

Submissions close on 2010-10-04 at 00:00. Programs will be marked in three categories: efficiency, conciseness, and quality of implementation.

Each winner is entitled to claim ownership of two (2) /prog/ memes of his choice, with the exception of SICP, HMA, and UMH.

Here's my completed entry.

Image: http://img706.imageshack.us/img706/751/testsl.png

type color = {r: int; g: int; b: int}
type vector = Vector of float * float * float
type solid = Sphere of vector * float * color
type interval = Interval of float * float * solid
type csg =
    Union of csg * csg
  | Difference of csg * csg
  | Intersection of csg * csg
  | Primitive of solid
type scene = Scene of vector * csg

(* boolean operations *)
let rec union a b =
  match (a, b) with
      (None, _) -> b
    | (_, None) -> a
    | (Some (Interval(a_l, a_h, a_s)), Some (Interval(b_l, b_h, b_s))) ->
      if a_l < b_l then a else b

let rec difference a b =
  match (a, b) with
      (None, _) -> b
    | (_, None) -> a
    | (Some (Interval(a_l, a_h, a_s)), (Some (Interval(b_l, b_h, b_s)))) ->
      if a_h < b_l then a else
      if b_h < a_l then b else
      if a_l < b_l then Some (Interval(a_l, b_l, a_s)) else
      if b_l < a_l then Some (Interval(b_l, a_l, b_s)) else
      None

let rec intersection a b =
  match (a, b) with
      (None, _) -> None
    | (_, None) -> None
    | (Some (Interval(a_l, a_h, a_s)), Some (Interval(b_l, b_h, b_s))) ->
      if a_h < b_l or b_h < a_l then None else
      let l = max a_l b_l in
      let h = min a_h b_h in
      if l == h then None else
      Some (Interval(l, h, if a_l <= b_l then b_s else a_s))

(* vectors *)
let vector_sub (Vector(a_x, a_y, a_z)) (Vector(b_x, b_y, b_z)) =
  Vector(a_x -. b_x, a_y -. b_y, a_z -. b_z)

let vector_dot (Vector(a_x, a_y, a_z)) (Vector(b_x, b_y, b_z)) =
  a_x *. b_x +. a_y *. b_y +. a_z *. b_z

let vector_mul (Vector(a_x, a_y, a_z)) k =
  Vector(a_x *. k, a_y *. k, a_z *. k) 

let vector_norm (Vector(a_x, a_y, a_z)) =
  sqrt(a_x *. a_x +. a_y *. a_y +. a_z *. a_z)

let vector_normalize v =
  vector_mul v (1. /. (vector_norm v))

(* line/sphere intersection *)
let line_sphere o i s =
  match s with Sphere(c, r, x) ->
    let c' = vector_sub c o in
    let i_c' = vector_dot i c' in
    let det = i_c' ** 2. -. (vector_dot c' c') +. r ** 2. in
    if det < 0. then None else
      let sqrt_det = sqrt det in
      Some (Interval(i_c' -. sqrt_det, i_c' +. sqrt_det, s))

(* normal calculation *)
let point_normal p (Sphere(c, _, _)) =
  vector_normalize (vector_sub p c)

(* ray / csg intersection *)
let clip i =
  match i with
      None -> None
    | Some (Interval(l, h, s)) ->
      if h < 0. then
        None
      else
        Some (Interval(max 0. l, h, s))

let primitive_ray o i s =
  match s with Sphere _ -> line_sphere o i s

let rec ray o i csg =
  match csg with
      Union(a, b) -> union (ray o i a) (ray o i b)
    | Difference(a, b) -> difference (ray o i a) (ray o i b)
    | Intersection(a, b) -> intersection (ray o i a) (ray o i b)
    | Primitive(s) -> clip (primitive_ray o i s)

(* pixel -> color *)
let view x y (Scene(light, csg)) =
  let o = Vector(x,y,0.) in
  let i = Vector(0.,0.,1.) in
  match ray o i csg with
      None -> {r=0;g=0;b=0}
    | Some Interval(z, _, s) ->
      let p = vector_sub o (vector_mul i (z *. -1.)) in
      let normal = point_normal p s in
      let r = vector_sub i (vector_mul normal (2. *. (vector_dot normal i))) in
      let spec_lum = max 0. (vector_dot light r) in
      let diff_lum = max 0. (vector_dot light normal) in
      let lum = spec_lum ** 4. *. 0.4 +. diff_lum ** 2. *. 0.6 in
      match s with Sphere(_d, _, c) ->
        let k = int_of_float (50. +. lum *. 205.) in
        {
          r = (c.r * k) / 255;
          g = (c.g * k) / 255;
          b = (c.b * k) / 255;
        }

(* csg -> ppm *)
let ppm w h scene =
  print_string "P6 ";
  print_int w;
  print_string " ";
  print_int h;
  print_string " 255";
  print_newline ();

  let w' = float_of_int w in
  let h' = float_of_int h in

  let rec yloop y =
    let rec xloop x =
      if x >= w' then yloop (y +. 1.) else
        let color = view ((x /. w') *. 3.) (y /. h') scene in
        print_char (char_of_int color.r);
        print_char (char_of_int color.g);
        print_char (char_of_int color.b);
        xloop (x +. 1.)
    in if y >= h' then () else xloop 0.
  in yloop 0.

(* Test scene *)
let sphere1 = Sphere(Vector(0.35,0.35,0.5), 0.3, {r = 255; g = 255; b = 0})
let sphere2 = Sphere(Vector(0.55,0.55,0.7), 0.4, {r = 0; g = 0; b = 255})
let sphere3 = Sphere(Vector(1.35,0.35,0.5), 0.3, {r = 255; g = 255; b = 0})
let sphere4 = Sphere(Vector(1.55,0.55,0.7), 0.4, {r = 0; g = 0; b = 255})
let sphere5 = Sphere(Vector(2.35,0.35,0.5), 0.3, {r = 255; g = 255; b = 0})
let sphere6 = Sphere(Vector(2.55,0.55,0.7), 0.4, {r = 0; g = 0; b = 255})
let csg1 = Union(Primitive(sphere1), Primitive(sphere2))
let csg2 = Difference(Primitive(sphere3), Primitive(sphere4))
let csg3 = Intersection(Primitive(sphere5), Primitive(sphere6))
let csg = Union(csg1, Union(csg2, csg3))
let light = vector_normalize (Vector(-1.0,-1.0,-2.0))
let scene = Scene(light, csg)

let _ = ppm 1500 500 scene