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]:
RdGy
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 [ ]: