plecto-green-waves.fs

module Plecto

let PHI = 1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374847 // https://oeis.org/A001622
let PI = System.Math.PI
let TAU = 2.*PI
let E = System.Math.E

let inline fma m a x = x*m + a
let inline mix x0 x1 x = fma (x1-x0) x0 x
let inline map x0 x1 y0 y1 =
  let m = (y1-y0)/(x1-x0)
  let a = y0 - x0*m
  fma m a
let inline norm x0 x1 =
  let m = 1./(x1-x0)
  let a = -x0*m
  fma m a

type V1 = float
type V2 = float * float

type F11 = V1 -> V1
type F22 = V2 -> V2

type CountingGrid(width:int, height:int) =
  let widthf = float width
  let heightf = float height
  let size = width*height
  let data = Array.create size 0u

  member this.Inc ((x, y):V2) : unit =
    let yi = y*heightf |> int
    let xi = x*widthf |> int
    if xi >= 0 && yi >= 0 && xi < width && yi < width then
      let i = yi*width + xi
      data.[i] <- data.[i] + 1u

  member this.SaveJpg (path:string, gamma:float) : unit =
    let ceiling = this.Ceiling() |> float
    use img = System.Drawing.Bitmap(width, height)
    let mutable i = 0
    for y = 0 to (height-1) do
      for x = 0 to (width-1) do
        let v = (((float data.[i])/ceiling)**gamma)*255. |> int
        let c = System.Drawing.Color.FromArgb(255-v, 255-v, 255-v)
        img.SetPixel(x, y, c)
        i <- i + 1
    img.Save(path, System.Drawing.Imaging.ImageFormat.Jpeg)

  member this.Ceiling () : uint32 = Array.max data

  member this.Sample (n:int) (f:F22) (rand:System.Random) =
    for i = 1 to n do
      (rand.NextDouble(), rand.NextDouble()) |> f |> this.Inc

  member this.ParallelSample (threads:int) (samples:int) (f:F22) =
    let rand = System.Random(0)
    [| for i in 1 .. threads -> this |]
    |> Array.map (fun x -> async { x.Sample (samples/threads) f (System.Random(rand.Next())) })
    |> Async.Parallel
    |> Async.RunSynchronously
    |> ignore

  member this.MergeOnto (dest:uint32[]) : unit =
    Array.iteri2 (fun i a b -> dest.[i] <- a + b) data dest

// http://www.flong.com/texts/code/shapers_bez/
let bzq (a:float, b:float) (x:float) : float =
  let a = if a = 0.5 then a + 1e-6 else a
  let om2a = 1. - 2.*a
  let t = ((sqrt (a*a + om2a*x)) - a)/om2a
  let y = (1.-2.*b)*(t*t) + (2.*b)*t
  y

let bzq' (a:float, b:float) : F11 =
  let a = if a = 0.5 then a + 1e-6 else a

  let m = fma -2. 1. a
  let mi = 1./m
  let linear1 = fma m (a*a)
  let linear2 = fma mi (-a*mi)

  let c1 = 2.*b
  let c2 = fma -2. 1. b
  let inline polynomial x = x*x*c2 + x*c1

  linear1 >> sqrt >> linear2 >> polynomial

let generate (width:int) (height:int) (samplesPerPixel:int) (fn:F22) =
  let samples = samplesPerPixel * width * height
  let grid = CountingGrid(width, height)

  let sw = System.Diagnostics.Stopwatch.StartNew()
  grid.ParallelSample System.Environment.ProcessorCount samples fn
  sw.Stop()
  printfn "%Ams" (int sw.ElapsedMilliseconds)

  let sw = System.Diagnostics.Stopwatch.StartNew()
  grid.SaveJpg("out.jpg", 1./E)
  sw.Stop()
  printfn "%Ams" (int sw.ElapsedMilliseconds)

[<EntryPoint>]
let main argv =
  let bz = bzq' (0.2, 0.9)

  let fn(x, y) =
    let x = bz x

    let ax = 0.0
    let ay = (x*x) + mix 0. 1. x
    let bx = 1.0
    let by = (sqrt y) + 0.02*(sin (y*PHI*4.*TAU)) + 0.05*(sin (x*PHI*8.*TAU)) + mix 0. 1. x

    let x', y' = fma (ax-bx) bx y, fma (ay-by) by y
    (x'+1.)%1., (y'+1.)%1.

  let width = 0x200
  let height = width
  let samplesPerPixel = 0x80
  generate width height samplesPerPixel fn
  0