The algorithms are presented using a C-like pseudocode that you should be able to adapt to the language of your choice. Even Python.
The code, as presented, treats complex numbers as individual (x,y) components and performs primitive arithmetic operations on those. Some bits might say simply "calculate |z|", but you should know how to compute the modulus from the individual components without having it explicitly spelled out. If you're using a language that has complex numbers as a primitive type, such as C99, then you should be able to implement many of these algorithms much more simply - you can divide complex numbers using the '/' operator and the compiler will generate all the boilerplate and presumably do so very efficiently. If, however, you're using a language where complex number support is via a library that gives you Complex objects with arithmetic methods called "add", "conjugate", "multiply", etc. that return new instances of the Complex class, steer well clear. Fractal generation typically requires billions of floating point calculations. This is not a problem - floating point performance these days is measured in gigaflops. However, creating billions of objects on the heap and then freeing them when they go out of scope is and will make your runtime hours rather than seconds. Using a built-in list type should be avoided as well; even if you're not reallocating Array objects, storing a value in an Array is orders of magnitude slower than a primitive arithmetic operation. Inside the iteration loops, you should also avoid function calls where possible, since the overhead of tens of millions of function calls is noticeable. That said, be sensible: it is much faster to calculate the exponent via a library call than to calculate the Taylor series approximation using basic arithmetic.
Things that are in ALL_CAPS are assumed to defined outside the functions (e.g. as object properties).
Smooth colouring requires the ability to obtain a colour in between two colours picked from a discrete palette:
function: getIntermediateColourValue
parameters: c1(R,G,B), c2(R',G',B'), m: 0 <= m <= 1
return: [R + m * (R' - R),
G + m * (G' - G),
B + m * (B' - B)]
Fractal images are typically coloured according to the iteration count after some bailout condition has been met. It is important not to use too few colours. Otherwise images of the fractal boundary will have a different, unrelated colour for every pixel, since the iteration counts tend to vary wildly. This will produce colour noise which will destroy detail and result in an ugly image.
define: paletteObject
properties: controlColours: Array(0..2N); RGB start points and end points of N colour channels
weightingPerGroup: Array(0..N); relative importance of each colour group; only used for the first cycle
coloursPerChannel: integer; number of colours in a group when the weighting is zero
Choose your control colours. If you're a colour theory expert, this should be a doddle. If not, don't use arithmetical relationships between R, G and B values since this will produce a really bad palette. What I do is use a decent colour picker that shows colours by their perceptual value. I then trace a diagonal line from a dark tone to a light tone and use the two end values as a start point and end point of a colour group. This provides lots of intermediate, perceptually related tones.
The weighting array is used to tune the palette for the fractal image. Typically, you'll want the first group transition to occur at the fractal boundary. You can use a negative weighting value to create a smaller colour group. You'll also probably want to assign a higher weighting to green tones than red / orange tones, since the eye is more sensitive to greens than reds and oranges.
A paletteObject is passed to a palette expansion function. The expansion function ensures two things:
function: makePalette
parameters: paletteDef: paletteObject; maxColours: integer >= 0 to limit the palette size with 0 signifying "unlimited"
let:
p = [],
toGenerate = maxColours > 0 : min(MAX_ITER, maxColours) : MAX_ITER,
controlColours = paletteDef.controlColours,
weightingPerGroup = paletteDef.weightingPerGroup,
coloursPerChannel = paletteDef.coloursPerChannel,
totalBias = 0,
weighting2 = [],
bias2 = 0,
groupCount = min(weightingPerGroup.length, controlColours.length / 2)
generated = 0,
thisRound = min(toGenerate, DEF_MAX_ITER);
// Calculate total bias and set up even weighting for subsequent rounds
for (i = 0..groupCount - 1) {
totalBias += weightingPerGroup[i];
weighting2[i] = 1;
++bias2;
}
while (generated < toGenerate) {
let:
coloursPerPointOfBias = ceil((thisRound - groupCount * coloursPerGroup) / totalBias);
for (i = 0..groupCount - 1) {
let:
idx = 2 * i,
c1 = controlColours[idx],
c2 = controlColours[idx + 1],
weighting = weightingPerGroup[i],
// The number of colours to create in this channel
colours = coloursPerGroup + coloursPerPointOfBias * weighting,
fInterval = 1 / colours;
// Add start point
push(p, c1);
// Intermediate points
for (j = 1..colours - 1) {
push(p, getIntermediateColourValue(c1, c2, j * fInterval);
}
// End point
push(p, c2);
}
generated += thisRound;
// Subsequent rounds use even weighting
weightingPerGroup = weighting2;
totalBias = bias2;
// Each subsequent round generates twice as many colours as the last
thisRound *= 2;
}
// White and black for set interior points
push(p, WHITE, BLACK);
return p;
Certain colouring operations - basins of attraction, distance shading - need to be done in the HSV domain. Therefore, we need to be able to convert RGB to HSV and HSV to RGB.
function rgb2hsv
parameters: r, g, b: integer in the range 0 - 255
// Convert to values between 0 and 1 such that 255 => 1.0
r = (r & 0xFF) / 255;
g = (g & 0xFF) / 255;
b = (b & 0xFF) / 255;
let:
Min = min(r, g, b),
Max = max(r, g, b);
if (Max == Min) {
// Greyscale
return [0, 0, Max];
}
let:
h = 0,
d = Max - Min,
s = d / Max,
v = Max;
// Get the sector in which the hue falls; r g b are 120° apart
when Max is
r:
h = (g - b) / d;
g:
h = (b - r) / d + 2;
b:
h = (r - g) / d + 4;
// Normalize h so that 0 <= h < 1
if (h < 0) {
h = h + 6;
}
return [h / 6, s, v];
function hsv2rgb
parameters: hue, saturation, value: real: normalized such that 0 <= v <= 1
if (hue == 1) {
// Hue is an angle around a colour wheel
hue = 0;
}
// Colour wheel is divided into six sectors
hue = hue * 6;
let:
C = saturation * value,
inter = hue % 2 - 1,
m = value - C;
if (inter < 0) {
inter = -inter;
}
let:
X = C * (1 - inter),
_r = 0,
_g = 0,
_b = 0,
sector = floor(hue);
when sector is
0:
_r = C;
_g = X;
1:
_r = X;
_g = C;
2:
_g = C;
_b = X;
3:
_g = X;
_b = C;
4:
_r = X;
_b = C;
5:
_r = C;
_b = X;
return [(_r + m) * 255, (_g + m) * 255, (_b + m) * 255];
If we are able to calculate a distance estimate to the fractal boundary we can use it to adjust the v-value in HSV colour space so that pixels closer to the boundary are darker. This allows us to pick out the boundary without having explicitly to colour pixels black. We use the boundaryFraction parameter to emphasize (small values) or de-emphasize the boundary.
function colourDistance
parameters: h, s, v: real: components of the HSV colour space;
distance: real: the distance estimate calculated during iteration;
pixelWidth: real: the width represented by a pixel on the screen;
boundaryFraction: real: value > 0 used to control prominence of the boundary
let:
dScale = log2(distance / pixelWidth),
factor = 0;
if (dScale > 0) {
factor = 1;
}
else if (dScale > -boundaryFraction) {
factor = (boundaryFraction + dScale) / boundaryFraction;
}
// Darken v when factor < 1
return hsv2rgb(h, s, v * factor);
The iteration function calculates an index into a palette which is generally real-valued allowing for continuous colouring and a distance estimation from the set boundary. The distance estimate may be null, such as when the Newton iteration function fails to converge. Given a palette (generated using the algorithm outlined above) we can use these two values to obtain a colour value for the pixel.
function lookupColour
parameters: iter: real: the iteration count before bailout (negative value signifies set interior);
distance: real: the distance estimate calculated during iteration;
p: array: RGB triplets
let:
// Last two colours are used for set interior
maxC = length(p) - 3,
r = 0, g = 0, b = 0, idx1 = 0, idx2 = 0, frac = 0;
if (iter >= 0) {
idx1 = floor(iter) % maxC,
idx2 = idx1 + 1;
if (idx2 > maxC) {
idx2 = 0;
}
frac = iter - floor(iter);
}
else {
// Set interior
idx1 = idx2 = length(p) + iter;
}
(r, g, b) = getIntermediateColourValue(p[idx1], p[idx2], frac);
if (distance) {
(h, s, v) = rgb2hsv(r, g, b);
(r, g, b) = colourDistance(h, s, v, distance, PIXEL_WIDTH, BOUNDARY_FRACTION);
}
// One of the colour options for a dark interior / light boundary
if (REVERSED_COLOURS) {
r = 255 - r;
g = 255 - g;
b = 255 - b;
}
return [r, g, b];
This takes an integer iteration value and turns it into a real value that can be fed into the lookupColour routine. It takes the escape value into account when the bailout condition has been reached. Identify the order of the polynomial, P. If P <= 0 we can't use this function because of the need to take a logarithm; in this case we simply use the integer value.
function: normalizeIterationCount_rp
parameters: r: the modulus of zn; logp: the natural logarithm of P; n: the iteration count at bailout
let:
inter = log2(r);
if (inter <= 0) {
return n;
}
let:
nu = 1 - log(inter) / logp;
return max(n + nu, 0);
To use it, calculate the iteration count, distance and zn as normal, then:
let:
r = sqrt(xn * xn + yn * yn),
niter = normalizeIterationCount_rp(r, log(P), iter),
return lookupColour(niter, distance, palette);
This is used as part of the basin colouring algorithm (see below). It is strictly correct only when the point converged on is a distinct root of the function. It can be used for general-purpose colouring but I personally prefer the results from the exponential smoothing algorithm.
Calculate the difference between two iteration values when the convergence criterion has been met. Then:
let:
q = log(log(diff) / log(CONVERGENCE_RADIUS)) / log(2)
This produces a value between 0 and 1 that should be subtracted from the iteration count.
This is the algorithm used for convergent fractals. The colour value needs to be calculated as part of the iteration loop:
let:
s = 0,
diff = 1,
distance = 0;
while (ITERATING) {
s += exp(-1 / diff);
Calculate diff = |(zk - zk-1)|;
Calculate distance estimate
}
return lookupColour(s, distance, PALETTE);
An option for Newtonian fractals is colouring according to which root of the function a point in the complex plane converges to. This algorithm allows for an unlimited number of basins to be coloured, although there is no guarantee that each basin will have a perceptually distinct colour. Colours are distinct with 6th order polynomials and distinguishable up to 10th order.
function makeBasinColour
parameters: rootIdx: the index of the root found by convergence; rootCount: the number of distinct roots of the function;
iter: the number of iterations required for convergence; diff: the absolute difference between the last two iterated values;
distance: the distance estimate calculated during iteration
let:
// Add to the hue angle to make the basin for the principal root an attractive colour
hueAngle = rootIdx / rootCount + 0.75;
// Ensure that adjacent roots don't have similar colours
if (rootIdx % 2) {
hueAngle = hueAngle + 0.5;
}
if (hueAngle >= 1) {
hueAngle = hueAngle - 1;
}
// Base value on normalized iteration count for a super attracting fixed point
let:
q = log(log(diff) / log(CONVERGENCE_RADIUS)) / log(2),
v = 1 - log(iter - q) / log(MAX_ITER);
return colourDistance(hueAngle, 1, v, distance, PIXEL_WIDTH, BOUNDARY_FRACTION);
The roots of the function are assumed to be sorted according to the complex number sorting algorithm given below.
We choose a fudge factor that makes the principal root (the one closest to the positive real axis) a colour that we find pleasing and go on from there. 0.75 corresponds to an angle of 3 * π / 2 and results in a nice shade of purple. A fudge factor of 0 gives pure red. Adjust to taste. Since the angles are evenly spaced around the colour wheel, a large number of roots will result in perceptually indistinct colouring for adjacent basins (since they'll all end up in the same π / 3 sector). Therefore, we make the angle of alternating basins opposite by adding π.
We only use two dimensions of a three-dimensional colour space as I've not been able to figure out a meaningful saturation calculation that produces good results. The value is based on iteration count - the fewer the brighter which makes the roots stand out as distinct bright spots. This also means the more the darker, so zooming into a singularity will make things look distinctly dim. Then again, zooming into a singularity is a fool's errand since there's nothing to see.
Passing the HSV co-ordinates and the distance estimate to the colourDistance routine picks out the boundary very nicely indeed.
Distance estimation is a method to calculate the distance of a point in the complex plane from the boundary of the fractal. We use this distance estimate to darken pixels close to the border. Boundary tracing controls noise, picks out detail and increases the subjective image quality.
We will only consider exterior distance estimation. Interior distance estimation for filled sets is possible, but has a prerequisite of identifying the period length which can be computationally prohibitive for points very close to the set boundary. Plus some seriously scary maths once we've done so.
This algorithm is for divergent fractals where infinity is a super-attracting fixed point and requires that we know the differential of the iteration function f. Since we're dealing with simple powers of z, this is straightforward. For powers of z >= 2, the differential of f is:
2z, 3z2, 4z3, 5z4, etc.
For powers of z <= -2, we have:
-2z-3, -3z-4, -4z-5, etc.
For intermediate powers of z, use the general rule:
d(zP) ----- = PzP-1 dz
The differential has to be calculated as part of the iteration loop. For Julia-style iterations, we are calculating the partial differential with respect to z while for Mandelbrot-style iterations the calculation is with respect to c. Fortunately, we can handle this difference with a simple change in initial values and a constant of addition. For illustration, I'll show the calculation for the z2 + c iteration function. Important: we have to calculate the differential using the results from the last iteration, otherwise we'll get things badly wrong.
Parameters: x, y, cX, cY
let:
x_ = x,
y_ = y,
iter = 0,
eV = 256,
rsquared = (x_ * x_ + y_ * y_),
dzx = ISJULIA ? 1 : 0,
dzy = 0,
dzPlus = ISMANDEL ? 1 : 0,
// Initialize distance to the null value
distance = null;
while (iter < MAX_ITER AND rsquared < eV) {
// Update dz
tmp = 2 * (dzx * x_ - dzy * y_) + dzplus;
dzy = 2 * (dzx * y_ + dzy * x_);
dzx = tmp;
// Update z
tmp = x_ * x_ - y_ * y_ + cX;
y_ = 2 * x_ * y_ + cY;
x_ = tmp;
rsquared = x_ * x_ + y_ * y_;
++iter;
}
At the end of the loop, we have calculated the differential. Now we can calculate the distance estimation:
if (iter < MAX_ITER) {
Let:
r = sqrt(rsquared),
dz = sqrt(dzx * dzx + dzy * dzy);
// This will be fed into lookupColour
distance = log(r * r) * r / dz;
}
Newtonian fractals are generated by applying the generalized Newton-Raphson root finding method to points in the complex plane:
f(zk)
zk+1 = zk - R·------
fʹ(zk)
z4 - (1+i)z3 - (1-i)z2 - (1-i)z - ionly has three distinct roots, despite being fourth order. It is the product of:
(z + 1)(z + i)(z - 1)2so that +1 is a root with a multiplicity of 2, making it attracting but not super-attracting (and makes Newton's method slow to boot).
For shorthand, let's call Newton's method g(z). To differentiate it we need to apply the quotient rule for differentiating a function of the form p/q:
pʹq - pqʹ --------- q2
Noting that q is fʹ(z) and that qʹ is fʹʹ(z) (the second differential of f), we can obtain the formula for differentiating g(z):
fʹ(z)2 - f(z)fʹʹ(z)
g'(z) = 1 - R·-------------------
fʹ(z)2
function dgz
parameters: fzx, fzy: the components of f(z); dfzx, dfzy: the components of fʹ(z);
d2fzx, d2fzy: the components of fʹʹ(z); rX, rY: the components of the relaxation parameter
let:
denomX = dfzx * dfzx - dfzy * dfzy,
denomY = 2 * dfzx * dfzy,
nomX = denomX - (fzx * d2fzx - fzy * d2fzy),
nomY = denomY - (fzx * d2fzy + fzy * d2fzx);
// Multiply by R
tmp = nomX * rX - nomY * rY;
nomY = nomY * rX + nomX * rY;
nomX = tmp;
// Multiply by conjugate of denominator
tmp = nomX * denomX + nomY * denomY;
nomY = nomY * denomX - nomX * denomY;
nomX = tmp;
denomX = denomX * denomX + denomY * denomY;
return [1 - nomX / denomX, nomY / denomX];
The cumulative differential of f(z) is calculated during the main iteration loop. Note that you have to keep track of the last value of the cumulative differential as well as the current one. Most references will indicate |dz| in the final distance estimate calculation, but this is several orders of magnitude out. The correct factor is |dz - dzlast| (credit: David Makin).
Parameters: x, y
Constants: rX, rY: the components of the relaxation parameter
let:
x_ = x,
y_ = y,
diff = 1,
dzx = 1,
dzy = 0,
dzxl = 0,
dzyl = 0,
eV = 1e-5,
iter = 0,
distance = null,
s = 0;
while (iter < MAX_ITER AND diff > eV) {
// Update colour index
s = s + exp(-1 / diff);
Compute [nomX, nomY] = f(z);
Compute [denomX, denomY] = fʹ(z);
Compute [d2zx, d2zy] = fʹʹ(z);
// Update differential
[_dzx, _dzy] = dgz(nomX, nomY, denomX, denomY, d2zx, d2zy, rX, rY);
dzxl = dzx;
dzyl = dzy;
tmp = dzx * _dzx - dzy * _dzy;
dzy = dzx * _dzy + dzy * _dzx;
dzx = tmp;
// Multiply by R
tmp = nomX * rX - nomY * rY;
nomY = nomY * rX + nomX * rY;
nomX = tmp;
// Multiply top and bottom by conjugate of denominator
tmp = nomX * denomX + nomY * denomY;
nomY = nomY * denomX - nomX * denomY;
nomX = tmp;
denomX = denomX * denomX + denomY * denomY;
if (!denomX) {
return BLACK;
}
// Compute g(z)
_x = x_ - nomX / denomX;
_y = y_ - nomY / denomX;
// Get the differences from last round
diffX = x_ - _x;
diffY = y_ - _y;
diff = sqrt(diffX * diffX + diffY * diffY);
// Set up the next round
x_ = _x;
y_ = _y;
++iter;
}
if (diff < eV) {
// Calculate distance estimate
let:
dz = sqrt((dzx - dzxl) * (dzx - dzxl) + (dzy - dzyl) * (dzy - dzyl));
distance = diff * -log(diff) / dz;
}
return lookupColour(s, distance, PALETTE);
Implementation note: for the convergent formulae, if the iteration loop ends without the escape value being reached, we don't automatically return the set interior colour (which basically means "no idea"). Injudicious values of R can prevent Newton's method converging though over/under relaxation. Nonetheless, the exponential smoothing algorithm still provides a usable colour value so we can produce (sometimes striking) images even when the iteration count is maximum. Assuming that we haven't screwed up Newton's method (R is sensible) and we know the roots that the iteration should converge to, we may still find that some starting values fail to converge on a root. There are certain combinations of coefficients and starting values where Newton's method ends up in a cycle. In this case, returning the set interior colour is the correct thing to do.
When finding roots of a function, we have no expectation of strict equality, since the last value of an iteration is only an approximation of the actual root within some limit of tolerance. Therefore, we need to be able to compare complex numbers for approximate equality.
function complexCompare
parameters: x1, y1, x2, y2: the values to compare; tolerance: the acceptable radius within which values are deemed to be equal
let:
_x = x1 - x2,
_y = y1 - y2,
delta = sqrt(_x * _x + _y * _y);
return (delta < tolerance);
This section assumes a parameterized sort method on a list where one of the parameters is a pointer or reference to a comparison function. It also assumes a list of points [x, y] representing values in the complex plane and that the semantics of comparing two values a, b are:
< 0: a sorts before b > 0: b sorts before a 0: a and b are equal
It also assumes an atan2 function that returns an angle theta:
Upper two quadrants of the complex plane: 0 <= theta <= π Lower two quadrants: 0 > theta > -π
The sort order for complex numbers is undefined. We cannot say which of -1-3i and -2-1i is more negative than the other. So we'll define our own sort order. Here's mine:
The basin colouring algorithm for the Newtonian fractal uses roots sorted using the following algorithm:
let:
theta1 = atan2(a.y, a.x),
theta2 = atan2(b.y, b.x),
r1 = sqrt(a.x * a.x + a.y * a.y),
r2 = sqrt(b.x * b.x + b.y * b.y);
if (theta1 < 0) {
theta1 += 2 * PI;
}
if (theta2 < 0) {
theta2 += 2 * PI;
}
let:
ret = theta1 - theta2;
if (0 == ret) {
ret = r1 - r2;
}
return ret;
Implementation note: compute the angles and moduli before sorting and sort objects with r and theta properties.
Basin colouring for Newtonian fractals requires us to know the roots of the function before we generate the fractal. We can also anticipate a future requirement of identifying critical points of a function for generating Mandelbrot sets. Therefore, we need a generic root finding algorithm. This algorithm uses - what else? - Newton's method.
It assumes the existence of a random function that returns a real value 0 <= x < 1.
definition rootObject
properties: x: the real component of z;
y: the imaginary component of z;
r: the modulus of z;
theta: the angle of z
function findRoots
parameters: f: a pointer to a function that takes a complex number, z, and computes f(z);
fprime: a pointer to a function that takes a complex number, z, and computes fʹ(z)
constants: XMAX, XMIN, YMAX, YMIN: the limits of the complex plane to be tested
let:
eV = 1e-8,
// For checking suspiciously small real / imaginary components
delta = eV * 100,
// For comparing roots found
tolerance = delta * 100,
rSpan = XMAX - XMIN,
iSpan = YMAX - YMIN,
// If we don't find a root after n iterations, give up and move on
iterations = 100,
// We'll test this many points in the x- and y- axes
points = 30,
rInt = rSpan / points,
iInt = iSpan / points,
roots = [];
// Test a reasonably large sample of points in the plane; assuming that we only have to
// do this when we regenerate a fractal, we can use a conservatively large value to minimize
// the chances of failing to find distinct roots that are close to each other
for (i = 0..points) {
for(j = 0..points) {
let:
x = XMIN + i * rInt,
y = YMIN + j * iInt,
diff = 1,
iter = 0;
while (iter < iterations AND diff > eV) {
// Get the numerator and denominator
[nomX, nomY] = f(x, y);
[denomX, denomY] = fprime(x, y);
// Divide f(z) by f'(z) - multiply through by conjugate of denominator
tmp = nomX * denomX + nomY * denomY;
nomY = nomY * denomX - nomX * denomY;
nomX = tmp;
denomX = denomX * denomX + denomY * denomY;
if (!denomX) {
// Try another number
break;
}
// Now we can divide
_x = x - nomX / denomX;
_y = y - nomY / denomX;
// Get the differences from last round
diffX = x - _x;
diffY = y - _y;
diff = sqrt(diffX * diffX + diffY * diffY);
// Set up the next guess
x = _x;
y = _y;
++iter;
}
if (diff < eV) {
let:
newRoot = TRUE;
// Zero suspiciously small real / imaginary components so that the roots can
// be sorted deterministically
if (abs(x) < delta) {
x = 0;
}
if (abs(y) < delta) {
y = 0;
}
for (k = 0..LENGTH(roots)) {
if (complexCompare(x, y, roots[k].x, roots[k].y, tolerance)) {
newRoot = FALSE;
break;
}
}
if (newRoot) {
push(roots, rootObject{ x : x, y : y, r : sqrt(x * x + y * y), theta : atan2(y, x) });
}
}
}
}
SORT roots;
return roots;
The idea behind the algorithm is that multiple independent threads of execution will find the same roots. The algorithm will treat distinct roots that differ by less than the tolerance threshold as multiple identical roots. Since any function for which this was true would be both stupid and deliberately constructed to have very close roots, we'll not care.
Back to the index