June 2018 – BitWise

Here’s another famous fractal that I wanted to try on Arduino: the Newton fractal. It’s named after Isaac Newton’s iterative method for finding roots of functions. The classic example applies the method to function f(z) = z³ – 1 , where z is a complex number. This function has three roots in the complex plane.

$classic newton fractal$

The sketch that produced this picture loops over the pixels of a 480×320 display, mapping each of them to a complex number z₀, that will serve as the initial value for Newton’s iteration process:

z_n+1 = z_n – f(z_n) / f'(z_n)

The basic color of a pixel (red, green or blue) depends on to which of the three roots the process converges. Its color intensity depends on the number of iterations that were needed to reach that root within a pre-defined precision. It’s just that simple!

Apart from playing with different color mappings (always essential for producing visually appealing fractals), I wanted to use modified versions of Newton’s method, as well as to apply them to different functions. The Arduino core has no support for complex calculus, and a library that I found didn’t even compile. So I wrote a couple of basic complex functions and put them in a functions.h file. There must better ways, but it works for me.

Once you have a basic sketch, the canvas of complexity is all yours!

f(z) = z⁴ – 1

This is my functions.h file. It must be in the same folder as the fractal sketch.

struct complex {
    double im;
    double re;
};

struct complex create(double real, double imag){
    struct complex t;
    t.re = real;
    t.im = imag;
    return t;
};

struct complex cpow2(complex z){
    struct complex t;
    t.re = z.re*z.re - z.im*z.im;
    t.im = 2.0*z.re*z.im;
    return t;
};

struct complex cpow3(complex z){
    struct complex t;
    t.re = pow(z.re,3) - 3.0*z.re*z.im*z.im;
    t.im = 3.0*z.re*z.re*z.im - pow(z.im,3);
    return t;
};

struct complex cpow4(complex z){
struct complex t;
t.re = pow(z.re,4) + pow(z.im,4) - 6.0*z.re*z.re*z.im*z.im;
t.im = 4.0*(pow(z.re,3)*z.im - pow(z.im,3)*z.re);
return t;
};

struct complex cadd(complex c, complex z) {
    struct complex t;
    t.re = c.re + z.re;
    t.im = c.im + z.im;
    return t;
}

struct complex csubs(complex c, complex z) {    //  c minus z
    struct complex t;
    t.re = c.re - z.re;
    t.im = c.im - z.im;
    return t;
}

struct complex cmult(complex c, complex z) {
    struct complex t;
    t.re = c.re*z.re - c.im*z.im;
    t.im = c.re*z.im + c.im*z.re;
    return t;
}

struct complex cdiv(complex c, complex z) {     // c divided by z
    struct complex t;
    t.re = (c.re*z.re + c.im*z.im)/(z.re*z.re + z.im*z.im);
    t.im = (c.im*z.re - c.re*z.im)/(z.re*z.re + z.im*z.im);
    return t;  
}

double cabs(complex z) {
  return sqrt(z.re*z.re + z.im*z.im);
}

struct complex {

double im;

double re;

};

struct complex create(double real, double imag){

struct complex t;

t.re = real;

t.im = imag;

return t;

};

struct complex cpow2(complex z){

struct complex t;

t.re = z.re*z.re - z.im*z.im;

t.im = 2.0*z.re*z.im;

return t;

};

struct complex cpow3(complex z){

struct complex t;

t.re = pow(z.re,3) - 3.0*z.re*z.im*z.im;

t.im = 3.0*z.re*z.re*z.im - pow(z.im,3);

return t;

};

struct complex cpow4(complex z){

struct complex t;

t.re = pow(z.re,4) + pow(z.im,4) - 6.0*z.re*z.re*z.im*z.im;

t.im = 4.0*(pow(z.re,3)*z.im - pow(z.im,3)*z.re);

return t;

};

struct complex cadd(complex c, complex z) {

struct complex t;

t.re = c.re + z.re;

t.im = c.im + z.im;

return t;

}

struct complex csubs(complex c, complex z) { // c minus z

struct complex t;

t.re = c.re - z.re;

t.im = c.im - z.im;

return t;

}

struct complex cmult(complex c, complex z) {

struct complex t;

t.re = c.re*z.re - c.im*z.im;

t.im = c.re*z.im + c.im*z.re;

return t;

}

struct complex cdiv(complex c, complex z) { // c divided by z

struct complex t;

t.re = (c.re*z.re + c.im*z.im)/(z.re*z.re + z.im*z.im);

t.im = (c.im*z.re - c.re*z.im)/(z.re*z.re + z.im*z.im);

return t;

}

double cabs(complex z) {

return sqrt(z.re*z.re + z.im*z.im);

}

And here’s the Newton fractal sketch. Note the #include “functions.h” on the first line.

#include "functions.h"
#include <SPI.h>
#include "Adafruit_GFX.h"
#include "Adafruit_HX8357.h"

#define TFT_CS D8
#define TFT_DC D1
#define TFT_RST -1  
Adafruit_HX8357 tft = Adafruit_HX8357(TFT_CS, TFT_DC, TFT_RST);

#define BLACK   0x0000
#define BLUE    0x001F
#define RED     0xF800
#define GREEN   0x07E0

// r1, r2 and r3 are the complex roots of z^3 - 1 = 0
struct complex r1 = create(1.0,0.0);
struct complex r2 = create(-0.5,0.8660254038);
struct complex r3 = create(-0.5,-0.8660254038);

struct complex real1 = create(1.0,0.0);
struct complex real3 = create(3.0,0.0);

int max_count = 128;
float precision = 0.0001;
int col_factor = 14;

struct complex z;

void setup() {
  tft.begin(HX8357D);
  tft.setRotation(3);  // landscape
  tft.fillScreen(BLACK);
  yield();
  for(int y=0;y<320;y++) {
    for(int x=0;x<480;x++) {
       int  count = 0;
       z.re = (x-240.0)/80.0; // -3 <= z.re < +3
       z.im = (y-160.0)/80.0; // -2 <= z.im < +2
       while( (count < max_count) && (cabs(csubs(z, r1)) >= precision) && (cabs(csubs(z, r2)) >= precision) && (cabs(csubs(z, r3)) >= precision) ) {
          if(cabs(z) > 0) {
            z = csubs(z,cdiv((csubs(cpow3(z), real1)),(cmult(cpow2(z), real3))));
          }
          yield();
          count++;
       }
       unsigned int kleur = max(0,255-col_factor*count); 
       if (cabs(csubs(z, r1)) < precision) tft.drawPixel(x,y,rgb24to16(0,0,kleur));
       if (cabs(csubs(z, r3)) < precision) tft.drawPixel(x,y,rgb24to16(kleur,0,0));   
       if (cabs(csubs(z, r2)) < precision) tft.drawPixel(x,y,rgb24to16(0,kleur,0)); 
       yield();
    }
    yield();
  }
}

void loop() {
  delay(100);
  yield();
}

unsigned int rgb24to16(byte r_value, byte g_value, byte b_value) {
  unsigned int tint = ((((r_value>>3)<<11) | ((g_value>>2)<<5) | (b_value>>3)));
  return tint;
}

#include "functions.h"

#include <SPI.h>

#include "Adafruit_GFX.h"

#include "Adafruit_HX8357.h"

#define TFT_CS D8

#define TFT_DC D1

#define TFT_RST -1

Adafruit_HX8357 tft = Adafruit_HX8357(TFT_CS, TFT_DC, TFT_RST);

#define BLACK 0x0000

#define BLUE 0x001F

#define RED 0xF800

#define GREEN 0x07E0

// r1, r2 and r3 are the complex roots of z^3 - 1 = 0

struct complex r1 = create(1.0,0.0);

struct complex r2 = create(-0.5,0.8660254038);

struct complex r3 = create(-0.5,-0.8660254038);

struct complex real1 = create(1.0,0.0);

struct complex real3 = create(3.0,0.0);

int max_count = 128;

float precision = 0.0001;

int col_factor = 14;

struct complex z;

void setup() {

tft.begin(HX8357D);

tft.setRotation(3); // landscape

tft.fillScreen(BLACK);

yield();

for(int y=0;y<320;y++) {

for(int x=0;x<480;x++) {

int count = 0;

z.re = (x-240.0)/80.0; // -3 <= z.re < +3

z.im = (y-160.0)/80.0; // -2 <= z.im < +2

while( (count < max_count) && (cabs(csubs(z, r1)) >= precision) && (cabs(csubs(z, r2)) >= precision) && (cabs(csubs(z, r3)) >= precision) ) {

if(cabs(z) > 0) {

z = csubs(z,cdiv((csubs(cpow3(z), real1)),(cmult(cpow2(z), real3))));

}

yield();

count++;

}

unsigned int kleur = max(0,255-col_factor*count);

if (cabs(csubs(z, r1)) < precision) tft.drawPixel(x,y,rgb24to16(0,0,kleur));

if (cabs(csubs(z, r3)) < precision) tft.drawPixel(x,y,rgb24to16(kleur,0,0));

if (cabs(csubs(z, r2)) < precision) tft.drawPixel(x,y,rgb24to16(0,kleur,0));

yield();

}

yield();

}

void loop() {

delay(100);

yield();

}

unsigned int rgb24to16(byte r_value, byte g_value, byte b_value) {

unsigned int tint = ((((r_value>>3)<<11) | ((g_value>>2)<<5) | (b_value>>3)));

return tint;

}

I haven’t been working on larger projects for a while, so here’s a couple of recent chips from the workbench.

In my ongoing quest for a faster HX8357D display driver, I was excited to read that the fast TFT_eSPI library now claims to support this beautiful but rather slow display. Unfortunately, the results were disappointing (unstable, readPixel freezes the screen). However, a nice bycatch of my experiments was the discovery that display updates become noticeably faster when the esp8266’s runs at 160MHz. Very useful for my fractal sketches, and for faster aircraft position updates on the Flight Radar‘s map.

Then I wondered how fast an esp8266 could actually drive my ili9341 driven 320×240 display. So I used my fast ‘offset-modulo’ Koch Snowflake algorithm for a snowflake animation speed test.

The result: an esp8266 at 160MHz can draw more than 100 (iteration depth 3) snowflakes per second on an ili9341 display, using the TFT_eSPI library. That’s about 5x faster than Adafruit’s HX8357 library can achieve on my 480×320 display. Hopefully, there will be a fast and reliable library for that display some day. Update: the TFT_eSPI library now supports driving the HX8457 in parallel mode and the results are impressive.

My renewed attention to the Koch Snowflake led to the discovery of a fractal that was new to me: the Plasma fractal. A commonly used algorithm to produce it is the Diamond-square algorithm. It’s apparently being used to generate nature-like backgrounds for video games.

I wrote a small demo sketch that produces random clouds (well, sort of…), but it can also be used to produce (non-musical) rock formations or landscapes.

6 samples of randomly generated clouds

(picture quality suffers from camera pixel interference with the TFT display)

The Diamond-square algorithm assigns initial values to the four corners of a (2ⁿ+1)x(2ⁿ+1) grid. Then it starts a process of calculating the remaining elements of the grid, based on neighbour values and a random noise of decreasing magnitude. The visual result depends on the four initial corner values, the random noise and the chosen color mapping.

The name Plasma fractal becomes clear when we show consecutive results in a loop. While still using random noise, the trick is to initialize the random generator with the same seed value for every new loop cycle. By slightly changing the initial corner values for each new cycle, the results will also be slightly different from the previous one, which is exactly what we want for a smooth animation.

The general formula that I used for initializing corner k with value h[k] is:

h[k] = d[k] + a[k]*sin( t*f[k] ) where t is a loop counter.

So the value of corner k will periodically swing around value d[k] with amplitude a[k]. The period of the swing is determined by f[k].

Now we can play with the d, a and f values for each corner, as well as with the color mapping, in order to find visually appealing animations. Here’s a first impression (esp8266 & ili9341 320×240 display). The sketch uses a 65×65 grid, with every point filling a corresponding 3×4 rectangle on the display with it’s color value. Since grid values are stored in 4-byte float variables, a 129×129 grid would become too large for the esp8266.

Used values:

const float a[4] = {70.0,80.0,50.0,120.0};
const float d[4] = {0.0,0.0,0.0,0.0};
const float f[4] = {0.053,0.101,0.167,0.211};

const float a[4] = {70.0,80.0,50.0,120.0};

const float d[4] = {0.0,0.0,0.0,0.0};

const float f[4] = {0.053,0.101,0.167,0.211};

Assignment:

h[0][0] = d[0] + a[0]*sin(t*f[0]);
h[64][0] = d[1] + a[1]*sin(t*f[1]);
h[64][64] = d[2] + a[2]*sin(t*f[2]);  
h[0][64] = d[3] + a[3]*sin(t*f[3]);

h[0][0] = d[0] + a[0]*sin(t*f[0]);

h[64][0] = d[1] + a[1]*sin(t*f[1]);

h[64][64] = d[2] + a[2]*sin(t*f[2]);

h[0][64] = d[3] + a[3]*sin(t*f[3]);

Color mapping:

// ci is an index for a 128 colors rainbow array
// h[c][r] is the value of grid element at row r and column c

ci = (128+(int)(h[c][r]))%128;

// ci is an index for a 128 colors rainbow array

// h[c][r] is the value of grid element at row r and column c

ci = (128+(int)(h[c][r]))%128;

Month: June 2018

Newton fractal

Snowflakes & Clouds