# Fractals 3 ways¶

This notebook illustrates three ways to obtain objects exhibiting fractal geometry:

1. Approaching a pointset with a fractal boundary, using deterministic iteration.
2. Explicit recursive subdivision.
3. Iterated function systems, using the random iteration algorithm.

I will use three different graphics libraries to do this: PyPlot (calling Python's matplotlib), Compose.jl (which is great for nested geometry), and Cairo (which can draw thousands of primitives in no time).

## Julia set¶

In [2]:
using PyPlot

In [3]:
# julia set
# (the familiar mandelbrot set is obtained by setting c==z initially)
function julia(z, c; maxiter=200)
for n = 1:maxiter
if abs2(z) > 4
return n-1
end
z = z*z + c
end
return maxiter
end

Out[3]:
# methods for generic function julia
julia(z,c) at In[3]:4
In [5]:
# varying the second argument to julia() tiny amounts results in a stunning variety of forms
@time m = [ uint8(julia(complex(r,i), complex(-.06,.67))) for i=1:-.002:-1, r=-1.5:.002:1.5 ];

elapsed time: 0.1899382 seconds (1502744 bytes allocated)

In [6]:
# the notebook is able to display ColorMaps
get_cmap("RdGy")

Out[6]:
In [7]:
imshow(m, cmap="RdGy", extent=[-1.5,1.5,-1,1])

Out[7]:
PyObject <matplotlib.image.AxesImage object at 0xc0649d0>

## Sierpinski triangle¶

In [83]:
using Compose

In [85]:
# drawing an SVG shape with Compose
draw(SVG(1inch,(sqrt(3)/2)inch), compose(canvas(), polygon((1,1),(0,1),(1/2,0))))

In [15]:
# directly from dcjones' Compose README
function sierpinski(n)
if n == 0
compose(canvas(), polygon((1,1), (0,1), (1/2, 0)))
else
t = sierpinski(n - 1)
compose(canvas(),
(canvas(1/4,   0, 1/2, 1/2), t),
(canvas(  0, 1/2, 1/2, 1/2), t),
(canvas(1/2, 1/2, 1/2, 1/2), t))
end
end

Out[15]:
# methods for generic function sierpinski
sierpinski(n) at In[15]:2
In [29]:
# get more detail by scaling up
draw(SVG(9inch,7inch), compose(sierpinski(8), fill(nothing), linewidth(0.1mm)))


## Iterated function systems¶

In [72]:
# draw from a discrete CDF
pick(p) = searchsortedfirst(p, rand())

Out[72]:
# methods for generic function pick
pick(p) at In[72]:2
In [72]:
# the iterated function system (IFS) "code" for Barnsley's fern.
# consists of four matrices and a probability for each.
fern_M = {[ 0    0    0; 0    .16 0;    0 0 1],
[ 0.85  .04 0; -.04 .85 1.6;  0 0 1],
[ 0.2  -.26 0;  .23 .22 1.6;  0 0 1],
[-0.15  .28 0;  .26 .24  .44; 0 0 1]}
fern_P = cumsum([.01, .85, .07, .07]);

In [59]:
# first, an implementation using Compose.
# This is a huge point cloud, which Compose is definitely not designed for, and runs quite slowly.

point(x,y) = polygon((x,y),(x+.001,y+.001))

function ifs(M, p, niter)
pt = [0.5,0.5,1]     # start at an arbitrary point
c = compose(canvas(), linewidth(0.3mm))
for i = 1:niter+10
pt = M[pick(p)]*pt
if i > 10        # wait 10 iterations to make sure we approach the attractor
c = compose(c, point((pt[1]+4)/14,1-pt[2]/10.2))
end
end
c
end

Out[59]:
# methods for generic function ifs
ifs(GC,M,p,niter) at In[15]:2
ifs(M,p,niter) at In[59]:4
In [157]:
@time draw(PNG(4.25inch,2.5inch), ifs(fern_M, fern_P, 16000))

elapsed time: 25.186542174 seconds (7203514612 bytes allocated)

In [8]:
using Cairo
using Color
using Base.Graphics

In [132]:
# next an implementation using Cairo directly to draw small circles.
# this is drastically faster (at least 100x)
function ifs(GC, M, p, niter, scale=0.01)
pt = [0.5,0.5,1]
for i = 1:niter+10
pt = M[pick(p)]*pt
if i > 10
circle(GC, pt[1], pt[2], scale)
stroke(GC)
end
end
end

Out[132]:
# methods for generic function ifs
ifs(GC,M,p,niter) at In[132]:5
ifs(M,p,niter) at In[59]:4
ifs(GC,M,p,niter,scale) at In[132]:5
In [133]:
# a convenient utility for setting up Cairo drawings with an initial fill color and coordinate system
function drawing(w,h; fill=nothing, coords=nothing, left=0, right=w, top=0, bottom=h)
s = CairoRGBSurface(w, h)
gc = creategc(s)
if coords != nothing
left, right, top, bottom = coords
end
Graphics.set_coords(gc, 0, 0, w, h, left, right, top, bottom)
if fill != nothing
set_source(gc, fill)
paint(gc)
end
set_source(gc, RGB(0,0,0))
gc
end

Out[133]:
# methods for generic function drawing
drawing(w,h) at In[133]:4
In [145]:
d = drawing(275, 250, coords=[-4, 4.5, 10.2, 0], fill=RGB(.94,.92,.9))
set_line_width(d,0.5)

In [146]:
@time ifs(d, fern_M, fern_P, 16000)

elapsed time: 0.204639029 seconds (2307208 bytes allocated)

In [147]:
d.surface

Out[147]:
In [144]:
# an IFS for the sierpinski triangle
# here you can see the relationship between the matrices and the explicit algorithm above
sier_M = {[0.5 0 1; 0 0.5 1; 0 0 1],
[0.5 0 1; 0 0.5 4; 0 0 1],
[0.5 0 4; 0 0.5 4; 0 0 1]}
sier_P = cumsum([1/3, 1/3, 1/3]);

In [148]:
d = drawing(275, 250, coords=[0, 9, 9, 1], fill=RGB(.94,.92,.9))
set_line_width(d,0.5)

In [149]:
ifs(d, sier_M, sier_P, 20000); d.surface

Out[149]:
In [154]:
tree_M = {[0     0     0;  0    0.5  0  ; 0 0 1],
[0.42 -0.42  0;  0.42 0.42 0.2; 0 0 1],
[0.42  0.42  0; -0.42 0.42 0.2; 0 0 1],
[0.1   0     0;  0    0.1  0.2; 0 0 1]}
tree_P = cumsum([.05,.4,.4,.15]);

In [155]:
d = drawing(275, 250, coords=[-0.4, 0.4, 0.5, 0], fill=RGB(.94,.92,.9))
set_line_width(d,0.5)

In [156]:
ifs(d, tree_M, tree_P, 10000, 0.001); d.surface

Out[156]:
In [ ]: