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

   Use globals.

   12/31/09
                                                                           *)

(* {{{ Data bits *)

(* Size of square output bitmap *)
let output_size = 300;;

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

(* Clockwise control points for flower petals *)
let outline1 =
[
   (15.24, -8.52);
   (39.81, -51.51);
   (-34.94, -52.33);
   (-9.32, -11.85)
];;
let fill1 =
[
   (14.42, -6.64);
   (37.08, -48.54);
   (-31.98, -48.81);
   (-7.94, -8.74)
];;

let outline2 =
[
   (3.29, -0.75);
   (19.61, -27.6);
   (-4.78, -24.2);
   (-1.8, -2.5)
];;
let fill2 =
[
   (2.46, 1.32);
   (17.7, -25.42);
   (-4.47, -21.44);
   (-1.92, -0.62)
];;

(* 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_size = output_size * 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);;

(* 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;;

(* Output buffer *)
let image =
   Array.make_matrix supersample_size supersample_size background_color;;

(* }}} *)

(* {{{ 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 =
   if value < 0 then
      supersample_size + (value mod supersample_size)
   else
      value mod supersample_size;;

(* Rasterize a single scanline *)
let rasterize_scanline path color 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) in
      let tx1 = wraparound (truncate x1) in

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

      (* Output scanline *)
      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 -> ();;

(* Rasterize a convex polygon *)
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;;

(* Rasterize a list of convex polygons *)
let rasterize_list paths color =
   List.iter (fun x -> rasterize x color) 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 pixel =
   let r, g, b = pixel in
   (r / sample_squared, g / sample_squared, b / sample_squared);;

(* Downsample a single scanline *)
let downsample_scanline scanline =
   let output = Array.make output_size (scanline.(0)) in
   for x = 0 to (output_size - 1) do
      (* Function to add consecutive samples *)
      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

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

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

   (* Combine scanlines *)
   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;

   (* Rescale component values *)
   Array.iteri
      (fun x pixel -> first_scanline.(x) <- scale_components pixel)
      first_scanline;;

(* Downsample image in-place *)
let downsample () =
   for y = 0 to (output_size - 1) do
      downsample_scanline_set y
   done;;

(* Print PPM image to stdout *)
let output_ppm () =
   Printf.printf "P3\n%d %d\n255\n" output_size output_size;
   Array.iteri
      (
         fun y scanline ->
            (* Output scanlines up to height pixels, ignoring the
               remaining scanlines.                               *)
            if y < output_size 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_size) /. 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_size) /. 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_size, y);;

(* Shift a row of points vertically for wraparound *)
let vertical_shift row =
   Array.map (fun (x, y) -> x, y +. float_of_int supersample_size) 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 path tx ty angle =
   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
      (
         interpolate_path
         (
            translate_path
               (rotate_path scaled_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_path fill_path fill_color tx ty angle =
   let outline = make_petals outline_path tx ty angle in
   let fill = make_petals fill_path tx ty angle in
   let background = List.map List.hd fill in
   rasterize background fill_color;
   rasterize_list outline outline_color;
   rasterize_list fill fill_color;;

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

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

(* Draw flowers at all positions *)
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;;

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

(* }}} *)

(* {{{ Program entry *)

Random.self_init();;

draw_all_flowers();;
downsample();;
output_ppm();;

(* }}} *)