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

   12/26/09
                                                                           *)

(* {{{ Data bits *)

(* Size of output bitmap *)
let output_width = 300;;
let output_height = 300;;

(* Sample this many pixels per output pixel *)
let supersample = 4;;
let supersample_scale = float_of_int supersample;;

(* Clockwise control points for flower petals.  All these are relative
   coordinates that describe a bezier curve starting at the origin.
   The paths are specified this way to make it easier to convert from
   the relative paths in SVG files.                                    *)
let outline1 =
[
   (24.565, -42.9843);
   (-50.1807, -43.8067);
   (-24.561, -3.32677)
];;
let fill1 =
[
   (22.6673, -41.9008);
   (-46.3969, -42.1701);
   (-22.3531, -2.10276)
];;

let outline2 =
[
   (16.3197, -26.85);
   (-8.07283, -23.4512);
   (-5.08858, -1.75512)
];;
let fill2 =
[
   (15.2366, -26.7315);
   (-6.92953, -22.7587);
   (-4.38189, -1.93425)
];;

(* Starting point offsets *)
let outline1_offset =
(
   23.3913 +. 24.561 -. 32.7075, 28.4083 +. 3.32677 -. 40.2591
);;
let fill1_offset =
(
   24.7717 +. 22.3531 -. 32.7075, 31.5165 +. 2.10276 -. 40.2591
);;
let outline2_offset =
(
   357.426 +. 5.08858 -. 359.222, 132.156 +. 1.75512 -. 134.66
);;
let fill2_offset =
(
   357.299 +. 4.38189 -. 359.222, 134.041 +. 1.93425 -. 134.66
);;

(* Inner circle radii for flower 1 *)
let circle1_outline_radius = 5.0 *. supersample_scale;;
let circle1_fill_radius = 3.5 *. supersample_scale;;

(* Spacing between flowers.  output_width and output_height must be some
   multiple of the unscaled number, otherwise there would be unnatural
   gaps between flowers in output.                                       *)
let flower_spacing = 75.0 *. supersample_scale;;

(* Internal dimensions *)
let supersample_width = output_width * supersample;;
let supersample_height = output_height * supersample;;

(* Colors *)
let outline_color = (0, 0, 0);;
let background_color = (124, 194, 169);;
let fill1_color = (233, 185, 185);;
let circle1_color = (221, 218, 167);;
let fill2_color = (234, 236, 237);;

(* Alternative set of colors.  These are closer to what's on the
   anime's official site, but they look too saturated on my monitor,
   especially when the output is tiled to fit the whole screen, so I
   prefer the colors above.
let background_color = (89, 160, 128);;
let fill1_color = (239, 166, 150);;
let circle1_color = (219, 203, 118);;
let fill2_color = (244, 244, 243);;
*)

(* Make a bezier path from starting offset and relative coordinates *)
let make_absolute_path start path =
   List.map
      (fun (x, y) -> x *. supersample_scale, y *. supersample_scale)
      (start :: List.map (fun (x, y) -> x +. fst start, y +. snd start) path);;

(* Generate clockwise circular path *)
let circle_path_steps = 10;;

let make_circle_path r =
   let step_to_radian =
      2.0 *. (atan2 0. (-1.0)) /. (float_of_int circle_path_steps) in
   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 circle_path_steps;;

(* }}} *)

(* {{{ Path management bits *)

(* Conversion factor from degrees to radian *)
let to_radian = atan2 0.0 (-1.0) /. 180.0;;

(* Rotate a path about the origin *)
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;;

(* Translate path *)
let translate_path path tx ty =
   List.map (fun (x, y) -> x +. tx, y +. ty) path;;

(* Interpolate bezier segment for one dimension *)
let interpolate_steps = 64;;
let interpolate a b c d =
   let rec interpolate_recurse i result =
      if i = interpolate_steps then
         List.rev result
      else
      (
         let t = (float_of_int i) /. (float_of_int interpolate_steps) in
         let ab = a +. t *. (b -. a) in
         let bc = b +. t *. (c -. b) in
         let cd = c +. t *. (d -. c) in
         let abc = ab +. t *. (bc -. ab) in
         let bcd = bc +. t *. (cd -. bc) in
         interpolate_recurse (i + 1) ((abc +. t *. (bcd -. abc)) :: result)
      ) in
   interpolate_recurse 0 [d];;

(* Interpolate a list of 4 elements *)
let interpolate_list x =
   match x with [a; b; c; d] -> interpolate a b c d
   | _ -> [];;

(* Interpolate a bezier path *)
let interpolate_path path =
   let xlist, ylist = List.split path in
   List.combine (interpolate_list xlist) (interpolate_list ylist);;

(* }}} *)

(* {{{ Rasterization bits *)

(* Given two vectors ab and ac, return true if ac lies clockwise from ab *)
let clockwise a b c =
   let abx = (fst b) -. (fst a) in
   let aby = (snd b) -. (snd a) in
   let acx = (fst c) -. (fst a) in
   let acy = (snd c) -. (snd a) in
   acx *. aby >= abx *. acy;;

(* Check if a point is inside a convex polygon with clockwise contour *)
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 [] ->
               (* Final point *)
               clockwise next_point first_point point
            | ia::ib ->
            (
               (* Intermediate points *)
               if not (clockwise next_point ia point) then
                  false
               else
                  inside_recurse ia ib
            ) in
         inside_recurse a b
   )
   | _ -> false;;

(* Check range and see if we should break recursion.
   This assumes that one unit equals one pixel in output. *)
let stop_recurse min_x max_x = max_x -. min_x < 0.5;;

(* Average two numbers *)
let avg a b = (a +. b) /. 2.0;;

(* Given a scanline and a convex polygon, find a point that is inside
   the polygon.  Raises Not_found if no such point exists.  Not very
   efficient when the polygon is leaning towards the right side of the
   bounding box.                                                       *)
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;;

(* Get minimum point where a scanline intersects a polygon.
   max_x must already be inside the polygon.                *)
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;;

(* Get maximum point where a scanline intersects a polygon.
   min_x must already be inside the polygon.                *)
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;;

(* Given a scanline, get the two ends where it intersects a polygon.
   Raises Not_found if scanline does not intersect polygon.          *)
let intersect_scanline path y min_x max_x =
   let start_x = find_point_on_scanline path y min_x max_x in
   (
      intersect_scanline_left path y min_x start_x,
      intersect_scanline_right path y start_x max_x
   );;

(* Get range of values from list elements.  Doesn't work if list is empty. *)
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;;

(* Wrap around a number such that it's within positive range *)
let wraparound value range_max =
   if value < 0 then
      range_max + (value mod range_max)
   else
      value mod range_max;;

(* Rasterize a single scanline *)
let rasterize_scanline path color output width height y min_x max_x =
   try
      (* Get scanline range *)
      let x0, x1 = intersect_scanline path y min_x max_x in

      (* Wraparound X *)
      let tx0 = wraparound (truncate x0) width in
      let tx1 = wraparound (truncate x1) width in

      (* Wraparound Y *)
      let sy = wraparound (truncate y) height in

      (* Output scanline *)
      if tx0 <= tx1 then
         Array.fill (output.(sy)) tx0 (tx1 - tx0) color
      else
      (
         Array.fill (output.(sy)) 0 (tx1 + 1) color;
         Array.fill (output.(sy)) tx0 (width - tx0) color
      )
   with Not_found -> ();;

(* Rasterize a convex polygon *)
let rasterize path color output width height =
   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 output width height
                         (float_of_int y) min_x max_x
   done;;

(* Rasterize a list of convex polygons *)
let rasterize_list paths color output width height =
   List.iter (fun x -> rasterize x color output width height) paths;;

(* }}} *)

(* {{{ Image manipulation bits *)

(* Add two pixels together *)
let add_pixel p1 p2 =
   let r1, g1, b1 = p1 in
   let r2, g2, b2 = p2 in
   (r1 + r2, g1 + g2, b1 + b2);;

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

(* Downsample a single scanline *)
let downsample_scanline scanline sample_size final_width =
   let output = Array.make final_width (scanline.(0)) in
   for x = 0 to (final_width - 1) do
      (* Function to add consecutive samples *)
      let rec add_sample index =
      (
         if index = 1 then
            scanline.(x * sample_size)
         else
            add_pixel (scanline.(x * sample_size + index - 1))
                      (add_sample (index - 1))
      ) in

      (* Add consecutive samples for each component *)
      output.(x) <- add_sample sample_size
   done;
   output;;

(* Downsample a set of scanlines in-place *)
let downsample_scanline_set image sample_size final_width output_y =
   (* Downsample scanlines in set *)
   for i = 0 to (sample_size - 1) do
      let input_y = output_y * sample_size + i in
      image.(input_y) <-
         downsample_scanline (image.(input_y)) sample_size final_width
   done;

   (* Combine scanlines *)
   image.(output_y) <- image.(output_y * sample_size);
   let first_scanline = image.(output_y) in
   for i = 1 to (sample_size - 1) do
      Array.iteri
         (fun x pixel ->
            first_scanline.(x) <- add_pixel (first_scanline.(x)) pixel)
         (image.(output_y * sample_size + i))
   done;

   (* Rescale component values *)
   let scale = scale_components (sample_size * sample_size) in
   Array.iteri
      (fun x pixel -> first_scanline.(x) <- scale pixel)
      first_scanline;;

(* Downsample image in-place *)
let downsample image sample_size =
   let final_height = (Array.length image) / sample_size in
   let final_width = (Array.length (image.(0))) / sample_size in
   for y = 0 to (final_height - 1) do
      downsample_scanline_set image sample_size final_width y
   done;;

(* Print PPM image to stdout *)
let output_ppm image width height =
   Printf.printf "P3\n%d %d\n255\n" width height;
   Array.iteri
      (
         fun y scanline ->
            (* Output scanlines up to height pixels, ignoring the
               remaining scanlines.                               *)
            if y < height then
               (* Number of pixels in the scanlines that will be printed
                  must match expected width.  This is guaranteed by
                  downsample_scanline.                                   *)
               Array.iter
                  (fun (r, g, b) -> Printf.printf "%d %d %d\n" r g b) scanline
            else
               ()
      )
      image;;

(* }}} *)

(* {{{ Position bits *)

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

(* Generate positions for a row of flowers *)
let flower_row_positions x0 y =
   let flower_count =
      truncate ((float_of_int supersample_width) /. flower_spacing) in
   Array.init
      flower_count
      (fun index ->
         (x0 +. (float_of_int index) *. flower_spacing +. random_jitter(),
          y +. random_jitter()));;

(* Generate positions for a grid of flowers *)
let flower_grid_positions () =
   let row_count =
      truncate ((float_of_int supersample_height) /. flower_spacing) in
   Array.init
      row_count
      (fun index ->
         flower_row_positions
            (random_jitter())
            ((float_of_int index) *. flower_spacing));;

(* Generate a single position that fills the gap between 4 points *)
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);;

(* Shift a point horizontally for wraparound *)
let horizontal_shift point =
   let x, y = point in (x +. float_of_int supersample_width, y);;

(* Shift a row of points vertically for wraparound *)
let vertical_shift row =
   Array.map
      (fun (x, y) -> x, y +. float_of_int supersample_height)
      row;;

(* Generate positions that fills the gap between two rows *)
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)]);;

(* Generate grid positions by filling in gaps *)
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)));;

(* }}} *)

(* {{{ Drawing bits *)

(* Make 5 rotated copies out of a single petal *)
let make_petals start path tx ty angle =
   let absolute_path = make_absolute_path start path in
   let rec rotate_and_duplicate_path index =
      if index = 0 then
         []
      else
      (
         interpolate_path
         (
            translate_path
               (rotate_path absolute_path
                            ((float_of_int index) *. 72.0 +. angle))
               tx ty
         )
      ) :: (rotate_and_duplicate_path (index - 1)) in
   rotate_and_duplicate_path 5;;

(* Rasterize flower petals *)
let draw_flower outline_offset outline_path
                fill_offset fill_path fill_color
                tx ty angle output =
   let outline = make_petals outline_offset outline_path tx ty angle in
   let fill = make_petals fill_offset fill_path tx ty angle in
   let background = List.map List.hd fill in
   rasterize background fill_color output supersample_width supersample_height;
   rasterize_list outline outline_color
                  output supersample_width supersample_height;
   rasterize_list fill fill_color output supersample_width supersample_height;;

(* Rasterize large flower *)
let large_flower tx ty angle output =
   draw_flower outline1_offset outline1
               fill1_offset fill1 fill1_color
               tx ty angle output;
   rasterize
      (translate_path (make_circle_path circle1_outline_radius) tx ty)
      outline_color output supersample_width supersample_height;
   rasterize
      (translate_path (make_circle_path circle1_fill_radius) tx ty)
      circle1_color output supersample_width supersample_height;;

(* Rasterize small flower *)
let small_flower tx ty angle output =
   draw_flower outline2_offset outline2
               fill2_offset fill2 fill2_color
               tx ty angle output;;

(* Draw flowers at all positions *)
let draw_flower_set flower_func positions output =
   Array.iter
      (fun row ->
         Array.iter
         (fun (x, y) -> flower_func x y (Random.float 360.0) output)
         row)
      positions;;

(* Draw large and small flowers *)
let draw_all_flowers output =
   let large_positions = flower_grid_positions() in
   let small_positions = average_grid_positions large_positions in
   draw_flower_set large_flower large_positions output;
   draw_flower_set small_flower small_positions output;;

(* }}} *)

(* {{{ Program entry *)

Random.self_init();;

let image =
   Array.make_matrix supersample_height supersample_width background_color;;

draw_all_flowers image;;
downsample image supersample;;
output_ppm image output_width output_height;;

(* }}} *)