Fractals/color mandelbrot

Coloring algorithms

In addition to plotting the set, a variety of algorithms have been developed to

• efficiently color the set in an aesthetically pleasing way
• show structures of the data (scientific visualisation)

general coloring algorithms: ( feature of the pixel -> color )

• discrete gradient:
  • last iteration (integer value) = Level Set Method = LSM
• continous gradient: takes normalised position (floating point number in [0.0, 1.0] range) as an input and gives color as an output, not related with fractals, but used in general computer graphic; $T:scalarvalue\to colorvalue$

Passing iterations into a color directly

One thing we may want to consider is avoiding having to deal with a palette or color blending at all. There are actually a handful of methods we can leverage to generate smooth, consistent coloring by constructing the color on the spot.

Here

• v refers to a normalized exponentially mapped cyclic iter count
• f(v) refers to the sRGB transfer function

A naive method for generating a color in this way is by directly scaling v to 255 and passing it into RGB as such

rgb = [v * 255, v * 255, v * 255]

One flaw with this is that RGB is non-linear due to gamma; consider linear sRGB instead. Going from RGB to sRGB uses an inverse companding function on the channels. This makes the gamma linear, and allows us to properly sum the colors for sampling.

srgb = [v * 255, v * 255, v * 255]

Continuous (smooth) coloring

This image was rendered with the escape time algorithm. There are very obvious "bands" of color

This image was rendered with the normalized iteration count algorithm. The bands of color have been replaced by a smooth gradient. Also, the colors take on the same pattern that would be observed if the escape time algorithm were used.

The escape time algorithm is popular for its simplicity. However, it creates bands of color, which, as a type of aliasing, can detract from an image's aesthetic value. This can be improved using an algorithm known as "normalized iteration count", which provides a smooth transition of colors between iterations. The algorithm associates a real number $\nu$ with each value of z by using the connection of the iteration number with the potential function. This function is given by

$\phi (z)=\lim _{n\to \infty }(\log |z_{n}|/P^{n}),$

where z_n is the value after n iterations and P is the power for which z is raised to in the Mandelbrot set equation (z_n+1 = z_n^P + c, P is generally 2).

If we choose a large bailout radius N (e.g., 10¹⁰⁰), we have that

$\log |z_{n}|/P^{n}=\log(N)/P^{\nu (z)}$

for some real number $\nu (z)$, and this is

$\nu (z)=n-\log _{P}(\log |z_{n}|/\log(N)),$

and as n is the first iteration number such that |z_n| > N, the number we subtract from n is in the interval [0, 1).

For the coloring we must have a cyclic scale of colors (constructed mathematically, for instance) and containing H colors numbered from 0 to H − 1 (H = 500, for instance). We multiply the real number $\nu (z)$ by a fixed real number determining the density of the colors in the picture, take the integral part of this number modulo H, and use it to look up the corresponding color in the color table.

For example, modifying the above pseudocode and also using the concept of w:linear interpolation would yield

for each pixel (Px, Py) on the screen do
    x0:= scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2.5, 1))
    y0:= scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-1, 1))
    x:= 0.0
    y:= 0.0
    iteration:= 0
    max_iteration:= 1000
    // Here N = 2^8 is chosen as a reasonable bailout radius.

    while x*x + y*y ≤ (1 << 16) and iteration < max_iteration do
        xtemp:= x*x - y*y + x0
        y:= 2*x*y + y0
        x:= xtemp
        iteration:= iteration + 1

    // Used to avoid floating point issues with points inside the set.
    if iteration < max_iteration then
        // sqrt of inner term removed using log simplification rules.
        log_zn:= log(x*x + y*y) / 2
        nu:= log(log_zn / log(2)) / log(2)
        // Rearranging the potential function.
        // Dividing log_zn by log(2) instead of log(N = 1<<8)
        // because we want the entire palette to range from the
        // center to radius 2, NOT our bailout radius.
        iteration:= iteration + 1 - nu

    color1:= palette[floor(iteration)]
    color2:= palette[floor(iteration) + 1]
    // iteration % 1 = fractional part of iteration.
    color:= linear_interpolate(color1, color2, iteration % 1)
    plot(Px, Py, color)

Exponentially mapped and cyclic iterations

Typically when we render a fractal, the range of where colors from a given palette appear along the fractal is static. If we desire to offset the location from the border of the fractal, or adjust their palette to cycle in a specific way, there are a few simple changes we can make when taking the final iteration count before passing it along to choose an item from our palette.

When we have obtained the iteration count, we can make the range of colors non-linear. Raising a value normalized to the range [0,1] to a power n, maps a linear range to an exponential range, which in our case can nudge the appearance of colors along the outside of the fractal, and allow us to bring out other colors, or push in the entire palette closer to the border.

$v=((\mathbf {i} /max_{i})^{\mathbf {S} }\mathbf {N} )^{1.5}{\bmod {\mathbf {N} }}$

where:

• i is our iteration count after bailout
• max_i is our iteration limit
• S is the exponent we are raising iters to
• N is the number of items in our palette.

This scales the iter count non-linearly and scales the palette to cycle approximately proportionally to the zoom.

We can then plug v into whatever algorithm we desire for generating a color.

HSV coloring

HSV Coloring can be accomplished by mapping iter count from [0,max_iter) to [0,360), taking it to the power of 1.5, and then modulo 360.

We can then simply take the exponentially mapped iter count into the value and return

hsv = [powf((i / max) * 360, 1.5) % 360, 100, (i / max) * 100]

This method applies to HSL as well, except we pass a saturation of 50% instead.

hsl = [powf((i / max) * 360, 1.5) % 360, 50, (i / max) * 100]

where:

$h=({\frac {i}{max_{i}}}360)^{1.5}\ mod\ 360$

LCH coloring

One of the most perceptually uniform coloring methods involves passing in the processed iter count into LCH. If we utilize the exponentially mapped and cyclic method above, we can take the result of that into the Luma and Chroma channels. We can also exponentially map the iter count and scale it to 360, and pass this modulo 360 into the hue.

${\textstyle {\begin{array}{lcl}x&\in &\mathbb {Q+} \\s_{i}&=&(i/max_{i})^{\mathbf {x} }\\v&=&1.0-cos^{2}(\pi s_{i})\\L&=&75-(75v)\\C&=&28+(75-75v)\\H&=&(360s_{i})^{1.5}{\bmod {3}}60\end{array}}}$

One issue we wish to avoid here is out-of-gamut colors. This can be achieved with a little trick based on the change in in-gamut colors relative to luma and chroma. As we increase luma, we need to decrease chroma to stay within gamut.

s = iters/max_i;
v = 1.0 - powf(cos(pi * s), 2.0);
LCH = [75 - (75 * v), 28 + (75 - (75 * v)), powf(360 * s, 1.5) % 360];

• 
  LCH gradient
• 
  Example of exponentially mapped cyclic LCH coloring without shading
• 
  Exponential Cyclic Coloring in LCH color space with shading

Histogram coloring

A more complex coloring method involves using a histogram which pairs each pixel with said pixel's maximum iteration count before escape/bailout. This method will equally distribute colors to the same overall area, and, importantly, is independent of the maximum number of iterations chosen.

This algorithm has four passes. The first pass involves calculating the iteration counts associated with each pixel (but without any pixels being plotted). These are stored in an array: IterationCounts[x][y], where x and y are the x and y coordinates of said pixel on the screen respectively.

The top row is a series of plots using the escape time algorithm for 10000, 1000 and 100 maximum iterations per pixel respectively. The bottom row uses the same maximum iteration values but utilizes the histogram coloring method. Notice how little the coloring changes per different maximum iteration counts for the histogram coloring method plots.

The first step of the second pass is to create an array of size n, which is the maximum iteration count: NumIterationsPerPixel. Next, one must iterate over the array of pixel-iteration count pairs, IterationCounts[][], and retrieve each pixel's saved iteration count, i, via e.g. i = IterationCounts[x][y]. After each pixel's iteration count i is retrieved, it is necessary to index the NumIterationsPerPixel by i and increment the indexed value (which is initially zero) -- e.g. NumIterationsPerPixel[i] = NumIterationsPerPixel[i] + 1 .

for (x = 0; x < width; x++) do
    for (y = 0; y < height; y++) do
        i:= IterationCounts[x][y]
        NumIterationsPerPixel[i]++

The third pass iterates through the NumIterationsPerPixel array and adds up all the stored values, saving them in total. The array index represents the number of pixels that reached that iteration count before bailout.

total: = 0
for (i = 0; i < max_iterations; i++) do
    total += NumIterationsPerPixel[i]

After this, the fourth pass begins and all the values in the IterationCounts array are indexed, and, for each iteration count i, associated with each pixel, the count is added to a global sum of all the iteration counts from 1 to i in the NumIterationsPerPixel array . This value is then normalized by dividing the sum by the total value computed earlier.

hue[][]:= 0.0
for (x = 0; x < width; x++) do
    for (y = 0; y < height; y++) do
        iteration:= IterationCounts[x][y]
        for (i = 0; i <= iteration; i++) do
            hue[x][y] += NumIterationsPerPixel[i] / total /* Must be floating-point division. */

...

color = palette[hue[m, n]]

...

Finally, the computed value is used, e.g. as an index to a color palette.

This method may be combined with the smooth coloring method below for more aesthetically pleasing images.

Muency

color From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2023

Java script version by Albert Lobo

// https://github.com/llop/mandelbrot-viewer-js/blob/master/mandelbrot-viewer.js
// http://axonflux.com/handy-rgb-to-hsl-and-rgb-to-hsv-color-model-c
  _hsvToRgb(h, s, v) {
    let r, g, b;
    const i = Math.floor(h * 6);
    const f = h * 6 - i;
    const p = v * (1 - s);
    const q = v * (1 - f * s);
    const t = v * (1 - (1 - f) * s);
    switch (i % 6) {
      case 0: r = v, g = t, b = p; break;
      case 1: r = q, g = v, b = p; break;
      case 2: r = p, g = v, b = t; break;
      case 3: r = p, g = q, b = v; break;
      case 4: r = t, g = p, b = v; break;
      case 5: r = v, g = p, b = q; break;
    }
    return [ r * 255, g * 255, b * 255 ];
  }
  
  
  // https://mrob.com/pub/muency/color.html
  _colorCheckered(k) {
    if (this.ns[k] >= this.maxN) return [ 255, 255, 255 ];
    
    const dwell = Math.floor(this.dwell[k]);
    const finalrad = this.dwell[k] - Math.floor(this.dwell[k]);
    const dscale = Math.log2(this.dist[k] / this.inc);
    
    let value = 0.0;
    if (dscale > 0.0) value = 1.0;
    else if (dscale > -10.0) value = (10.0 + dscale) / 10.0;
    
    let p = Math.log(dwell) / Mandelbrot.LOG_BIG_NUM;
    let angle = 0.0;
    let radius = 0.0;
    if (p < 0.5) {
      p = 1.0 - 1.5 * p;
      angle = 1.0 - p;
      radius = Math.sqrt(p);
    } else {
      p = 1.5 * p - 0.5;
      angle = p;
      radius = Math.sqrt(p);
    }
    
    if (dwell % 2) {
      value *= 0.85;
      radius *= 0.667;
    }
    
    if (this.finalang[k] > 0.0) {
      angle += 0.02;
    }
    angle += 0.0001 * finalrad;
    
    let hue = angle * 10.0;
    hue -= Math.floor(hue);
    let saturation = radius - Math.floor(radius);
    
    return this._hsvToRgb(hue, saturation, value);
  }
  
  // b&w version of checkerboard
  _colorCheckeredBlackAndWhite(k) {
    const color = this._colorCheckered(k);
    const gray = Math.round(color[0] * 0.299 + color[1] * 0.587 + color[2] * 0.114);
    return [ gray, gray, gray ];
  }
  
  
  // color pixel k
  _color(k) {
    if (this.ns[k] == 0) return [ 0, 0, 0 ];  // C not yet processed: color black
    this.pix[k] = 0;                          // mark pixel as colored
    return this.colorFunc(k);                 // return color for pixel
  }
  
  
  // paint what we have on the canvas
  render() {
    if (!this.scanDone) {
      let offset = 0;
      for (let k = 0; k < this.imgSize; ++k, offset += 4) {
        if (this.pix[k]) {
          const color = this._color(k);
          this.image.data[offset] = color[0];
          this.image.data[offset + 1] = color[1];
          this.image.data[offset + 2] = color[2];
          this.image.data[offset + 3] = 255;
        }
      }
      if (!this.scanning) this.scanDone = true;
    }
    this.context.clearRect(0, 0, this.imgWidth, this.imgHeight);
    this.context.putImageData(this.image, 0, 0);
  }
  
  
  // sleep function. use 'await this.sleep()' in async functions
  _sleep() { return new Promise(requestAnimationFrame); }
  
}

See also

• computing the Mandelbrot Set Exterior
• coloring the dynamic plane and th Julia and the Fatou sets
• general coloring algorithms: ( feature of the pixel -> color ): $T:scalarvalue\to colorvalue$