This posting was inspired by the ideas in ‘Composing Fractals’ by Mark P Jones.

Fractals are generated by applying a function to the result of applying that function to the result of applying that function…

Haskell has a very simple way of expressing this iterative function application:

1 |
iterate :: (a -> a) -> a -> [a] |

i.e. iterate takes a function, that takes an ‘a’ and produces an ‘a’, along with an initial ‘a’ – the seed value. The function is then repeatedly applied ‘forever’ giving an infinite list of [a].

1 2 |
-- here x is the seed iterate f x == [x, f x, f (f x), ...] |

Of course we can’t really use all of an infinite list so we apply Haskell’s ‘take’ function to obtain however many items from the list.

For example in ghci:

1 2 3 4 5 6 7 8 9 |
Prelude λ-> xs = iterate (+1) 0 Prelude λ-> take 5 xs [0,1,2,3,4] Prelude λ-> take 20 xs [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] Prelude |

Two of the most well known fractals are the Mandelbrot set and the Julia set and both are generated in similar but subtly different ways. With the Mandelbrot set the sequence is calculated with a seed of 0 and for the Julia set the seed is the current point from the (x, y) grid being considered. In either case the values in the sequence could become larger and larger but if they don’t then they are counted as ‘being in the set’.

These ideas about sequences converging, seeds etc are explained in a very accessible way at https://plus.maths.org/content/unveiling-mandelbrot-set

So let’s see some code!

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 |
{-# LANGUAGE OverloadedStrings #-} {-# OPTIONS -Wall -fwarn-tabs -fno-warn-type-defaults -fno-warn-unused-do-bind #-} -- The package juicypixels needs to have been installed. import Data.Complex import Codec.Picture type C = Complex Double type Pnt = (Double, Double) type Grid a = [[a]] w :: Int w = 500 h :: Int h = 500 -- Complex constants for the Julia set rendering. fixC1, fixC2, fixC3, fixC4 :: C fixC1 = (-1.037) :+ 0.17 fixC2 = (-0.52) :+ 0.57 fixC3 = 0.295 :+ 0.55 fixC4 = (-0.624) :+ 0.435 juliaConstants :: [(Int, C)] juliaConstants =[ (1, fixC1), (2, fixC2), (3, fixC3), (4, fixC4)] --iterate f x == [x, f x, f (f x), ...] -- in mandel the seed is fixed at 0 mandel :: C -> [C] mandel c = iterate (\z -> z * z + c) (0.0 :+ 0.0) --iterate f x == [x, f x, f (f x), ...] -- in julia, for a given fixed constant, the seed varies julia :: C -> C -> [C] julia fixc = iterate (\z -> z * z + fixc) closeToZero :: C -> Bool closeToZero z = magnitude z < 5 -- get the pallete value corresponding to the 'depth' at which the value fails the test chooseColor :: [color] -> ( C -> [C]) -> C -> color chooseColor colors fz z = (colors !!) . length . take n . takeWhile closeToZero $ fz z where n = length colors - 1 -- A grid of points that is defined by its bottom left, top right and the number of subdivision per row and per column grid :: Int -> Int -> Pnt -> Pnt -> Grid Pnt grid c r (xmin, ymin) (xmax, ymax ) = [[(x, y) | x <- for c xmin xmax ] | y <- for r ymin ymax ] where for :: Int ->Double ->Double ->[Double] for n mn mx = take n [mn, mn + delta ..] where delta = (mx - mn) / fromIntegral (n - 1) -- a grid grid1 :: Grid Pnt grid1 = grid w h (-1.5, -1.5) (1.5, 1.5) -- calculate sequences over the grid and for each sequence get an 'end' colour mapOverGrid :: Grid Pnt -> (C -> [C]) -> [color] -> [[color]] mapOverGrid grd fz colors = [ [chooseColor colors fz (x :+ y) | (x,y) <- row] | row <- grd ] -- used by juicy pixels to get a pixel value at x, y gridFunc ::Grid Pnt -> (C -> [C]) -> Int -> Int -> PixelRGB8 gridFunc grd fz x y = ((mapOverGrid grd fz rgb8) !! y) !! x generateImg :: (Int -> Int -> PixelRGB8) -> DynamicImage generateImg gf = ImageRGB8 (generateImage gf w h) -- Juicy pixels type rgb8 :: [PixelRGB8] rgb8 = foldr (\x ac -> PixelRGB8 x x x : ac) [] [255, 254..0] allJulia :: IO() allJulia = mapM_ makeSaveImage juliaConstants where makeSaveImage :: (Int, C) -> IO() makeSaveImage (n, c) = savePngImage name $ generateImg . gridFunc grid1 . julia $ c where name = "imageJ" ++ show n ++ ".png" main :: IO () main = do allJulia savePngImage "imageM.png" $ generateImg . gridFunc grid1 $ mandel print "Done." |

I think the code is fairly self explanatory and quite succinct given what it does.

The images below have been generated by the above code and they are the Mandelbrot and the four Julia images corresponding to ‘fixC1’ to ‘fixC4’. (The background colour has been manually changed to contrast with the page background).

Qué?

Pretty pictures 😉

Hi Kate, 🙂 Yes very pretty. With more effort in rendering you can get things like this!

https://commons.wikimedia.org/wiki/File:Mandel_zoom_08_satellite_antenna.jpg