"I write descriptions of natural forms on the walls, scratching them on the tile surface with a diamond."
–Donald Barthelme, "Game"

About the Pictures



I’ve had several requests to explain the mathematics behind the pictures I’ve generated for this blog. I’m flattered by your interest and hope I can both slake and pique your curiosity with this post, which I may refer back to and update in the future.

(The formatting for this may not translate well to your dashboard; on mobile it will look disastrous. I’d suggest viewing it on the full blog.  If you’re viewing this on the dashboard, click that folded-over upper right-hand corner.)

On a language and environment level, I’m programming in Objective-C using Xcode. I still think in terms of BASIC, though, having grown up with a hand-me-down Commodore VIC-20 as old as I was. In college, in the last days of DOS compatibility, I started exploring these fractal spaces with  (obsolete but approachable) QBASIC. Really, all you have to understand is pixels and for-loops. It’s great if you can write directly to bitmaps, but you can (I did) still get by drawing to the screen and then using screen-capture. You can even do some of this in Excel with conditional color-formatting of cells. I’ll include some examples here in pseudocode, so you can translate it to whatever environment and language you’re familiar with.

So there are two types of bitmap-building algorithms I use. We’ll work up to them from mathematical examples.

First, let’s talk about something called the Sierpinski Carpet, which is typically defined using a substitution system. You start with a square. Let’s denote a filled-in square with a “1” and an empty square space with a “0”:


    Filled square:   1

    Empty square:    0


Now we can “zoom in” to either of these and replace them with a 3x3 grid:


    Zoom in:    1   ->   1 1 1        0   ->   0 0 0
                         1 1 1                 0 0 0
                         1 1 1                 0 0 0


To make the Sierpinski Carpet, we remove the center filled-in square when we zoom in:


    Sierpinksi Carpet:

                1   ->   1 1 1        0   ->   0 0 0
                         1 0 1                 0 0 0
                         1 1 1                 0 0 0


Then you repeat the square-to-eight replacement indefinitely.


    1   ->   1 1 1   ->   1 1 1 1 1 1 1 1 1
             1 0 1        1 0 1 1 0 1 1 0 1
             1 1 1        1 1 1 1 1 1 1 1 1
                          1 1 1 0 0 0 1 1 1
                          1 0 1 0 0 0 1 0 1
                          1 1 1 0 0 0 1 1 1
                          1 1 1 1 1 1 1 1 1
                          1 0 1 1 0 1 1 0 1
                          1 1 1 1 1 1 1 1 1


    ->   1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1
         1 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 0 1
         1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1
         1 0 1 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 0 1
         1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1
         1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1
         1 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 1
         1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1
         1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1
         1 0 1 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 0 1
         1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1
         1 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 0 1 1 0 1 0 0 0 1 0 1
         1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 1 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
         1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1 1 0 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1


    -> etc.


Here it is again with non-numeric ASCII formatting:


    *   ->   * * *   ->   * * * * * * * * * 
             *   *        *   * *   * *   * 
             * * *        * * * * * * * * * 
                          * * *       * * * 
                          *   *       *   * 
                          * * *       * * * 
                          * * * * * * * * * 
                          *   * *   * *   * 
                          * * * * * * * * * 


    ->   * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         *   * *   * *   * *   * *   * *   * *   * *   * *   * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         * * *       * * * * * *       * * * * * *       * * * 
         *   *       *   * *   *       *   * *   *       *   * 
         * * *       * * * * * *       * * * * * *       * * * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         *   * *   * *   * *   * *   * *   * *   * *   * *   * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         * * * * * * * * *                   * * * * * * * * * 
         *   * *   * *   *                   *   * *   * *   * 
         * * * * * * * * *                   * * * * * * * * * 
         * * *       * * *                   * * *       * * * 
         *   *       *   *                   *   *       *   * 
         * * *       * * *                   * * *       * * * 
         * * * * * * * * *                   * * * * * * * * * 
         *   * *   * *   *                   *   * *   * *   * 
         * * * * * * * * *                   * * * * * * * * * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         *   * *   * *   * *   * *   * *   * *   * *   * *   * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         * * *       * * * * * *       * * * * * *       * * * 
         *   *       *   * *   *       *   * *   *       *   * 
         * * *       * * * * * *       * * * * * *       * * * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 
         *   * *   * *   * *   * *   * *   * *   * *   * *   * 
         * * * * * * * * * * * * * * * * * * * * * * * * * * * 


    -> etc.


Mathematically, the fractal is defined as the set of points which are never removed from the initial square (the infinite intersection of all iterations, with a “zoom in” or lit-square-to-eight replacement considered as a clarification of the previous, less-detailed stage). Of course we can’t draw an infinitely detailed picture, but we can produce iteration “N” (starting with N=0 as the solid square) on a canvas of 3^N x 3^N pixels. In fact, we can do that pixel-by-pixel without having to compute every intermediary stage.

As a convention, we’ll index the pixels starting at 0, so the range of pixels in each direction for iteration N is 0 to 3^N - 1. Here’s an algorithm in pseudocode which will generate the picture we want:


    int x, y, iteration    // indices
    int n = 5              // for example
    bool writePixel

    For x = 0 to 3^n - 1
       For y = 0 to 3^n - 1

          writePixel = TRUE    // so each pixel will be written unless …

          For iteration = 0 to n

             // … unless at any iteration x and y both fall within the middle third of that square:

             If
               x\(3^iteration) mod 3 = 1 And y\(3^iteration) mod 3 = 1
             Then
               writePixel = FALSE
             End If

             // mathematics unpacked below

          Next iteration

          If writePixel = TRUE Then setPixel(color, x, y)    // draw "color" at coordinates x, y

       Next y
    Next x


At the heart of the program, we have an expression including “x\(3^iteration) mod 3” — what does that mean? Well, the backwards slash is integer division, which discards the fractional remainder of division: a\b = int(a/b). And mod refers to modular arithmetic. These operations “stretch out” and “stripe out” the number to evaluate where it falls within each iteration of thirds-dividing. Here, take a look:


    x            0  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
    x mod 3      0  1  2  0  1  2  0  1  2  0  1  2  0  1  2  0  1  2  0  1  2  0  1  2  0  1  2

    x\3          0  0  0  1  1  1  2  2  2  3  3  3  4  4  4  5  5  5  6  6  6  7  7  7  8  8  8
    x\3 mod 3    0  0  0  1  1  1  2  2  2  0  0  0  1  1  1  2  2  2  0  0  0  1  1  1  2  2  2

    x\9          0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2


So using these operations, the program numerically tests the coordinates to see if they fall in the middle of any square at any scale.

If you can implement this, then you can explore the beautiful variations. For one thing, you can turn writePixel into another integer variable and link that to the pixel color — instead of white and black, you can have a gradient showing how well each pixel conforms to the condition you’ve established. See an example here for which the pixel is brighter the more scales at which the Sierpinski Carpet criteria (i.e., non-middle-square residence) are satisfied. You can consider that picture the sum of 1-and-0 grids from various iterations.


     1 1 1 1 1 1 1 1 1         1 1 1 1 1 1 1 1 1         2 2 2 2 2 2 2 2 2
     1 1 1 1 1 1 1 1 1         1 0 1 1 0 1 1 0 1         2 1 2 2 1 2 2 1 2
     1 1 1 1 1 1 1 1 1         1 1 1 1 1 1 1 1 1         2 2 2 2 2 2 2 2 2
     1 1 1 0 0 0 1 1 1         1 1 1 1 1 1 1 1 1         2 2 2 1 1 1 2 2 2
     1 1 1 0 0 0 1 1 1    +    1 0 1 1 0 1 1 0 1    =    2 1 2 1 0 1 2 1 2
     1 1 1 0 0 0 1 1 1         1 1 1 1 1 1 1 1 1         2 2 2 1 1 1 2 2 2
     1 1 1 1 1 1 1 1 1         1 1 1 1 1 1 1 1 1         2 2 2 2 2 2 2 2 2
     1 1 1 1 1 1 1 1 1         1 0 1 1 0 1 1 0 1         2 1 2 2 1 2 2 1 2
     1 1 1 1 1 1 1 1 1         1 1 1 1 1 1 1 1 1         2 2 2 2 2 2 2 2 2


Furthermore, why stick with just this one shape? By defining an array and linking your conditionals to that, you can build a grid of whatever size and configuration you like. A 2x2 grid (that illustrates bitwise XOR) makes a particularly pleasing pattern, and when you stack the layers:


    1 1 1 1 0 0 0 0       1 1 0 0 1 1 0 0       1 0 1 0 1 0 1 0       3 2 2 1 2 1 1 0
    1 1 1 1 0 0 0 0       1 1 0 0 1 1 0 0       0 1 0 1 0 1 0 1       2 3 1 2 1 2 0 1
    1 1 1 1 0 0 0 0       0 0 1 1 0 0 1 1       1 0 1 0 1 0 1 0       2 1 3 2 1 0 2 1
    1 1 1 1 0 0 0 0   +   0 0 1 1 0 0 1 1   +   0 1 0 1 0 1 0 1   =   1 2 2 3 0 1 1 2
    0 0 0 0 1 1 1 1       1 1 0 0 1 1 0 0       1 0 1 0 1 0 1 0       2 1 1 0 3 2 2 1
    0 0 0 0 1 1 1 1       1 1 0 0 1 1 0 0       0 1 0 1 0 1 0 1       1 2 0 1 2 3 1 2
    0 0 0 0 1 1 1 1       0 0 1 1 0 0 1 1       1 0 1 0 1 0 1 0       1 0 2 1 2 1 3 2
    0 0 0 0 1 1 1 1       0 0 1 1 0 0 1 1       0 1 0 1 0 1 0 1       0 1 1 2 1 2 2 3


You get something like this, the first picture I posted to this blog. Just to be clear: an algorithm nested like the one I have above will calculate the sum of all the layers (i.e., iterations) pixel by pixel, not layer-by-layer. This is easier on system resources than producing a bitmap or integer array for each layer and then adding each equivalently placed pixel together into a final layer. But the result is the same.

(Programming side note: some languages don’t have an integer exponentiation operator built in — for example, Objective-C, in which I actually write my programs.  In my case, I define an additional variable called “scale” which, in the example above, is identical to 3^iteration. That is, it starts as “1” where writePixel is defined and then is multiplied by 3 after each iteration. For different fractals, replace “3” with the appropriate magnification factor.)

All of the colored fractals I’ve produced are variations on this type of algorithm as well. In that case, your writePixel entity must contain three values, one each for red, green, and blue (or however your color space is set up). These values are calculated by reference to three substitution grids, again one for each color. You can also vary the substitution grids at each step: reversing, rotating, randomizing, whatever your imagination seizes upon. The variations are endless and beautiful.

So that’s the first type of grid-based fractal. Now let’s talk about recursion.

Actually, let’s step back a second. Here’s another grid replacement rule to consider.


    1   ->   1 1
             1 0
 


It is left as an exercise to the reader to show that this produces a right-angled Sierpinski Triangle in the basic black-and-white substitution scheme outlined before. However, there’s another way to build this exact same pattern with a different algorithm. Let’s talk about recursion.

So what we’re going to recreate here is something called “Pascal’s Triangle modulo 2.” That’s all Wikipedia stuff, though. Right now, all you need to consider is a different way of determining pixels. This will again be a simple two-color scheme using 1s and 0s, but instead of performing mathematics on the numerical coordinates of each pixel, we’re going to (1) first establish some initial conditions — i.e., assign states to pixels beforehand, then (2) sweep across the area under consideration using a recursion scheme — i.e., a rule for which the calculated state of each pixel will depend on the states of the pixels nearby.

The recursion rule in this first case says, to determine the subsequent state of any given pixel, consider that pixel itself, the pixel directly above, and the pixel directly to the left. Add these three values modulo 2.


    *  *      use the rule      x0 x1      which we summarize as    0 1
    *  ?                        x1 x1                               1 1
                                + mod 2                             modulus = 2


If a neighboring pixel is outside of the area under consideration (i.e., beyond the edge of the bitmap) then we will ignore it.

You may be intuiting generalizations already. This is good.

The initial conditions are: Pixel(0, 0) = 1, everything else = 0. Just the one nonzero point at the origin. This is the sole initial condition I’ve used for all the bitmaps I’ve posted, although you can investigate all sorts of interesting alternative patterns. For example, if you populate the upper and left edges of your bitmap under this scheme with random 0- and 1-pixels, then the rest of the bitmap is split about 50-50 between the two states, which as you’ll see is distinct from the pattern that emerges with just the single initial nonzero point.

So let’s carry out this recursion scheme and see how it goes. Remember, we start with a single nonzero point at the origin:


    1


This point at the origin won’t change since there are no visible cells that could alter it. We proceed outward (to the right and downward) by having each new cell correspond to the mod-2 addition of the cell to the left and the cell above, if either is present.

An area of 2x2 appears like so:


    1 1
    1 0


Note how the original pattern, the single “1,” is repeated to the right and below the original with a void expanse in the lower right quadrant. Now let’s carry it out to an area of 4x4:


    1 1 1 1
    1 0 1 0
    1 1 0 0
    1 0 0 0


The pattern repeats.


    1 1 1 1 1 1 1 1
    1 0 1 0 1 0 1 0
    1 1 0 0 1 1 0 0
    1 0 0 0 1 0 0 0
    1 1 1 1 0 0 0 0
    1 0 1 0 0 0 0 0
    1 1 0 0 0 0 0 0
    1 0 0 0 0 0 0 0


And again further:


    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
    1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0
    1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0
    1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0
    1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0
    1 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0
    1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0
    1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
    1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
    1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0
    1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0
    1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
    1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
    1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
    1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0


In fact, it repeats likewise at every doubled expanse: the smaller pattern is reproduced adjacently to the right and below, with the lower right quadrant empty. This corresponds exactly with the Iterated Function System definition of the right-angled Sierpinski Triangle, and the pattern also matches the substitution algorithm mentioned before. There are other ways to prove the correspondence as well if you refer back to the “Pascal’s Triangle mod 2” interpretation of this arithmetical procedure. Many roads lead to this pattern. (I like to think of it as a Triforce of Triforces.)

Here’s some pseudocode to express what we just did.


    int x, y, xOffset, yOffset    // indices
    int pixelState
    int recursionGrid[2][2] = {{0, 1}, {1, 1}}
      //  recursion grid as summarized earlier; indexing starts at 0
    int modulus = 2

    setPixel(1, 0, 0)    // initial conditions: point(0, 0) has color 1

    For x = 0 to 255     // an example size; x increases rightward, y downward
       For y = 0 to 255

          pixelState = 0    // reset pixelState before computing new value

          For xOffset = -1 to 0
             For yOffset = -1 to 0

                If x + xOffset >= 0 and y + yOffset >= 0 Then
                   pixelState =
                     pixelState + recursionGrid[xOffset + 1]⁠[yOffset + 1]
                     * getPixel(x + xOffset, y + yOffset)
                   // treat previous three lines as one; mathematics unpacked below
                End If

             Next yOffset
          Next xOffset

          pixelState = pixelState mod modulus
          setPixel(pixelState, x, y)

       Next y
    Next x


Okay, again we have a relatively simple algorithm with a mathematical knot at the center. Let’s untangle that. Basically, what we’re doing is applying those multiplicative weights — stored in the recursionGrid array — to the nearby pixel values, accessed by getPixel, provided that they’re not off the grid. We add each appropriate weight-times-value to the variable pixelState. After all the values have been added, we adjust the total per the modulus. Then we move on to the next pixel in our sweep.

Note that this only sweeps across the bitmap once. To get an animation, you could nest the bitmap-sweep loops (x and y) inside an another loop to see how it proceeds further. (Note to self: do this.) Strictly speaking, this isn’t a cellular automaton in the Conway’s Game of Life sense, since the order in which the cells are computed makes a difference. But it falls within the broader definition of cellular automata and we could always add an extra state to the cells describing which are changeable (all the initial zeroes) and which are permanent (the 1-cell at the origin and any computed value left- or below-adjacent to a permanent cell) to force it more in line with conventions. I also suspect we could do some coordinate transformation to make it equivalent to a Wolfram-style one-dimensional recursion as well, but I haven’t bothered to shore up that intuition with proof yet. Regardless, I’m still filing these under the “cellular automata” tag since they partake of the same visual aesthetic; react as you will.

Here’s one variation left as another exercise for the reader: take a look at the pattern generated by the following recursion grid:


    1 1
    1 1
    modulus = 3


You may be pleasantly surprised.

Now, for the final trick, we can combine the monochromatic map generated by recursion (mod 2) with the zoom-layering as described for the grid-based fractals. This makes a lovely gradient picture. For extra beauty (with little mathematical import, but great aesthetic charm), you can arrange four copies of the gradient in a symmetrical pattern, origins at the outer corners or at the center. Beautiful! I have a backlog of about a thousand (2^10) of these, several hundred queued.

Furthermore, and relevant to that backlog story, once you’ve figured out any particular design you can generalize and iterate over the rule-space surrounding it. Here’s how: figure out (count) how many possible recursion/substitution grids there are for any given conceptual framework, then devise a way to translate each of those integers into a specific rule (probably with integer division and modular arithmetic). Place that translation within another for-loop, set up some memory management, and you’ll generate hundreds or thousands of pictures at one go. For things like color spaces in which the variations are effectively infinite, randomize and assemble your favorites into a bouquet.

It’s a beautiful world out there in the abstract spaces. Let’s go exploring.



  1. boombadoom-he-got-that-superbass reblogged this from neutrinoprism
  2. nebulaprocella reblogged this from neutrinoprism
  3. quondam reblogged this from neutrinoprism
  4. 1cetec reblogged this from neutrinoprism
  5. nosuch-thing reblogged this from neutrinoprism
  6. mooooooosings-on-the-whiteboard reblogged this from neutrinoprism
  7. invisioned said: Omg
  8. neutrinoprism posted this