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

   Obscure tweaks.

   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;;
exception Not_Intersect;;

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 =
    float_of_int0 supersample,
    supersample * supersample,
    output_size * supersample;;

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

let rec make_circle_path r step =
   if step = 0 then []
   else
   (
      let angle = 0.2 *. pi *. float_of_int0 step in
      (r *. cos angle, r *. sin angle) :: (make_circle_path r (step - 1))
   );;

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

let interpolate_list x =
   match x with [a; b; c; d] ->
      let rec interpolate 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 (i + 1) ((step_i abc bcd) :: result)
         ) in
      interpolate 0 [d]
   | _ -> [];;

let dx a b = (fst b) -. (fst a);;
let dy a b = (snd b) -. (snd a);;
let not_clockwise a b c = (dx a c) *. (dy a b) < (dx a b) *. (dy a c);;

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

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

let wraparound_int value =
   let mod_value = (truncate0 value) mod supersample_size in
      if mod_value < 0 then
         supersample_size + mod_value
      else
         mod_value;;

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_list path_list color =
   L.iter
      (
         fun path ->
            let xlist, ylist = list_split path in
            let (path_min_x, path_max_x), (path_min_y, path_max_y) =
               (get_range xlist), get_range ylist in
            for iy = (truncate0 path_min_y) to
                     (truncate0 (path_max_y +. 0.5)) do
               let y = float_of_int0 iy in
               let inside x =
                  match path with first_point::a::b ->
                  (
                     if not_clockwise first_point a (x, y) then
                        false
                     else
                        let rec inside_recurse next_point tail_points =
                           match tail_points
                           with [] ->
                              not_clockwise next_point (x, y) first_point
                           | ia::ib ->
                           (
                              if not_clockwise next_point ia (x, y) then
                                 false
                              else
                                 inside_recurse ia ib
                           ) in
                        inside_recurse a b
                  )
                  | _ -> false in

               let rec find_point_on_scanline min_x max_x =
                  if stop_recurse min_x max_x then
                     raise Not_Intersect
                  else
                     let x = avg min_x max_x in
                     if inside x then
                        x
                     else
                        try find_point_on_scanline min_x x
                        with Not_Intersect -> find_point_on_scanline x max_x in

               let rec intersect_scanline_left 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 x then
                        intersect_scanline_left min_x x
                     else
                        intersect_scanline_left x max_x in

               let rec intersect_scanline_right 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 x then
                        intersect_scanline_right x max_x
                     else
                        intersect_scanline_right min_x x in

               try
                  let start_x = find_point_on_scanline path_min_x path_max_x in
                  let tx0, tx1, sy = 
                     wraparound_int
                        (intersect_scanline_left path_min_x start_x),
                     wraparound_int
                        (intersect_scanline_right start_x path_max_x),
                     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_Intersect -> ()
            done
      )
      path_list;;

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

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

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

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 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 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
               (
                  let angle =
                     ((float_of_int0 index) *. 72.0 +. angle) *. pi /. 180.0 in
                  let ca, sa = cos angle, sin angle in
                  list_map
                     (fun (x, y) -> ca *. x -. sa *. y, sa *. x +. ca *. y)
                     scaled_path
               )
               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, fill =
      make_petals_fixed outline_path,
      make_petals_fixed fill_path in
   let background = list_map L.hd fill in
   rasterize_list [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_list
      [(translate_path
         (make_circle_path (5.0 *. supersample_scale) 10)
         tx ty)]
      outline_color;
   rasterize_list
      [(translate_path
         (make_circle_path (3.5 *. supersample_scale) 10)
         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 =
   array_init
      row_count
      (fun index ->
         flower_row_positions
            (random_jitter())
            ((float_of_int0 index) *. flower_spacing));;

let small_positions =
   array_init
      row_count
      (fun index ->
         if index = 0 then
            average_positions
               (large_positions.(row_count - 1))
               (
                  A.map (fun (x, y) -> x, y +. float_of_int0 supersample_size)
                        (large_positions.(0))
               )
         else
            average_positions (large_positions.(index))
                              (large_positions.(index - 1)));;

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) <-
      (
         let scanline = (image.(input_y)) in
         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
      )
   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
)
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 / sample_squared)
                                               (g / sample_squared)
                                               (b / sample_squared))
               scanline
         else
            ()
   )
   image;;