jdh_ray_tracer.ml

module K = struct
  let ( * ) = ( *. )
  let ( + ) = ( +. )
  let ( / ) = ( /. )
  let ( - ) = ( -. )
  let ( ~- ) = ( ~-. )
end

module V = struct
let zero = 0., 0., 0.
let ( * ) s (rx, ry, rz) = K.( s * rx, s * ry, s * rz )
let ( + ) (ax, ay, az) (bx, by, bz) = K.( ax + bx, ay + by, az + bz )
let ( - ) (ax, ay, az) (bx, by, bz) = K.( ax - bx, ay - by, az - bz )
let ( / )  (rx, ry, rz) s  = K.( rx / s, ry / s, rz / s )
let ( |*| ) (ax, ay, az) (bx, by, bz) = K.( ax * bx + ay * by + az * bz )
let unitise r =  r / sqrt (r |*| r )
end

let rec intersect o d (l, _ as hit) (c, r, s) =
  let open V in
  let v = c - o in
  let b = ( v |*| d ) in
  let disc =  sqrt K.( b * b - ( v |*| v ) + r * r ) in
  let t1 = K.( b - disc ) and t2 = K.( b + disc ) in
  let l' = if t2>0. then if t1>0. then t1 else t2 else infinity in
  if l' >= l then hit else match s with
    [] -> l', unitise (o + l' * d - c)
  | ss -> List.fold_left (intersect o d) hit ss

let light = V.unitise (1., 3., -2.) and ss = 4

let rec create level c r =
  let obj = c, r, [] in
  if level = 1 then obj else
    let a = K.( 3. * r / sqrt 12. ) in
    let aux x' z' = create (level - 1) V.(c + (x', a, z')) K.(0.5 * r) in
    c, K.(3. * r ), K.[obj; aux (-a) (-a); aux a (-a); aux (-a) a; aux a a]

let level, n =
  try int_of_string Sys.argv.(1), int_of_string Sys.argv.(2) with _ -> 9, 512

let scene = create level (0., -1., 4.) 1.

let rec ray_trace dir =
  let open V in
  let l, n = intersect zero dir (infinity, zero) scene in
  let g = ( n |*| light ) in
  if g <= 0. then 0. else
    let p = V.(l * dir + sqrt epsilon_float * n) in
    if fst(intersect p light (infinity, zero) scene) < infinity then 0. else g

let aux x d = K.( float x - float n / 2. + float d / float ss );;
Printf.printf "P5\n%d %d\n255\n%!" n n;
let ss2 = ss * ss in
for y = n - 1 downto 0 do
  for x = 0 to n - 1 do
    let g = ref 0. in
    for d = 0 to ss2 - 1 do
      g := !g +. ray_trace (V.unitise(aux x (d mod ss), aux y (d/ss), float n))
    done;
    let g = K.( 0.5 + 255. * !g / float ss2 ) in
    print_char @@ char_of_int @@ int_of_float g
  done;
done