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

   Removed comments.

   12/31/09
                                                                           *)

let output_size = 300;;

let supersample = 4;;
let supersample_scale = float_of_int supersample;;
let sample_squared = supersample * supersample;;

let flower_spacing = 75.0 *. supersample_scale;;

let supersample_size = output_size * supersample;;

let outline_color = (0, 0, 0);;

let pi = (atan2 0. (-1.0));;
let step_to_radian = 0.2 *. pi;;

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

let image =
   Array.make_matrix supersample_size supersample_size (124, 194, 169);;

let to_radian = pi /. 180.0;;

let rotate_path path angle =
   let ca = cos (angle *. to_radian) in
   let sa = 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
      List.rev result
   else
      let step_i x y = x +. (float_of_int i) *. (y -. x) /. 64.0 in
      let ab = step_i a b in
      let bc = step_i b c in
      let cd = step_i c d in
      let abc = step_i ab bc in
      let bcd = 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 (truncate 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 =
   List.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 = get_range xlist in
   let min_y, max_y = get_range ylist in
   for y = (truncate min_y) to (truncate (max_y +. 0.5)) do
      rasterize_scanline path color (float_of_int y) min_x max_x
   done;;

let rasterize_list paths color =
   List.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 = Array.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 downsample_scanline_set output_y =
   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;;

let downsample () =
   for y = 0 to (output_size - 1) do
      downsample_scanline_set y
   done;;

let output_ppm () =
   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;;

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

let flower_row_positions x0 y =
   let flower_count =
      truncate ((float_of_int supersample_size) /. flower_spacing) in
   Array.init
      flower_count
      (fun index ->
         (x0 +. (float_of_int index) *. flower_spacing +. random_jitter(),
          y +. random_jitter()));;

let row_count = truncate ((float_of_int supersample_size) /. flower_spacing);;
let flower_grid_positions () =
   Array.init
      row_count
      (fun index ->
         flower_row_positions
            (random_jitter())
            ((float_of_int 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_int supersample_size, y);;

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

let average_positions row1 row2 =
   let count = Array.length row1 in
   Array.init
      count
      (fun index ->
         if index = 0 then
            fill_gap [row1.(count - 1); horizontal_shift (row1.(0));
                      row2.(count - 1); horizontal_shift (row2.(0))]
         else
            fill_gap [row1.(index); row1.(index - 1);
                      row2.(index); row2.(index - 1)]);;

let average_grid_positions positions =
   let count = Array.length positions in
   Array.init
      count
      (fun index ->
         if index = 0 then
            average_positions (positions.(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_int index) *. 72.0 +. angle))
               tx ty
         ) in
         List.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 List.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;;

Random.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;;

downsample();;
output_ppm();;