(* kazari6.ml - Don Yang (uguu.org)

   Use tuple assignment.

   01/01/10
                                                                           *)

let output_size,
    supersample,
    pi,
    outline_color,
    float_of_int0,
    truncate0 =
    300,
    4,
    (atan2 0. (-1.0)),
    (0, 0, 0),
    float_of_int,
    truncate;;

module L = List;;
module A = Array;;
module R = Random;;

(* These functions can't be combined with other let statements due to
   http://caml.inria.fr/mantis/view.php?id=4951                       *)
let list_map,
    list_split,
    array_fill,
    array_init,
    array_iter,
    array_iteri,
    random_float,
    printf_printf =
    L.map,
    L.split,
    A.fill,
    A.init,
    A.iter,
    A.iteri,
    R.float,
    Printf.printf;;

let supersample_scale,
    sample_squared,
    supersample_size,
    to_radian =
    float_of_int0 supersample,
    supersample * supersample,
    output_size * supersample,
    pi /. 180.0;;

let flower_spacing,
    image =
    75.0 *. supersample_scale,
    A.make_matrix supersample_size supersample_size (124, 194, 169);;

let step_to_radian x = 0.2 *. pi *. float_of_int0 x;;

let make_circle_path r =
   let rec make_path_recurse step =
      if step = 0 then []
      else
      (
         (r *. cos (step_to_radian step),
          r *. sin (step_to_radian step)) ::
         make_path_recurse (step - 1)
      ) in
   make_path_recurse 10;;

let rotate_path path angle =
   let ca, sa = cos (angle *. to_radian), sin (angle *. to_radian) in
   list_map (fun (x, y) -> ca *. x -. sa *. y, sa *. x +. ca *. y) path;;

let translate_path path tx ty =
   list_map (fun (x, y) -> x +. tx, y +. ty) path;;

let rec interpolate a b c d i result =
   if i = 64 then
      L.rev result
   else
      let step_i x y = x +. (float_of_int0 i) *. (y -. x) /. 64.0 in
      let ab, bc, cd = step_i a b, step_i b c, step_i c d in
      let abc, bcd = step_i ab bc, step_i bc cd in
      interpolate a b c d (i + 1) ((step_i abc bcd) :: result);;

let interpolate_list x =
   match x with [a; b; c; d] -> interpolate a b c d 0 [d]
   | _ -> [];;

let dx a b = (fst b) -. (fst a);;
let dy a b = (snd b) -. (snd a);;
let clockwise a b c = (dx a c) *. (dy a b) >= (dx a b) *. (dy a c);;

let inside path point =
   match path with first_point::a::b ->
   (
      if not (clockwise first_point a point) then
         false
      else
         let rec inside_recurse next_point tail_points =
            match tail_points
            with [] ->
               clockwise next_point first_point point
            | ia::ib ->
            (
               if not (clockwise next_point ia point) then
                  false
               else
                  inside_recurse ia ib
            ) in
         inside_recurse a b
   )
   | _ -> false;;

let stop_recurse min_x max_x = max_x -. min_x < 0.5;;

let avg a b = (a +. b) /. 2.0;;

let rec find_point_on_scanline path y min_x max_x =
   if stop_recurse min_x max_x then
      raise Not_found
   else
      let x = avg min_x max_x in
      if inside path (x, y) then
         x
      else
         try find_point_on_scanline path y min_x x
         with Not_found -> find_point_on_scanline path y x max_x;;

let rec intersect_scanline_left path y min_x max_x =
   if stop_recurse min_x max_x then
      max_x
   else
      let x = avg min_x max_x in
      if inside path (x, y) then
         intersect_scanline_left path y min_x x
      else
         intersect_scanline_left path y x max_x;;

let rec intersect_scanline_right path y min_x max_x =
   if stop_recurse min_x max_x then
      min_x
   else
      let x = avg min_x max_x in
      if inside path (x, y) then
         intersect_scanline_right path y x max_x
      else
         intersect_scanline_right path y min_x x;;

let wraparound value =
   if value < 0 then
      supersample_size + (value mod supersample_size)
   else
      value mod supersample_size;;

let wraparound_int value =
   wraparound (truncate0 value);;

let intersect_scanline path y min_x max_x =
   let start_x = find_point_on_scanline path y min_x max_x in
   (
      wraparound_int (intersect_scanline_left path y min_x start_x),
      wraparound_int (intersect_scanline_right path y start_x max_x)
   );;

let get_range values =
   L.fold_left (fun (min_x, max_x) x -> (min x min_x, max x max_x))
               (infinity, neg_infinity)
               values;;

let rasterize_scanline path color y min_x max_x =
   try
      let tx0, tx1 = intersect_scanline path y min_x max_x in

      let sy = wraparound_int y in

      if tx0 <= tx1 then
         array_fill (image.(sy)) tx0 (tx1 - tx0) color
      else
      (
         array_fill (image.(sy)) 0 (tx1 + 1) color;
         array_fill (image.(sy)) tx0 (supersample_size - tx0) color
      )
   with Not_found -> ();;

let rasterize path color =
   let xlist, ylist = list_split path in
   let (min_x, max_x), (min_y, max_y) = (get_range xlist), get_range ylist in
   for y = (truncate0 min_y) to (truncate0 (max_y +. 0.5)) do
      rasterize_scanline path color (float_of_int0 y) min_x max_x
   done;;

let rasterize_list paths color =
   L.iter (fun x -> rasterize x color) paths;;

let add_pixel p1 p2 =
   let r1, g1, b1 = p1 in
   let r2, g2, b2 = p2 in
   (r1 + r2, g1 + g2, b1 + b2);;

let scale_components pixel =
   let r, g, b = pixel in
   (r / sample_squared, g / sample_squared, b / sample_squared);;

let downsample_scanline scanline =
   let output = A.make output_size (scanline.(0)) in
   for x = 0 to (output_size - 1) do
      let rec add_sample index =
      (
         if index = 1 then
            scanline.(x * supersample)
         else
            add_pixel (scanline.(x * supersample + index - 1))
                      (add_sample (index - 1))
      ) in

      output.(x) <- add_sample supersample
   done;
   output;;

let random_jitter () =
    (random_float (flower_spacing *. 0.4)) -. (flower_spacing *. 0.2);;

let flower_row_positions x0 y =
   let flower_count =
      truncate0 ((float_of_int0 supersample_size) /. flower_spacing) in
   array_init
      flower_count
      (fun index ->
         (x0 +. (float_of_int0 index) *. flower_spacing +. random_jitter(),
          y +. random_jitter()));;

let row_count = truncate0 ((float_of_int0 supersample_size) /. flower_spacing);;
let flower_grid_positions () =
   array_init
      row_count
      (fun index ->
         flower_row_positions
            (random_jitter())
            ((float_of_int0 index) *. flower_spacing));;

let fill_gap points =
   let rec fill_gap_component values sum =
      match values with a::b -> fill_gap_component b (sum +. a)
      | [] -> (sum /. 4.0) +. random_jitter() in
   let xlist, ylist = list_split points in
   (fill_gap_component xlist 0.0, fill_gap_component ylist 0.0);;

let horizontal_shift point =
   let x, y = point in (x +. float_of_int0 supersample_size, y);;

let vertical_shift row =
   A.map (fun (x, y) -> x, y +. float_of_int0 supersample_size) row;;

let average_positions row1 row2 =
   array_init
      row_count
      (fun index ->
         if index = 0 then
            fill_gap [row1.(row_count - 1); horizontal_shift (row1.(0));
                      row2.(row_count - 1); horizontal_shift (row2.(0))]
         else
            fill_gap [row1.(index); row1.(index - 1);
                      row2.(index); row2.(index - 1)]);;

let average_grid_positions positions =
   array_init
      row_count
      (fun index ->
         if index = 0 then
            average_positions (positions.(row_count - 1))
                              (vertical_shift (positions.(0)))
         else
            average_positions (positions.(index))
                              (positions.(index - 1)));;

let make_petals tx ty angle path =
   let scaled_path =
      list_map (fun (x, y) -> (x *. supersample_scale,
                               y *. supersample_scale))
               path in
   let rec rotate_and_duplicate_path index =
      if index = 0 then
         []
      else
      (
         let xlist, ylist = list_split 
         (
            translate_path
               (rotate_path scaled_path
                            ((float_of_int0 index) *. 72.0 +. angle))
               tx ty
         ) in
         L.combine (interpolate_list xlist) (interpolate_list ylist)
      ) :: (rotate_and_duplicate_path (index - 1)) in
   rotate_and_duplicate_path 5;;

let draw_flower outline_path fill_path fill_color tx ty angle =
   let make_petals_fixed = make_petals tx ty angle in
   let outline = make_petals_fixed outline_path in
   let fill = make_petals_fixed fill_path in
   let background = list_map L.hd fill in
   rasterize background fill_color;
   rasterize_list outline outline_color;
   rasterize_list fill fill_color;;

let large_flower tx ty angle =
   draw_flower
      [
         (15.24, -8.52);
         (39.81, -51.51);
         (-34.94, -52.33);
         (-9.32, -11.85)
      ]
      [
         (14.42, -6.64);
         (37.08, -48.54);
         (-31.98, -48.81);
         (-7.94, -8.74)
      ]
      (233, 185, 185)
      tx ty angle;
   rasterize
      (translate_path
         (make_circle_path (5.0 *. supersample_scale))
         tx ty)
      outline_color;
   rasterize
      (translate_path
         (make_circle_path (3.5 *. supersample_scale))
         tx ty)
      (221, 218, 167)
      ;;

let small_flower tx ty angle =
   draw_flower
      [
         (3.29, -0.75);
         (19.61, -27.6);
         (-4.78, -24.2);
         (-1.8, -2.5)
      ]
      [
         (2.46, 1.32);
         (17.7, -25.42);
         (-4.47, -21.44);
         (-1.92, -0.62)
      ]
      (234, 236, 237)
      tx ty angle;;

let draw_flower_set flower_func positions =
   array_iter
      (fun row ->
         array_iter
         (fun (x, y) -> flower_func x y (random_float 360.0))
         row)
      positions;;

R.self_init();;

let large_positions = flower_grid_positions();;
let small_positions = average_grid_positions large_positions;;

draw_flower_set large_flower large_positions;;
draw_flower_set small_flower small_positions;;

for output_y = 0 to (output_size - 1) do
(
   for i = 0 to (supersample - 1) do
      let input_y = output_y * supersample + i in
      image.(input_y) <- downsample_scanline (image.(input_y))
   done;

   image.(output_y) <- image.(output_y * supersample);
   let first_scanline = image.(output_y) in
   for i = 1 to (supersample - 1) do
      array_iteri
         (fun x pixel ->
            first_scanline.(x) <- add_pixel (first_scanline.(x)) pixel)
         (image.(output_y * supersample + i))
   done;

   array_iteri
      (fun x pixel -> first_scanline.(x) <- scale_components pixel)
      first_scanline
)
done;;

printf_printf "P3\n%d %d\n255\n" output_size output_size;
array_iteri
   (
      fun y scanline ->
         if y < output_size then
            array_iter
               (fun (r, g, b) -> printf_printf "%d %d %d\n" r g b) scanline
         else
            ()
   )
   image;;