# Iterated Function Systems

A few quick notes before we start:

1. For this chapter, we will be following the methodology set by the plotting chapter. That is to say that the code presented in this chapter will output another file that can be easily plotted by an external plotter. If you like to use a plotter provided by your language of choice, please modify the code provided to do so.

2. This chapter is currently a subsection to the plotting chapter, but we may extend the algorithm archive in the future with other fractal generation methods, which would require created a new section on fractals, in particular. This would include a chapter with more rigorous definitions on fractals, which is largely missing from the following discussion. Please let us know if you are interested!

In this chapter, we will show you how to make one of the most famous fractals via Iterated Function Systems (IFSs): the Sierpinski triangle. We will also introduce a number of interesting concepts for further exploration, such as chaos games, Hutchinson operators, and attractors.

## The Sierpinski Triangle

To begin the discussion of Iterated Function Systems (IFSs), we will first discuss what might be one of the most famous fractals currently known, the Sierpinski triangle, shown below: This image is clearly a set of triangles embedded in a larger triangle in such a way that it can be continually cut into three identical pieces and still retain its internal structure. This idea is known as self-similarity , and it is usually the first aspect of fractals that will catch an audience's attention. In fact, there are plenty of uses of fractals and their mathematical underpinnings, such as estimating the coastline of Britain , identifying fingerprints , and image compression . In many more rigorous definitions, a fractal can be described as any system that has a non-integer Hausdorff dimension . Though this is an incredibly interesting concept, the discussion of this chapter will focus primarily on methods to generate fractal patterns and will refrain from delving into this discussion for now and instead focus on the methods by which fractals can be generated through IFSs.

To start, imagine creating a triangle from three points, , , and . These points can be arbitrarily chosen, but for this conversation, we will constrict them to the vertices of an equilateral triangle, as shown below: Now let's create three separate functions that can act on a 2-dimensional space:

Each function will input a particular location in space (here, ) and output a new location that is the midpoint between the input location and , , or for functions 1, 2, and 3 respectively. The union of all of these functions (the set of all possible functions available for use) is often notated as the Hutchinson operator , and for this case it would look like this:

By iteratively using this operator, we can traverse through all possible movements in the set. For example, let's generate 3 new points that are halfway between and , and , and and , which will be called , , and respectively. This is shown below: From here, each new point (, , and ) will spawn 3 children, and each child will move according to one of the three possible functions in the Hutchinson operator, as shown below:

Here, all red children come from , green children come from and blue children come from . From here, the children will then spawn 3 more children, each of which will move according to a different function. Those children will then spawn more children, who act accordingly. As this process continues on and on, we begin to see an interesting pattern form:

This is the Sierpinski triangle. At first, it might seem like mathematical magic that a simple set of 3 functions can create such a pattern. After all, why aren't any of the children migrating to the empty spaces in the structure? This will require some thought, but the simplest answer is that no function within the Hutchinson operator allows for children to enter those spaces.

When I learned about this for the first time, I began to wonder about other shapes. Could we create fractal squares? Hexagons? Circles? Such shapes seem like natural extensions to the triangular Hutchinson operator provided above, but there's a bit of a hitch...

First, let's take 4 points, , , , and , this time located at the four vertices of a square, like so: In a similar fashion, we'll create 4 functions with , and such that:

If we then create 5 initial points located between all the vertices and allow these points to continually spawn children like before, something peculiar happens:

Well, this is essentially a square of squares. What happened to the self-similar structure we were getting before? Why isn't this more interesting?

The best answer I have for now is that some Hutchinson operators are interesting and some are not. On the other hand, this square is a bit more interesting than it first appears, but I think this is easiest to understand when we use the Hutchinson operator in a slightly different way.

## Chaos games and attractors

Until now, we have been using the Hutchinson operator in a computationally costly way. Every iteration, we generate 3 or 4 new children per child per step of the simulation. This scales exponentially and means that we will quickly have millions of children to keep track of!

A much more computationally feasible method to use the Hutchinson operator comes in the form of the chaos game . Here, instead of tracking children of children, we track a single individual that chooses randomly between the Hutchinson functions, as shown here:

function chaos_game(n::Int, shape_points)

# Initializing the output array and the initial point
output_points = zeros(n,2)
point = [rand(), rand()]

for i = 1:n
output_points[i,:] .= point
point = 0.5*(rand(shape_points) .+ point)
end

return output_points

end

// This is a function to simulate a "chaos game"
PointVector chaosGame(int numOutputPoints, const PointVector& inputPoints) {
// Choose first point randomly
Point curPoint = {drand(), drand()};

// For each output point, compute midpoint to random input point
PointVector outputPoints(numOutputPoints);
for (auto& outPoint : outputPoints) {
outPoint = curPoint;
curPoint = 0.5 * (curPoint + choose(inputPoints));
}

return outputPoints;
}

def chaos_game(n, shape_points):
# Initialize the starting point
point = [random(), random()]

for _ in range(n):
# Update the point position and yield the result
point = [(p + s) / 2 for p, s in zip(point, choice(shape_points))]
yield point

void chaos_game(struct point *in, size_t in_n, struct point *out,
size_t out_n) {

struct point cur_point = {drand(), drand()};

for (int i = 0; i < out_n; ++i) {
out[i] = cur_point;
struct point tmp = random_element(in, in_n);
cur_point.x = 0.5 * (cur_point.x + tmp.x);
cur_point.y = 0.5 * (cur_point.y + tmp.y);
}
}

data point(x=0, y=0):
return point(self.x + other.x, self.y + other.y)

def __rmul__(self, other):
return point(self.x * other, self.y * other)

def chaos_game(n, shape_points):
p = point(random(), random())

for _ in range(n):
p = (1/2) * (p + choice(shape_points))
yield p


If we set the initial points to the on the equilateral triangle we saw before, we can see the Sierpinski triangle again after a few thousand iterations, as shown below:

That said, there is something peculiar about the way the chaos game starts. Essentially, our lone child begins their journey at a random location. What if that location is off the triangle? Well, let's test that!

If we start the point at the center of the triangle, it will eventually find its way onto the triangle, and then all subsequent iterations will be on the fractal. Here, I have plotted the first 20 steps where the wanderer is still looking for the correct shape:

It would seem that the wanderer is attracted to the Sierpinski triangle with this set of functions. That is actually the case. The truth is that the word attractor is a very loaded term in the literature, but for the purposes of our discussion here, the attractor is any shape defined by the iteration through Hutchinson operator functions.

So let's go back to the square, which seemed like a somewhat random distribution of points: Even with the chaos game, this will not change; however, we now know that the random distribution isn't truly random. Rather, it's an attractive plane where our lone wanderer can exist happily within. That is to say, the 2-dimensional square is, itself, the attractor for that Hutchinson operator, and if we start our journey person off of the square, they will eventually find themselves within it, similar to the triangle before.

The next obvious question is whether a square can create any more interesting fractally patterns, and the answer is "yes, but only if we restrict the movement a bit." Which brings us to another topic entirely: restricted chaos games. Discussing restricted chaos games in more detail is a chapter in its own right, so I will forego the discussion here. If you are interested, please let me know and I will be more than willing to add the chapter in the future!

## Video Explanation

Here is a video describing iterated function systems:

## Example Code

For the code in this chapter, we have decided to write it specifically for the Chaos game, not the hutchinson animations shown at the start of the chapter. This is because that animation is slightly tricky to create and distracts from the overall purpose of this chapter. In addition, we have written the chaos game code to take in a set of points so that it is not hard-coded for the Sierpinski triangle and can be easily extended to other shapes like the square or restricted chaos games, as we mentioned before!

using DelimitedFiles

# This is a function to simulate a "chaos game"
function chaos_game(n::Int, shape_points)

# Initializing the output array and the initial point
output_points = zeros(n,2)
point = [rand(), rand()]

for i = 1:n
output_points[i,:] .= point
point = 0.5*(rand(shape_points) .+ point)
end

return output_points

end

# This will generate a Sierpinski triangle with a chaos game of n points for an
# initial triangle with three points on the vertices of an equilateral triangle:
#     A = (0.0, 0.0)
#     B = (0.5, sqrt(0.75))
#     C = (1.0, 0.0)
# It will output the file sierpinski.dat, which can be plotted after
shape_points = [[0.0, 0.0],
[0.5, sqrt(0.75)],
[1.0, 0.0]]
output_points = chaos_game(10000, shape_points)
writedlm("out.dat", output_points)

#include <cmath>
#include <fstream>
#include <random>
#include <vector>

// Simple X-Y point structure, along with some operators
struct Point {
double x, y;
};

Point operator+(Point lhs, Point rhs) { return {lhs.x + rhs.x, lhs.y + rhs.y}; }
Point operator*(double k, Point pt) { return {k * pt.x, k * pt.y}; }
Point operator*(Point pt, double k) { return k * pt; }

using PointVector = std::vector<Point>;

// Returns a pseudo-random number generator
std::default_random_engine& rng() {
// Initialize static pseudo-random engine with non-deterministic random seed
static std::default_random_engine randEngine(std::random_device{}());
return randEngine;
}

// Returns a random double in [0, 1)
double drand() {
return std::uniform_real_distribution<double>(0.0, 1.0)(rng());
}

// Returns a random integer in [0, numElems-1]
std::size_t randrange(std::size_t numElems) {
return std::uniform_int_distribution<std::size_t>(0, numElems - 1)(rng());
}

// Return a random point from the non-empty PointVector
Point choose(const PointVector& points) {
return points[randrange(points.size())];
}

// This is a function to simulate a "chaos game"
PointVector chaosGame(int numOutputPoints, const PointVector& inputPoints) {
// Choose first point randomly
Point curPoint = {drand(), drand()};

// For each output point, compute midpoint to random input point
PointVector outputPoints(numOutputPoints);
for (auto& outPoint : outputPoints) {
outPoint = curPoint;
curPoint = 0.5 * (curPoint + choose(inputPoints));
}

return outputPoints;
}

int main() {
// This will generate a Sierpinski triangle with a chaos game of n points for
// an initial triangle with three points on the vertices of an equilateral
// triangle.
PointVector inputPoints = {{0.0, 0.0}, {0.5, std::sqrt(0.75)}, {1.0, 0.0}};
auto outputPoints = chaosGame(10000, inputPoints);

// It will output the file sierpinski.dat, which can be plotted after
std::ofstream ofs("sierpinski.dat");
for (auto pt : outputPoints)
ofs << pt.x << '\t' << pt.y << '\n';
}

from random import random, choice
from math import sqrt

# This generator simulates a "chaos game"
def chaos_game(n, shape_points):
# Initialize the starting point
point = [random(), random()]

for _ in range(n):
# Update the point position and yield the result
point = [(p + s) / 2 for p, s in zip(point, choice(shape_points))]
yield point

# This will generate a Sierpinski triangle with a chaos game of n points for an
# initial triangle with three points on the vertices of an equilateral triangle:
#     A = (0.0, 0.0)
#     B = (0.5, sqrt(0.75))
#     C = (1.0, 0.0)
# It will output the file sierpinski.dat, which can be plotted after
shape_points = [[0.0, 0.0],
[0.5, sqrt(0.75)],
[1.0, 0.0]]
with open("sierpinski.dat", "w") as f:
for point in chaos_game(10000, shape_points):
f.write("{0}\t{1}\n".format(*point))

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

struct point {
double x, y;
};

double drand() {
return ((double) rand() / (RAND_MAX));
}

struct point random_element(struct point *array, size_t n) {
return array[rand() % (int)n];
}

void chaos_game(struct point *in, size_t in_n, struct point *out,
size_t out_n) {

struct point cur_point = {drand(), drand()};

for (int i = 0; i < out_n; ++i) {
out[i] = cur_point;
struct point tmp = random_element(in, in_n);
cur_point.x = 0.5 * (cur_point.x + tmp.x);
cur_point.y = 0.5 * (cur_point.y + tmp.y);
}
}

int main() {
struct point shape_points  = {{0.0,0.0}, {0.5,sqrt(0.75)}, {1.0,0.0}};
struct point out_points;

srand(time(NULL));

chaos_game(shape_points, 3, out_points, 1000);

FILE *fp = fopen("out.dat", "w+");

for (int i = 0; i < 1000; ++i) {
fprintf(fp, "%f\t%f\n", out_points[i].x, out_points[i].y);
}

fclose(fp);

return 0;
}

from math import sqrt
from random import random, choice

data point(x=0, y=0):
return point(self.x + other.x, self.y + other.y)

def __rmul__(self, other):
return point(self.x * other, self.y * other)

def chaos_game(n, shape_points):
p = point(random(), random())

for _ in range(n):
p = (1/2) * (p + choice(shape_points))
yield p

# This will generate a Sierpinski triangle with a chaos game of n points for an
# initial triangle with three points on the vertices of an equilateral triangle:
#     A = (0.0, 0.0)
#     B = (0.5, sqrt(0.75))
#     C = (1.0, 0.0)
# It will output the file sierpinski.dat, which can be plotted after
shape_points = [point(0.0, 0.0),
point(0.5, sqrt(0.75)),
point(1.0, 0.0)]
with open("sierpinski.dat", "w") as f:
for p in chaos_game(10000, shape_points):
f.write("{0}\t{1}\n".format(p.x, p.y))


### Bibliography

 1 Wikipedia: Self-similarity, 2019. 2 Mandelbrot, Benoit, How long is the coast of Britain? Statistical self-similarity and fractional dimension, American Association for the Advancement of Science, 1967. 3 Jampour, Mahdi et al., A new fast technique for fingerprint identification with fractal and chaos game theory, World Scientific, 2010. 4 Wikipedia: Fractal Compression, 2019. 5 Saupe, Dietmar and Hamzaoui, Raouf, A review of the fractal image compression literature, ACM, 1994. 6 Sanderson, G, 3blue1brown: Fractals are typically not self-similar, 2017. 7 Wikipedia: Hausdorff dimension, 2019. 8 Gneiting, Tilmann and Ševčíková, Hana and Percival, Donald B, Estimators of fractal dimension: Assessing the roughness of time series and spatial data, JSTOR, 2012. 9 Wikipedia: Hutchinson Operator, 2019. 10 Hutchinson, John E, Fractals and self similarity, JSTOR, 1981. 11 Wikipedia: Chaos Game, 2019. 12 Wolfram: Chaos Game, 2019.

##### Code Examples 