Simple Fractals in Haskell.

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:

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

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].

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

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:

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

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!

{-# 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."

{-# 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).