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
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 zn is the value after n iterations and P is the power for which z is raised to in the Mandelbrot set equation (zn+1 = znP + c, P is generally 2).

If we choose a large bailout radius N (e.g., 10100), 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 |zn| > 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.

$\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 360 \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];

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