Kaleidoscopic IFS Fractals

By: Gijs Bellaard
Created: 29 September 2019
Last modified: 29 September 2020

The best way to understand how Kaleidoscopic Iterated Function System(KIFS) Fractals work is to just start building one naively and keep working on it until the general form reveals itself. Rendering KIFS fractals is based on Distance fields and I therefor would suggest familiarizing yourself with the concept of Distance Fields before continuing. The code is written in GLSL because almost all real-time rendering of these fractals is done through fragment shaders.

Recreating a simple fractal

We will recreate the following fractal:

This fractal is made by repeatedly placing squares on the sides of all the squares, that are three times smaller, starting with a unit square. So let's start with a unit square centered at the origin:

float distanceFractal(vec2 p){
    float dis = INFINITY;
    
    dis = min(dis, distanceUnitSquare(p));
    
    return dis;
}

We continue by placing the four squares on the sides:

float distanceFractal(vec2 p){
    float dis = INFINITY;
    
    dis = min(dis, distanceUnitSquare(p));
    
    dis = min(dis, distanceUnitSquare(p*3.+vec2(4,0))/3.);
    dis = min(dis, distanceUnitSquare(p*3.-vec2(4,0))/3.);
    dis = min(dis, distanceUnitSquare(p*3.+vec3(0,4))/3.);
    dis = min(dis, distanceUnitSquare(p*3.-vec3(0,4))/3.);
    
    return dis;
}

Because we are working with distance fields this can be drastically optimized using some folds. We will fold the space towards the cube that is placed at positive x.

float distanceFractal(vec2 p){
    float dis = INFINITY;
    
    dis = min(dis, distanceUnitSquare(p));

    if(p.y>p.x)  p = vec2(p.y,p.x);
    if(p.y<-p.x) p = vec2(-p.y,p.x);
    
    dis = min(dis, distanceUnitSquare(p*3.-vec2(4,0))/3.);
    
    return dis;
}

This gives us the exact same distance field but we have shaved of three distance evaluations. The next step is to realise we don't actually want to place cubes on the side of the unit cube, but the fractal itself:

float distanceFractal(vec2 p){
    float dis = INFINITY;
    
    dis = min(dis, distanceUnitSquare(p));

    if(p.y>p.x)  p = vec2(p.y,p.x);
    if(p.y<-p.x) p = vec2(-p.y,p.x);
    
    dis = min(dis, distanceFractal(p*3.-vec2(4,0))/3.);
    
    return dis;
}

Mathematically this would be correct, however in a practical sense this is total bogus: the function never halts. Also, even worse, in GLSL recursion is not allowed! Because the recursion does not ''branch'' it is easily unrolled. We also hardcode a limit, i.e. the amount of iterations. All in all our final code is:

float distanceFractal(vec2 p){
    float dis = INFINITY;
    float s = 1.;
    
    for(int i=0;i<ITERATIONS;i++){
        dis = min(dis, distanceUnitSquare(p)*s);

        if(p.y>p.x)  p = vec2(p.y,p.x);
        if(p.y<-p.x) p = vec2(-p.y,p.x);

        p   *= 3.;
        s   /= 3.;
        p.x -= 4.;
    }
    
    return dis;
}

This form reveals why we call these iterated function system fractals: we have a function that performs some operations (folds, scalings and offsets), that is repeatedly (for loop) applied on the input (p).