The Night of Fractals

Algorithms

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!”

Fractal

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….

Fractal

Fractal

Fractal

Let’s take look at some formulas. Suppose we have a finite set of transformations of plane Ti: R2 →R2, i = 1, …, k. For a random E set, we define that T(E) = T1(E) ∪… ∪ Tk(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 Ti 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.

Sierpinski triangle

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 z0 point in the plane, then randomly choose the Ti1 transformation and calculate z1= Ti1(z0), and then randomly choose Ti2 and calculate z1= Ti2(z0), 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)
    }
}

Fractal

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 z2 + c, we’ll obtain a Julia set fractal as the result of this algorithm.

Fractal Julia Set

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

In addition, we assume that we can partition the F set, so that Ti(F) ∩ Tj(F) = ∅, i != j (this assumption is not always satisfied). Partition the whole plane R2 into pieces A1, …., Ak, so that Ti(F) is the subset Ai for all i. Now, let’s define the f(z) function piecewise: put f(z) = Tk−1(z) on the Ai 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).

Fractal. Sierpinski Triangle.

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:

Fractal

We’ve obtained a lovely Sierpinski triangle!

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

Fractal

Fractal

Fractal

Fractal

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:

Fractal

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))
    }
}

Fractal

A few more pictures.

Fractal

Fractal

Fractal

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) =z2 + 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.

Comments

  1. Indeed beautiful.
940

Ropes — Fast Strings

Most of us work with strings one way or another. There’s no way to avoid them — when writing code, you’re doomed to concatinate strings every day, split them into parts and access certain characters by index. We are used to the fact that strings are fixed-length arrays of characters, which leads to certain limitations when working with them. For instance, we cannot quickly concatenate two strings. To do this, we will at first need to allocate the required amount of memory, and then copy there the data from the concatenated strings.