(* A comand line program to generate png images of fractals

By Daniel Thomas an extension from ML Tick6* Foundations Of Computer Science at Cambridge University

Copyright (C) 2008  Daniel Thomas

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You can get a copy of the GNU General Public License from 
http://www.gnu.org/licenses/old-licenses/gpl-2.0.html or, write to the 
Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*)


(*to compile: mosmlc -standalone ./fractals.sml
(do this in the right directory)
(on the pwfs the code compiled by mosmlc using normal methods expects camlrunm to be somewhere it isn't and so the code breaks - thus we specify standalone and get a much larger file - but it works - also it should be more portable)*)


exception Fail of string (* something has gone wrong in a terminal fashion *)
(* PROBLEM 4. Write a function mapi whidh applies a function to every pixel in a Gdimage.image *)

(*f is the colouring function, img is the image to which the function is being applied *)

(* a function should only return one thing but with functions that return () what can you do with them... therefore this little function *)

fun u((),()) = ();
infix u;

fun mapi f img = (* f is the colouring function int*int->int*int*int *)
    let
        val s = Gdimage.size img (*s is the size of the image *)
        fun dP (f, img, (x,y))= Gdimage.drawPixel img (Gdimage.Color((Gdimage.color img (f(x,y))))) (x,y) (* dP is a shorthand for this long call *) (*Color *)
        fun ix (f, img, (0,y)) = dP(f, img, (0,y)) (* ix iterates through the x values *)
          | ix (f, img, (x,y)) = dP(f, img, (x,y)) u ix(f, img, (x - 1, y))
        fun iy (f, img, (x,0)) = ix(f, img, (x,0))
          | iy (f, img, (x,y)) = ix(f, img, (x, y)) u iy(f, img, (x, y - 1)) (* iy iterates through the y values calling ix for each y value *)
    in
        iy(f, img, s)
    end;

(*mapi returns unit because it modifies an existing image rather than creating a new one*)


(* PROBLEM 6. Add the following function to your source code.*)

fun mandelbrot maxIter (x,y) =
    let
        fun solve (a,b) c = 
            if c = maxIter then 1.0
            else
                if (a*a + b*b <= 4.0) then
                    solve (a*a - b*b + x, 2.0*a*b + y) (c+1)
                else (real c)/(real maxIter)
    in
        solve (x, y) 0
    end;


(* PROBLEM 7. Add the following function to your source code. This function selects a RGB colour based on the number of iterations returned by the mandelbrot function.*)

fun chooseColour n =
    let
        val r = round ((Math.cos n) * 256.0)
        val g = round ((Math.cos n) * 256.0)
        val b = round ((Math.sin n) * 256.0)
    in
        (r, g, b)
    end;


(* PROBLEM 8. Write a function to convert pixel values to positions on the real number plane as this is what the mandelbrot function requires. rescale (w, h) (cx, cy, s) (x, y) should return the tuple (p, q) which is the real number position on the number plane which corresponds to the pixel value (x, y). The tuple (w, h) gives the size of the image in pixels, and the tuple (cx, cy, s) specifies the central point of the image and the size of the windows of interest.*)

fun rescale (w, h) (cx, cy, s) (x, y) = ((real(x)/w) * s - s/2.0 + cx, (real(y)/h) * s - s/2.0 + cy);


(* PROBLEM 9. Write a function mcompute which combines your rescale function with mandelbrot and chooseColour in order to compute the Mandelbrot set of a particular region and save the result to a file on disk (mandelbrot.png). Your function mcompute should have the following type real*real*real -> unit and have structure mcompute (cx, cy, s)*)

fun mcompute (cx, cy, s, file, maxIter, size) = 
    let
        val (sx,sy) = size (* for when we need to the x and y components separately*)
        val img = Gdimage.image size (255,255,255)
        val r = rescale (real(sx), real(sy)) (cx, cy, s) (* pass rescale its first two arguments so that we only have to pass it the (x,y) later*)
        val m = mandelbrot maxIter  (*pass mandelbrot the maximum number of iterations now so that we don't have to later *)
        val nothing = mapi (fn x => chooseColour (m (r x))) img
    in
        Gdimage.toPng img file
    end;

(* test with:   mcompute (~0.74364990, 0.13188204, 0.00073801); *)

(* Koch curve generating function *)
fun koch(dL, 0, x1, y1, x2, y2) = dL((round(x1), round(y1)), (round(x2), round(y2))) (* draw the line between the two points *)
  | koch(dL, i, x1, y1, x2, y2) = 
    let
        val x3 = x1 + (x2 - x1)/3.0       (*x coordinate of point 1/3 of the way along the line *)
        val y3 = y1 + (y2 - y1)/3.0
        val x4 = x1 + (x2 - x1)*(2.0/3.0)
        val y4 = y1 + (y2 - y1)*(2.0/3.0)
        val alpha = Math.atan2((y2-y1),(x2-x1)) (* angle of original line to horizontal *)
        val l = Math.sqrt((x4 - x3) * (x4 - x3) + (y4 - y3) * (y4 - y3))    (* length of the central line segment *)
        val theta = Math.pi/3.0 + alpha       (* total angle of new line from horizontal *)
        val x5 = x3 + l * Math.cos(theta)   (* x coordinate of new point at peak of triangle *)
        val y5 = y3 + l * Math.sin(theta)
    in
        koch(dL, i-1, x1, y1, x3, y3) u koch(dL, i-1, x3, y3, x5, y5) u koch(dL, i-1, x5, y5, x4, y4) u koch(dL, i-1, x4, y4, x2, y2)
    end;

(* takes arguments for the koch curve and processes them for passing to koch *)
fun kcompute file i x1 y1 x2 y2 sx sy =
    let
        val img = Gdimage.image (sx,sy) (255,255,255)
        val dL = Gdimage.drawLine img (Gdimage.Color((Gdimage.color img (0,0,0))))
        val nothing = koch(dL, i, x1, y1, x2, y2)
    in
        Gdimage.toPng img file
    end;

(* takes arguments for the koch snowflake and processes them for passing to koch *)
fun kscompute file i x1 y1 x2 y2 x3 y3 sx sy = 
    let
        val img = Gdimage.image (sx,sy) (255,255,255)
        val dL = Gdimage.drawLine img (Gdimage.Color((Gdimage.color img (0,0,0))))
        val nothing = koch(dL, i, x1, y1, x2, y2) u koch(dL, i, x2, y2, x3, y3) u koch(dL, i, x3, y3, x1, y1)
    in
        Gdimage.toPng img file
    end;


(* Calculates the sierpinski triangle *)
fun sTriangle(dT, 0, x1, y1, x2, y2, x3, y3) = dT #[(round(x1), round(y1)), (round(x2), round(y2)), (round(x3), round(y3))]
  | sTriangle(dT, i, x1, y1, x2, y2, x3, y3) = 
    let
        val x4 = (x1 + x2)/2.0
        val y4 = (y1 + y2)/2.0
        val x5 = (x1 + x3)/2.0
        val y5 = (y1 + y3)/2.0
        val x6 = (x2 + x3)/2.0
        val y6 = (y2 + y3)/2.0
    in
        sTriangle(dT, i-1, x1, y1, x4, y4, x5, y5) u sTriangle(dT, i-1, x4, y4, x2, y2, x6, y6) u sTriangle(dT, i-1, x5, y5, x6, y6, x3, y3)
    end;

fun stcompute file i x1 y1 x2 y2 x3 y3 sx sy = 
    let
        val img = Gdimage.image (sx,sy) (255,255,255)
        val dT = Gdimage.fillPolygon img (Gdimage.Color((Gdimage.color img (0,0,0))))
        val nothing = sTriangle(dT, i, x1, y1, x2, y2, x3, y3)
    in
        Gdimage.toPng img file
    end;


(* Calculates the sierpinski carpet *)
fun sCarpet(dS, 0, x1, y1, x2, y2, x3, y3, x4, y4) = dS #[(round(x1), round(y1)), (round(x2), round(y2)), (round(x3), round(y3)), (round(x4), round(y4))]
  | sCarpet(dS, i, x1, y1, x2, y2, x3, y3, x4, y4) = 
    let
        fun newcalc (xa,ya) (xb,yb) = ( ((xa + (xb-xa)/3.0), (ya + (yb-ya)/3.0)), ((xa + (xb-xa)*(2.0/3.0)), (ya + (yb-ya)*(2.0/3.0) )) )
        val ((x5,y5),(x6,y6)) = newcalc (x1,y1) (x2,y2)
        val ((x7,y7),(x8,y8)) = newcalc (x2,y2) (x3,y3)
        val ((x9,y9),(x10,y10)) = newcalc (x3,y3) (x4,y4)
        val ((x11,y11),(x12,y12)) = newcalc (x4,y4) (x1,y1)
        val ((x13,y13),(x14,y14)) = newcalc (x12,y12) (x7,y7)
        val ((x15,y15),(x16,y16)) = newcalc (x8,y8) (x11,y11)
    in
        sCarpet(dS, i-1, x1, y1, x5, y5, x13, y13, x12, y12) u sCarpet(dS, i-1, x5, y5, x6, y6, x14, y14, x13, y13) u sCarpet(dS, i-1, x6, y6, x2, y2, x7, y7, x14, y14) u sCarpet(dS, i-1, x14, y14, x7, y7, x8, y8, x15, y15) u sCarpet(dS, i-1, x15, y15, x8, y8, x3, y3, x9, y9) u sCarpet(dS, i-1, x16, y16, x15, y15, x9, y9, x10, y10) u sCarpet(dS, i-1, x11, y11, x16, y16, x10, y10, x4, y4) u sCarpet(dS, i-1, x12, y12, x13, y13, x16, y16, x11, y11)
    end;

fun sccompute file i x1 y1 x2 y2 x3 y3 x4 y4 sx sy = 
    let
        val img = Gdimage.image (sx,sy) (255,255,255)
        val dS = Gdimage.fillPolygon img (Gdimage.Color((Gdimage.color img (0,0,0))))
        val nothing = sCarpet(dS, i, x1, y1, x2, y2, x3, y3, x4, y4)
    in
        Gdimage.toPng img file
    end;


(* Calculates a brownian curve *)
fun brownian 0 dL x y r gen = ()
  | brownian i dL x y r gen = 
    let
        val x1 = x + 1.0
        val y1 = y + (Random.random gen - 0.5)*r
    in
        (dL ((round(x),round(y)),(round(x1),round(y1)))) u (brownian (i-1) dL x1 y1 r gen)
    end;

fun bcompute file i x y r sx sy= 
    let
        val gen = Random.newgen ()
        val img = Gdimage.image (sx,sy) (255,255,255)
        val dL = Gdimage.drawLine img (Gdimage.Color((Gdimage.color img (0,0,0))))
        val nothing = brownian i dL x y r gen
    in
        Gdimage.toPng img file
    end;

(* Calculates a random walk *)
fun randomwalk 0 dL x y r gen = ()
  | randomwalk i dL x y r gen = 
    let
        val x1 = x + (Random.random gen - 0.5)*r
        val y1 = y + (Random.random gen - 0.5)*r
    in
        (dL ((round(x),round(y)),(round(x1),round(y1)))) u (randomwalk (i-1) dL x1 y1 r gen)
    end;

fun rwcompute file i x y r sx sy= 
    let
        val gen = Random.newgen ()
        val img = Gdimage.image (sx,sy) (255,255,255)
        val dL = Gdimage.drawLine img (Gdimage.Color((Gdimage.color img (0,0,0))))
        val nothing = randomwalk i dL x y r gen
    in
        Gdimage.toPng img file
    end;


(* Argument collection and processing *)

exception Arguments of string; (*invalid arguments*)

fun main () = 
    let
        val args = CommandLine.arguments()
        val choice = List.nth(args, 0)
    in
        if choice = "mandelbrot" orelse choice = "m" then 
            let
                val cx = Option.valOf(Real.fromString(List.nth(args, 1)))
                val cy = Option.valOf(Real.fromString(List.nth(args, 2)))
                val s = Option.valOf(Real.fromString(List.nth(args, 3)))
                val file = List.nth(args, 4)
                val amaxIter = 150 + floor(1.0/(Math.sqrt(s))) (* automatic maximum number of iterations calculation *)
                val maxIter = Option.valOf(Int.fromString(List.nth(args, 5))) handle Subscript => if amaxIter < 8192 then amaxIter else 8192 (* The deeper you go the more iterations you need to see clearly - not sure that this is the correct relation - but it might help*)
                val sx = Option.valOf(Int.fromString(List.nth(args, 6))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 7))) handle Subscript => sx (*square looks better*)
                val size = (sx,sy)
            in
                mcompute (cx, cy, s, file, maxIter, size)
            end
        else if choice = "koch" orelse choice = "k" then 
            let
                val file = List.nth(args, 1)
                val i = Option.valOf(Int.fromString(List.nth(args, 2)))
                val x1 = Option.valOf(Real.fromString(List.nth(args, 3))) handle Subscript => 10.0
                val y1 = Option.valOf(Real.fromString(List.nth(args, 4))) handle Subscript => 320.0
                val x2 = Option.valOf(Real.fromString(List.nth(args, 5))) handle Subscript => 630.0
                val y2 = Option.valOf(Real.fromString(List.nth(args, 6))) handle Subscript => 320.0
                val sx = Option.valOf(Int.fromString(List.nth(args, 7))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 8))) handle Subscript => sx (*square looks better*)
            in
                kcompute file i x1 y1 x2 y2 sx sy
            end
        else if choice = "kochsnowflake" orelse choice = "ks" then
            let
                val file = List.nth(args, 1)
                val i = Option.valOf(Int.fromString(List.nth(args, 2)))
                val x1 = Option.valOf(Real.fromString(List.nth(args, 3))) handle Subscript => 320.0
                val y1 = Option.valOf(Real.fromString(List.nth(args, 4))) handle Subscript => 40.0
                val x2 = Option.valOf(Real.fromString(List.nth(args, 5))) handle Subscript => 80.0
                val y2 = Option.valOf(Real.fromString(List.nth(args, 6))) handle Subscript => 501.0
                val x3 = Option.valOf(Real.fromString(List.nth(args, 7))) handle Subscript => 560.0
                val y3 = Option.valOf(Real.fromString(List.nth(args, 8))) handle Subscript => 501.0
                val sx = Option.valOf(Int.fromString(List.nth(args, 9))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 10))) handle Subscript => sx (*square looks better*)
            in
                kscompute file i x1 y1 x2 y2 x3 y3 sx sy
            end
        else if choice = "sierpinski" orelse choice = "s" orelse choice = "sierpiński" then 
            let
                val file = List.nth(args, 1)
                val i = Option.valOf(Int.fromString(List.nth(args, 2)))
                val x1 = Option.valOf(Real.fromString(List.nth(args, 3))) handle Subscript => 320.0
                val y1 = Option.valOf(Real.fromString(List.nth(args, 4))) handle Subscript => 40.0
                val x2 = Option.valOf(Real.fromString(List.nth(args, 5))) handle Subscript => 80.0
                val y2 = Option.valOf(Real.fromString(List.nth(args, 6))) handle Subscript => 501.0
                val x3 = Option.valOf(Real.fromString(List.nth(args, 7))) handle Subscript => 560.0
                val y3 = Option.valOf(Real.fromString(List.nth(args, 8))) handle Subscript => 501.0
                val sx = Option.valOf(Int.fromString(List.nth(args, 9))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 10))) handle Subscript => sx (*square looks better*)
            in
                stcompute file i x1 y1 x2 y2 x3 y3 sx sy
            end
        else if choice = "carpet" orelse choice = "c" then
            let
                val file = List.nth(args, 1)
                val i = Option.valOf(Int.fromString(List.nth(args, 2)))
                val x1 = Option.valOf(Real.fromString(List.nth(args, 3))) handle Subscript => 20.0
                val y1 = Option.valOf(Real.fromString(List.nth(args, 4))) handle Subscript => 20.0
                val x2 = Option.valOf(Real.fromString(List.nth(args, 5))) handle Subscript => 620.0
                val y2 = Option.valOf(Real.fromString(List.nth(args, 6))) handle Subscript => 20.0
                val x3 = Option.valOf(Real.fromString(List.nth(args, 7))) handle Subscript => 620.0
                val y3 = Option.valOf(Real.fromString(List.nth(args, 8))) handle Subscript => 620.0
                val x4 = Option.valOf(Real.fromString(List.nth(args, 9))) handle Subscript => 20.0
                val y4 = Option.valOf(Real.fromString(List.nth(args, 10))) handle Subscript => 620.0
                val sx = Option.valOf(Int.fromString(List.nth(args, 11))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 12))) handle Subscript => sx (*square looks better*)
            in
                sccompute file i x1 y1 x2 y2 x3 y3 x4 y4 sx sy
            end
        else if choice = "brownian" orelse choice = "b" then 
            let
                val file = List.nth(args, 1)
                val i = Option.valOf(Int.fromString(List.nth(args, 2)))
                val x = Option.valOf(Real.fromString(List.nth(args, 3))) handle Subscript => 20.0
                val y = Option.valOf(Real.fromString(List.nth(args, 4))) handle Subscript => 320.0
                val r = Option.valOf(Real.fromString(List.nth(args, 5))) handle Subscript => 5.0
                val sx = Option.valOf(Int.fromString(List.nth(args, 6))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 7))) handle Subscript => sx (*square looks better*)
            in
                bcompute file i x y r sx sy
            end
        else if choice = "randomwalk" orelse choice = "rw" then 
            let
                val file = List.nth(args, 1)
                val i = Option.valOf(Int.fromString(List.nth(args, 2)))
                val x = Option.valOf(Real.fromString(List.nth(args, 3))) handle Subscript => 320.0
                val y = Option.valOf(Real.fromString(List.nth(args, 4))) handle Subscript => 320.0
                val r = Option.valOf(Real.fromString(List.nth(args, 5))) handle Subscript => 5.0
                val sx = Option.valOf(Int.fromString(List.nth(args, 6))) handle Subscript => 640
                val sy = Option.valOf(Int.fromString(List.nth(args, 7))) handle Subscript => sx (*square looks better*)
            in
                rwcompute file i x y r sx sy
            end
        else raise Arguments "Invalid choice of fractal"
    end;


(*              _DOCUMENTATION_
_Orbit Fractals_
mandelbrot/m
    cx:real cy:real s:real file:string [maxIter:int [sx:int [sy:int ]]]
    cx is central x coordinate
    cy is central y coordinate
    s is zoom/size smaller is closer
    file is the file name for the image to be out put to - in the current directory
    [maxIter is the maximum number of iterations in the mandebrot set = 150 + floor(1.0/(s*10.0))
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]
   test with: fractals m -0.74364990 0.13188204 0.00007801 fimage.png


_Iterated Function Fractals_
koch/k
    file:string i:int [x1:real [y1:real [x2:real [y2:real [sx:int [sy:int]]]]]]
    file is the file name for the image to be out put to - in the current directory
    i is iterations
    [x1 is the initial x coordinate
    [y1 is the initial y coordinate
    [x2 is the final x coordinate
    [y2 is the final y coordinate
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]]]]

kochsnowflake/ks
    file:string i:int [x1:real [y1:real [x2:real [y2:real [x3:real [y3:real [sx:int [sy:int]]]]]]]]
    file is the file name for the image to be out put to - in the current directory
    i is iterations
    [x1 is the initial x coordinate
    [y1 is the initial y coordinate
    [x2 is the second x coordinate
    [y2 is the second y coordinate
    [x3 is the final x coordinate
    [y3 is the final y coordinate
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]]]]]]


sierpinski/s
    file:string i:int [x1:real [y1:real [x2:real [y2:real [x3:real [y3:real [sx:int [sy:int]]]]]]]]
    file is the file name for the image to be out put to - in the current directory
    i is iterations
    [x1 is the initial x coordinate
    [y1 is the initial y coordinate
    [x2 is the second x coordinate
    [y2 is the second y coordinate
    [x3 is the final x coordinate
    [y3 is the final y coordinate
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]]]]]]

carpet/c (sierpinski)
    file:string i:int [x1:real [y1:real [x2:real [y2:real [x3:real [y3:real [x4:real [y4:real [sx:int [sy:int]]]]]]]]]]
    file is the file name for the image to be out put to - in the current directory
    i is iterations
    [x1 is the initial x coordinate
    [y1 is the initial y coordinate
    [x2 is the second x coordinate
    [y2 is the second y coordinate
    [x3 is the third x coordinate
    [y3 is the third y coordinate
    [x4 is the fourth x coordinate
    [y4 is the fourth y coordinate
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]]]]]]]]

_Random Fractals_
brownian/b
    file:string i:int [x:real [y:real [r:real [sx:int [sy:int]]]]]
    i is iterations
    [x is the initial x coordinate
    [y is the initial y coordinate
    [r is the range of the random movement
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]]]

randomwalk/rw
    file:string i:int [x:real [y:real [r:real [sx:int [sy:int]]]]]
    i is iterations
    [x is the initial x coordinate
    [y is the initial y coordinate
    [r is the range of the random movement
    [sx is the width of the png file produced in px = 640
    [sy is the height of the png file produced in px = sx]]]]]



                _TESTING_
fractals mandelbrot -0.74364990 0.13188204 0.00007801 fmandelbrot.png
fractals koch fkoch.png 5
fractals kochsnowflake fkochsnowflake.png 5
fractals sierpinski fsierpinski.png 7
fractals carpet fcarpet.png 5
fractals brownian fbrownian.png 600
fractals randomwalk frandomwalk.png 20000

example output for these commands can be found at: http://www.srcf.ucam.org/~drt24/fractals/1.0/example/
NOTE: Random Fractals will probably look different every time.

*)

val _ = main();

(* Initiate program *)