The Night of Fractals

It was a late Sunday night, and I was thinking to go to sleep, when a friend of mine sent me a picture from some website and wrote: “Beautiful!”. I drew such pictures five years ago by using the so-called escape-time algorithm. But to apply this algorithm, it was necessary to know how to partition the plane into regions for a given set of transformations. I couldn’t figure it out, and did not get back to this algorithm anymore. But now I knew exactly what should I do and wrote him back: “Random IFS first, then kNN, and then the Escape-Time Algorithm!”

I had an old netbook at hand, as my laptop was under repair. My friend was saying something to me, and I was answering him, but in my mind I was already writing code, and was searching for at least some compiler on the netbook. And I found the C++ Builder 6! Then I realized that I would meet the morning with the Borland compiler. In five hours, I sent new pictures to my friend, but as any normal person, he was sleeping….

Let’s take look at some formulas. Suppose we have a finite set of transformations of plane T_i: R² →R², i = 1, …, k. For a random E set, we define that T(E) = T₁(E) ∪… ∪ T_k(E), that is to say that we will affect the E set with each transformation, and then combine the results. We can prove that if T_i representation was compressing, the E, T(E), T(T(E)),… sequence will reduce to the F set for any nonempty compact E. This construction is known as the iterated function system.

For example, if we take a smiley as a nonempty compact and consider three transformations, each of which is a composition of compressing and the moving to the i-th vertex of a regular triangle, the first iterations will look like the following, and we’ll obtain the Sierpinski triangle at the limit.

As a rule, instead of directly calculating the E, T(E), T(T(E)),… sequence we use the so-called Chaos game for creating fractals. The game is as follows. Choose a random z₀ point in the plane, then randomly choose the T_i1 transformation and calculate z₁= T_i1(z₀), and then randomly choose T_i2 and calculate z₁= T_i2(z₀), etc. We can show that everything will be fine, and the set of obtained points will, in some sense, approach the F set defined above. I will mention this algorithm bellow as Random IFS.

z = (0, 0)
for (i = 0; i < maxIterNum; ++i) {
    cl = random([p1, ..., pk]) // pi – the probability, with which we choose the Ti transformation.
    z = T[cl](z)
    if (i > skipIterNum) { // The first few iterations may be far enough from the attractor.
        draw(z)
    }
}

Escape-Time Algorithm

It’s high time to go to the description of the escape-time algorithm. Suppose we have one transformation of the f plane. For each z point of the plane, we will start calculating the f(z),f(f(z)), f(f(f(z)),… sequence, either till the number of iterations exceeds a predetermined N number, or the norm of the z number becomes greater than that of the B. After that, we choose the color of the point, according to the number of performed iterations.

for (y = y1; y < y2; y += dy) {
    for (x = x1; x < x2; x += dx) {
        z = (x, y);
        iter = 0;
        while (iter < N && ||z|| < B) {
             z = f(z)
             iter += 1;
        }
        drawPixel((x, y), color(iter))
    }
}

If we assume that our plane is complex, and the f(z) transformation is equal to z² + c, we’ll obtain a Julia set fractal as the result of this algorithm.

Suppose now that we have a system of iterated functions, defined by a set of reversible compressing transformations of the T₁, …, T_k plane. Assume that F is the attractor of this system.

In addition, we assume that we can partition the F set, so that T_i(F) ∩ T_j(F) = ∅, i != j (this assumption is not always satisfied). Partition the whole plane R₂ into pieces A₁, …., A_k, so that T_i(F) is the subset A_i for all i. Now, let’s define the f(z) function piecewise: put f(z) = T_k⁻¹(z) on the A_i set for all i.

For example, we’ll consider the following partition for the Sierpinski triangle (there are some problems with three points, but let’s turn a blind eye to them).

The most important question is what will happen if we apply the escape-time algorithm to the function f built in such way?

Let’s see:

We’ve obtained a lovely Sierpinski triangle!

Turns out, it was no accident. Here are a few more examples:

In these examples, we can easily define the corresponding partition analytically, with the help of Boolean combinations of circles and half-planes, using the method of gazing. But we often fail to guess simple conditions. So, instead of guessing, we will teach a computer to determine a partition by itself.

The k-nearest neighbors algorithm will help us here.

First, we generate several thousand points, using Random IFS. For each point, we memorize the number of transformation, with the help of which it has been obtained. During the operation of the Escape-Time Algorithm, for each pixel, we define the area where it gets with the help of 1NN.

For instance, 1NN gives the following partition into 4 pieces for this star:

Putting it together, we obtain:

points = RandomIFS(Ts)
classifier = kNN(points);
for (y = y1; y < y2; y += dy) {
    for (x = x1; x < x2; x += dx) {
        z = (x, y)
        iter = 0
        while (iter < maxIterNum && ||z|| < sqrdBound) {
            cl = classifier.getClass(z);
            z = T[cl].applyInvert(z);
            iter += 1;
        }
        draw((x, y), color(iter))
    }
}

A few more pictures.

That’s all. Finally, two remarks.

Firstly, the attentive reader can ask that if we can use the escape-time algorithm to build the fractals that are built with the help of Random IFS, is it possible to build the Julia set via Random IFS? Turns out, you can. You just need to invert the f(z) =z² + c mapping, remembering how to extract the square root of a complex number. (However, there are great difficulties when applying this method to make Julia set images).

x = z0.re
y = z0.im
for (i = 0; i < N; ++i) {
    x -= c.re;
    y -= c.im;
    len = sqrt(x*x + y*y);
    ang = atan2(y, x);
    if (rand() % 2 == 0) { // We need something more interesting here.
        x = sqrt(len) * cos(ang / 2);
        y = sqrt(len) * sin(ang / 2);
    } else {
        x = sqrt(len) * cos(ang / 2 + M_PI);
        y = sqrt(len) * sin(ang / 2 + M_PI);
    }
    draw(x, y)
}

Secondly, the article has told you what happens. If you want to know why this happens, I recommend the “Fractals Everywhere” by M. Barnsley.

Source code is available on github.