PostScript Fractals

In memory of Robert Chrismas
Pretty? It's a screenshot of the output from, which is a PostScript program for drawing fractals using random iteration.


Roughly speaking, any structure is said to be fractal if it is self-similar; i.e. if examining part of the structure yields something similar to the whole. Of course, this does all depend on how you define "similar". Look at the word "chaos" in the picture above; you should be able to see that every line in every letter of the word is in fact the word "chaos" again. Every line in those little "chaos" is another "chaos", and so on. The triangle (it's known as Sierpinski's Triangle) is fractal too - it's made of three smaller triangles, each of which is contains three smaller triangles.

It turns out that knowing exactly how the picture is self-similar is the only thing you need in order to reproduce it. We can describe the self-similarity of a fractal as a set of transformations, each of which maps the the image onto part of itself, such that combining all of the transformed copies produces the original image again.

The Mathsy Bit

We could represent transformations as 2x2 matrices, but we want to be able to translate too, so we express the transformations as 3x3 matrices with 6 nontrivial elements:

Under this scheme, Sierpinski's Triangle is represented by:

a  b  c  d  e  f
.5  0  0 .5  0  0
.5  0  0 .5 .5  0
.5  0  0 .5  0 .5
Briefly, each copy of the image is scaled to half size, keeping the bottom left corner (the origin of the coordinate space) fixed. One copy stays there, one is translated up by half the original size, one is translated right by that same amount.

We don't like infinite lists of transformations. So, while we could describe a fern as an infinitely many smaller ferns placed in two lines, we don't. Instead, we break down the fern's self-similarity into just four transformations. Note that we1 can describe the fern in terms of a finite number of copies of itself; one for the first leaf on each side, one squashed to a line for the stem, and one shrunk a little for the whole of the rest of the fern.
Fractal decomposition of a fern
This image is in the public domain.

To produce a plot of a fractal from a set of transformations, we use the random iterations algorithm [Barnsley]. In short, we repeatedly transform a single point, each time using a transformation randomly chosen from the list. If we plot the points visited, we get an approximation to the fractal; proof of this is well beyond the scope of this page.


PostScript is a stack-based language in which programs and data are described in Reverse Polish Notation. This is somewhat counterintuitive at first, but it makes it very easy to write a fast and robust parser, and for graphics programs to use it as an output language.

I am no expert on PS, having only started coding in it earlier this week - I make no representations as to the quality, reliability or conformity of my program. Specifically, bear in mind that if you send this file to a PostScript printer, the algorithm will be executed on the printer, not on your PC. What takes your machine a second to preview may take your printer several minutes to produce - you have been warned! [ 1.9K ]

My Code...

The first part of the program is the data for the three fractals. Each line specifies a transformation and a probability associated with it; these are used to weight the choice of transformation, so that the image produced appears to have an equal density of dots for all areas that are on the fractal.

The makeptable function takes a fractal specification and produces a table of cumulative probabilities.
    /makeptable {
        dup length array 0
        0 1 3 index length 1 sub {
            4 copy 4 -1 roll exch
            get 0 get add 
            3 -1 roll exch
            4 -1 roll pop dup 4 1 roll
        } for
        pop exch pop
        dup dup length 1 sub 1.0 put
    } bind def
In retrospect, makeptable was unnecessary, as you can produce a deviate from a given distribution without a cumulative probability table, just by repeatedly subtracting and comparing. Nonetheless, randidx uses the table to produce a random index into the list of transformations, weighted by the probabilities.
    /randidx {
        rand 2147483648.0 div
            1 index 3 index 2 index get le {
                3 1 roll pop pop exit
            } if
                1 add
        } loop
    } bind def
The real work is done in randfractal; this is an implementation of Barnsley's random iteration algorithm, described above. The .0025 constants are the width and height of the dot drawn on the page. Using larger dots results in the fractal requiring fewer iterations but produces a coarser image.
    /randfractal {
        dup makeptable /pmap exch def
        /eqns exch def
        /n exch def
        0 0
        1 1 n {
            10 gt { 2 copy .0025 .0025 rectfill} if
            eqns pmap randidx get 1 get transform
        } for
    pop pop
    } bind def
A quick function to conveniently scale and translate each rendered fractal...
    /render {
        matrix currentmatrix
        7 1 roll
        translate scale randfractal
    } bind def
... and finally the fractals can be plotted. Each line consists of a number of iterations, the name of the fractal, a position to translate to, a size to scale to, and a call to the render function.
    40000 Sierpinski 200 200 100 300 render
    200000 Fern 40 40 350 350 render
    100000 Chaos 25 25 50 150 render

PS can be small...

Note the size of the file: it's around 2K. That's less than a tenth the size of the screenshot at the top of this page (which is just under 30K). Unlike that screenshot, the PS file can be rendered at any resolution, and (if the number of iterations and dot size are adjusted) will produce as much detail as the output device can handle.

As an experiment, I decided to see exactly how small I could get the code. This was also a good way to figure out exactly what was and was not legal in the language. My favourite trick was deliberately nesting the '[]' and '{}' brackets "incorrectly". In most languages, this would be invalid code; in PS, the left bracket is just a marker, and may be manipulated by other instructions. I used this to help compress the data, and to implement stack frames which could be popped more efficiently.

/B{bind def}def/r{roll}B/i{index}B/N{0 0 5 2 r]}B/R{0 5 1 r 0 3 1 r]}B/H{20 div .2}B/V{20 div -.2}B/f{matrix currentmatrix[9 2 r translate scale 0 0 1 1 5 -1 r{10 gt{2 copy .0025 dup rectfill}if 3 i[4 i rand 21474836.48 div 0{1 i 3 i 2 i get le{4 1 r]pop exit}if 1 add}loop get transform}for]pop setmatrix}B[[.5 .5 0 0 N[.5 .5 .5 0 N[.5 .5 0 .5 N][33 66 100]40000 200 200 100 300 f[[0 .16 0 0 N[.85 -.04 .04 .85 0 1.6][.2 .23 -.26 .22 0 1.6][-.15 .26 .28 .24 0 .44]][1 86 93 100]200000 40 40 350 350 f[[5 V 1 0 R[5 V 5 0 R[5 V 7 0 R[4 V 11 1 R[3 V 13 1 R[3 V 16 1 R[2 V 9 0 R[1 H 5 2 N[1 H 17 3 N[1 H 19 1 N[2 H 1 0 N[2 H 1 4 N[2 H 8 4 N[2 H 8 2 N[2 H 9 0 N[3 H 17 4 N[3 H 17 2 N[3 H 17 0 N[4 H 12 0 N[4 H 12 4 N][8 16 24 31 36 41 45 47 49 51 55 59 63 67 71 76 81 86 93 100]100000 25 25 50 150 f
The best I managed (without changing the output at all) was 809 bytes :-)

Further information



[1] Actually, the fern is (like Sierpinski's Triangle) a pretty standard example; the first to use it was (probably) Barnsley


[Barnsley] Michael Barnsley, Fractals Everywhere (ISBN 0-12-079062-9)


Creative Commons License
© 2006 Will Benfold (email: will AT benfold DOT com)
This work is licensed under a Creative Commons Attribution 2.5 License.